text:
\section{Introduction}
The fundamental nature of dark matter (DM) has puzzled scientists for decades. While we have long observed the gravitational pull it exerts on regular baryonic matter \cite{1933AcHPh...6..110Z,1937ApJ....86..217Z,1980ApJ...238..471R}, no conclusive hint of the particle physics governing DM has so far shown up in laboratory experiments (see Refs.~\cite{Aalseth:2010vx,Aalseth:2011wp,Angloher:2011uu,Agnese:2013rvf,Bernabei:2013xsa} for tentative signals). This does not necessarily imply that the physics of DM is trivial or uninteresting; it merely tells us that it couples very weakly to the visible sector, thus allowing it to hide its potentially rich physics. To make progress, one can turn to astrophysical observations of DM dominated objects such as dwarf spheroidals \cite{Walker:2011zu,2012MNRAS.420.2034S,BoylanKolchin:2011de,BoylanKolchin:2011dk,Walker:2012td,Zavala:2012us,Laporte:2013fwa,Amorisco:2013uwa}, galaxies \cite{2004MNRAS.351..903G,deNaray:2009xj,deNaray:2011hy,Vogelsberger:2012ku}, and merging galaxy clusters \cite{Markevitch:2003at,Randall:2007ph,Merten:2011wj,Dawson:2011kf,Dawson:2012fx}.  These objects are however highly non-linear and one cannot in general neglect the impact of baryon physics (see e.g.~Refs~\cite{Oh:2010mc,Pontzen:2011ty,Governato:2012fa,Zolotov:2012xd,Brooks:2012vi,Shen:2013wva}) on their evolution. Numerical simulations are therefore necessary to assess the impact of non-minimal DM physics on these objects.

Since DM dominates the matter density on cosmological scales, it is natural to ask whether cosmological data can shed light on the fundamental physics of DM. While the cold DM paradigm \cite{Blumenthal:1984bp,Davis:1985rj} provides a very good fit to data on large cosmological scales, it is possible that a subdominant DM component could display very different properties. For instance, there could be new dark forces \cite{Foot:2004pa,ArkaniHamed:2008qn,Ackerman:2008gi,Feng:2009mn,Baldi:2012ua,Tulin:2012wi,Hooper:2012cw,Aarssen:2012fx,Tulin:2013teo} that couple only to a fraction of the DM particles or a portion of the DM could be warm \cite{Bode:2000gq,Dalcanton:2000hn,Zentner:2003yd,Smith:2011ev}.

A particularly interesting case is one in which a fraction of the DM can interact with or via a massless (or nearly massless) particle. In this scenario, the interacting DM component is prohibited from forming gravitationally bound structures until it kinematically decouples from the light state. If this decoupling happens at relatively late times (that is, close to or after matter-radiation equality), an imprint similar in many ways to the baryon acoustic oscillation feature should be left on the matter density field on cosmological scales \cite{CyrRacine:2012fz}. This dark acoustic oscillation (DAO) feature generically arises in any model where DM is coupled to relativistic particles until relatively late times. For instance, they occur if DM couples to neutrinos \cite{Mangano:2006mp,Serra:2009uu,Aarssen:2012fx,Diamanti:2012tg} or photons \cite{Boehm:2001hm,McDermott:2010pa,Wilkinson:2013kia}, if DM interacts with a dark $U(1)_D$ gauge boson \cite{Ackerman:2008gi,Feng:2009mn,Kaplan:2009de,Kaplan:2011yj, CyrRacine:2012fz}, if DM couples to a light scalar field \cite{1992ApJ...398..407G}, or for the so-called ``cannibal'' DM models \cite{1992ApJ...398...43C,1994ApJ...431...41M,1995ApJ...452..495D}. We emphasize for the reader unfamiliar with the above body of work that if the DM sector couples purely gravitationally to the visible sector, some of these scenarios are surprisingly unconstrained. For instance, even a model in which the totality of the DM interacts via a massless $U(1)_D$ gauge boson \cite{Ackerman:2008gi,Feng:2009mn} has a large allowed parameter space.

Intriguingly, a subset of these scenarios for which a fraction of the DM couples directly to the light state (as opposed to coupling to the light state via a massive mediator) also has the potential to impact structure formation on galactic scales. Indeed, since a fraction of the DM can lose energy by radiating a light mediator particle, its distribution inside galaxies may differ significantly from that of the standard cold DM component due to the formation of a DM disk \cite{Fan:2013tia,Fan:2013yva,McCullough:2013jma}. As it also generically predicts a DAO feature on large cosmological scales, this Double-disk DM scenario has the unique feature of affecting the DM density field on a large range of scales, hence providing us with a multi-pronged approach to constrain its properties. This contrasts with usually-considered self-interacting DM scenarios where only the smallest astrophysical scales are affected (see e.g.~Refs.~\cite{Rocha:2012jg,Peter:2012jh}).

Moreover, scenarios in which DM couples to a relativistic component naturally incorporates some form of dark radiation (DR), which might be necessary to reconcile the current cosmic microwave background (CMB) data \cite{Hinshaw:2012aka,planckXVI} with recent measurements of the local Hubble expansion rate \cite{Riess:2011yx,Freedman:2012ny}. We emphasize however that the DR in these models is different than that usually parametrized via the quantity $\Delta N_{\rm eff}$, which effectively describes the number of fermionic free-streaming species that have the same temperature as the cosmic neutrino background assuming instant decoupling. Indeed, since DR couples to DM at early times, DR is not free-streaming at all times and therefore does not behave like a neutrino. While this subtle change leaves unchanged the background expansion history of Universe, it \emph{does} affect the way density fluctuations evolve, implying that we cannot straightforwardly interpret the DR in terms of the tacit (and sometimes obscure, see e.g.~Ref.~\cite{Steigman:2013yua}) parameter $\Delta N_{\rm eff}$. As we will see below, the cosmological constraints on the abundance of a DR component that couples to DM until late times can be vastly different than the current limits on $\Delta N_{\rm eff}$, implying that one must exert caution in interpreting bounds on this latter parameter.

In this paper, we use cosmological data from the CMB, baryon acoustic oscillations (BAO), and large-scale structure to constrain the possible fraction of interacting DM as well as the strength of its interaction with light relativistic particles. Since, for the models we consider, interacting DM only affects the cosmological observables through its gravitational interaction, our bounds on interacting DM are very general and apply to any hidden-sector model in which large-scale DAO arises. We will make this correspondence explicit below. As we will see, cosmological data alone place strong limits on both the possible fraction of interacting DM as well as on the strength of its interaction with the light state.

The paper is organized as follow. In section \ref{PIDM}, we present the details of our partially-interacting DM scenario. We then introduce in section \ref{CosmoPert} the physics of DAOs and review the perturbation equations necessary to describe the cosmological evolution of this model. In section \ref{observables}, we describe the consequence of our scenario on the cosmological observables, including galaxy clustering, the CMB, and CMB lensing. The details of the data used for our analysis are given in section \ref{data}. Our main results are presented in section \ref{results} and we discuss their implications for the formation of a dark disk inside galaxies in section \ref{galform}. We discuss the implications of our results in section \ref{sec:conc}.

\section{Partially-Interacting Dark Matter}\label{PIDM}
We consider a hybrid DM sector comprised of particles interacting with a light relativistic state and of a standard non-interacting cold DM component. We shall refer to such a scenario as partially-interacting DM (PIDM).  For definitiveness and simplicity, we take the interacting component to be made of dark atoms \cite{Goldberg:1986nk,Kaplan:2009de,Kaplan:2011yj,Behbahani:2010xa,Petraki:2011mv,Cline:2012is,CyrRacine:2012fz}, a well-studied model that allows us to make exact quantitative predictions. Atomic DM naturally encompasses hidden-charged DM models \cite{Ackerman:2008gi,Feng:2009mn} and can mimic the behavior of CDM in some limits. Importantly, this model readily incorporates a dark radiation (DR) component. Since its impact on DM fluctuations (acoustic oscillations, damping) are very generic, atomic DM can be viewed as a simple toy model to parametrize deviations from a pure CDM scenario. As we discussed in the introduction, the limits we obtain in this work are very general and apply to a variety of models. We briefly review the physics of dark atoms in section \ref{adm} below. We refer the reader to Ref.~\cite{CyrRacine:2012fz} for more details.

Throughout this work, we assume that the relic abundances of both the CDM and the interacting components are set by some unspecified UV physics. Examples of such UV completion are given in Refs.~\cite{Kaplan:2011yj,Petraki:2011mv}. We do not expect the details of the UV completion to affect the low-energy interactions responsible for modifying the growth of DM fluctuations on small scales. To retain generality, we will refer to the fraction of DM made of dark atoms as the ``interacting DM'' component in order to distinguish it from the standard collionsionless CDM. We further assume that the CDM and the interacting component only interact through gravity.

\subsection{Atomic Dark Matter}\label{adm}
In the atomic DM scenario, two oppositely-charged massive fermions can interact through a new unbroken $U(1)_D$ gauge force to form hydrogen-like bound states. In the early Universe, the dark atoms are all ionized by the hot thermal bath of DR (that is, the $U(1)_D$ gauge boson) and form a plasma similar to that of standard baryons and photons. When the temperature of the DR falls significantly below their binding energy, dark atoms are allowed to recombine into neutral bound states if the recombination rate is larger than the expansion rate of the Universe. The DR eventually decouples from the atomic DM and begin to free-stream across the Universe. We note that the order and the dynamics of the different important transitions of the dark plasma (recombination, onset of DR free-streaming, atomic DM drag epoch, DM thermal decoupling, etc.) can be very different than in the standard baryonic case. We refer the reader to Ref.~\cite{CyrRacine:2012fz} for more details.

To retain generality and emphasize that the PIDM scenario we are considering is quite general, we shall refer to the massless $U(1)_D$ ``dark photons'' simply as DR. For simplicity, we also denote the lightest fermion as  ``dark electron''  (mass $m_{\bf e}$) while the heaviest fermion is referred to as  ``dark proton'' (mass $m_{\bf p}$). We assume that these two oppositely-charged components come in equal number such that the dark sector is overall neutral under the $U(1)_D$ interaction. This model is characterized by five parameters which are the mass of the dark atoms $m_D$, the dark fine-structure constant $\alpha_D$, the binding energy of the dark atoms $B_D$, the present-day ratio of the DR temperature ($T_D$) to the cosmic microwave background temperature $\xi\equiv(T_D/T_{\rm CMB})|_{z=0}$, and the fraction of the overall DM density contained in interacting DM (here, dark atoms), $f_{\rm int}\equiv\rho_{\rm int}/\rho_{\rm DM}$, where $\rho_{\rm DM}=\rho_{\rm int}+\rho_{\rm CDM}$ and where $\rho_{\rm int}$ is the energy density of the interacting DM component. These parameters are subject to the consistency condition $m_D/B_D \geq 8/\alpha_D^2-1$, which ensures that the relationship $m_{\bf e}+m_{\bf p}-B_D=m_D$ is satisfied. We note that if the visible and dark sectors were coupled above the electroweak scale, we naturally expect $\xi\sim0.5$ \cite{Fan:2013yva}. A smaller value would either require new degrees of freedom in the visible sector or that the two sectors were never in thermal equilibrium in the first place.

The evolution of the dark plasma is largely governed by the opacity $\tau_D^{-1}$ of the medium to DR.  For the model we considered, the main contributions\footnote{In this work, we neglect the small contribution to the opacity from photoionization processes.} to this opacity are Compton scatterings of DR off charged dark fermions and Rayleigh scatterings off neutral dark atoms, that is,
\begin{equation}\label{opacity}
\tau_D^{-1}=\tau_{\rm Compton}^{-1}+\tau_{\rm R}^{-1},
\end{equation}
where
\begin{equation}
\tau_{\rm Compton}^{-1}=an_{\rm ADM}x_D\sigma_{{\rm T},D}\left[1+\left(\frac{m_{\bf e}}{m_{\bf p}}\right)^2\right],
\end{equation}
and
\begin{eqnarray}\label{rayleigh}
\tau_{\rm R}^{-1}&=&an_{\rm ADM}(1-x_D)\langle\sigma_{\rm R}\rangle\nonumber\\
&\simeq&32\pi^4an_{\rm ADM}(1-x_D)\sigma_{{\rm T},D}\left(\frac{T_D}{B_D}\right)^4.
\end{eqnarray}
Here, $\sigma_{{\rm T},D}\equiv 8\pi\alpha_D^2/(3 m_{\bf e}^2)$ is the dark Thomson cross section, $a$ is the scale factor describing the expansion of the Universe, $x_D$ is the ionized fraction of the dark plasma, $n_{\rm ADM}$ is the number density of dark atoms, $\sigma_{\rm R}$ is the Rayleigh scattering cross section, and where the angular bracket denotes thermal averaging. We note that the second line of Eq.~\ref{rayleigh} is only valid if $T_D<B_D$. It is out of the scope of this paper to discuss in detail the evolution of the ionized fraction and of the DM temperature. We refer the reader to Ref.~\cite{CyrRacine:2012fz} for a thorough investigation of dark atom recombination and thermal history.
\subsection{$\xi$ vs $\Delta N_{\rm eff}$}
We note that, as far as the \emph{background cosmological expansion} is concerned, varying the temperature of the DR is equivalent to changing the effective number of relativistic species $N_{\rm eff}$ according to the mapping
\begin{equation}\label{delta_n_eff}
\Delta N_{\rm eff} \leftrightarrow \frac{8}{7}\left(\frac{11}{4}\right)^{4/3}\xi^4.
\end{equation}
\emph{However}, since $\xi$ affects the evolution of cosmological fluctuations in a different way than $\Delta N_{\rm eff}$ (because the DR couples to DM and is not always free-streaming), we emphasize that one cannot blindly translate the constraints on $\Delta N_{\rm eff}$ from, say, Planck \cite{planckXVI} to a bound on $\xi$. In fact, as we discuss below, the bounds on $\xi$ are generally much more stringent than the equivalent  limit on $\Delta N_{\rm eff}$.

\section{Cosmological Evolution}\label{CosmoPert}

\subsection{Dark Acoustic Oscillation Scale}
Since a fraction of the DM forms a tightly-coupled plasma in the early Universe, the evolution of cosmological fluctuations in the PIDM model departs significantly from that of a standard $\Lambda$CDM Universe. Indeed, as Fourier modes enter the causal horizon, the DR pressure provides a restoring force opposing the gravitational growth of over densities, leading to the propagation of dark acoustic oscillations (DAO) in the plasma. These acoustic waves propagate until DR kinematically decouples from the interacting DM component. Similar to the baryon case, the scale corresponding to the sound horizon of the dark plasma at kinetic decoupling remains imprinted on the matter field at late times. This so-called DAO scale is given by
\begin{equation}\label{rDAO}
r_{\rm DAO}\equiv\int_0^{\eta_{D}}c_D(\eta)d\eta,
\end{equation}
where $c_D$ is the sound speed of the dark plasma, $\eta$ is the conformal time, and $\eta_D$ denotes the conformal time at the epoch at which atomic DM kinematically decouples from the DR bath. The DAO scale is a key quantity of cosmologically-interesting interacting DM models. Indeed, much like the free-streaming length of warm DM models, the DAO scale divides the modes that are strongly affected by the DM interactions (through damping and oscillations) from those that behave mostly like in the CDM paradigm. We note however that, in contrast to warm DM models, the suppression of small-scale fluctuations in the PIDM scenario is mostly due to acoustic (also known as collisional) damping \cite{Boehm:2000gq,CyrRacine:2012fz}, while residual free-streaming after kinematic decoupling can play a minor role.

In the tight-coupling limit of the dark plasma, the sound speed takes the form $c_D=1/\sqrt{3(1+R_D^{-1})}$, where $R_D\equiv4\rho_{\tilde{\gamma}}/3\rho_{\rm int}$. Here, $\rho_{\tilde{\gamma}}$ stands for the the energy density of the DR. In a matter-radiation Universe, the integral of Eq.~\ref{rDAO} can be performed analytically
\begin{eqnarray}\label{r_DAO_theory}
r_{\rm DAO} &=& \frac{4 \xi^2\sqrt{\Omega_{\gamma}}}{3 H_0 \sqrt{f_{\rm int}\Omega_{\rm DM}\Omega_{\rm m}}}\times\\
&&\ln{\left[\frac{\sqrt{\gamma_{\rm int}}\sqrt{\Omega_{\rm r}+\Omega_{\rm m} a_D}+\sqrt{\Omega_{\rm m}+\gamma_{\rm int}a_D}}{\sqrt{\gamma_{\rm int}\Omega_{\rm r}}+\sqrt{\Omega_{\rm m}}}\right]},\nonumber
\end{eqnarray}
where we have defined
\begin{equation}
\gamma_{\rm int} \equiv \frac{3f_{\rm int}\Omega_{\rm DM}}{4\xi^4\Omega_{\gamma}},
\end{equation}
$a_D$ is the scale factor at the epoch of atomic DM kinematic decoupling, and $H_0$ is the present-day Hubble constant. $\Omega_{\gamma}$, $\Omega_{\rm r}$, and $\Omega_{\rm m}$ stand for the energy density in photons, radiation (including neutrinos and DR), and non-relativistic matter, respectively, all in units of the critical density of the Universe. We observe that the DAO scale depends most strongly on the ratio $\xi^2/\sqrt{f_{\rm int}}$ and that the details of the interacting DM microphysics only enter through a logarithmic dependence on $a_D$. The scale factor at the epoch of dark decoupling can be estimated from the criterion $n_{\rm ADM}x_D\sigma_{{\rm T},D}\simeq H$, since Thomson scattering is the dominant mechanism responsible for the opacity of the dark plasma. Here, $H$ is the Hubble parameter. We give an expression for $a_D$ in terms of the dark parameters in Appendix \ref{app:a_D}.

The scale factor at the epoch of dark decoupling (and, consequently, $r_{\rm DAO}$) is largely determined by the following combination of parameters
\begin{equation}\label{eq:pidm_mod}
\Sigma_{\rm DAO}\equiv \alpha_D \left(\frac{B_D}{\rm eV}\right)^{-1}\left(\frac{m_D}{\rm GeV}\right)^{-1/6}.
\end{equation}
\begin{figure}[t!]
\begin{centering}
\includegraphics[width=0.5\textwidth]{rDAO_f_ADM.pdf}
\includegraphics[width=0.5\textwidth]{rDAO_xi.pdf}
\caption{Comoving DAO scale as a function of the parameter $\Sigma_{\rm DAO}$ for strongly-coupled atomic DM models ($\alpha_D>0.025$). In the upper panel, we fix $\xi=0.5$ and vary the fraction of interacting DM. In the lower panel, we fix $f_{\rm int}=5\%$ and let $\xi$ vary. Here, take $H_0=69.57$ km/s/Mpc, $\Omega_{\rm m}=0.3048$, $\Omega_{\rm DM}h^2=0.1198$, and $N_{\rm eff}=3.046$.}
\label{fig:rDAO}
\end{centering}
\end{figure}
This quantity is directly proportional to the scattering rate between DR and interacting DM. Its non-trivial dependence on the dark parameters $\alpha_D$, $B_D$, and $m_D$ is caused by $x_D$ which itself depends on these dark parameters (see Appendix \ref{app:a_D}). To give a sense of scale, we note that for regular baryonic hydrogen we have $\Sigma_{\rm BAO}\simeq5.4\times10^{-4}$. We emphasize that, while the definition given in Eq.~\ref{eq:pidm_mod} is very specific to the atomic DM model considered, $\Sigma_{\rm DAO}$ is a simple proxy for the cross section between DM and DR at the epoch of kinematic decoupling ($\sigma_{\rm DM-DR}(a_D)$) over the DM mass. Explicitly, the relation between $\Sigma_{\rm DAO}$ and the DM-DR cross section over the DM mass is
\begin{eqnarray}
\left(\frac{\sigma_{\rm DM-DR}(a_D)}{m_D}\right)&= &1.9\times10^{-4}\left(\frac{\xi}{0.5}\right)\left(\frac{\Sigma_{\rm DAO}}{10^{-3}}\right)\nonumber\\
&&\,\,\,\times\left(\frac{f_{\rm int}\Omega_{\rm DM}h^2}{0.12}\right)^{-1}\frac{\rm cm^2}{\rm g}.
\end{eqnarray}
It should be clear from the above expression that any constraints we put on $\xi$ and $\Sigma_{\rm DAO}$ can be directly translated to model-independent limits on the DM-DR cross section over the DM mass\footnote{We note that this cross section between DM and DR should not be confused with DM self-interaction cross section relevant for small-scale astrophysical objects.} at kinematic decoupling. In the remainder of this work, we shall parametrize the strength of the DM-DR interaction in terms of $\Sigma_{\rm DAO}$ but the reader should keep in mind that our results can be translated to any model in which a fraction of the DM interacts with a bath of DR.

We plot the DAO scale as a function of $\Sigma_{\rm DAO}$  in Fig.~\ref{fig:rDAO} for different values of the interacting DM fraction (top panel) and for different values of $\xi$ (lower panel). We observe that for $\Sigma_{\rm DAO} > 10^{-3}$ and $\xi>0.2$, the $r_{\rm DAO}$ lies in the range of scales currently probed by galaxy surveys and CMB experiments. Looking ahead, we thus expect these datasets to severely constrain PIDM models lying in this region of parameter space. In particular, a typical double-disk DM model \cite{Fan:2013tia,Fan:2013yva} with $\Sigma_{\rm DAO}\sim10^{-2.5}$, $\xi\sim0.5$, and $f_{\rm int}\sim5\%$, has $r_{\rm DAO}\sim80h^{-1}$Mpc, which is ruled out by current data as we will discuss in section \ref{results}.

\subsection{Cosmological Fluctuations}
We now turn our attention to the evolution of cosmological perturbations in the PIDM scenario. The equations describing the evolution of interacting DM density and velocity fluctuations are
\begin{equation}\label{delta_D}
\dd{\delta}_{D}+\theta_{D}-3\dd{\phi}=0,
\end{equation}
\begin{equation}\label{theta_D}
\dd{\theta}_{D}+\frac{\dd{a}}{a}\theta_D-c_D^2k^2\delta_{D}-k^2\psi= \frac{R_D}{\tau_D}(\theta_{\tilde{\gamma}}-\theta_D),
\end{equation}
where we closely followed the notation of Ref.~\cite{Ma:1995ey} in conformal Newtonian gauge. Here, a dot atop a quantity denotes a derivative with respect to conformal time, $\delta_D$ is the interacting DM density contrast, $\theta_D$ and $\theta_{\tilde{\gamma}}$  are the divergence of the interacting DM and DR velocity, respectively; $\phi$  and $\psi$ are the gravitational scalar potentials, and $k$ is the wavenumber of the mode. Here, the subscript $\tilde{\gamma}$ always refers to the DR. The right-hand side of Eq.~(\ref{theta_D}) represents the collision term between the atomic DM and the DR. At early times, we generally have $R_D\gg1$ and $\tau_D\ll H^{-1}$, implying that the interacting DM component is effectively dragged along by the DR.
The latter evolves according to the following Boltzmann equations:
\begin{equation}\label{delta_gD}
\dd{\delta}_{\tilde{\gamma}}+\frac{4}{3}\theta_{\tilde{\gamma}}-4\dd{\phi}=0;
\end{equation}
\begin{equation}\label{theta_gD}
\dd{\theta}_{\tilde{\gamma}}-k^2(\frac{1}{4}\delta_{\tilde{\gamma}}-\frac{F_{\tilde{\gamma}2}}{2})-k^2\psi= \frac{1}{\tau_D}(\theta_D-\theta_{\tilde{\gamma}});
\end{equation}
\begin{equation}\label{FgD2}
\dd{F}_{\tilde{\gamma}2} = \frac{8}{15}\theta_{\tilde{\gamma}}-\frac{3}{5}kF_{\tilde{\gamma}3}-\frac{9}{10\tau_D}F_{\tilde{\gamma}2};
\end{equation}
\begin{equation}\label{FgDl}
\dd{F}_{\tilde{\gamma}l}=\frac{k}{2l+1}\left[l F_{\tilde{\gamma}(l-1)}-(l+1)F_{\tilde{\gamma}(l+1)}\right]-\frac{1}{\tau_D}F_{\tilde{\gamma}l}.
\end{equation}
Eqs.~(\ref{delta_gD}) and (\ref{theta_gD}) describe the evolution of the DR over-densities ($\delta_{\tilde{\gamma}}$) and of the DR velocity, respectively. It is also necessary to solve for the hierarchy of DR multipoles (Eqs.~(\ref{FgD2}) and (\ref{FgDl})) to properly account for DR diffusion and its impact on interacting DM perturbations. We note that the opacity of the dark plasma can be written in terms of $\Sigma_{\rm DAO}$.

We solve these equations numerically together with those describing the evolution of CDM, baryon, photon, and neutrino fluctuations using a modified version of the code \texttt{CAMB} \cite{Lewis:1999bs}. We first precompute the evolution of the dark plasma opacity as described in Ref.~\cite{CyrRacine:2012fz}. We assume purely adiabatic initial conditions
\begin{equation}
\delta_D(z_{\rm i}) =\delta_{\rm c}(z_{\rm i})\qquad \delta_{\tilde{\gamma}}(z_{\rm i})=\delta_{\gamma}(z_{\rm i}),
\end{equation}
\begin{equation}
\theta_D(z_{\rm i})=\theta_{\tilde{\gamma}}(z_{\rm i})=\theta_{\gamma}(z_{\rm i}),
\end{equation}
\begin{equation}
F_{\tilde{\gamma}l} = 0,\quad l\geq2.
\end{equation}
where $z_{\rm i}$ is the initial redshift which is determined such that all modes of interest are superhorizon at early times, $k\tau(z_{\rm i})\ll1$. Here, the subscripts ``c"  and ``$\gamma$" refer to CDM and regular photon, respectively. At early times when $k\tau_D\ll1$ and $\tau_D/\tau\ll1$, Eqs.~(\ref{theta_D}) and (\ref{theta_gD}) are very stiff and we use a second-order tight-coupling scheme similar to that used for the baryon-photon plasma at early times \cite{2011PhRvD..83j3521C,Pitrou:2010ai}. Due to the presence of tightly-coupled DR at early times, the neutrino initial conditions need to be modified to take into account the different free-streaming fraction.
\section{Cosmological Observables}\label{observables}
Since it modifies the growth of DM fluctuations on a variety of scales, PIDM can imprint its signatures on cosmological observables such as the CMB, CMB lensing, and the matter power spectrum. These observables are mostly sensitive to the momentum transfer rate between the interacting DM and the DR, which itself determines the kinetic decoupling epoch. This rate is largely determined by the parameter $\Sigma_{\rm DAO}$ defined in Eq.~\ref{eq:pidm_mod}. For a relatively strongly-coupled ($\alpha_D > 0.025$) dark sector, any changes to $\alpha_D$, $B_D$, and $m_D$ that leaves $\Sigma_{\rm DAO}$ invariant lead to the same cosmological observables. Indeed, the recombination process for these models is well-described by the Saha approximation until the epoch of dark decoupling, which itself is determined entirely by $\Sigma_{\rm DAO}$ (see Eq.~\ref{eq:sol_a_D}). For smaller values of the dark fine-structure constant however, the details of the dark recombination process becomes important and the observables develop a small explicit dependence on $\alpha_D$. In the following subsections, we focus our attention on strongly-coupled models (unless otherwise stated) which are fully characterized by $\Sigma_{\rm DAO}$, but will address the constraints on weakly-coupled models in section \ref{results}. Unless otherwise noted, the cosmological observables plotted in this section assume $100\Omega_{\rm b}h^2=2.22$, $\Omega_{\rm DM}h^2=0.1253$, $H_0=69.57$ km/s/Mpc, $10^9A_{\rm s}=2.21$, $n_{\rm s}=0.969$, and $\tau=0.092$. In the following, we often compare our PIDM observables with those from a standard $\Lambda$CDM cosmology with a ``equivalent'' number of neutrinos (denoted $\Lambda$CDM$+\nu$) to ensure an identical cosmological expansion history. This equivalent number of neutrinos is given by $N_{\rm eff,equiv}=3.046+(8/7)(11/4)^{4/3}\xi^4$.

\subsection{Galaxy Clustering}
The matter power spectrum describing the clustering of matter in the Universe depend on a variety of cosmological parameters, and for this reason it has been used (together with its Fourier transform, the correlation function) to set constraints on, among others, dark energy parameters~\cite{Samushia:2013vn}, models of gravity~\cite{Raccanelli:2012fk}, neutrino mass~\cite{Putter:2012kx, Zhao:2012ly}, the growth of structures~\cite{Samushia:2011uq, Reid:2012zr}, and non-Gaussianity~\cite{Ross:2012ys3}. Since PIDM models can generally have a large impact on the clustering of matter in the Universe, we expect that recent measurements of the galaxy power spectrum and correlation function \cite{Anderson:2012ve, Samushia:2013vn, Reid:2012zr, Sanchez:2013qf} can provide useful limits on scales where non-linearities can be neglected\footnote{Since PIDM models generally predict a different shape and amplitude for the small-scale power spectrum as compared to a standard $\Lambda$CDM model, one cannot use tools such as Halofit \cite{Smith:2002dz} to model non-linearities.}.
\begin{figure}[t!]
\begin{centering}
\includegraphics[width=0.5\textwidth]{PIDM_correlation_function.pdf}
\includegraphics[width=0.5\textwidth]{PIDM_correlation_function_f_ADM.pdf}
\caption{Angle averaged galaxy correlation function $\tilde{\xi}_0(r)$ for different PIDM models. In the upper panel, we take $f_{\rm int}=5\%$, $\xi=0.5$ and vary $\Sigma_{\rm DAO}$ and $\alpha_D$. In the lower panel, we fix $\Sigma_{\rm DAO}=10^{-3}$, $\alpha_D=0.01$ and $\xi=0.5$, but let the fraction of interacting DM vary. We set the galaxy bias to $b=2.2$ and the dilation scale to $\alpha=1.016$. We compare theoretical predictions with BOSS-DR9 measurements from Ref.~\cite{Anderson:2012ve}, and we also show a standard $\Lambda$CDM model with an equivalent number of effective neutrinos. In this work, we focus uniquely on linear scales, which lie to the right of the dashed vertical line on the plot.}
\label{fig:correlation}
\end{centering}
\end{figure}

Similar to the baryon acoustic oscillation signature in the galaxy correlation function of the standard $\Lambda$CDM model, the DAO feature appears as a local enhancement of the correlation function at $r_{\rm DAO}$.  While the location of the DAO feature depends mostly on $\xi=(T_D/T_{\rm CMB})|_{z=0}$, $f_{\rm int}$, and $\Sigma_{\rm DAO}$ (it also depends somewhat on $\alpha_D$ through $\epsilon_D$, see Eq.~\ref{eq:S_D_w_epsilon}), the shape on the DAO feature does depend on the ``microphysical'' dark parameters such as $\alpha_D$, $m_D$, and $B_D$, with the feature being generally sharper for small values of $\alpha_D$. Moreover, since fluctuations on scales smaller than the DAO scale are suppressed by acoustic damping, the correlation function (and the matter power spectrum) will generally be damped on small scales as compared to a standard $\Lambda$CDM cosmology without a DAO feature. If it affects the linear cosmological scales, we expect this damping to play a major role in our constraints on PIDM.

We illustrate in Fig.~\ref{fig:correlation} the predicted galaxy correlation function for different PIDM models. We plot the galaxy linear correlation function, that is computed from the linear matter correlation as
\begin{equation}
\tilde{\xi}_g(r) =b^2 \tilde{\xi}_m(\alpha r) \, ,
\end{equation}
where $b$ is the galaxy bias and $\alpha$ is the scale dilation parameter compensating for the difference between the fiducial cosmology used to compute the correlation function from the data and the actual cosmology (see~\cite{Xu:2012uq} for more details). The matter correlation function is computed from the linear matter power spectrum, $P_m(k)$, via the relation
\begin{equation}
\tilde{\xi}_m (r) = \int \frac{k^2}{2 \pi^2} P_m(k) \, j_0 (kr)  \, dk.
\end{equation}

In this work, we focus exclusively on the linear cosmological scales (corresponding to comoving $k\leq0.12 h/{\rm Mpc}$, see Section~\ref{sec:pk_data} for more details). Nevertheless, we also plot in Fig.~\ref{fig:correlation} the predictions for smaller scales to highlight the considerable damping at small scales for PIDM models. This shows how important it would be, in order to further reduce the allowed parameter space, to be able to model the quasi-linear regime.

In the upper panel of Fig.~\ref{fig:correlation}, we vary $\Sigma_{\rm DAO}$ and $\alpha_D$ for a fixed fraction of interacting DM $f_{\rm int}=5\%$. In all cases considered the DAO feature is clearly visible, providing a characteristic signature for these models. We observe that for the majority of the models shown ($\Sigma_{\rm DAO}> 10^{-4}$), even such a small fraction of interacting DM is in tension with measurements of the galaxy correlation function from the BOSS survey. In the lower panel, we fix $\Sigma_{\rm DAO}= 10^{-3}$ and instead vary the fraction of interacting DM between $2\%$ and $20\%$. We observe the scaling of the DAO scale with $f_{\rm int}$, $r_{\rm DAO}\propto1/\sqrt{f_{\rm int}}$, and also that a $\sim2\%$ fraction of strongly interacting DM seems to be compatible with current data. As we discuss in section \ref{results}, these qualitative observations will turn out to be supported by quantitative analyses.

For quantitative statistical analyses, it is usually computationally easier to consider the matter power spectrum directly. The signatures of interacting DM on the matter power spectrum has been extensively studied in Ref.~\cite{CyrRacine:2012fz} and we only review them briefly here. First, the presence of the DAO scale in PIDM models generally appears as extra oscillations in the matter power spectrum on scales with $k> k_{\rm DAO}\sim\pi/r_{\rm DAO}$ . Second, just as the correlation function is suppressed on small-scales due to acoustic damping in the dark plasma, the matter power spectrum displays less power at large wave numbers as compared with an equivalent $\Lambda$CDM model.

In Fig.~\ref{fig:matterpower}, we show the linear galaxy power spectrum for different PIDM models, along with the measured power spectrum from the BOSS-CMASS data~\cite{Anderson:2012ve}. In Section~\ref{sec:pk_data}, we explain how we convert theoretical PIDM matter power spectra to the shown galaxy power spectra, and give more details on the measurement of the BOSS-CMASS power spectrum and the computation of its errors. The upper panel of Fig.~\ref{fig:matterpower} displays how the power spectrum varies as $\Sigma_{\rm DAO}$ changes for the case of only $5\%$ of interacting DM. The lower panel illustrates the variations in the power spectrum as $f_{\rm int}$ changes from $2\%$ to $20\%$ for a fixed $\Sigma_{\rm DAO}=10^{-3}$ (these are the same models as those plotted in the lower panel of Fig.~\ref{fig:correlation}). The most obvious signature of PIDM in these plots is the damping of small-scale power. The actual acoustic oscillations are only clearly visible for models with $f_{\rm int}\gtrsim10\%$, indicating that dark oscillations are probably better illustrated through the DAO scale in the correlation function for models with small interacting DM fraction (see Fig.~\ref{fig:correlation}).

For the purpose of this work we limit our analysis to linear scales and avoid modeling the small and very large scales, where the galaxy clustering needs to include corrections due to non-linearities~\cite{Taruya:2009bs, Taruya:2010ij, Kwan:2011fu, Jennings:2012fv, Torre:2012dz} and large-scale effects~\cite{Raccanelli:2010fk, Yoo:2010cr, Challinor:2011nx, Bonvin:2011oq, Bertacca:2012hc, Yoo:2012dq, Yoo:2012bh, Montanari:2012kl, Raccanelli:2013tg}, respectively.
Our constraints on PIDM models using measurements of the galaxy power spectrum will be presented in section \ref{results}.
\begin{figure}[t!]
\begin{centering}
\includegraphics[width=0.5\textwidth]{PIDM_BOSS_5.pdf}
\includegraphics[width=0.5\textwidth]{PIDM_BOSS_f_ADM.pdf}
\caption{Linear galaxy power spectra for different PIDM models. In the upper panel, we fix $f_{\rm int}=5\%$, $\xi=0.5$ and vary $\Sigma_{\rm DAO}$. The lower panel uses $\Sigma_{\rm DAO}=10^{-3}$ and $\xi=0.5$ but let the fraction of interacting DM vary. To compare with galaxy power spectrum from the CMASS catalogue, we have convolved our linear matter power spectra with the BOSS window function and multiplied the results by a scale-independent galaxy bias $b=2.01$ (see section \ref{sec:pk_data} for more details).  For comparison, we also show a standard $\Lambda$CDM model with an equivalent number of effective neutrinos. In this work, we focus uniquely on linear scales, which lie to the left of the dashed vertical line on the plot.}
\label{fig:matterpower}
\end{centering}
\end{figure}
\subsection{Cosmic Microwave Background}\label{section:CMB}
The CMB probes cosmological fluctuations 380,000 years after the big bang. At that epoch, DM accounts for about $65\%$ of the energy budget of the Universe, hence making the CMB a particularly good probe of nonstandard DM physics. The PIDM scenario affects the CMB in three different ways. First, the presence of extra DR mimics the presence of extra neutrino species and affects the expansion history of the Universe, possibly modifying the epoch of matter-radiation equality, the CMB Silk damping tail, and the early Integrated Sachs-Wolfe effect. However, unlike standard free-streaming neutrinos, the DR forms a tightly-coupled fluid at early times, leading to distinct signatures on CMB fluctuations (see e.g. Ref.~\cite{Cyr-Racine:2013jua}). Second, the DR pressure prohibits the growth of interacting DM fluctuations on length scales entering the causal horizon before the epoch of DM kinematic decoupling. This weakens the depth of gravitational potential fluctuations on these scales, hence affecting the source term of CMB temperature fluctuations. Finally, as discussed in the previous subsection, the modified matter clustering in the Universe due to nonstandard DM properties will affect the lensing of the CMB as it travel from the last-scattering surface to us. We briefly review these signatures below but refer the reader to Ref.~\cite{CyrRacine:2012fz} for more detail. CMB lensing will be covered in the next subsection.

The impact of extra free-streaming relativistic species on the CMB has been well-studied in the literature \cite{Bashinsky:2003tk,Hou:2011ec}. In the following, we focus on the CMB signatures that arise due the modified evolution of cosmological fluctuations. These are sensitive to the actual physical properties of the relativistic species, hence allowing one to discriminate between, say, tightly-coupled DR or extra free-streaming neutrinos. In the standard $\Lambda$CDM universe, free-streaming neutrinos can establish gravitational potential perturbations beyond the sound horizon of the photon-baryon plasma, leading to a suppression of photon fluctuations as well as inducing a shift to their phase shorty after horizon entry. As shown in \cite{Bashinsky:2003tk}, these suppressions and phase shifts are directly proportional to the free-streaming fraction of the total amount of radiation in the Universe. In the PIDM scenario, this fraction is altered by the presence of a DR component that is tightly-coupled to the interacting DM at early times. This results in the CMB displaying a higher amplitude and having its peak structure shifted toward smaller scales (higher $l$) for multipoles that receive contributions from scales entering the horizon before the epoch of dark kinematic decoupling. On the other hand, CMB multipoles receiving most of their contributions from scales entering the horizon well after dark kinetic decoupling should be essentially indistinguishable from a $\Lambda$CDM model having an equivalent effective number of free-streaming species.

In addition to horizon-entry effects, the delayed onset of DR free-streaming as compared to neutrinos generally results in a modified evolution of sub-horizon shear stress perturbations. The modified anisotropic stress will in turn affect the amplitude of the photon quadrupole moment ($F_{\gamma2}$), leading to changes in the CMB temperature and polarization spectra. However, this effect can only be important if the impact on the photon quadrupole is significant near the peak of the CMB visibility function. We thus expect it to be largest for PIDM models with a kinematic decoupling epoch near or after the time of CMB last scattering. Indeed, in theses scenarios, the sudden growth of anisotropic stress due to the onset of DR free-streaming will modify the photon quadrupole just as the primary anisotropies are imprinted on the CMB. Similarly, the relative absence (as compared with an equivalent free-streaming neutrino model) of anisotropic stress near the peak of the visibility function in models with late kinetic decoupling will also result in a modified photon quadrupole, hence leading to different CMB anisotropies. The magnitude and sign of this effect depends on the relative phase between the photon quadrupole and the shear stress perturbation. Since the CMB polarization signal is very sensitive to the photon quadrupole moment near the epoch of last scattering, we expect this effect to be most important for the  polarization power spectrum.
\begin{figure}[t!]
\begin{centering}
\includegraphics[width=0.5\textwidth]{psi.pdf}
\includegraphics[width=0.5\textwidth]{delta_gam+psi.pdf}
\caption{{\it Upper panel:} Time-evolution of the gravitational potential $\psi$ for different values of $\Sigma_{\rm DAO}$. Here, $\psi$ is given in units of $2\zeta/3$, where $\zeta$ is the curvature perturbation on constant density hypersurfaces, which is conserved on super-horizon scale for pure adiabatic fluctuations. Before horizon entry, we have $\psi\rightarrow (2\zeta/3)(1+\frac{4}{15}R_{\rm f-s})^{-1}$, where $R_{\rm f-s}$ is the free-streaming fraction of the total radiation content of the Universe. {\it Lower panel:} Monopole source term for the CMB temperature anisotropies.  For both panels, we have taken $k = 0.15$ Mpc$^{-1}$, $\xi=0.5$, and $f_{\rm int}=100\%$. For comparison, we also show a standard $\Lambda$CDM model with an equivalent number of effective neutrinos. }
\label{fig:delta_gam+psi}
\end{centering}
\end{figure}

Before interacting DM kinematically decouples from the DR, fluctuations in the former cannot grow, leading to substantially shallower gravitational potentials in the matter-dominated era on sub-horizon scales (see upper panel of Fig.~\ref{fig:delta_gam+psi}). Physically, since perturbations in the photon-baryon plasma essentially obey a harmonic oscillator equation driven by the force of gravity, this amounts to a severe damping of the driving term\footnote{Mathematically, this can be understood as a damping of the particular solution to the differential equation governing perturbations in the photon-baryon plasma, while leaving unchanged its homogeneous solution.}. This has for consequences of weakening the gravitationally-induced compression phase of the acoustic oscillations (corresponding to odd $C_l^{\rm TT}$ peaks) while strengthening the pressure-supported expansion phase (corresponding to even $C_l^{\rm TT}$ peaks) of the oscillations (see Ref.~\cite{Foot:2012ai} for a similar effect in a different context). In lower panel of Fig.~\ref{fig:delta_gam+psi}, we illustrate the time-evolution of the monopole source term of CMB temperature fluctuations $(\delta_{\gamma}+\psi)^2$ for different values of $\Sigma_{\rm DAO}$. The enhancement of the expansion peaks and the suppression of the compression peaks are clearly visible for models with $\Sigma_{\rm DAO}\gtrsim10^{-3}$. The effect is also clearly visible at $l\gtrsim600$ in the $C_l^{\rm TT}$ spectra shown in Fig.~\ref{fig:cmb_nolens} for $\Sigma_{\rm DAO}=10^{-3}$.
\begin{figure}[]
\begin{centering}
\includegraphics[width=0.5\textwidth]{delta_ClsTT_ISW.pdf}
\includegraphics[width=0.5\textwidth]{delta_ClsTT_f_ADM_ISW.pdf}
\caption{Difference between PIDM models and a $\Lambda$CDM model with an equivalent number of neutrinos for the pure ISW contribution to $C_l^{\rm TT}$, $\Delta C_l^{\rm ISW} = C_l^{\rm ISW,PIDM}-C_l^{\rm ISW,CDM}$. In the upper panel, we take $\xi=0.5$, $f_{\rm int}=100\%$, and vary $\Sigma_{\rm DAO}$. In the lower panel, we take $\xi=0.5$, $\Sigma_{\rm DAO}=10^{-2}$, but let $f_{\rm int}$ vary.}
\label{fig:ISW}
\end{centering}
\end{figure}

For PIDM models that kinematically decouple after the last scattering of CMB photons, the continuous decay of the gravitational potential due to the DR pressure leads also to an enhanced early integrated Sachs-Wolfe effect (ISW). This is illustrated in Fig.~\ref{fig:ISW} where we show the difference in the pure early ISW contribution (that is, the ISW-ISW autocorrelation) to $C_l^{\rm TT}$ between a few PIDM models and a $\Lambda$CDM model with an equivalent number of neutrinos. We see that the temperature anisotropies of models with $\Sigma_{\rm DAO} > 10^{-3}$ can receive a large positive contribution from the early ISW effect for $200 \lesssim l\lesssim1200$. At larger multipoles however, the early ISW effect in PIDM models can actually be weaker than that of an equivalent $\Lambda$CDM model, although this effect is much less relevant since the ISW effect is subdominant on those scales.  We note that the actual impact on the temperature anisotropies can be larger than illustrated here since the ISW contribution adds in phase with the photon monopole source term. In the lower panel of Fig.~\ref{fig:ISW}, we observe that the impact of PIDM on the early ISW effect can be fairly large even if the fraction of interacting DM is subdominant.

Putting all these effects together, we obtain the unlensed temperature and polarization power spectra shown in Fig.~\ref{fig:cmb_nolens} where we have taken $f_{\rm int}=100\%$ to magnify the impact of PIDM on the CMB. First, we note that models with $\Sigma_{\rm DAO}\lesssim10^{-5}$ are essentially undistinguishable from a $\Lambda$CDM model with an equivalent number of neutrinos. As $\Sigma_{\rm DAO}$ is increased to $10^{-4}$, both temperature and polarization spectra begin to display the rise in amplitude and the phase shift associated with the DR being tightly-coupled at early times. If $\Sigma_{\rm DAO}$ is further increased to $10^{-3}$, the damping of the gravitational potential perturbations on small scales leads to the suppression of the odd peaks and the enhancement of the even peaks for the temperature anisotropies, which is clearly visible for $l\gtrsim600$. The early ISW also enhances the first acoustic peak of the temperature spectrum for this model. For the polarization spectrum, the impact of the DR shear stress on the photon quadrupole moment leads to a fairly complex variation, with some peak being enhanced, while some are suppressed as compared with an equivalent $\Lambda$CDM model.  Finally, as $\Sigma_{\rm DAO}$ is increased to $10^{-2}$, all the physical effects discussed above are present, with the early ISW leading to strong enhancement of the third temperature peak, which is almost sufficient to offset the suppression caused by the weak gravitational driving force. The early ISW also enhances the first temperature peak while the second peak is somewhat suppressed due the DR shear stress, which also affect the polarization spectrum in a non-trivial way.

In summary, as illustrated in Fig.~\ref{fig:cmb_nolens}, PIDM models can lead to rich CMB signatures that cannot be easily mimicked by varying other standard cosmological parameters. This makes the CMB an excellent probe of interacting DM and DR physics. As we will see in Section \ref{results}, the latest CMB data does indeed provide strong constraints on PIDM models.
\begin{figure}[t!]
\begin{centering}
\includegraphics[width=0.49\textwidth]{ClsTT_GamDAO_nolens.pdf}
\includegraphics[width=0.49\textwidth]{CMB_EE_GamDAO_nolens.pdf}
\caption{CMB unlensed temperature (upper panel) and E polarization (lower panel) power spectra for four different PIDM models with $f_{\rm int}=100\%$. We have taken $\xi=0.5$. For comparison, we also show a standard $\Lambda$CDM model with an equivalent number of effective neutrinos. }
\label{fig:cmb_nolens}
\end{centering}
\end{figure}
\subsection{CMB Lensing}
As the CMB photons free-stream from the last-scattering surface to us, they encounter large DM structures which can deflect their path and rotate their polarization state. This CMB lensing (see \cite{LewChal06} for a review) by foreground matter structures has now been detected at high statistical significance ($\sim25\sigma$, \cite{planck17_2013}) and can be used to study the distribution of matter throughout the Universe. Since PIDM models generally predict a modified matter distribution as compared to a pure CDM model, CMB lensing can by itself provide useful constraints on interacting DM scenarios.

The gravitational deflection potential $\phi$, of which the gradient gives the lensing displacement vector on the sky,
is related to the gravitational potential perturbation $\psi$ projected along the line of sight in the ${\bf \hat{n}}$ direction, via
\begin{equation}
\phi({\bf \hat{n}}) = -2\int_0^{\chi_*}d\chi\,\psi(\chi {\bf \hat{n}};\eta_0-\chi)\frac{\chi_*-\chi}{\chi\chi_*},
\end{equation}
where $\chi_*$ is the comoving distance to the last scattering surface and $\eta_0$ is the comoving size of the causal horizon today. The lensing potential power spectrum can be written as
\begin{equation}
C_l^{\phi\phi} = 16\pi \int\frac{dk}{k}P_{\mathcal{R}}(k) |\Delta_{\psi}(k)|^2,
\end{equation}
where
\begin{equation}
\Delta_{\psi}(k)= \int_0^{\chi_*}d\chi T_{\psi}(k;\eta_0-\chi))j_l(k\chi)\frac{\chi_*-\chi}{\chi\chi_*},
\end{equation}
and where $P_{\mathcal{R}}(k)$ is the primordial spectrum of comoving curvature fluctuations. The transfer function $T_{\psi}(k,\eta)$ is defined by $\psi(k,\eta) = T_{\psi}(k,\eta)\mathcal{R}(k)$, where $\mathcal{R}(k)$ stands for the comoving curvature fluctuation.

We show in Fig.~\ref{fig:lensing} the CMB lensing power spectrum for different PIDM models. In the upper panel, we display the spectra for increasing values of $\Sigma_{\rm DAO}$. It should be clear from this plot that the most extreme models with $\Sigma_{\rm DAO}\gtrsim10^{-3}$ are ruled out by current data if interacting DM forms the totality of the DM. In the lower panel of Fig.~\ref{fig:lensing}, we fix  $\Sigma_{\rm DAO}=10^{-3}$ but instead vary the fraction of interacting DM. We observe that even a fraction as small as $5\%$ can have a sizable effect on the lensing power spectrum. This indicates that current and future CMB lensing measurements could potentially be very sensitive probes of nonstandard DM physics.
\begin{figure}[t!]
\begin{centering}
\includegraphics[width=0.5\textwidth]{PIDM_lensing_100.pdf}
\includegraphics[width=0.5\textwidth]{PIDM_lensing_f_ADM.pdf}
\caption{CMB lensing power spectrum for different PIDM models. For both panels we use $\xi=0.5$. In the upper panel, we vary $\Sigma_{\rm DAO}$ while leaving $f_{\rm int}=100\%$ fixed. The model with $\Sigma_{\rm DAO}=10^{-5}$ is essentially undistinguishable from the $\Lambda$CDM+$\nu$ model. In the lower panel, we vary $f_{\rm int}$ but leave $\Sigma_{\rm DAO}=10^{-3}$ fixed. We show the eight band powers used in the Planck lensing likelihood.  For comparison, we also show a $\Lambda$CDM model with an equivalent number of neutrinos.}
\label{fig:lensing}
\end{centering}
\end{figure}

Lensing by foreground matter structure also distorts the CMB temperature and polarization power spectra presented in Fig.~\ref{fig:cmb_nolens} above. Essentially, lensing acts to smooth out the oscillatory structure of the spectra, filling in the troughs and damping the peaks. As we discussed above, since PIDM models generally predict different amount of lensing, the associated smoothing of the CMB spectra provides yet another handle (albeit correlated with other CMB signatures) to constrain interacting DM. We illustrate lensed CMB spectra in Figs.~\ref{fig:cmbEE} and \ref{fig:cmbTT} for increasing values of $\Sigma_{\rm DAO}$ and for $f_{\rm int}=1$. Besides the PIDM signatures discussed in section \ref{section:CMB}, we observe that the TT and EE spectra display sharper peaks and troughs in the damping tail as $\Sigma_{\rm DAO}$ is increased, which is in line with our expectations that these models should be \emph{less} affected by gravitational lensing. We also note that the lensing signatures can obscure some of the effects discussed in section \ref{section:CMB}, especially the enhancement of the even acoustic peaks in the damping tail of the temperature spectrum.

Taken as a whole, it is clear that the CMB and its lensing by foreground matter structures provide an exquisite probe of DM physics and of its possible interaction with new relativistic species. Having described the CMB signatures predicted by PIDM models and the physics behind them, we now turn our attention to the data and what they can tell us about the physics of DM.
\begin{figure}[b!]
\begin{centering}
\includegraphics[width=0.5\textwidth]{Cls_EE_GamDAO.pdf}
\caption{CMB E polarization power spectrum for four different PIDM models with $f_{\rm int}=100\%$. and $\xi=0.5$.  For comparison, we also show a standard $\Lambda$CDM model with an equivalent number of effective neutrinos. }
\label{fig:cmbEE}
\end{centering}
\end{figure}

\begin{figure*}[]
\begin{centering}
\includegraphics[width=\textwidth]{ClsTT_GamDAO.pdf}
\caption{CMB temperature power spectrum for four different PIDM models with $f_{\rm int}=100\%$. and $\xi=0.5$. We also show the Planck \cite{planckXVI}, SPT \cite{Story:2012wx}  and ACT  \cite{Dunkley:2013vu,Das:2013zf} band powers.  For comparison, we also show a standard $\Lambda$CDM model with an equivalent number of effective neutrinos.  }
\label{fig:cmbTT}
\end{centering}
\end{figure*}

\section{Data}\label{data}

\subsection{Cosmic Microwave Background and Lensing}
To constrain interacting DM, we use the CMB data from the Planck satellite \cite{planckXVI}. We utilize both the low-multipole and high-multipole temperature data, incorporating the required ``nuisance'' parameters describing foregrounds and instrumental effects, and also include the WMAP low-$l$ polarization data.  We refer to this dataset as ``Planck+WP''. We also incorporate the high-resolution temperature data from the South Pole Telescope (SPT) and the Atacama Cosmology Telescope (ACT). As in the Planck analysis, we only include the ACT $148\!\times\!148$ spectra for $l\!\geq\!1000$, the ACT $148\!\times\!218$ and $218\!\times\!218$ spectra for $l\!\geq\!1500$ \cite{Dunkley:2013vu,Das:2013zf}, and the SPT data described in \cite{Story:2012wx} for $l\!\geq\!2000$. We fully incorporate the nuisance parameters describing foregrounds and calibration uncertainties for both SPT and ACT. We collectively refer to this dataset as ``High-$l$''.

Our likelihood also makes use of the measurement of the CMB lensing potential power spectrum
by the Planck collaboration \cite{planck17_2013}. This data set
consists of eight bandpower estimates of the lensing power spectrum,
covering the multipole range $l = 40 - 400$, as shown in Fig.~\ref{fig:lensing}.
For a fixed realization of the lensing potential,
CMB lensing causes correlations between different Fourier modes of the CMB temperature field.
The lensing potential can therefore be reconstructed
by averaging over products of pairs of temperature modes (see e.g.~Ref.~\cite{okahu03}).
The measurement of the power spectrum of this reconstructed lensing potential
is thus a temperature {\it trispectrum} measurement and is non-zero because when averaged over realizations
of the lensing potential, the lensed temperature field is non-Gaussian. We refer to the Planck lensing dataset as ``Lens''.
\subsection{Baryon Acoustic Oscillation}

We also include in our analysis baryon acoustic oscillation (BAO) data from a reanalysis
of the Sloan Digital Sky Survey DR7 \cite{Padmanabhan:2012hf}, from the 6-degree Field
survey \cite{Beutler:2011hx}, and from the Baryon Oscillation Spectroscopic Survey \cite{Anderson:2012ve}.
To measure the BAO scale from the galaxy power spectrum (or correlation function),
typically a fit is performed to a template spectrum based on a $\Lambda$CDM fiducial
cosmology. The template is allowed to shift in the horizontal direction (corresponding
to a dilation of distances) and the best fit dilation parameter
gives a measurement of the ratio of the true BAO scale over the true distance to the galaxy sample, relative
to this ratio in the fiducial cosmology. To account for modeling systematics,
modifications of the broadband shape and amplitude of the power spectrum are allowed, and parametrized by nuisance parameters,
which are then marginalized over. This procedure is based on the assumption
that modifications to the template power spectrum are slowly varying, i.e.~that the deviation
of the true spectrum from the model does not have any sharp, or rapidly varying features.

While this assumption is very reasonable when the main systematics come from non-linear corrections, etc.,
to a $\Lambda$CDM power spectrum, one might worry that it breaks down in the presence of
DAO wiggles. For example, if the DAO scale is close to the BAO scale, the BAO measurement procedure may be biased or might even inadvertently pick out the DAO scale instead.
However, it turns out that most of the region of PIDM parameter space where the DAO scale is large and might cause such confusion is ruled out already by the CMB data alone, so that for the parameter space where the BAO measurement has any effect, the standard BAO measurement procedure is appropriate. Only for very small fractions of interacting DM can the DAO scale be close to the BAO scale and not be in tension with CMB data. For these cases, we confirm our results by using the full shape of the galaxy power spectrum which is free from the assumptions made in order to reconstruct the BAO feature. We can thus safely use the BAO measurements given in the literature as a simple set of low-redshift distance priors.

\subsection{BOSS Galaxy Power Spectrum}
\label{sec:pk_data}

We use the galaxy power spectrum from the Baryon Oscillation Spectroscopic
Survey (BOSS \cite{dawsonetal13}), which is a component of the Sloan Digital Sky Survey (SDSS \cite{yorketal00}),
specifically SDSS-III \cite{sdss3}.
The data set used here is the CMASS sample of luminous galaxies (see e.g.~\cite{whiteetal11,Anderson:2012ve}), released as part of
Data Release 9 \cite{DR9}.
This sample consists of $\sim 264,000$ galaxies in the redshift range $z = 0.43 - 0.7$ (effective redshift
$z_{\rm eff} = 0.57$), covers 3275 deg$^2$ of the sky, and has an effective volume
$V_{\rm eff} =  2.2$ Gpc$^3$.
We quantify galaxy clustering by the angle-averaged power spectrum of this sample,
which has also been used to obtain the strongest measurement to date of the baryon acoustic oscillation (BAO)
scale in \cite{Anderson:2012ve}, to constrain neutrino mass in \cite{Zhao:2012ly},
and to study primordial non-Gaussianity \cite{Ross:2012ys3}.

The power spectrum is measured using the Feldman-Kaiser-Peacock (FKP \cite{FKP}) estimator,
see \cite{reidetal10,Anderson:2012ve} for details.
To calculate constraints on the PIDM model, we use only the wave vector
range $k = 0.03 - 0.12 h/$Mpc. The upper limit $k_{\rm max} = 0.12 h/$Mpc serves to ensure
that the power spectrum can be modeled using linear perturbation theory (this will be discussed in more detail below).

While the lower limit $k_{\rm min} = 0.03 h/$Mpc minimizes its importance,
we include a template to subtract
any spurious clustering signal on large scales
due to systematics that have not been accounted for \cite{arossetal12,Ross:2012ys3},
\begin{equation}
\label{eq:sys}
P_{\rm meas}(k) = P_{\rm meas, w}(k) - S \, P_{\rm sys}(k).
\end{equation}
Here, $P_{\rm meas, w}(k)$ is the observed power spectrum described in \cite{Anderson:2012ve},
using the standard weights that take into account known systematics (including stellar density).
The template $P_{\rm sys}(k)$ is equal to the contribution to the galaxy power spectrum
due to the correlation of observed galaxy density with stellar density (the dominant known systematic on large scales),
if stellar density weights had not been included. It is argued in \cite{arossetal12,Ross:2012ys3} that
the $k$-dependence of the contribution from other systematics will be the same as that of $P_{\rm sys}(k)$,
so that, following these works, we model unknown systematics
by simply treating the amplitude $S$ in Eq.~(\ref{eq:sys}) as a free parameter and marginalizing over it. We restrict $S$ to be in the range from -1 to 1.

The likelihood for a given PIDM cosmology is obtained by comparing the measured power spectrum $P_{\rm meas}(k)$
to a model power spectrum, $P_{\rm model}(k)$, which we describe below.
To start, we model the true galaxy power spectrum by
\begin{equation}
\label{eq:Pg}
P_g(k) = b^2 P_m(k) + P_0.
\end{equation}
Here, $P_m(k)$ is the {\it linear} matter power spectrum, computed using a modified version
of CAMB (see Section \ref{CosmoPert}), and $b$ is a scale-independent, linear galaxy bias.
The first term on the right hand side of Eq.~(\ref{eq:Pg})
describes the galaxy power spectrum in the linear regime.
While our choice $k_{\rm max} = 0.12 h/$Mpc limits the importance of non-linear corrections,
we also include a nuisance parameter $P_0$ to model possible deviations from the linear description
due to scale-dependent galaxy bias, and/or due to imperfect shot-noise subtraction.
This simple galaxy bias model is motivated by the halo model \cite{seljak2000,seljak01,schulzwhite06,guziketal07}
and local bias \cite{scherrwein98,colesetal99,saitoetal09} approaches to galaxy clustering.

The power spectrum model also needs to take into account how the galaxy power spectrum is estimated from data.
We follow closely the approach of, e.g., \cite{Anderson:2012ve,Zhao:2012ly} here.
First of all, to estimate the power spectrum, cosmic distances were calculated assuming
a fixed fiducial cosmology. To take this into account, we
dilate the wave vector $k$ appearing in the theory power spectrum by
a factor $\alpha = D_{V, {\rm fid}}/D_V$ (see e.g.~\cite{Tegetal2006}),
$k \to \alpha k$,
where $D_{V, {\rm fid}}$ is the volume weighted distance measure to $z_{\rm eff} = 0.57$
in the fiducial cosmology, and $D_V$ is the same quantity in the cosmology in which the
galaxy power spectum is modeled.
The amplitude of the power spectrum should also be rescaled by a factor $\alpha^3$, but we
absorb this shift into the galaxy bias parameters.
Secondly, we account for the effect of the survey geometry
by convolving the model power spectrum with the Fourier transform of the survey window
function \cite{percetal07,arossetal12}.
Schematically, we thus have the following model power spectrum,
\begin{equation}
P_{\rm model}(k) = W * P_g(\alpha k),
\end{equation}
where $P_g$ was defined in Eq.~(\ref{eq:Pg}).

To calculate the likelihood, we assume the power spectrum estimator follows a Gaussian distribution
and use a covariance matrix based on 600 CMASS mock catalogs \cite{maneraetal13}.
Since the galaxy bias parameter $b^2$ corresponds to an overall
scaling of the model power spectrum (after a trivial redefinition $P_0 \to P_0/b^2$),
we marginalize over it analytically, thus reducing by one the number of free parameters to sample. We shall refer to this dataset as ``DR9''.

\section{Results}\label{results}
We determine constraints on PIDM models using the publicly available Markov Chain Monte Carlo code \texttt{CosmoMC} \cite{Lewis:2002ah}. In all our chains, we let the six $\Lambda$CDM cosmological parameter vary ($\Omega_{\rm b}h^2$, $\Omega_{\rm DM}h^2$, $\theta_{\rm MC}$, $\tau$, $\ln(10^{10} A_{\rm s})$, and $n_{\rm s}$) in addition to some combinations of the dark parameters specified below. We use uniform priors on the standard cosmological parameters and nuisance parameters similar to those described in \cite{planckXVI}. We assume a flat Universe with massless neutrinos and fix $N_{\rm eff}=3.046$. We compute the helium abundance consistently \cite{Simha:2008zj}, taking into account the baryon density and the extra contributions to the radiation density from DR and dark electrons if the latter are relativistic at big bang nucleosynthesis. To explore the $\Sigma_{\rm DAO}$ constraints and ensure maximum chain mobility, we vary $\log_{10}(\Sigma_{\rm DAO})$ instead of $\Sigma_{\rm DAO}$ itself. The ranges of the uniform priors we use for the dark parameters are listed in Table \ref{table_prior}. To ensure chain convergence, we run 8 independent chains for each dataset combinations and make sure that the Gelman and Rubin criterion is $R-1 \leq 0.02$.
\begin{table}[t!]
\begin{center}
\begin{tabular}{|c|c|}
\hline
Parameters & Prior Range \\
\hline
\hline
$\log_{10}(\Sigma_{\rm DAO})$ & [-6,1] \\
$f_{\rm int}$ & [0,1] \\
$\xi$ & [$10^{-3}$,1] \\
\hline
\end{tabular}
\caption{\label{table_prior} Prior range of the dark parameters used in our analysis.\vspace{-0.7cm}}
\end{center}
\end{table}

Before looking at the quantitative results, it is important to note that PIDM models reduce to a standard CDM model in a few different limits. First, as $\xi$ is driven to small values, interacting DM looses the pressure support that prevents it from clustering and basically behaves like CDM. Second, as $\Sigma_{\rm DAO}$ is reduced, the epoch of DM kinematic decoupling is pushed toward earlier times, hence restricting the impact of interacting DM to small non-linear scales, unobservables with the CMB and (linear) galaxy clustering. Finally, the signatures of PIDM obviously become less and less important as $f_{\rm int}$ is reduced toward zero. The caveat here is that a smaller interacting DM fraction also implies a larger DAO scale, hence potentially bringing the small signatures of PIDM to potentially observable scales. Of course, as $f_{\rm int}\rightarrow0$, we recover the standard CDM observables. Since there are many ways to effectively ``hide'' the interacting DM component, it usually makes more sense to consider the cosmological constraints on specific slices of the DM parameter space. We discuss the most interesting slices in the next few subsections.

In the following subsections, we vary different subsets of the dark parameters to study their joint probability density function.

\subsection{ $\Sigma_{\rm DAO}$ + $\xi$ + $f_{\rm int}$}
\begin{figure}[b]
\begin{centering}
\includegraphics[width=0.5\textwidth]{loggamDAO_xi_2D_mpk_marg_f_ADM.pdf}
\caption{Marginalized constraints on $\xi$ and $\Sigma_{\rm DAO}$ for three combination of datasets. Here, the fraction of interacting DM has been marginalized over. We display the $68\%$ and $95\%$ confidence regions.}
\label{GamDAOvsxi_1}
\end{centering}
\end{figure}

For these constraints, we let the trio of dark parameters $\{\Sigma_{\rm DAO},f_{\rm int},\xi\}$ freely vary within the prior ranges listed in Table \ref{table_prior}. Here, we focus on strongly-coupled models and fix $\alpha_D=0.05$ and $m_D=10$ GeV. Constraints for different values of $\alpha_D$ and $m_D$ can be obtained by appropriately rescaling $\Sigma_{\rm DAO}$, as long as $\alpha_D > 0.025$. Constraints for smaller values of the dark fine-structure constant will be discussed in Section \ref{sec:var_alphaD}. In Fig.~\ref{GamDAOvsxi_1}, we illustrate the 2D marginalized posterior in the $\xi-\log_{10}(\Sigma_{\rm DAO})$ plane for three different combinations of datasets.  We note that the two combinations of datasets that make use of galaxy clustering data (here identified as BAO and DR9) yield very similar limits, indicating that most of the constraining power from galaxy power spectrum data comes from the BAO feature. The red contours display the constraints if no galaxy clustering data is included in the analysis. We observe that for $\Sigma_{\rm DAO} > 10^{-3}$ most of the constraining power comes from the CMB data alone (including lensing), while the large-scale structure data act to strengthen the constraints for lower values of $\Sigma_{\rm DAO}$.

In all cases, we observe a sharp regime change around $\Sigma_{\rm DAO}=10^{-4.5}$. For $\Sigma_{\rm DAO}\lesssim10^{-4.5}$, the bound on $\xi$ reads $\xi \lesssim 0.6$ ($95\%$ C.L., Planck+WP+High-$l$+BAO+Lens) which is equivalent to current cosmological constraints on $N_{\rm eff}$ for the datasets shown. This indicates that in this region of parameter space, PIDM essentially behaves like CDM and DR is indistinguishable from extra species of free-streaming neutrinos, in agreement with our discussion from section \ref{observables}. On the other hand, the allowed range of $\xi$ values rapidly shrinks as $\Sigma_{\rm DAO}$ is increased beyond $10^{-4.5}$. The strong constraints on $\xi$ in this region of parameter space have very different origins than the standard constraints on extra relativistic species. Indeed, the suppressed DM fluctuations and the absence of early DR free-streaming for these values of $\Sigma_{\rm DAO}$ have a significant impact on the CMB and the matter power spectrum  which forces $\xi$ to take small values in order to hide these effects. The constraints shown in Fig.~\ref{GamDAOvsxi_1} implies that we must have $\xi<0.2$ for $\Sigma_{\rm DAO} > 10^{-3}$. To give a sense of scale, we note that $\xi=0.2$ corresponds, according to Eq.~\ref{delta_n_eff}, to $\Delta N_{\rm eff}\simeq0.007$, indicating that even such a small amount of DR can dramatically affect the evolution of cosmological perturbations if it couples strongly enough to DM. Since the PIDM models lying inside the $95\%$ confidence region are, for all practical purposes, undistinguishable from a $\Lambda$CDM universe, the fraction of interacting DM is largely unconstrained in the allowed region. Constraints on models lying on the edge of the $95\%$ confidence region will however depend on the value of $f_{\rm int}$ and we explore these limits in the next section.

\subsection{$\Sigma_{\rm DAO}$ + $\xi$}
\begin{figure}[b!]
\begin{centering}
\includegraphics[width=0.5\textwidth]{loggamDAO_xi_2D_f_pidm.pdf}
\caption{Marginalized constraints on $\xi$ and $\Sigma_{\rm DAO}$ for three fixed values of $f_{\rm int}$. We display the $68\%$ and $95\%$ confidence regions for the dataset ``Planck+WP+High-$l$+BAO+Lens''.}
\label{GamDAOvsxi_2}
\end{centering}
\end{figure}
In this section, we keep the fraction of interacting DM fixed while letting $\xi$ and $\log_{10}(\Sigma_{\rm DAO})$ vary, allowing to determine how the constraint contours change as a function of the interacting DM fraction. As before, we fix $\alpha_D=0.05$ and $m_D=10$ GeV. We show in Fig.~\ref{GamDAOvsxi_2} the marginalized constraints for three different values of $f_{\rm int}$, using the dataset Planck+WP+High-$l$+BAO+Lens. While we observe the constraints becoming progressively weaker as $f_{\rm int}$ is reduced, the difference between the $f_{\rm int}=50\%$ and $f_{\rm int}=5\%$ limits is surprisingly modest. This indicates that our constraints are robust to changes in the interaction DM fraction (for $f_{\rm int}\gtrsim5\%$). It also show that it matters little if $f_{\rm int}=5\%$, $50\%$, or $100\%$ in the ruled out regions: there, PIDM affects the cosmological observables in a way that is incompatible with the current data and lowering the fraction of interacting DM only slowly improves the fit. This is in agreement with our discussion of sections \ref{CosmoPert} and \ref{observables} where we showed that shrinking the fraction of interacting DM does reduce the impact on the cosmological observables but at the price of \emph{increasing} the DAO scale and bringing it to linear observable scales.

\subsection{Varying the Dark Fine-structure Constant}\label{sec:var_alphaD}
\begin{figure}[t!]
\begin{centering}
\includegraphics[width=0.5\textwidth]{loggamDAO_xi_2D_alpha.pdf}
\caption{Marginalized constraints on $\xi$ and $\Sigma_{\rm DAO}$ for three fixed values of $\alpha_D$. Here, we have fixed $f_{\rm int}=10\%$ and $m_D=10$ GeV. We display the $68\%$ and $95\%$ confidence regions for the dataset ``Planck+WP+High-$l$+BAO+Lens''.}
\label{GamDAOvsxi_3}
\end{centering}
\end{figure}
In this section, we study the effect of varying $\alpha_D$ on the cosmological limits on the $\xi$ and $\Sigma_{\rm DAO}$ parameters. We display the constraints on these two dark parameters in Fig.~\ref{GamDAOvsxi_3} for three values of the dark fine-structure constant. Here, we fix $f_{\rm int}=10\%$ which yields constraints representative of a broad range of interacting DM fraction (see previous section). For $\Sigma_{\rm DAO}> 10^{-2.5}$, we observe that the constraint on $\xi$ is largely independent of $\alpha_D$, indicating that our limits are robust to changes in the dark sector microphysics in that region of parameter space. At smaller values of $\Sigma_{\rm DAO}$, the constraints become \emph{stronger} as $\alpha_D$ is reduced. This somewhat counterintuitive result is a consequence of the definition of $\Sigma_{\rm DAO}$: at fixed $\Sigma_{\rm DAO}$, lowering $\alpha_D$ leads to a smaller values of the atomic binding energy, hence bringing dark recombination and kinematic decoupling closer to the last scattering surface of CMB photons and leading to a larger effects on the cosmological observables. From Fig.~\ref{GamDAOvsxi_3}, we observe that the main impact of varying $\alpha_D$ is to modify the shape of the $\xi$ constraint in the transition region delimiting the parameter space where the limit is similar to the $N_{\rm eff}$ bound (small $\Sigma_{\rm DAO}$) and where it is dominated by PIDM effects (large $\Sigma_{\rm DAO}$).

\subsection{$\Sigma_{\rm DAO}$ + $f_{\rm int}$}\label{DAO+f}
\begin{figure}[]
\begin{centering}
\includegraphics[width=0.5\textwidth]{loggamDAO_f_ADM_2D.pdf}
\caption{Marginalized constraints on $f_{\rm int}$ and $\Sigma_{\rm DAO}$ for three values of $\xi$. Here, we have fixed $\alpha_D=0.05$ and $m_D=10$ GeV. We display the $68\%$ and $95\%$ confidence regions for the dataset ``Planck+WP+High-$l$+BAO+Lens''.}
\label{GamDAOvsxi_4}
\end{centering}
\end{figure}

In this section, we explore the constraints in the $f_{\rm int}-\log_{10}(\Sigma_{\rm DAO})$ plane for fixed values of $\xi$. This analysis allows to determine \emph{how much} interacting DM is allowed by the current data given an interaction strength and a certain density of DR. We displays the constraints in Fig.~\ref{GamDAOvsxi_4} for three values of $\xi$. As in the bounds presented in the previous subsections, there is a sharp transition around $\Sigma_{\rm DAO}\sim10^{-3}-10^{-4}$ for which the fraction of interacting DM goes from being largely unconstrained to being tightly bounded with$f_{\rm int}\lesssim5\%$. The exact constraint depends somewhat on the value of $\Sigma_{\rm DAO}$, with larger values of the latter leading to smaller allowed interacting DM fraction, as we intuitively expect. For $\xi\sim0.4-0.5$, there is an intriguing high-probability region near $\Sigma_{\rm DAO}\sim10^{-2.5}$ for small values of $f_{\rm int}$. To obtain a better picture of what is going on there and to determine a rigorous constraint on the allowed fraction of very strongly interaction DM (that is, interacting with strength equal or stronger than regular baryons), we fix $\Sigma_{\rm DAO}=10^{-2.5}$ and let $f_{\rm int}$ vary freely for the three values of $\xi$ shown in Fig.~\ref{GamDAOvsxi_4}. The resulting marginalized posteriors are shown in Fig.~\ref{allowed_frac} for both the ``Planck+WP+High-$l$+BAO+Lens'' and ``Planck+WP+High-$l$+DR9+Lens'' datasets.

We indeed observe that both datasets display a mild preference for a nonzero fraction of interacting DM for $\xi=0.5$, while this preference largely goes away as $\xi$ is decreased. It is interesting that this preference becomes stronger when the full shape of the BOSS DR9 galaxy power spectrum is included in the analysis. The penchant for $f_{\rm int}=(2.0\pm1.6)\%$ ($95\%$ C.L.) when $\xi=0.5$ and $\Sigma_{\rm DAO}=10^{-2.5}$ can be understood by looking at the DAO scale for these models. Indeed, in this corner of parameter space the DAO scale lies very close to the standard BAO scale (see Fig.~\ref{fig:rDAO}) leading to substantial overlap and interaction between the DAO and BAO features. The results is effectively a modified amplitude and shape of the BAO feature which, incidentally, improve the fit to the data. However, given the uncertainties in modeling and reconstructing the BAO bump in the galaxy correlation function, it is possible that this preference for a non-vanishing interacting DM fraction is purely coincidental. It is nonetheless intriguing that a $\sim 2\%$ DM fraction interacting very strongly with a DR bath at a temperature $T_D\sim0.5 T_{\rm CMB,0}$ provides an excellent fit to current cosmological data. We emphasize that this allowed PIDM model is particularly interesting since it departs significantly from a pure CDM scenario, in contrast with other allowed regions of the PIDM parameter space which have low values of $\Sigma_{\rm DAO}$ and/or $\xi$, and thus for which PIDM is for all practical purposes indistinguishable from CDM on linear cosmological scales.

We list in Table~\ref{table_f_ADM} the $95\%$ confidence limits for the three values of $\xi$ shown in Figs.~\ref{GamDAOvsxi_4} and \ref{allowed_frac}. While the exact numbers somewhat vary depending on the dataset considered, we observe that current data bound the deviation from a pure CDM scenario to be at most $\sim5\%$ of the overall cosmological DM density for a DM candidate that strongly interacts  with a cosmologically significant amount\footnote{We note that $\xi=0.3$ corresponds to $\Delta N_{\rm eff}\simeq0.036$ which is only marginally cosmologically significant.} of DR. This is a key result of our paper. Since it is based purely on the gravitational effects of interacting DM with the rest of the cosmological constituents, this result is very general and model-independent. It only relies on having a fraction of the DM interacting with a cosmologically significant amount of DR (which can be any type of relativistic particles) and assumes no particle coupling between the visible and dark sector.

\begin{figure}[t!]
\begin{centering}
\includegraphics[width=0.5\textwidth]{allowed_frac.pdf}
\includegraphics[width=0.5\textwidth]{allowed_frac_mpk.pdf}
\caption{Marginalized constraints on the fraction of very strongly interacting DM for three values of $\xi$. Here, we have fixed $\alpha_D=0.05$, $\Sigma_{\rm DAO}=10^{-2.5}$ and $m_D=10$ GeV. The top panel displays the constraints for the ``Planck+WP+High-$l$+BAO+Lens'' dataset, while the lower panel substitutes the BAO measurements for the BOSS DR9 CMASS galaxy power spectrum \cite{Anderson:2012ve}.}
\label{allowed_frac}
\end{centering}
\end{figure}
\begin{table}[]
\begin{center}
\begin{tabular}{|c|c|c|}
\hline
$\xi$ &CMB+DR9&CMB+BAO  \\
\hline
\hline
$0.3$ & $< 0.085$ & $< 0.049$ \\
\hline
$0.4$ & $< 0.018$ & $<0.035$\\
\hline
$0.5$ & $0.02\pm0.016$ & $< 0.038$\\
\hline
\end{tabular}
\caption{\label{table_f_ADM} Limits on $f_{\rm int}$ for the ``Planck+WP+High-$l$+DR9+Lens'' (middle column) and  ``Planck+WP+High-$l$+BAO+Lens'' (right column) dataset. Here, we have fixed $\alpha_D=0.05$, $\Sigma_{\rm DAO}=10^{-2.5}$ and $m_D=10$ GeV. We display the $95\%$ confidence limits. The severe constraint for $\xi=0.4$ arising when the shape of the BOSS DR9 galaxy power spectrum is taken into account is caused by the DAO scale being just below the BAO scale, leading to a gravitational damping of the BAO bump which is incompatible with data.
\vspace{-0.7cm}}
\end{center}
\end{table}
\subsection{Limits on the DAO Scale}

We can translate our constraints on the dark parameters to an upper limit on the size of the sound horizon of interacting DM when it kinematically decouples from the DR bath.  Fixing the value of the interacting DM fraction, we obtain the following limits on the comoving value $r_{\rm DAO}$
\begin{equation}
r_{\rm DAO} < 3.7 h^{-1} {\rm Mpc}\quad(f_{\rm int}=100\%),
\end{equation}
\begin{equation}
r_{\rm DAO} < 5.3 h^{-1} {\rm Mpc}\quad(f_{\rm int}=50\%),
\end{equation}
\begin{equation}
r_{\rm DAO} < 15.2 h^{-1} {\rm Mpc}\quad(f_{\rm int}=10\%),\\
\end{equation}
\begin{equation}
r_{\rm DAO} < 27.9 h^{-1} {\rm Mpc}\quad(f_{\rm int}=5\%),
\end{equation}
where we are giving the $95\%$ confidence limits for the ``Planck+WP+High-$l$+BAO+Lens'' dataset and where we used Eq.~\ref{r_DAO_theory} to compute the DAO scale. These constraints are valid for $\alpha_D>0.025$, but only become slightly stronger for lower values of the dark fine-structure constant. For $f_{\rm int}\geq5\%$, our constraints imply that the DAO scale must lie on relatively small scales where nonlinear effects can be important. Improving upon these constraints will therefore necessitate a prescription to model small-scale nonlinearities in PIDM models. One might worry that some of the bounds listed above are on comoving scales smaller than those probed by the data used in our analysis. The resolution to this apparent paradox lies in the shape and width of the DAO feature, which can be quite broad and affect a large range of scales (see Fig.~\ref{fig:correlation}). Therefore, even if the peak of the DAO feature (that is, $r_{\rm DAO}$) is outside the reach of the data considered, the tail of the DAO bump can have an effect on observable scales, hence the above limits.  For  $f_{\rm int}\lesssim2\%$, the DAO feature becomes very small and the DAO scale is thus largely unconstrained. If instead of fixing $f_{\rm int}$ we marginalize over it, the constraint reads
\begin{equation}
r_{\rm DAO} < 8.5h^{-1} {\rm Mpc}\quad(f_{\rm int} \text{ marginalized}),
\end{equation}
where again we are displaying the $95\%$ confidence limits for the ``Planck+WP+High-$l$+BAO+Lens'' dataset.

\section{Impact on Galaxy Formation: Double-Disk Dark Matter}\label{galform}
In this section, we turn our attention to the late-time consequences of having a fraction of the dark matter interacting directly with a massless gauge boson.  Much like the baryonic case, the coupling between DM and DR allows the interacting DM to cool via the emission of DR. This cooling can have a large impact on the structure of DM halos around galaxies and can possibly lead to the formation of a DM disk, as was shown in Refs.~\cite{Fan:2013tia,Fan:2013yva}. Here, we point out that the parameters required to obtain the large amount of cooling necessary to form a dark disk can also lead to large-scale cosmological signatures that might be incompatible with current data.
\begin{figure*}[t]
\begin{centering}
\subfigure{\includegraphics[width=0.49\textwidth]{DDDM_const_xi_05.pdf}}
\subfigure{\includegraphics[width=0.49\textwidth]{DDDM_const_xi_03.pdf}}\\
\subfigure{\includegraphics[width=0.49\textwidth]{DDDM_const_fADM_2perc_xi_05.pdf}}
\subfigure{\includegraphics[width=0.49\textwidth]{DDDM_const_fADM_2perc_xi_03.pdf}}\\
\caption{PIDM parameter space where a galactic dark disk is likely to form superimposed on the cosmological constraints from the ``Planck+WP+High-$l$+BAO+Lens'' dataset. The red regions below the black dot-dashed line are ruled out at $95\%$ confidence level, while the gray regions delimitated by the black long-dashed line are ruled out at $68\%$ confidence level. As indicated, interacting DM can cool in equilibrium in the green regions while only out-of-equilibrium cooling is possible in the yellow regions. Along the solid purple line, the cooling timescale ($t_{\rm cool}$, the minimum of either bremsstrahlung or Compton scattering) is equal to the age of the Universe ($t_{\rm Univ}$), while the solid blue line denotes the parameters for which the Coulomb equilibration time ($t_{\rm eq}$) is equal to the cooling time. The cyan solid line denotes the parameters for which the virial temperature ($T_{\rm vir}$) of the halo is equal to the dark atomic binding energy. The short-dashed black line shows the value of $\alpha_D$ for which the thermal relic abundance of atomic DM is equal to $f_{\rm int}\Omega_{\rm DM}$. Throughout, we take $n_{\bf e}=n_{\bf p}=7.3\times10^{-3}(f_{\rm int}/0.05)(m_D/{\rm GeV})^{-1}$ cm$^{-3}$ and assume a DM halo mass $M_{\rm DM}=10^{12}M_{\odot}$.}
\label{DDDM1}
\end{centering}
\end{figure*}

The formation of a dark disk in the PIDM scenario requires DM to cool on a fast enough time scale through either the emission of soft dark photons (bremsstrahlung) or through Compton scattering off colder dark photons\footnote{If dark atoms behave at all similarly to regular baryonic hydrogen, we note that molecular cooling could also be very important for the formation of the dark disk.}. Both mechanisms are only effective if the dark atoms are ionized which generally requires the virial temperature of the halo to be larger than the dark atomic binding energy. Since bremsstrahlung and Compton scattering are mostly effective at cooling the lighter dark electron, the timescale for dark proton and dark electron to equilibrate through Coulomb scattering must also be folded into the analysis. As long as the Coulomb scattering rate is faster than the overall cooling rate (through either bremsstrahlung or Compton scattering), the interacting DM cools in equilibrium and a dark disk generally arises. In the opposite scenario, the dark electrons cool faster than they can equilibrate with the heavy dark protons. Ref.~\cite{Fan:2013yva} argues that the overall interacting DM sector could also cool in this regime, but actual simulations will be required to determine the exact outcome. The expressions for the bremsstrahlung cooling, Compton cooling, and Coulomb equilibration times are respectively given in Eqs.~[23], [24], and [27] of Ref.~\cite{Fan:2013yva}.

We illustrate in Fig.~\ref{DDDM1} both the PIDM parameter space where the formation of a dark disk is plausible and our cosmological constraints on the model derived in the previous section. To ease the comparison with the work of Ref.~\cite{Fan:2013yva}, we display the constraints in the $m_{\bf e}-\alpha_D$ plane, where $m_{\bf e}$ is the mass of the dark electron\footnote{In terms of these variables, $\Sigma_{\rm DAO}\simeq\frac{2}{\alpha_D}\left(\frac{m_{\bf e}}{\rm eV}\right)^{-1}\left(\frac{m_D}{\rm GeV}\right)^{-1/6}$.}. We illustrate the double-disk DM parameter space for the optimistic case that DM is distributed according to a Navarro-Frenk-White profile with a characteristic scale $R_s=20$ kpc, leading to central number density $n_{\bf e}=n_{\bf p}=7.3\times10^{-3}(f_{\rm int}/0.05)(m_D/{\rm GeV})^{-1}$ cm$^{-3}$. Note that we have taken the interacting DM fraction inside the halo to be equal to the cosmological mean value, but in general the former might differ from the latter. The green and yellow (appearing as orange when superimposed on the ruled out red regions) regions show the parameter space where equilibrium and non-equilibrium cooling happen, respectively. In each panel, the red region below the black dot-dashed line shows the parameter space that is ruled out at $95\%$ confidence level by the ``Planck+WP+High-$l$+BAO+Lens'' dataset, while the gray regions denote PIDM models ruled out at the $68\%$ confidence level.

We first note that all of the double-disk parameter space with $\xi=0.5$ and $f_{\rm int}=5\%$ (top-left panel of Fig.~\ref{DDDM1}) is ruled out by cosmological data with high confidence. Models with larger values of the interacting DM fraction or of $\xi$ would result is even \emph{stronger} constraints and are therefore also ruled out. As $\xi$ is decreased, the cosmological constraints slowly weaken, opening up some interesting parameter space where a dark disk could form. The top-right panel of Fig.~\ref{DDDM1} illustrates the case with $\xi=0.3$ where we see that a model with $\alpha_D\sim0.01$ and $B_D\sim50$ eV has the right parameters to form a dark disk while lying within the two-sigma contour of current cosmological data. For $\xi\lesssim0.2$ and $f_{\rm int}=5\%$, the double-disk DM parameter space is largely unconstrained by cosmological data, as indicated by the green contours in Fig.~\ref{GamDAOvsxi_2}.

Another avenue to weaken the bounds on double-disk DM is to further reduce the interacting DM fraction below $5\%$. The lower panels of Fig.~\ref{DDDM1} show the constraints on the double-disk parameter space for $f_{\rm int}=2\%$. For $\xi=0.5$ (left bottom panel), we observe that a large swath of the parameter space where a dark disk could form is in good agreement with current cosmological data. As we discussed in Section~\ref{DAO+f}, the CMB and galaxy clustering data are well fitted by a model with $\xi=0.5$, $f_{\rm int}=2\%$, and $\Sigma_{\rm DAO}\sim10^{-2.5}$, which explain the large allowed region overlapping with the double-disk DM parameter space. As $\xi$ is decreased, this preferred region closes up, but most parameter values where DM can cool and form a disk remain within the allowed $95\%$ confidence region.

In summary, most PIDM models that could lead to a dark matter disk within galaxies are ruled out if $f_{\rm int}\gtrsim5\%$ or $\xi\gtrsim0.2$. Nevertheless, some interesting models remain viable if $f_{\rm int}\sim 2\%$. However, given the small fraction of interacting DM in these scenarios, it remains to be seen if such models could lead to a significant impact on galactic dynamics and on direct and indirect DM searches. Simulations will be necessary to assess the relevance of these allowed models on galactic scales. In any case, our results compellingly highlight the complementarity between the largest cosmological scales and the much smaller galactic scales in pinpointing the nature of dark matter.

\section{Discussion}\label{sec:conc}
In this paper we have shown that if all or a fraction of the DM were coupled to a bath of DR in the early Universe we expect the combined DM-DR system to give rise to acoustic oscillations of the dark matter until it decouples from the DR.  Much like the standard baryon acoustic oscillations, these DAO imprint a characteristic scale, the sound horizon of dark matter, on the matter distribution in the Universe. We have seen that a having such a fraction of interacting DM can lead to potentially unique signatures on the CMB and large-scale structure data. Although we have modeled the interacting DM and DR system as dark atoms coupled to a bath of dark photons, our results can be straightforwardly applied to a broad class of models that couple DM particles to various light relativistic species. These include, for instance, models where dark matter is coupled to light scalar states or models where the dark sector couples to light states via heavy mediators (analogous to the neutrinos coupling via the weak force).

We have determined that PIDM models with $\Sigma_{\rm DAO}\gtrsim 10^{-3}$, $\xi\gtrsim 0.2$, and $f_{\rm int}\gtrsim 5\%$ are generally in severe tension with the most recent cosmological data. For much lower values of $\Sigma_{\rm DAO}$, the fraction of interacting DM becomes largely unconstrained while the bounds on $\xi$ reflects the current limits on the effective number of relativistic species, $N_{\rm eff}$. It is particularly interesting that the transition between this last regime and the regime where $\xi$ is severely constrained happens for values of $\Sigma_{\rm DAO}$ similar to that of standard baryons (remember that for baryons, $\Sigma_{\rm BAO}\sim10^{-3.3}$). This is not a coincidence. For $\Sigma_{\rm DAO}\ll\Sigma_{\rm BAO}$, the kinematic decoupling of interacting DM happens much before the epoch of CMB last scattering and any change to the matter power spectrum is limited to scales smaller than the BAO scale. On the other hand, for $\Sigma_{\rm DAO}\gg\Sigma_{\rm BAO}$, interacting DM stays coupled to the DR bath after the epoch of CMB last scattering and the clustering of matter is affected on large cosmological scales, leading to severe constraints on these PIDM models.

For PIDM models with interaction strength equal or greater than that of baryons, we have determined that at most $\sim5\%$  of the DM could be interacting with a cosmologically-significant ($\xi\gtrsim0.3$) DR bath. For $\xi\gtrsim0.4$, the constraint is even more restrictive with $f_{\rm int}\lesssim4\%$. To our knowledge, this is the first time that the allowed deviation from a pure collisionless CDM scenario on large scales is rigorously quantified. The surprise here is that there exists a class of models with $\xi\sim0.5$, $\Sigma_{\rm DAO}\sim10^{-2.5}$, and $f_{\rm int}\sim2\%$ which provides a very good fit to the data, although the improvement to the fit is marginal compared to a simple $\Lambda$CDM model.  This class of models is nonetheless interesting since $\xi\sim0.5$ is the ``natural'' value that we expect if the visible and dark sectors were coupled above the weak scale. Moreover, these are the only PIDM models where the data actually prefers a non-vanishing value of the interacting DM, albeit only at the $\sim2\sigma$ level. Finally, these models are generally expected to impact galactic and possibly cluster dynamics due to the expected cooling of interacting DM via the emission of DR, implying that these models could be probed on a large range of scales.

We have also determined that current cosmological data allow a large fraction of DM interacting with strength less than that of standard baryons. For $\Sigma_{\rm DAO}\lesssim10^{-4.5}$, the fraction of interacting DM is largely unconstrained (see Fig.~\ref{GamDAOvsxi_4}) and the latter could therefore form all of the DM. Improving upon these constraints will require a prescription to model the small-scale non-linearities in PIDM scenarios. Since these models generally predict a different overall shape for the matter power spectrum, it will likely be necessary to run $N$-body simulations to determine how non-linear structures form and evolve. Depending on the exact PIDM model considered, these simulations could be a lot more involved than standard CDM simulations because of the interacting nature of a fraction (or all) of the DM. Even if we restrict the parameter space to regions where radiative processes such as cooling are inefficient, PIDM will generally be self-interacting at some level inside halos, which can affect their central density profiles. Encouragingly, recent work on $N$-body simulations \cite{Vogelsberger:2012ku,Rocha:2012jg,Peter:2012jh,Zavala:2012us} have started exploring self-interacting DM for some simple cases, and hence are building up the knowledge necessary to eventually conduct and interpret realistic PIDM simulations.
\begin{figure}[t!]
\begin{centering}
\includegraphics[width=0.5\textwidth]{CMB_BB_GamDAO.pdf}
\includegraphics[width=0.5\textwidth]{CMB_BB_f_ADM.pdf}
\caption{Lensed CMB B-mode spectra for various PIDM models. The top panel fixes $\xi=0.5$ and $f_{\rm int}=5\%$ while letting $\Sigma_{\rm DAO}$ varies. The bottom panel fixes $\xi=0.5$ and $\Sigma_{\rm DAO}=10^{-3}$ while varying the fraction of interacting DM. We also show the data from Ref.~\cite{Hanson:2013hsb} obtained by combining data from SPTpol \cite{Austermann:2012ga} and Herschel \cite{Pilbratt:2010mv} (more specifically, we illustrate the $(\hat{E}^{150}\hat{\phi}^{\rm CIB})\times\hat{B}^{150}$ cross correlation). }
\label{fig:b-modes}
\end{centering}
\end{figure}

A promising avenue to improve the constraints on PIDM is CMB lensing. Upcoming CMB polarization data from Planck and from ground-based telescopes will dramatically improve the reconstruction of the lensing potential, which itself depends strongly on the matter power spectrum. Since it can probe the lensing potential down to smaller angular scales \cite{Smith:2006nk}, the conversion of polarization E-modes into B-modes via gravitational lensing may well provide the strongest bounds on PIDM scenarios.  Moreover, since lensed B-modes are more sensitive to the matter distribution at redshifts higher than those probed by current and near-future galaxy surveys \cite{Smith:2006nk}, it is less affected by non-linearities and therefore has the potential to provide more reliable limits on PIDM. We illustrate a few lensed B-mode spectra in Fig.~\ref{fig:b-modes} as well as the first-detection data from Ref.~\cite{Hanson:2013hsb}. We observe that even a $2\%$ fraction of interacting DM can have a sizable effect on the B-mode spectrum.

One topic that we have not touched upon in this work is whether PIDM can address the apparent discrepancies between state-of-the-art CMB data and the Hubble parameter inferred from local measurements and supernovae on the one hand, and between the CMB data and the Sunyaev-Zel'dovich (SZ) cluster count on the other hand (see Refs.~\cite{planckXVI,Ade:2013lmv}). Since PIDM models naturally contain a DR component, they tend to prefer a higher value of the Hubble parameter which could ease the tension between the CMB and other probes of the local expansion rate. Moreover, since PIDM generically predicts a damping of small-scale power, it could potentially reconcile the SZ cluster count with the CMB data. We leave such an analysis for future work.

On that note, we point out that it would be interesting to explore the possible degeneracy between PIDM and models with massive neutrinos (active or sterile). While PIDM can indeed mimic some of the signatures of massive neutrinos like the damping of the matter power spectrum on small scales, the CMB could provide enough discriminatory power to distinguish the two scenarios. For instance, in the PIDM case, the DR radiation transitions from being tightly-coupled to the DM to a free-streaming state at dark kinematic decoupling while in the massive neutrino case, the neutrinos transition from a free-streaming state to a cold non-relativistic state when their temperature falls below their mass. Since these two types of transition impact the CMB differently, we expect the PIDM and the massive neutrino scenarios to be be distinguishable to some degree. We leave the exploration of these degeneracies to future work.

\section{Conclusion}\label{sec:conc}

In this work we have found that if at least $5\%$ of the dark matter was coupled to a bath of dark radiation in the early Universe, its sound horizon must lie on small non-linear scales. For a smaller interacting dark matter fraction, the DAO scale becomes progressively unconstrained as the DAO feature shrinks in importance.
Our results imply that more than $95\%$ of the dark matter must behave like collisionless CDM, long decoupled from any radiation component, on large cosmological scales.
Like most knowledge we have gleaned about dark matter since inferring its existence \cite{1933AcHPh...6..110Z,1937ApJ....86..217Z,1980ApJ...238..471R}, this constraint rests on the  gravitational influence of dark matter (or dark radiation) on the observable Universe.
Provided gravity is universal the pull of dark matter betrays its distribution even though it remains otherwise elusive and invisible.
Our conclusions further rely on the extraordinarily detailed description of the physics of cosmological perturbations in the early Universe that are now required to enable cosmologists to make precise predictions for the CMB and LSS observables --- and reveal ever more about the nature of our Universe.\\

{\noindent \bf Acknowledgements---} We thank  JiJi Fan, Andrey Katz, Lisa Randall, Matthew Reece, and Manoj Kaplinghat for useful discussions. We further thank Abhilash Mishra, Marius Millea, and Matthew Reece for insightful comments on an earlier version of this draft. We are also grateful to Olga Mena and Shun Saito for generously providing an initial version of the code for the BOSS galaxy power spectrum likelihood. This work was performed in part at the California Institute of Technology for the Keck Institute for Space Studies, which is funded by the W. M. Keck Foundation. The research of KS is supported in part by a National Science and Engineering Research Council (NSERC) of Canada Discovery Grant.  FYCR and KS thank the Aspen Center for Physics, where part of this work was completed, for their hospitality.  This research was supported in part by the National Science Foundation under Grant No. NSF PHY11-25915.  Part of the research described in this paper was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. Part of this work is supported by NASA ATP grant 11-ATP-090.

redpajama_set_name: RedPajamaArXiv