\section{Purcell effect and Umklapp scattering}
The significance of compact Wannier orbitals can be understood from a comparison to photon emission in atomic physics, which is directly proportional to the local density of optical states~\cite{Scully1997}. Replacing  electron orbitals in atoms with Wannier orbitals---and photons with  phonons---we expect the phonon emission rate to be modulated by a Wannier orbital formfactor and enhanced when the orbitals are  localized. As an illustration, we consider orbitals that are well-localized on the SL scale, such that the orbital radius $\xi$ is much smaller than the SL period $a$: $\xi \ll a$ (see schematic in Fig.~\ref{fig:model}a). Since Wannier orbitals define a formfactor that modulates the el-ph transition matrix element~\cite{Ziman1960,Koshino2019b}, the typical momentum transfer in the el-ph scattering processes will be of the order $q\sim 1/\xi$, i.e. much greater than the SL reciprocal vector. As a result, the el-ph scattering can ``outcouple" electrons in the flat bands to phonons with momenta far outside the SL Brillouin zone. Scattering by such large-momentum phonons cannot occur through the usual el-ph transitions in which phonon momentum transfer is $q\lesssim 2k_F$. The large-$q$ phonon emission enabled by compact Wannier orbitals is expected to produce a Purcell-like enhancement of the scattering rate. Indeed, theories of flat bands in moir\'e graphene predict Wannier orbitals that are well-localized on the SL  scale~\cite{Laissardiere2010,Kang2018,Koshino2018,Carr2019}. As we will see below, such compact Wannier orbitals can produce a considerable Purcell-like enhancement of el-ph scattering for phonons of wavelength $\lambda\lesssim a$.

\begin{figure}[t]
  \includegraphics[width=\columnwidth]{figA.pdf}
  \caption{{\bf Electron-phonon scattering and cooling power in a superlattice.} {\bf a} The schematic figure of the model we consider. The gray balls are the atoms and the yellow cloud is the Wannier orbital of electrons. Electrons hop between the Wannier orbitals. Black lines show the unit cell of the superlattice (SL) and the two arrows are the primitive vectors of the SL. The Brillouin zone and electron-phonon scattering processes of the model are shown in {\bf b}. The dashed line in the center is the first Brillouin zone of SL and the blue circle inside is the Fermi surface. Besides the normal scattering (purple arrow), the Umklapp processes (yellow and orange arrows, see panel {\bf e}.) contribute to the scattering. The temperature dependence of cooling power (Eq.~\eqref{eq:JG}) is in {\bf c.} and {\bf d.} The two figures are for different radii of Wannier orbitals $\xi=a/6$ and $\xi=a/3$, respectively. The color of the dots corresponds to the contribution from different scattering processes in {\bf b}. The results are for $t=4.67$ meV, $a=11.7$ nm, $s=1.5\times10^4$ ms$^{-1}$, $g^2=9.37\,{\rm{eV}\AA^{2}}$, and $n=0.25\times10^{12}$ cm$^{-2}$ (see Eqs.~\eqref{eq:He}, \eqref{eq:Hp}, and \eqref{eq:g_q}). Note that the bandwidth is $W=3t=14.01$~meV. {\bf e.} is a schematic picture of the dispersion in the SL Brillouin zone. The acoustic band is folded in the 1st Brillouin zone of SL. We define the normal process as the scattering that involves the phonons in the lowest-band phonons (purple line) and the rest as the $n$th-order Umklapp process.}
  \label{fig:model}
\end{figure}

It is  instructive to note an analogy with the Purcell-like enhancement of plasmon emission observed in experiments on nanowire systems~\cite{Akimov2007}, in which electronic orbitals confined to quantum dots were employed to boost the plasmon emission rate. Unlike the electronic states in this experiment, that were fixed-size, the Wannier orbitals in moir\'e graphene can be reshaped by changing the twist angle. This offers an appealing possibility to make Purcell enhancement tunable.

The Purcell enhancement can also be interpreted in a simple way in terms of electron Umklapp scattering, the process in which electron momentum changes by a SL reciprocal vector as illustrated in Fig.~\ref{fig:model}. Crucially, unlike  Umklapp processes in pristine graphene monolayer, which do not occur below room $T$, the small SL period in moir\'e graphene makes these processes relevant even at low $T$. Furthermore, since the flat bands result from strong Bragg scattering of electronic waves by the SL potential, the amplitudes of electronic Umklapp scattering in a SL are expected to be relatively strong.  This mechanism will activate scattering processes that are prohibited in pristine graphene; the amplitudes of these processes are given by the Fourier harmonics of the Wannier orbital, in agreement with the discussion above.

There are many different Umklapp processes arising in this manner, some of which are illustrated in Fig.~\ref{fig:model}b. In Fig.~\ref{fig:model}e, the acoustic-phonon band is shown folded in the SL Brillouin zone; different parts of this band are accessed through different Umklapp pathways as indicated by dotted lines. The purple arrow in Fig.~\ref{fig:model}b is the momentum-conserving ${\vec G}=0$ scattering by the low-energy phonons in Fig.~\ref{fig:model}e. The processes shown in orange and yellow involve higher-energy parts of the phonon band and violate momentum conservation, i.e., ${\vec G}\ne0$. To delineate the contributions of different processes we will refer to the former processes as normal scattering (shown in purple in Fig.~\ref{fig:model}e), and those by the higher parts of the phonon band as order-$n$ Umklapp scattering ($n=1,2,3,\cdots$).

We note parenthetically that in addition to the electron Umklapps one may also consider phonon Umklapps. However, since phonon Bragg scattering is considerably weaker than electron Bragg scattering~\cite{Koshino2019}, we will focus on the electronic Umklapps and ignore the phonon Umklapps. For the sake of simplicity, we consider the longitudinal acoustic phonon band only; the electron interaction with transverse acoustic phonons are subject to a similar  Purcell enhancement that will be discussed elsewhere.

\section{Model}
To explore the effect of compact Wannier orbitals on el-ph coupling, we employ a toy model in which the nearly-flat band is described as a tight-binding problem on a honeycomb lattice with the periodicity and symmetry matching those of moir\'e SLs; Wannier orbitals are modeled as gaussian functions $w({\vec r})\sim e^{-r^2/2\xi^2}$ with the radius $\xi$ adjusted to match the realistic orbital size, as pictured in Fig.~\ref{fig:model}a. The electron Hamiltonian reads
\begin{nalign}\label{eq:He}
H_{\rm e}=&\sum_{\bf k}
\left(\begin{array}{c}
\psi_{{\vec k},1}^\dagger \\ \psi_{{\vec k},2}^\dagger
\end{array}\right)^{\rm T}
\left(\begin{array}{cc}
0 & h_{\vec k} \\
h_{\vec k}^\ast & 0
\end{array}\right)
\left(\begin{array}{c}
\psi_{{\vec k},1} \\
\psi_{{\vec k},2}
\end{array}\right),\\
h_{\vec k}=&t\left(1+e^{-i\vec k\cdot\vec a_1}+e^{-i\vec k\cdot\vec a_2}\right)
\end{nalign}
where $t$ is the nearest-neighbor hopping, ${\vec a}_{1,2}=a\left(\pm\frac12,\frac{\sqrt3}2\right)$, $|\vec a_{1,2}|=a$, are primitive vectors and ${\vec k}=(k_x,k_y)$ is electron momentum. The phonon field is the displacement of background gray atoms placed with a distance $a_{\rm cc}\ll a$. Since temperatures of interest are much lower than graphene monolayer Debye's temperature, we focus on the $a_{\rm cc}\to0$ limit; $a_{\rm cc}$ does not appear in our calculation in this limit. The phonon Hamiltonian reads
\begin{nalign}\label{eq:Hp}
H_{\rm p}=&\sum_{\vec q, \vec G}\hbar\omega_{\vec q+\vec G}\hat a_{\vec q+\vec G}^\dagger\hat a_{\vec q+\vec G},
\quad
\omega_{\vec q+\vec G}=s|{\vec q+\vec G}|
,
\end{nalign}
where $\hbar\omega_{\vec q+\vec G}$ is the energy of the phonon with momentum $\vec q + \vec G$. The sum runs over the momenta $\vec q$ in the SL Brillouin zone and the reciprocal lattice vectors ${\vec G}=n_1{\vec G}_1+n_2{\vec G}_2$, $n_{1,2}=0,\pm 1,\pm 2...$, ${\vec G}_{1,2}=\frac{2\pi}a\left(\pm1,\frac1{\sqrt3}\right)$. Recent study of phonons in moir\'e graphene finds a dispersion that, at low energies, is similar to that of the monolayer graphene~\cite{Koshino2019,Koshino2019b}. Therefore, we used the sound velocity $s=1.5\times10^4\;{\rm ms^{-1}}$, a typical number observed in graphene~\cite{Koshino2019,Cong2019}. The electrons and phonons interact through deformation coupling,
\begin{nalign}
H_{\rm ep} = -\sum_{{\vec q},{\vec k},{\vec G},n} \frac{g_{\vec q+\vec G}}{\sqrt{V}} &\sqrt{\hbar\omega_{\vec q+\vec G}} \hat\psi^\dagger_{{ \overline{\vec k+\vec q}},n}\hat\psi_{{\vec k},n}\\
& \times (\hat a_{\vec q+\vec G}+\hat a_{-{\vec q-\vec G}}^\dagger).\label{eq:Heph}
\end{nalign}
Here, $V\equiv\sqrt3a^2N/2$ is the area of the system with $N$ being the number of the SL unit cells. We will use the notation  $\overline{\vec k+\vec q}$ to denote momenta $\vec k+\vec q+\vec G$ shifted to the first SL Brillouin zone by adjusting $\vec G$. The modulation of the electron-phonon coupling due to the Wannier orbital is described by a formfactor
\begin{equation}\label{eq:g_q}
g_{\vec q}=g\,e^{-q^2\xi^2/4},
\end{equation}
where we assumed $g^2=9.37\,{\rm eV\AA^2}$; the value of $g$ corresponds to the deformation potential $D=20\,{\rm eV}$~\cite{Chen2008,Efetov2010} and mass density $1.52\times10^{-10}\,{\rm kg\,cm^{-2}}$.

The ${\vec G}\ne{\vec 0}$ terms in Eq.~\eqref{eq:Heph}, which embody the enhancement of the effective el-ph coupling  due to Wannier orbital localization, appear as a consequence of the electron Bragg scattering by SL potential. In the absence of the Bragg scattering, the ${\vec G}\ne{\vec 0}$ terms are prohibited by the momentum conservation. On the other hand, when electron Bragg scattering creates well-localized Wannier function, as it does in moir\'e graphene~\cite{Kang2018,Koshino2018,Carr2019}, the momentum non-conserving $\vec G\ne\vec0$ el-ph processes become allowed. These scattering processes are illustrated in Fig.~\ref{fig:model}c. The solid circles represent a multisheet Fermi surface in the SL Brillouin zone (gray hexagon), pictured for carrier density such that the Fermi level is away from the Dirac points. The purple arrow in the SL Brillouin zone (gray hexagon) is the normal process that involves phonons with $\bf q$ in the first Brillouin zone. Besides the normal process, the Umklapp processes shown by orange and yellow arrows scatter the electron to the same final state as the normal one. The Umklapp processes contribute to the physics at  temperatures $T\gtrsim T^\ast=\hbar\omega_{\vec G}/2k_B$ ($k_B$ is Boltzmann constant) at which higher-energy phonons can be thermally activated (Fig.~\ref{fig:model}c). The quantity $T^\ast$ plays a role similar to that of the Debye temperature, however its value is much smaller than the crabin lattice Debye's temperature. For moir\'e graphene, it is $T^\ast={\cal O}(10K)$, which is close to the Bloch-Gruneisen $T$ of graphene monolayer~\cite{Koshino2019}. Therefore, the Umklapp processes contribute significantly to the electron scattering at not-too-low temperatures and, as we will see, are relevant for the temperatures of interest.

\begin{figure}[tb]
  \includegraphics[width=\columnwidth]{figC.pdf}
  \caption{{\bf Electron cooling rate and heat capacity.} {\bf a.} Temperature dependence of electron cooling rate (Eq.~\eqref{eq:J'})
  for Wannier-function radius $\xi=a_0/6$ and $\xi=a_0/3$. {\bf b.} Temperature dependence of electron heat capacity (Eq.~\eqref{eq:C_Te}). The error bars in the two figures are the estimated errors of numerical integral. The results are for $t=4.67$ meV, $a=11.7$ nm, $s=1.5\times10^4$ ms$^{-1}$, $g^2=9.37\,{\rm eV\AA^{2}}$, $n=0.25\times10^{12}$ cm$^{-2}$, and $T_{\rm ph}=10$ K. Note that the bandwidth is $W=3t=14.01$ meV.
  }
  \label{fig:coolingrate}
\end{figure}

\section{Electron cooling}
One quantity that is highly sensitive to the nature of el-ph scattering is the electron cooling. We evaluate the cooling power using the two-temperature model~\cite{Bistritzer2009,Song2012}. This model assumes the electrons and phonons are in equilibrium with different temperatures $T_{\rm e}$ and $T_{\rm ph}$, respectively. Using Fermi's golden rule~\cite{Bistritzer2009,Song2012} and the interaction Hamiltonian $H_{\rm ep}$ from Eq.~\eqref{eq:Heph}, the cooling power per a unit area reads $J=\sum_{\vec G}J_{\vec G}$, where
\begin{align}
J_{\vec G}\equiv\sum_{n,n'}\int\frac{dk^2}{(2\pi)^2}\frac{dk'^2}{(2\pi)^2}&(\varepsilon_{n\vec k}-\varepsilon_{n\vec k'})W_{n{\vec k}',n{\vec k}}^{\vec G}f_{n\vec k}(1-f_{n\vec k'}),\label{eq:JG}
\end{align}
is the contribution from the $\vec G$ phonon band with
\begin{nalign}\label{eq:summary:Wkk}
W_{n'{\vec k}',n{\vec k}}^{\vec G}&=\delta_{n'n}2\pi|g_{\overline{\vec k'-\vec k}+{\vec G}}|^2\omega_{\overline{\vec k'-\vec k}+{\vec G}}\\
&\times\left\{N_{\overline{\vec k'-\vec k}+{\vec G}}\delta(\varepsilon_{n'\vec k'}-\varepsilon_{n\vec k}-\omega_{\overline{\vec k'-\vec k}+{\vec G}})\right.\\
&\quad+\left.(N_{\overline{\vec k'-\vec k}+{\vec G}}+1)\delta(\varepsilon_{n'\vec k'}-\varepsilon_{n\vec k}+\omega_{\overline{\vec k'-\vec k}+{\vec G}})\right\}.
\end{nalign}
Here, $N_{\vec q}=\frac1{e^{\beta_{\rm ph}sq}-1}$ is the phonon Bose-Einstein distribution and $\beta=1/T$ is the inverse temperature. The total cooling power is the sum of the contributions from different scattering processes.

Figures~\ref{fig:model}c~and~\ref{fig:model}d show the numerical result of total cooling power $J$ (solid line) and $J_{\vec G}$ for $|{\vec G}|=0$ (normal), $|{\vec G}|=\frac{4\pi}{a\sqrt3}$ (first-order Umklapp), and $|{\vec G}|=\frac{4\pi}{a}$ (second-order Umklapp). In Fig.~\ref{fig:model}c, the contribution from $J_{{\vec G}\ne{\bf0}}$ is larger than that of the normal process $J_{\vec 0}$. Indeed, $J$ is an order of magnitude larger than $J_{\vec 0}$. This trend is robust against the change of $\xi$ as seen from the $\xi=a/3$ result in Fig.~\ref{fig:model}d. The significant contribution from $J_{{\vec G}\ne{\bf0}}$ is a consequence of two factors. First, each $J_{\vec G\ne\vec0}$ gives a contribution comparable to the normal process $J_{{\vec G}={\bf0}}$ when $T_e\gtrsim T^\ast$. Second, there are six bands for each order of Umklapp processes, which enhances the contribution by almost an order of magnitude. The low frequency of the phonons and the multiplicity of $|{\vec G}|\ne0$ bands results in a large contribution.

Experimentally, the cooling power affects the electron cooling rate, i.e., the decrease of electron temperature per a unit time. The cooling rate reads
\begin{align}\label{eq:J'}
J'(T_{\rm e},T_{\rm ph})=\frac{J(T_{\rm e},T_{\rm ph})}{C(T_{\rm e})\Delta T},
\end{align}
where $\Delta T\equiv T_e-T_{ph}$, $C(T_e)$ is the heat capacity, and $T_{\rm e}$ and $T_{\rm ph}$ are respectively the temperatures of electrons and phonons. We calculate the heat capacity using
\begin{align}\label{eq:C_Te}
C(T_{\rm e})=&\sum_n\int_{1BZ}\frac{dk^2}{(2\pi)^2}\frac{4k_B\beta_{\rm e}^2\varepsilon_{n{\vec k}}^2}{\cosh^2(\beta_{\rm e}\varepsilon_{n\vec k}/2)},
\end{align}
because the temperature is comparable to the bandwidth. Figure~\ref{fig:coolingrate} shows the cooling rate [Fig.~\ref{fig:coolingrate}a] and the electron heat capacity [Fig.~\ref{fig:coolingrate}b] between 1 K $\le T_{\rm e}\le$ 350 K. The heat capacity $C(T_{\rm e})$ shows a Schottky peak at $T_{\rm e}\sim20$ K, which manifests a crossover from the low-temperature ($k_BT_{\rm e}\ll t$) to high-temperature ($k_BT_{\rm e}\gg t$) regime. Interestingly, the cooling rate $J'(T_{\rm e},T_{\rm ph})\sim J/C$, pictured in Fig.~\ref{fig:coolingrate}a, increases monotonically and approximately linearly with $T_{\rm e}$ showing no peak, a behavior that reflect the decrease of $C(T_{\rm e})$ at the high temperature.

The $T_{\rm e}$-linear dependence of $J'(T_{\rm e},T_{\rm ph})$ in the high-$T_e$ is a manifestation of the asymptotic behavior in the high-$T_{\rm e}$ limit. In the limit $k_{\rm B}T_{\rm e}\gg k_{\rm B}T_{\rm ph}, \hbar\omega_{\vec G}, t$, the cooling power reads
\begin{align}
J(T_{\rm e},T_{\rm ph})\sim\frac{2\pi}\hbar\sum_{n}\int_{1BZ}\frac{dk^2}{(2\pi)^2}\eta_{n{\vec k},{\vec G}},
\end{align}
where
\begin{nalign}
\eta_{n{\vec k}}\equiv&\sum_{\vec G}\int_{1BZ}\frac{dk'{}^2}{(2\pi)^2}(\hbar\omega_{\overline{{\vec k}'-{\vec k}}+{\vec G}})^2|g_{n,\overline{{\vec k}'-{\vec k}}+{\vec G}}|^2\\
&\hspace{2cm}\times\delta\left(\varepsilon_{n{\vec k}'}-\varepsilon_{n{\vec k}}-\hbar\omega_{\overline{{\vec k}'-{\vec k}}+{\vec G}}\right),
\end{nalign}
and
\begin{align}
C(T_e)=&\frac1{k_{\rm B}T_{\rm e}^2}\sum_n\int_{1BZ}\frac{dk^2}{(2\pi)^2}\varepsilon_{n{\vec k}}^2.
\end{align}
The approximately linear $T$-dependence $J'(T_{\rm e},T_{\rm ph})\propto T_{\rm e}$, which is valid even in the high-$T_{\rm e}$ limit, arises because
$J\propto{\rm const.}$ and $C\propto T_{\rm e}^{-2}$ at high $T$.

\begin{figure}[tb]
  \centering
  \includegraphics[width=\columnwidth]{figD.pdf}
  \caption{{\bf Temperature dependence of resistivity $\rho(T)$ and the contributions from different scattering processes.} The two figures are for different Wannier orbital radius $\xi$: {\bf a.} $\xi=a/6$ and {\bf b.} $\xi=a_0/3$. The solid line shows the resistivity and the points are the contribution from normal, first-order Umklapp, and second-order Umklapp processes (Eq. \eqref{eq:rho_G}). The results are for $t=4.67$ meV, $a=11.7$ nm, $s=1.5\times10^4$ ms$^{-1}$, $g^2=9.37\,{\rm eV\AA^{2}}$ and $n=0.25\times10^{12}$ cm$^{-2}$. Note that the bandwidth is $W=3t=14.01$ meV.}
  \label{fig:rho}
\end{figure}

\section{Resistivity}
Umklapp processes also contribute to the electrical resistivity. The total resistivity $\rho$ of a system is related to the entropy production by the electron-phonon scattering. We consider only the longitudinal current as the transverse current is prohibited by symmetry. Within the variational method assuming the modulation of electron distribution $\delta f_{n\vec k}\equiv f_{n\vec k}-f_{n\vec k}^0=\tau e\vec E\cdot\vec v_{n\vec k}\beta f^0_{n{\vec k}}(1-f^0_{n{\vec k}})$~\cite{Ziman1960}, $\rho(T)$ follows Matthiessen's rule $\rho(T)=\sum_{\vec G}\rho^{\vec G}(T)$ where
\begin{align}\label{eq:rho_G}
\rho^{\vec G}(T)=&\frac{\sum_{n,n'}\int\frac{dk^2}{(2\pi)^2}\frac{dk'^2}{(2\pi)^2}\beta P^{\vec G}_{n{\vec k},{\vec k'-\vec k};n'{\vec k}'}\left(v^x_{n{\vec k}}-v^x_{n'{\vec k}'}\right)^2}{\left[e\sum_{n}\int\frac{dk^2}{(2\pi)^2}(v_{n{\vec k}}^x)^2\beta f^0_{n{\vec k}}(1-f^0_{n{\vec k}})\right]^2},
\end{align}
are the contributions of the phonon $\vec G$ bands with
\begin{align}
P_{n{\vec k},{\vec q};n'{\vec k'}}^{\vec G}=&2\pi\delta_{n,n'}|g_{\vec q+\vec G}|^2\omega_{\vec q+\vec G} N_{\vec q+\vec G}f^0_{n{\vec k}}\left(1-f^0_{n'{\vec k'}}\right)\nonumber\\
&\qquad\times\delta(\varepsilon_{n{\vec k}}+\hbar\omega_{\vec q+\vec G}-\varepsilon_{n'{\vec k}'}).
\end{align}
Here, $\vec E$ is the external electric field, and $\tau$ is the variational parameter. We assume $T=T_{\rm e}=T_{\rm ph}$. The temperature dependence of $\rho$ and $\rho^{\bf G}$ is shown in Fig.~\ref{fig:rho}. The main contribution comes from the first-order Umklapp process at temperatures $T\gtrsim T^\ast$, which is larger than that of the normal process.

\section{Electron-phonon coupling in moir\'e graphene}
To validate the conclusions drawn from the toy model, we made a detailed comparison to the microscopic continuum model of moir\'e graphene~\cite{Bistritzer2011} (See Supplementary Information for technical details). The results for different Umklapp processes, shown in Fig.~\ref{fig:gq}, are in good agreement with the toy model that links the Umklapp scattering rates to the Wannier function formfactor $g_q$. As Fig.~\ref{fig:gq} illustrates, the best fit to the microscopic calculation is indeed provided by tightly localized Wannier orbitals as in Eq.~\eqref{eq:g_q} with $\xi/a\sim1/8-1/5$, supporting the Purcell-like enhancement of el-ph scattering rates due to wannier orbitals localization. Furthermore, this best-fit value appears to be independent of the twist angle, adding an air of confidence to the toy model approach.

\begin{figure}[t]
  \centering
  \includegraphics[width=\columnwidth]{figE.pdf}
  \caption{{\bf The electron-phonon coupling form factor $g_{\bf q}$ for moir\'e graphene (Eq.~\eqref{eq:g_q}).} The disks and diamonds show the form factor for twist angles $\theta=1.05^\circ$ and $1.20^\circ$, respectively. They are calculated from the Wannier functions of the continuum model (See Supplementary Information). The gray lines are the gaussian form factor $g_{\bf q}\propto\exp(-\xi^2q^2/4)$ with different Wannier function radius $\xi$.}
  \label{fig:gq}
\end{figure}

\section*{Discussion}
We see that reshaping of Wannier orbitals impacts the observables such as the cooling power and resistivity. The enhancement in the el-ph scattering rates due to Wannier orbital compression originates from the same physics as the Purcell effect in quantum optics. The enhancement has a nontrivial dependence on phonon momentum, being particularly important for large-$q$ phonons. As the analysis above indicates, these phonons are responsible for a substantial enhancement of the cooling power and the linear-$T$ resistivity. The exact value for such an enhancement will depend on system details on the SL scale;  for realistic parameter values it can be as high as tenfold, or more.

A key parameter that controls the enhancement is the Wannier orbital radius $\xi$ which suppresses the electron-phonon scattering for $q\xi\gtrsim1$ (Eqs.~\ref{eq:Heph},~\ref{eq:g_q}).  Shrinking of Wannier orbitals, described by a reduction in $\xi$, leads to a stronger effective el-ph coupling arising because the large-$q$ phonons with $|q|\gtrsim\pi/a$ can now be activated. Indeed,   the el-ph formfactor $g_q\sim e^{-q^2\xi^2}$, Eq.~\ref{eq:g_q}, is small for large-$q$ phonons $|q|\gtrsim\pi/a$ when the Wannier orbitals extend throughout the SL unit cell, i.e., $\xi \sim a$, but increases as the Wannier orbitals shrink, $\xi\ll a$. This is exactly what happens in moir\'e graphene due to its localized Wannier orbitals~\cite{Laissardiere2010,Kang2018,Koshino2018,Carr2019}.

The tunability of the el-ph interactions of moir\'e materials through reshaping the Wannier orbitals has many interesting implications. The dramatic enhancement of the cooling and momentum relaxation rates, the two effects that  are straightforward to measure experimentally, will offer an independent confirmation of the Purcell enhancement mechanism. The connection between enhanced interactions and reshaped Wannier orbitals will provide unique means to control interactions in moir\'e bands. As a new direction in the study of moir\'e materials, it will facilitate the microscopic understanding of superconductivity and other strongly-correlated phenomena.

We thank Andrey Shytov for useful discussions. HI was supported by JSPS KAKENHI (Grant Numbers JP18H03676, JP18H04222, and JP19K14649) and UTokyo Global Activity Support Program for Young Researchers. LL acknowledges support from the Science and Technology Center for Integrated Quantum Materials, NSF Grant No. DMR-1231319; and Army Research Office Grant W911NF-18-1-0116. FG was supported by funding from the European Commission, under the Graphene Flagship, Core 3, grant number 881603, and by the grants NMAT2D (Comunidad de Madrid, Spain),  SprQuMat and SEV-2016-0686, (Ministerio de Ciencia e Innovaci\'on, Spain).

\section{Electron-phonon coupling for continuum model}

The electron-phonon coupling in graphene includes a scalar contribution, described by the deformation potential, and a gauge contribution, due to the modulation of the nearest neighbor electron hoppings induced by the phonons. We consider first the scalar coupling, discussed in the main text. In the twisted bilayer, the coupling between states in the narrow bands with momenta $\vec{k}$ and $\vec{k} + \vec{q}$ through a longitudinal acoustic phonon of momentum $\vec{q} + \vec{G}$, where $\vec{G}$ is a reciprocal lattice vector,  can be expressed as
\begin{align}
g_{\vec{k}, \vec{k} + \vec{q}} &= D \left\langle \psi_{\vec{k} + \vec{q}} \left| {\cal I}_\sigma {\cal I}_\tau e^{i ( \vec{q} + \vec{G} ) \vec{r}}  \right| \psi_{\vec{k}} \right\rangle
\label{mat}
\end{align}
where the wavefunction $\left|\psi_\vec{k} \right>$ is a four component spinor, whose components are labeled by the sublattice, $\sigma$, and layer, $\tau$. Note that the full electron-phonon coupling includes a factor proportional to $| \vec{q} + \vec{G} |$. The matrices ${\cal I}_\sigma$ and $ {\cal I}_\tau$ are $2 \times 2$ unit matrices in sublattice and layer indices. The use of the matrix $ {\cal I}_\tau$ implies that we consider phonons with displacements which are equal in both layers.  In order to estimate the relative strength of Umklapp scattering processes, we analyze the square of the matrix elements, which is the quantity which appears in the transport coefficients. We calculate the dimensionless quantity.
\begin{align}
M_{\vec{K} , \vec{G} } &=  \left| \left\langle \psi_{\vec{K}} \left| {\cal I}_\sigma {\cal I}_\tau e^{i  \vec{G}  \vec{r}}  \right| \psi_{\vec{K}} \right\rangle \right|^2
\label{mat2}
\end{align}
where, in addition, we have used the states at the corners of the Brillouin Zone, $\vec{K}$, as representatives of the $\vec{G}$ dependence of the matrix elements.  By construction, $M_{K , \vec{G} = 0} = 1$. For $\vec{G} \ne 0$, we average over equivalent vectors. The calculations are using the parameters in\cite{Koshino2018s} and an angle $\theta =1.015^\circ$.Results are shown in Fig.~\ref{fig:mat}.

\begin{figure}[b]
\begin{tabular}{cc}
\includegraphics[width=3.in]{umklapp_matrix_2_1p015.pdf} &
\includegraphics[width=3.in]{umklapp_matrix_2_1p2.pdf}
\end{tabular}
\caption{Matrix elements squared, see Eq.~\eqref{mat2}, for Umklapp processes involving different reciprocal lattice vectors. The reciprocal lattice vectors are labeled by two integers, $\{ m , n \}$ such that $\vec{G} = m \vec{G}_1 + n \vec{G}_2$, where $\vec{G}_1$ and$\vec{G}_2$ are the unit vectors of the reciprocal lattice. It is assumed that the atomic displacements induced by the phonon are equal in both layers.}
\label{fig:mat}
\end{figure}

\begin{figure}[h]
\begin{tabular}{cc}
\includegraphics[width=3.in]{umklapp_matrix_2_1p015_2.pdf} &
\includegraphics[width=3.in]{umklapp_matrix_2_1p2_2.pdf}
\end{tabular}
\caption{As in Fig.~\ref{fig:mat}, but assuming that the atomic displacements have opposite signs in the two layers.}
\label{fig:mat_2}
\end{figure}

The calculation can be extended to phonons whose displacements are opposite in the two layers. We have checked that the corresponding couplings squared are, at least, one order of magnitude smaller than the ones discussed previously, see Fig.~\ref{fig:mat_2}.

\begin{figure}[h]
\begin{tabular}{cc}
\includegraphics[width=3.in]{umklapp_matrix_2_1p015_vx.pdf} &
\includegraphics[width=3.in]{umklapp_matrix_2_1p2_vx.pdf}
\end{tabular}
\caption{As in Fig.~\ref{fig:mat}, but replacing the operator ${\cal I}_\sigma$ by $\sigma_x$ in Eq.~\eqref{mat} (gauge coupling, and assuming that the atomic displacements have the same sign in the two layers.}
\label{fig:mat_3}
\end{figure}

\begin{figure}[h]
\begin{tabular}{cc}
\includegraphics[width=3.in]{umklapp_matrix_2_1p015_vx_2.pdf} &
\includegraphics[width=3.in]{umklapp_matrix_2_1p2_vx_2.pdf}
\end{tabular}
\caption{As in Fig.~\ref{fig:mat}, but replacing the operator ${\cal I}_\sigma$ by $\sigma_x$ in Eq.~\eqref{mat} (gauge coupling, and assuming that the atomic displacements have opposite sign in the two layers.}
\label{fig:mat_4}
\end{figure}

For completeness, we have also considered the coupling to transverse acoustic phonons, where the coupling is proportional to the in-layer velocity operators, $v_x , v_y \propto \sigma_x , \sigma_y$. For a twisted graphene bilayer, the electronic states in the narrow bands are a superposition of Bloch states in both layers with almost opposite velocities\cite{Bistritzer2011s}. Hence, the coupling to transverse phonons with equal displacements in the two layers is suppressed by a factor proportional to the reduction in the Fermi velocity. This suppression does not take place if the displacements have opposite signs in the two layers. These antisymmetric phonons can contribute to the resistivity and the electron cooling. As in the case of longitudinal phonons, a Purcell like enhancement exists, and Umklapp processes have a significant contribution, see Figs.~\ref{fig:mat_3} and \ref{fig:mat_4}.

\FloatBarrier


\section{Introduction}
Ever since the first inference of the existence of dark matter almost a century ago, the nature of it has remained elusive.
Significant progress has been made in establishing its properties. For example it is known that dark matter can only interact weakly with standard model particles (perhaps only through gravitation), and that it must be \textit{cold} in the sense of having sufficiently small thermal velocity so that small scale structure can form.
The standard $\Lambda$CDM model of cosmology is based on dark matter with exactly these properties and in general provides an excellent fit to observational data.

However, despite many advances in both theory and observations much still remains unknown about the nature of dark matter. For example dark matter could have significant interactions, either with itself or with other particles in a dark sector, separate from the standard model.
Such models have been invoked as possible explanations of a variety of anomalies in cosmology, such as the missing satellites problem, the cusp-core problem, and the Hubble tension.

In this paper we investigate the possibility that dark matter is cold and consists of massive particles, but that these particles are unstable and decay.
The decay product cannot simply be electromagnetic radiation, since that would be detectable, so instead it is assumed that the decaying cold dark matter (DCDM) decays into massless particles in a dark sector. In line with previous work on the topic we shall refer to the decay product as \textit{dark radiation} (DR), so that we can sketch the process as
\begin{align*}
    X_{_{\mathrm{dcdm}}}\xrightarrow{}\gamma_{_{\mathrm{dr}}}\,.
\end{align*}
When comparing a model where all dark matter is DCDM with our observational data (Planck, etc.), we need a very large lifetime which renders our DCDM basically stable. Because of that, we add another degree of freedom to our model which is the initial fraction of DCDM to the total amount of cold dark matter (CDM),
\begin{equation}
    f_{\mathrm{dcdm}}=\frac{\Omega^{\rm ini}_{\mathrm{dcdm}}}{\Omega^{\rm ini}_{\mathrm{dcdm}}+\Omega_{\mathrm{scdm}}}\,,
\end{equation}
where SCDM is the stable CDM component. The physical interpretation of having this fractional decay could be~\cite{Poulin_16}:
\begin{enumerate}[label={\arabic*)}]
    \item CDM is multi-component such that only a part of the dark matter decays,
    \item all dark matter decays, but the decay product is both dark radiation and a stable CDM component with fraction $(1-f_{\mathrm{dcdm}})$.
\end{enumerate}
Since the amount of dark matter decreases, we define $\Omega^{\mathrm{ini}}_{\mathrm{dcdm}}$ as the density parameter of DCDM today as if none of it had decayed. We otherwise adopt the same notation as in Ref.~\cite{MaBert}, which is now standard.

\subsection{Previous work and constraints}\label{sec:previous}
The field of decaying cold dark matter has been growing ever since the first analysis by Ref.~\cite{Ichiki_2004} in 2004, where the simple model with dark matter decaying into dark radiation was also used. A lot has happened in the last two decades with the field being much wider now with numerous models of varying sophistication. The different regimes of the decay rate, $\Gamma_{\mathrm{dcdm}}$, is investigated in Ref.~\cite{Poulin_16} and the fractional amount of DCDM is constrained to $f_{\mathrm{dcdm}}<3.8\times10^{-2}$ using Cosmic Microwave Background (CMB) data from Planck-2015. This is in agreement with the results found in Ref.~\cite{Chudaykin_2016} where the same CMB data was used. In this paper we update the results and constraints from Ref.~\cite{Poulin_16} using the newest Planck-2018 data as well as Baryonic Acoustic Oscillations (BAO) data.

It is well known that allowing dark matter to decay can potentially relieve the Hubble tension as proposed in Ref.~\cite{Berezhiani_2015}, and this has been demonstrated numerous times in the literature. This includes Ref.~\cite{Pandey_2020} where both the Hubble tension and the milder $S_8$ tension of matter fluctuations between $\Lambda$CDM-model and the model-independent measurements in the local Universe were relieved using CMB data from Planck-2018 and the same DCDM model as we are using in this paper. Including other data sets such as BAO data and intermediate-redshift data, however, leads to only a slight decrease in both tensions. Another attempt at relieving the Hubble tension is made in Ref.~\cite{Chudaykin_2018} where the tension is reduced by $1.5\sigma$ using Planck-2015 data, BAO data, Redshift Space Distortion (RSD) measurements and a DCDM model with decay time after recombination. By allowing the CMB lensing power amplitude to be a free fitting parameter, the tension is reduced by $3.3\sigma$. They furthermore found an upper limit on the fractional amount of DCDM at the level $f_{\rm dcdm}\lesssim5\%$ for DCDM decaying after recombination. The fractional amount of short-lived DCDM is also constrained in Ref.~\cite{Xiao_2020} where CMB data from Planck-2015, BAO data, and RSD measurements are used to infer the upper limit $f_{\rm dcdm}\lesssim 2.73\%$, but this only leads to a slight reduction of the Hubble tension. They furthermore show the impacts of DCDM on the kinetic Sunyaev-Zel'dovich effect and note that this could lead to further constraints on DCDM in the future.

Better and more data has allowed the constraints to be more refined, with different studies using different combinations of data sets. An upper bound on the decay rate of late time DCDM at $\Gamma_{\mathrm{dcdm}} \leq (175 \,\mathrm{Gyr})^{-1}$ is inferred in Ref.~\cite{Enqvist_2020} by using CMB data and cluster counts from Planck-2015 along with weak lensing data from KiDS450, while a more tight constraint of $\Gamma_{\rm eff}<9.1\times 10^{-9}\,{\rm Gyr^{-1}}$ is found in Ref.~\cite{Oldengott_2016} to the effective decay rate by using Planck-2015 data and analysing cosmic reionisation and dark matter decay simultaneously. A different approach is taken in Ref.~\cite{Blanco_2019} where they assume that the decay product is detectable, i.e. a pair production of an elementary particle from the Standard Model (quarks, leptons or bosons), which leads to lower limits on the lifetime in the range $\tau\sim(1-5)\times 10^{28}\,\mathrm{s}$ using data from the isotropic gamma-ray background.

There are many studies beyond the simple DCDM model, where the decay product is only DR, and a more agnostic approach is found in Ref.~\cite{Bringmann_2018} where the conversion of dark matter to dark radiation is treated in a more general aspect without assuming the transition being caused by a decay. They then find the impacts of the conversion on the CMB spectrum as well as constraints on their general model parameters using \textsc{camb}\, and \textsc{CosmoMC}, which lead to a reduction in the Hubble tension. They then treat a conversion model in detail where dark matter particles interact via a light mediator particle leading to Sommerfeld-enhanced self-annihilation.

Another, more general analysis is done in Ref.~\cite{Dienes_2017} where they investigate Dynamical Dark Matter in which the dark sector is comprised of a large ensemble of particles with different decay widths. Their analysis shows that the constraints allow energy scales ranging from GeV scale to the Planck scale, but the cosmological abundance of the dark sector today must be spread across an increasing number of states in the ensemble as the energy scale is decreased down to GeV scale from the Planck scale. A third continuation of the DCDM model is found in Ref.~\cite{Raveri_2017} which features Partially Acoustic Dark Matter, where a subdominant part of the dark matter is strongly coupled to the DR fluid, and this combined fluid undergoes acoustic oscillations below the effective dark sound horizon. This is used to reduce the tensions in the Hubble constant and the lensing amplitude between CMB data and direct measurements. They find that even though the model can be used to reduce the tensions separately, it cannot accommodate both at once, since additional CDM is required by the CMB data to preserve the shape of the acoustic peaks.

It has recently been increasingly popular to assume a slightly warm component of the decaying dark matter sector, and an analysis of the background equations can be found in Ref.~\cite{Vattis_2019}, where the cold dark matter decays to both DR and a warm daughter particle. They show that this is can potentially relieve the Hubble tension, which makes the model very interesting to future work in the field of cosmology. More thorough analyses of the two-body decay model are found in Refs.~\cite{Blackadder:2014wpa,Clark:2020miy,Haridasu:2020xaa}, and the same model has also been shown to potentially relieve the $S_8$ tension in various studies, including Refs.~\cite{Wang:2012eka,Abellan:2020pmw,Abellan:2021bpx}. This is due to the recoil velocity of the massive daughter particle inducing a free-streaming suppression of matter fluctuations. The model has also shown promise when it comes to the problems of small-scale structure, and treatments of this using N-body simulations can be found in Refs.~\cite{Borzumati:2008zz,Peter:2010jy,Peter:2010sz,Wang:2014ina}. More general studies, where both daughter particles can have arbitrary masses, are treated in Refs.~\cite{Aoyama:2011ba,Aoyama:2014tga}. A slightly different model, where the decaying component is warm and the decay product is only DR, is treated thoroughly in Ref.~\cite{Blinov_2020}. The idea of this model is to not modify the evolution of the gravitational potentials, which otherwise leads to inconsistencies with data as in the case of decaying cold dark matter. They find that this can significantly reduce the Hubble tension as well.

\subsection{Outline of this paper}
In this paper, we are updating the constraints on the DCDM model parameters from Ref.~\cite{Poulin_16} using the Markov Chain Monte Carlo (MCMC) sampler \textsc{MontePython}~\cite{montepython} and the Einstein--Boltzmann code \textsc{class}~\cite{class}. To describe the cosmological framework we will use the following set of parameters:
\begin{equation}
    \Theta = \{\Omega_{\rm b} h^2, \Omega_{\rm cdm} h^2, h^2, A_s, n_s, \tau_{\rm reio} \},
\end{equation}
in addition to the decay parameter vector $\{ f_{\rm dcdm}, \Gamma_{\rm dcdm} \}$.
As was also done in Ref.~\cite{Poulin_16}, we will split our MCMC runs into two categories, named \textit{short-lived} and \textit{long-lived}. The reason for this is that the likelihood contour in the full parameter space has a shape which makes convergence exceedingly slow. When the lifetime becomes very short, essentially all DCDM has decayed well before matter-radiation equality. This makes it impossible to constrain $\Gamma_{\rm dcdm}$, and the only constrainable quantities are then $f_{\rm dcdm}$ and the product $\Gamma f_{\rm dcdm}$. This leads to a funnel-like likelihood surface stretching towards very large values of $\Gamma_{\rm dcdm}$. This particular part of the parameter space is best probed using a logarithmic prior on $\Gamma_{\rm dcdm}$, whereas for the long-lived regime, where only a fraction of the dark matter has decayed before the present, a prior which is flat in $\Gamma_{\rm dcdm}$ is more suitable. We will elaborate on the technicalities of this split in section~\ref{sec:MCMC}.

We will furthermore look at an analogous scheme to the short-lived DCDM, i.e. a model with an increase in $N_{\rm eff}$ instead of a decaying dark matter component. This is treated both analytically and numerically using the \textsc{class}\, code.

The structure of this paper is as follows. We introduce the formal theoretical framework of DCDM in the synchronous gauge in section~\ref{sec:boltzmann}, where both the background equations and the perturbation equations are presented. In section~\ref{sec:Neff_mapping} we will treat the mapping between a correction to $N_{\rm eff}$ and short-lived DCDM by deriving an analytical expression for $\Delta N_{\rm eff}$ corresponding to a short-lived DCDM component and comparing this to numerical results from \textsc{class}. In section~\ref{sec:MCMC} we will present our results from the MCMC sampler \textsc{MontePython} along with improved constraints in the long-lived and short-lived regimes, and in section~\ref{sec:tensions} we will investigate the impact of our model on the Hubble and $\sigma_8$ tensions. Lastly, we will conclude and summarise in section~\ref{sec:conclusion}.

\subsection{Cosmological data}
As our data sets we use the newest CMB data from Planck-2018~\cite{Planck2018} which has a higher quality in the polarisation data than its predecessor from 2015~\cite{Planck2015}. In all of our computations using Planck-2018, we use both polarisation and temperature likelihoods for both high-$\ell$ and low-$\ell$ as well as the lensing likelihood. When comparing to Planck-2015 data, we of course use the corresponding likelihoods from this data release, which is the same combination used in Ref.~\cite{Poulin_16}.

We furthermore use the BAO data from BOSS data release 12 (DR-12)~\cite{bao}. While the BAO data presumably has little impact in parameter constraints in the short-lived regime, due to the decay happening much earlier than the formation of the BAO, it is quite important in the limit of very long-lived DCDM.

\section{Boltzmann equations for DCDM and DR}\label{sec:boltzmann}
The behaviour of DCDM and DR can be calculated using the Boltzmann equation, which takes the following form for DCDM~\cite{Poulin_16}:

\begin{equation}\label{eq:boltzmann}
    \dv{f}{\tau}=\pdv{f}{\tau}+\pdv{f}{x^i}\dv{x^i}{\tau}+\pdv{f}{p}\dv{p}{\tau}+\pdv{f}{\hat{p}^i}\dv{\hat{p}^i}{\tau}=\pm a\Gamma_{\mathrm{dcdm}}f_{\mathrm{dcdm}}\,,
\end{equation}
with $-$ and $+$ for DCDM and DR respectively.

\subsection{Background equations}
The zeroth moment of the Boltzmann equation, eq.~\eqref{eq:boltzmann}, i.e. integrating it over phase-space and keeping only terms of zeroth order, leads to the continuity equation, which is different for DCDM and DR than for the homogeneous universe in having source terms dependent of the decay rate $\Gamma_{\mathrm{dcdm}}$ with respect to proper time~\cite{Poulin_16}:
\begin{equation}\label{eq:dcdm_fluid}
\begin{aligned}
    \rho'_{\mathrm{dcdm}}&=-3\frac{a'}{a}\rho_{\mathrm{dcdm}}-a\Gamma_{\mathrm{dcdm}}\rho_{\mathrm{dcdm}}\,,\\
    \rho'_{\mathrm{dr}}&=-4\frac{a'}{a}\rho_{\mathrm{dr}}+a\Gamma_{\mathrm{dcdm}}\rho_{\mathrm{dcdm}}\,,
\end{aligned}
\end{equation}
where the prime denotes derivatives with respect to conformal time, $\tau$.

\subsection{Perturbation equations}
One way of obtaining the relevant perturbation equations is again through the Boltzmann equation, eq.~\eqref{eq:boltzmann}. We are doing our calculations in the comoving synchronous gauge, so we need to express the perturbation equations in this gauge. Perturbations in the synchronous gauge can be expressed in the following way using the space-time interval~\cite{MaBert}:
\begin{equation}
    \d s^2=\,a^2(\tau)\left(-\d\tau^2+\left[\delta_{ij}+h_{ij}(\vec{x},\tau)\right]\d x^i\d x^j\right)\,,
\end{equation}
where the scalar part of the perturbation $h_{ij}(\vec{x},\tau)$ can be expressed through its Fourier transform
\begin{equation}
h_{ij}(\vec{x},\tau)=\mathop{\mathlarger{\mathlarger{\mathlarger{\int}}}}\d^3k\,\mathrm{e}^{i\vec{k}\cdot\vec{x}}\left(h(\vec{k},\tau)\hat{k}_i\hat{k}_j +6\eta(\vec{k},\tau)\left[\hat{k}_i\hat{k}_j-\frac{1}{3}\delta_{ij}\right]\right)\,,
\end{equation}
with $h(\vec{k},\tau)$ and $\eta(\vec{k},\tau)$ being the metric perturbations in the synchronous gauge. Integrating the Boltzmann equation, eq.~\eqref{eq:boltzmann}, over phase-space and keeping the first order terms leads to an expression for the evolution of the density perturbation, $\delta_{\mathrm{dcdm}}=\rho_{\mathrm{dcdm}}/\bar{\rho}_{\mathrm{dcdm}}-1$, where the bar denotes the average value as in a homogeneous universe. We can furthermore find the evolution of the divergence, $\theta_{\mathrm{dcdm}}$, of the fluid velocity by taking the divergence of the first moment of eq.~\eqref{eq:boltzmann}, which is found by multiplying the equation by $\overrightarrow{p}/E$ (with $\overrightarrow{p}$ and $E$ being the 3-momentum and the energy of a particle, respectively) and again integrating over phase-space. The resulting equations are
\begin{align}
    \delta'_{\mathrm{dcdm}}&=-\frac{h'}{2}\,,\\
    \theta'_{\mathrm{dcdm}}&=-\mathcal{H}\theta_{\mathrm{dcdm}}=0\,.
\end{align}
These two equations are enough to describe the evolution of DCDM since it per definition is cold, which means that we have neglected all second (or higher) order terms of the momentum in our calculations~\cite{dodelson}. All higher moments would therefore be zero.

The equations turn out a bit more complicated for DR, since this species is not cold and we therefore cannot neglect higher orders of momentum. Following the same procedure as for DCDM in the synchronous gauge, we get the following equations for the evolution of the density perturbations and the velocity divergence:
\begin{align}
    \delta'_{\mathrm{dr}}&=-\frac{2}{3}h'+a\Gamma_{\mathrm{dcdm}}\frac{\rho_{\mathrm{dcdm}}}{\rho_{\mathrm{dr}}}(\delta_{\mathrm{dcdm}}-\delta_{\mathrm{dr}})\,,\\
    \theta'_{\mathrm{dr}}&=\frac{k^2}{4}\delta_{\mathrm{dr}}-k^2\sigma_{\mathrm{dr}}-a\Gamma_{\mathrm{dcdm}}\frac{\rho_{\mathrm{dcdm}}}{\rho_{\mathrm{dr}}}\theta_{\mathrm{dr}}\,.
\end{align}
These certainly do not look as simple as those for DCDM, and we notice the parameter $\sigma_{\mathrm{dr}}$ in the bottom equation which is the next moment of the Boltzmann equation - the shear stress. Generally the evolution of a moment, $l$, will depend on the next moment, $l+1$, which leads to an infinite Boltzmann hierarchy for the moments containing the same information as the momentum-dependent Boltzmann equation itself. We still need to truncate the hierarchy at some large $l$-value to work with it numerically. In the \textsc{class}\, code, the cut-off $l$-value can be user-specified and is $l_{\mathrm{max}}=17$ by default.

\section{Mapping short-lived DCDM to $N_{\rm eff}$}\label{sec:Neff_mapping}

In the very short-lived regime, all DCDM has decayed well before matter-radiation equality. In this regime the primary effect of DCDM should be to enhance $N_{\rm eff}$ through the decay product, DR.
We define the correction, $\Delta N_{\rm eff}$, as the ratio between the energy densities of the additional radiation (which is DR) and a single massless neutrino. However, we need to evaluate this ratio after the DCDM has fully decayed, where no more DR is produced and the energy density scales as $a^{-4}$ (like any type of radiation),
\begin{equation}\label{eq:dN_eff_def}
\Delta N_{\rm eff} = \left.\frac{\rho_{\rm dr}}{\rho^{N=1}_{\nu}}\right|_{t\gg t_d}\,,
\end{equation}
where $\rho^{N=1}_{\nu}=7/8\,(4/11)^{4/3}\rho_{\gamma}$ represents only a single massless neutrino~\cite{Lesgourgues:2018oca}, and $t_d$ is the time of the decay.

We can estimate how $\Delta N_{\rm eff}$ should scale by assuming an instant decay, so the energy density of DR just after the decay equals that of DCDM just before the decay. We can thus evaluate the energy density of DCDM at the time of decay instead of that of DR. We can then scale that to the current time, introducing $a_d$ as the scale factor at the time of decay
\begin{equation}
\Delta N_{\rm eff} \approx \left.\frac{\rho_{\rm dcdm}}{\rho^{N=1}_{\nu}}\right|_{t=t_d} \sim a_d\cdot \frac{f_{\rm dcdm}}{1-f_{\rm dcdm}} \cdot\frac{\Omega_{\rm dm,0}}{\Omega^{N=1}_{\nu,0}}\,.
\label{eq:neff}
\end{equation}
Of course this is a very simplistic calculation, but it shows that, as $a_d \to 0$ (and thus $\Gamma_{\rm dcdm} \to \infty$) the effect of DCDM vanishes and the model becomes identical to standard $\Lambda$CDM.

\subsection{Analytic solution of Boltzmann equations}\label{sec:Neff_analytical}

We can derive an analytic approximation to the Boltzmann equations fairly easily.
First, we define
\begin{equation}
Y = a^3\,\rho_{\rm dcdm}/\rho_{\nu,0}\,, \qquad X = a^4\,\rho_{\rm dr}/\rho_{\nu,0}\,,
\end{equation}
i.e. $X$ and $Y$ are constant in the absence of decays. Here $\rho_{\nu,0}$ is the energy density of all neutrinos today as described by the effective number $N_{\rm eff}$. The Boltzmann equations~\eqref{eq:dcdm_fluid}, can then be expressed in proper time as
\begin{equation}\label{eq:N_eff_boltzmann}
    \dv{Y}{t} = -\Gamma Y(t)\,, \qquad    \dv{X}{t} = a\Gamma Y(t)\,,
\end{equation}
where we have dropped the subscript on the decay rate $\Gamma$. The Boltzmann equation for the decaying component then has the solution
\begin{equation}
    Y(t) =Y_i e^{-\Gamma t}\,,
\end{equation}
\begin{equation}\label{eq:Yi}
    Y_i = a_i^3\frac{(\rho_{\rm dcdm})_i}{\rho_{\nu,0}}=\frac{\Omega^{\rm ini}_{\rm dcdm}}{\Omega_{\nu,0}}=\frac{\Omega_{\rm dm,0}}{\Omega_{\nu,0}}\cdot\frac{f_{\rm dcdm}}{1-f_{\rm dcdm}}\,,
\end{equation}
where the subscript $i$ represents some initial starting point before the decay, and the last equal sign assumes that all DCDM has decayed today, thus making $\Omega_{\rm dm,0}$ both the current density parameter of all dark matter and that of only stable dark matter, since no decaying component is left.

The Boltzmann equation for DR can now be solved using the solution to that of DCDM. We switch coordinates to the scale factor divided by the initial scale factor at the starting point, $\tilde a \equiv a/a_i$, and assume that the Universe is radiation dominated throughout the decay so that $a \propto t^{1/2}$. We furthermore define
\begin{equation}
    \alpha \equiv \frac{H_i}{\Gamma} = \left(\frac{8\pi G}{3}\rho_{\rm r,0}\,a_i^{-4}\right)^{1/2}\cdot\frac{1}{\Gamma}\,,
\end{equation}
where $\alpha \gg 1$. We then integrate the differential equation from the starting point $\tilde a = 1$ to find the solution for $X$ (a similar result can be found in Ref.~\cite{Scherrer:1984fd})
\begin{equation}
    X(\tilde a) = X_i + \left[\sqrt{\frac{\pi}{2}}{\rm erf}\left(     \frac{\tilde a}{\sqrt{2 \alpha}}  \right)- \frac{\tilde a e^{-\tilde a^2/(2\alpha)}}{\sqrt{\alpha}}\right]\left(\frac{8 \pi G }{3}\rho_{\rm r,0}\right)^{1/4} \frac{\Omega_{\rm dm,0}}{\Omega_{\nu,0}}\cdot\frac{f_{\rm dcdm}}{1-f_{\rm dcdm}}\Gamma^{-1/2}\,,
\end{equation}
where $X_i$ is the initial value of $X(\tilde a)$ proportional to the initial DR component (here we have $X_i=0$ since there is no pre-existing DR component). In the asymptotic limit ($\tilde a \to \infty$) we find
\begin{equation}
    X(\tilde a \to \infty) = X_i + \sqrt{\frac{\pi}{2}}\left(\frac{8 \pi G \rho_{\rm r,0}}{3}\right)^{1/4} \frac{\Omega_{\rm dm,0}}{\Omega_{\nu,0}} \cdot\frac{f_{\rm dcdm}}{1-f_{\rm dcdm}}\Gamma^{-1/2}\,,
    \label{eq:neff2}
\end{equation}
and we finally note that $\Delta N_{\rm eff} = N_{\rm eff} X(\tilde a \to \infty)$ because of our definition of $X$.
Here we notice the scaling relation $\Delta N_{\rm eff} \propto \Gamma^{-1/2} f_{\rm dcdm}/(1-f_{\rm dcdm})$. Inserting numbers in eq.~\eqref{eq:neff2}, we arrive at an approximate relation of the form
\begin{equation}
    \Delta N_{\rm eff} \sim 5.74 \, \Omega_{\rm dm,0} \, h^2 \, \frac{f_{\rm dcdm}}{1-f_{\rm dcdm}}\Gamma_6^{-1/2}\,,
\end{equation}
where $\Gamma_6 = \Gamma/(10^6 \, {\rm Gyr}^{-1})$.

\subsection{Numerical analysis using \textsc{CLASS}}\label{sec:N_eff_class}

Using the \textsc{class}\, code, we can get an exact numerical correction to $N_{\rm eff}$ caused by short-lived DCDM. The values of input parameters which are held constant in all of the calculations can be seen in table~\ref{tab:class_input} ($N_{\rm eff}$ is modified by $\Delta N_{\rm eff}$ in figure~\ref{fig:CMB}).

\begin{table}[t]
\centering
\setlength\extrarowheight{4pt}
\begin{tabular}{| c c c c c c c  |}
\hline
\multicolumn{1}{| c }{$\Omega_{\rm b}h^2$}&\multicolumn{1}{ c }{$\Omega_{\rm cdm}h^2$}&\multicolumn{1}{ c
}{$h^2$}&\multicolumn{1}{ c
}{$A_s \times 10^{9}$}&\multicolumn{1}{ c
}{$n_s$}&\multicolumn{1}{ c
}{$\tau_{\rm reio}$}&\multicolumn{1}{ c |}{$N_{\rm eff}$}\\
\hline
0.022032  &
0.11933  &
0.67556  &
2.215  &
0.9619  &
0.0925  &
3.046 \\
\hline
\end{tabular}
\caption{Table of input parameters in \textsc{class}\, which are held constant through all calculations in section~\ref{sec:N_eff_class}. $N_{\rm eff}$ is, however, modified with $\Delta N_{\rm eff}$ in figure~\ref{fig:CMB}.}
\label{tab:class_input}
\end{table}

Running \textsc{class}\, with a high $\Gamma_{\rm dcdm}$ ($>10^3\,\mathrm{Gyr}^{-1}$) ensures that all DCDM has decayed well before today, which means that the energy density of the resulting DR scales like normal radiation. DR and neutrinos therefore scale similarly, which makes the ratio of their energy densities constant after the decay. This ratio is used to define the correction $\Delta N_{\rm eff}$ cf. eq.~\eqref{eq:dN_eff_def}. However, because we do not only have a single neutrino, we need to divide the neutrino density with $N_{\rm eff}$. The numerical value of $\Delta N_{\rm eff}$ can in fact be evaluated at any time after the decay as long as the energy density of DCDM is no longer of any significant size. We, however, evaluate at the current time in order to get a density of DCDM as low as possible,
\begin{equation}\label{eq:dN_eff_class}
\Delta N^{\rm (numerical)}_{\rm eff} = \left. N_{\rm eff}\frac{\rho_{\rm dr}}{\rho_{\nu}}\right|_{t=t_0},
\end{equation}
where $\rho_{\nu}$ is the total energy density of all massless neutrino species.

At first, we are interested in how well our analytical approach fits the numerical results. Using eq.~\eqref{eq:dN_eff_class}, we calculate $\Delta N_{\rm eff}$ for different decay rates, $\Gamma_{\rm dcdm}$, and fractional amounts of DCDM, $f_{\rm dcdm}$, and plot the results as functions hereof. This can be seen in the figures~\ref{fig:dNeff_Gamma} and \ref{fig:dNeff_f} respectively along with the analytical results and the ratio between the two. The numerical and analytical results agree with increasing precision for higher values of $\Gamma_{\rm dcdm}$, which is also supported by their ratios approaching unity. Here we also note that both results approach zero for increasing decay rate, which supports that we should recover the $\Lambda$CDM model for $\Gamma_{\rm dcdm}\to\infty$ as we argued in the beginning of section~\ref{sec:Neff_mapping}. From figure~\ref{fig:dNeff_f} we, however, see that the analytical and numerical results agree with decreasing precision for higher values of $f_{\rm dcdm}$, which makes sense according to eq.~\eqref{eq:neff2} due to the divergence at $f_{\rm dcdm}=1$. It should of course behave this way because we cannot have that short-lived DCDM comprises all of the dark matter, since that would correspond to no dark matter at all after recombination, which contradicts some assumptions made in the derivation in section~\ref{sec:Neff_analytical}. Note that the value of $\Delta N_{\rm eff}$ also approaches zero for $f_{\rm dcdm}\to 0$ which again leads to the $\Lambda$CDM model. We also note that the ratio between the analytical and numerical results in figure~\ref{fig:dNeff_f} does not approach unity exactly, but rather a slightly higher value, since the decay rate is fixed at a finite value ($\Gamma_{\rm dcdm}=10^{7}\,{\rm Gyr}^{-1}$) and the ratio only approaches unity for $\Gamma_{\rm dcdm}\to\infty$, according to figure~\ref{fig:dNeff_Gamma}.

\begin{figure}[t]
\centering
\subfloat[\label{fig:dNeff_Gamma}]{\hspace{.57\textwidth}}
\subfloat[\label{fig:dNeff_f}]{\hspace{.34\textwidth}}\\
\includegraphics[width=\textwidth]{Figures/dN_eff_full.pdf}
\caption{$\Delta N_{\rm eff}$ as a function of different parameters for both the numerical result calculated using eq.~\eqref{eq:dN_eff_class} and the analytical result calculated using eq.~\eqref{eq:neff2}. The lower panels show the ratio between the numerical and analytical $\Delta N_{\rm eff}$ as a function of the same parameters. The additional parameters used for the simulations are found in table~\ref{tab:class_input}. (a) $\Delta N_{\rm eff}$ as a function of the decay rate, $\Gamma_{\rm dcdm}$, with the fractional amount of DCDM fixed to $f_{\rm dcdm}=0.3$. (b) $\Delta N_{\rm eff}$ as a function of the fractional amount of DCDM, $f_{\rm dcdm}$, with the decay rate fixed to $\Gamma_{\rm dcdm}=10^{7}[{\rm Gyr^{-1}}]$.\label{fig:dN_eff}}
\end{figure}

\begin{figure}[t]
    \centering
    \includegraphics{Figures/dN_eff_fig.pdf}
    \caption{CMB TT power spectra for different values of $\Gamma_{\rm dcdm}$ and the corresponding $\Delta N^{\rm numerical}_{\rm eff}$ (calculated used eq.~\eqref{eq:dN_eff_class}) relative to the power spectrum of the $\Lambda$CDM model with no correction to $N_{\rm eff}$. The additional parameters used for the simulations are found in table~\ref{tab:class_input}, and the fractional amount of DCDM is fixed to $f_{\rm dcdm}=0.2$.}
    \label{fig:CMB}
\end{figure}

We can now run \textsc{class}\, again without DCDM, but with the corresponding $\Delta N^{\rm (numerical)}_{\rm eff}$ from the DCDM computation instead. According to the theory, we would expect the cosmology of such a run to be similar to that of very short-lived DCDM. Figure~\ref{fig:CMB} shows the CMB TT power spectra of the simulations with DCDM and $\Delta N_{\rm eff}=\Delta N^{\rm (numerical)}_{\rm eff}$ for two values of the decay rate, $\Gamma_{\rm dcdm}\in\{10^5\,{\rm Gyr}^{-1},\; 10^8\,{\rm Gyr}^{-1}\}$, and a fraction of initial DCDM to all dark matter of $f_{\rm dcdm}=0.2$. From this figure it is very clear that the mapping to $\Delta N_{\rm eff}$ becomes more accurate for higher values of $\Gamma_{\rm dcdm}$, as predicted by the theory. For the lower value of $\Gamma_{\rm dcdm}=10^5\,{\rm Gyr}^{-1}$, we still see the same behaviour in the power spectra, but they do not overlap as well due to the decay happening much closer to recombination, so a more significant amount of DR production is still ongoing at the time of the primary CMB signal formation. We also note that the magnitude of the relative power spectra is decreasing for larger $\Gamma_{\rm dcdm}$ as we would also expect from the theory.

\section{Current constraints on decay parameters}\label{sec:MCMC}

In order to obtain constraints on decay parameters we have used the publicly available code \textsc{MontePython}~\cite{montepython}.
This Markov Chain Monte Carlo (MCMC) sampler uses the Metropolis-Hastings algorithm to sample the probability distribution assuming flat priors. The code is run with four or seven sample chains on a computer cluster until convergence of the chains. Our Gelman-Rubin criterion for convergence of the chains is $R-1\lesssim0.01$, which means that the largest $R-1$ value of any parameter should be around or smaller than $0.01$. The prior of the $\Gamma_{\mathrm{dcdm}}$ parameter is set with different upper and lower bounds depending on which regime we want to analyse. These bounds and the number of sample chains will be stated when necessary.

\subsection{Long-lived DCDM regime}
The long-lived regime with a decay rate in the interval $\Gamma_{\mathrm{dcdm}}\in [0,10^3]\,\mathrm{Gyr}^{-1}$ corresponds to anything in between the DCDM starting to decay around the time of recombination (and thus having fully decayed by now) and an infinite lifetime where none of the DCDM has decayed yet. The common thing here is that the decay has not finished before recombination which leads to an ongoing DR production after the emission of the CMB. The late decay will thus have an effect on e.g. the angular diameter distance and other background parameters.

Using \textsc{MontePython} we search the parameter space with four sample chains using the likelihoods of Planck-2015 and Planck-2018 to see the differences of the two data sets. The result of this is found in figure~\ref{fig:15-18-long}.
\begin{figure}[t]
    \centering
    \includegraphics{Figures/Planck2018-VS-Planck2015_long.pdf}
    \caption{Triangle plot of posteriors in the long-lived regime using the two Planck data sets from 2015 and 2018.}
    \label{fig:15-18-long}
\end{figure}
We have also tried to further constrain the parameters using a combination of Planck-2018 and Baryonic Acoustic Oscillation data (BAO) from BOSS DR12. These results are plotted along with that of only Planck-2018 in figure~\ref{fig:BAO-long}.
\begin{figure}[t]
    \centering
    \includegraphics{Figures/Planck-VS-BAO_long.pdf}
    \caption{Triangle plot of posteriors in the long-lived regime using Planck-2018 data with and without BAO data from BOSS DR12.}
    \label{fig:BAO-long}
\end{figure}

From figure~\ref{fig:15-18-long} we see that the Planck-2018 data constrains the amount of DCDM through the parameter $f_{\rm dcdm}$ much more than the Planck-2015 data does, and our results using Planck-2015 is in great agreement with those found in Ref.~\cite{Poulin_16}. The same is true regarding the constraints on the decay rate, $\Gamma_{\rm dcdm}$, which means that the Planck-2018 data favours very long-lived DCDM even more than Planck-2015 data does. The plateau for higher values of $\Gamma_{\rm dcdm}$ is due to the fact that the data cannot distinguish between the models corresponding to this plateau given that they all imply a very late decay of DCDM (much later than recombination). The plateau continues all the way to the short-lived regime, but we only show the lowest values here. The disagreement in the parameter $\omega^{\rm ini}_{\rm cdm} \equiv (\Omega^{\rm ini}_{\rm dcdm}+\Omega_{\rm cdm})h^2$ is mainly due to the disagreement in $\Omega_{\rm cdm}$ of the two Planck data sets~\cite{Planck2015,Planck2018}, since the very long-lived DCDM behaves like a stable component and we thus would not expect a significant change to the total dark matter component from that of the $\Lambda$CDM model.

In figure~\ref{fig:BAO-long} we see that the addition of BAO data relaxes the amount of DCDM through $f_{\rm dcdm}$ a little, and we see a slight change in the posterior of the decay rate as well, raising the plateau for higher values while faintly narrowing the peak towards lower values. The parameter $\omega^{\rm ini}_{\rm dcdm}$ is lowered by BAO as well,
and this is in agreement with a narrower peak in the decay rate when including BAO, i.e. if the decay rate is smaller, the lifetime is longer, and we therefore need a smaller amount of initial dark matter since less of it has time to decay before the present.

The relevant best-fit parameters and constraints in the long-lived regime can be seen in the middle panel of table~\ref{tab:constraints} where we see a much tighter constraint on $f_{\rm dcdm}$ inferred by data from Planck-2018 as opposed to data from Planck-2015, and the same is true for the parameter $\Gamma f_{\rm dcdm}$. We also notice that the inclusion of BAO data loosens the constraints of the parameters while slightly modifying the best-fit values of $H_0$ and $\Omega^{\rm ini}_{\rm cdm}$ as well.

\subsubsection{Very long-lived DCDM regime}
In order to compare with the previous results of Ref.~\cite{Poulin_16}, we have also investigated a \textit{very long-lived} sub-regime with $\Gamma_{\mathrm{dcdm}}\lesssim H_0$. This yields slightly different results, since the long-lived runs are strongly undersampled in the very long-lived limit, but the tendencies are the same with tighter constraints from the newest Planck-2018 data. The inclusion of BAO data actually further tightens the constraint of the parameter $\Gamma f_{\rm dcdm}$ even though this was opposite in the long-lived regime.

We find that this regime allows for a much larger fractional amount of DCDM ($f_{\rm dcdm}\to 1$), which is clear from figure~\ref{fig:very_longlived}, since the lifetime is larger than the age of the Universe and the fractional amount thus has only little impact on current observables and no impact on early-time observables such as the CMB. Constraints from this regime can be found in the upper panel of table~\ref{tab:constraints}. We notice that we are only able to constrain the parameter $f_{\rm dcdm}$ when including the BAO data, and even then it is only possible to $1\sigma$. This is due to the degenerate nature of the parameter in this regime, where all values are allowed by the data. The degeneracy is lifted when including BAO data since late-time measurements can rule out too high values of $f_{\rm dcdm}$, given that these would feel the effects of the decay as opposed to the early-time measurements. This is also apparent from figure~\ref{fig:very_longlived}.

\begin{figure}[t]
    \centering
    \includegraphics{Figures/1D_posteriors_very_longlived.pdf}
    \caption{1D posteriors of DCDM parameters in the very long-lived sub-regime using Planck-2015 data along with Planck-2018 data with and without BAO data.}
    \label{fig:very_longlived}
\end{figure}

\subsection{Short-lived DCDM regime}
In the short-lived regime, the decay rate is rather high ($\Gamma_{\mathrm{dcdm}}>10^3\,[\mathrm{Gyr}^{-1}]$) resulting in most of the DCDM decaying well before recombination. We now search the parameter space using seven sample chains, and in order to search more efficiently, we implement the logarithm of the decay rate, $\mathrm{log}_{10}(\Gamma_{\mathrm{dcdm}})$, as a parameter in \textsc{class}, which allows us to use this parameter in \textsc{MontePython} as well. The parameter space is effectively reduced in size which makes the sampling much faster. We cut the prior at the upper bound $\Gamma_{\mathrm{dcdm}}<10^6\,[\mathrm{Gyr}^{-1}]$ (similar to Ref.~\cite{Poulin_16}) for convergence reasons, but we have checked that the posterior for the decay rate flattens out for increasing values, thus creating another plateau where the data cannot distinguish between the models. To further increase the efficiency we use the flag \texttt{-T=2.0} (default \texttt{1.0}) when running the code, which corresponds to increasing the statistical temperature of the chains so they are more likely to jump further in the parameter space.
\begin{figure}[t]
    \centering
    \includegraphics{Figures/Planck2018-VS-Planck2015_short.pdf}
    \caption{Triangle plot of posteriors in the short-lived regime using the two Planck data sets from 2015 and 2018.}
    \label{fig:15-18-short}
\end{figure}

The figures~\ref{fig:15-18-short} and \ref{fig:BAO-short} show the posterior distributions in the short-lived regime of Planck-2018 data and either Planck-2015 data or including BAO data, respectively. The 1D-posteriors of the parameter $\log_{10}(\Gamma_{\rm dcdm}[{\rm Gyr}^{-1}])$ are not shown in their entirety, but rather a closeup of the interesting region is presented by only showing the vertical interval $[0,0.2]$ (the posteriors are normalised so highest value equals unity).

\begin{figure}[t]
    \centering
    \includegraphics{Figures/Planck-VS-BAO_short.pdf}
    \caption{Triangle plot of posteriors in the short-lived regime using Planck-2018 data with and without BAO data from BOSS DR12.}
    \label{fig:BAO-short}
\end{figure}

Figure~\ref{fig:15-18-short} shows that Planck-2018 data allows for a slightly higher initial dark matter density, $\omega^{\rm ini}_{\rm cdm}$, than Planck-2015 data does, which is apparent from the broader peak and longer exponential tail towards higher values in the 1D-posterior. The amount of DCDM allowed by Planck-2018 is also slightly higher which is apparent from the 1D-posterior of $f_{\rm dcdm}$ being slightly broader than that of Planck-2015. The posterior of the decay rate consists of two sub-regimes as explained in Ref.~\cite{Poulin_16}, where the small peak (plateau) in the 1D-posterior of Planck-2015 (2018) represents a decay happening mostly in between matter-radiation equality and recombination while the increase in the posterior at larger values of the decay rate represents a decay taking place mostly before matter-radiation equality. Our 1D-posterior of the decay rate for the Planck-2015 data is reminiscent of the one presented in Ref.~\cite{Poulin_16}, but with a much less significant peak even though the same data has been used. The difference here could possibly be due to the Gelman-Rubin convergence criterion of the MCMC sampling, where our criterion is an order of magnitude lower than the one used in Ref.~\cite{Poulin_16}. The conclusion, however, does not change by letting the chains converge further, but we do get that the region between the small peak and the subsequent increase is populated more frequently with a stronger convergence criterion, thus diminishing the profoundness of the peak and making it more like a plateau. Using Planck-2018 data we completely transform the peak into a plateau instead.

\begin{table}[t]
\centering
\begin{tabular}{? Sc | Sc  Sc  Sc  Sc ?}
\specialrule{.12em}{0em}{0em}
\bf Very long-lived &
${\Omega^{\rm ini}_{\rm cdm} h^2\times 10^2}$ &
${H_0}$ &
${f_{\rm dcdm} \times 10^2}$  & $\Gamma f_{\rm dcdm}\times 10^{3}$\\
\bf regime &
& $[{\rm km\,s^{-1}Mpc^{-1}}]$
&
& $[{\rm Gyr^{-1}}]$\\
\hline
Planck-2015  & $11.9^{+0.3}_{-0.3}$ &  $67.7^{+1.3}_{-1.2}$  & \,\,---  & $<7.12$\\
Planck-2018  & $11.9^{+0.2}_{-0.2}$  &  $67.5^{+1.2}_{-1.2}$  & \,\,---  & $<4.01$\\
Planck-2018 + BAO  & $11.9^{+0.2}_{-0.2}$ &  $67.6^{+0.9}_{-0.9}$   & $<8.05^*$  & $<3.72$ \\
\specialrule{.12em}{0em}{0em}
\bf Long-lived regime &
${\Omega^{\rm ini}_{\rm cdm} h^2\times 10^2}$ &
${H_0}$ &
${f_{\rm dcdm} \times 10^2}$  & $\Gamma f_{\rm dcdm}\times 10^{2}$\\
&
& $[{\rm km\,s^{-1}Mpc^{-1}}]$
&
& $[{\rm Gyr^{-1}}]$\\
\hline
Planck-2015  & $12.0^{+0.3}_{-0.3}$ &  $67.3^{+1.6}_{-1.5}$  & $<4.14$   & $<10.99$\\
Planck-2018  & $12.1^{+0.3}_{-0.3}$  &  $67.1^{+1.2}_{-1.2}$  & $<2.44$  & $<\,\,\,5.18$\\
Planck-2018 + BAO  & $11.9^{+0.2}_{-0.2}$ &  $67.7^{+1.0}_{-0.9}$   & $<2.62$  & $<\,\,\,5.84$ \\
\specialrule{.12em}{0em}{0em}
\bf Short-lived regime &
${\Omega^{\rm ini}_{\rm cdm} h^2\times 10^2}$ &
${H_0}$ &
${f_{\rm dcdm}}\times 10^{1}$ &  $\Gamma f_{\rm dcdm}\times 10^{-4}$\\
&
& $[{\rm km\,s^{-1}Mpc^{-1}}]$
&
& $[{\rm Gyr^{-1}}]$\\
\hline
Planck-2015  & $12.6^{+1.5}_{-0.5}$ &  $67.7^{+1.6}_{-1.6}$  & $<0.98$ & $<2.35$\\
Planck-2018  & $12.7^{+1.6}_{-0.8}$  &  $67.8^{+1.4}_{-1.5}$  & $<1.31$ & $<3.01$\\
Planck-2018 + BAO  & $13.1^{+1.8}_{-1.0}$ &  $68.6^{+1.2}_{-1.4}$  & $<1.49$ & $<3.78$ \\
\specialrule{.12em}{0em}{0em}
\end{tabular}
\caption{Table of parameter constraints in the long-lived and short-lived regimes as well as in the very long-lived sub-regime. We present 2$\sigma$ constraints corresponding to a 95\% confidence level. The $*$ refers to 1$\sigma$ constraints only.}
\label{tab:constraints}
\end{table}

From figure~\ref{fig:BAO-short} it is clear that the inclusion of BAO data does not have a significant impact in the short-lived regime, as we previously argued. We see almost the exact same features in the posteriors with the exception that the 1D-posteriors of the initial dark matter density and the fractional amount of DCDM seem broader when including BAO data as opposed to just using Planck-2018 data. This is, however, not significant and we should not draw any conclusions of it, since it could very well be due to our chains not being converged enough. Increasing the amount of sample points or chains could possibly make the difference between the posteriors vanish.

The lower panel of table~\ref{tab:constraints} shows the relevant best-fit values and constraints in the short-lived regime. Here we see more loose constraints of the DCDM parameters from Planck-2018 data than from Planck-2015 data, and even more so when including BAO data. The DCDM parameter constraints are, however, not affected as much in the short-lived regime as in the very long-lived limit, where the inclusion of BAO data lifts the degeneracy of the $f_{\rm dcdm}$ parameter, and this agrees with the idea that BAO data only has little impact in short-lived regime.

We note that at first it might seem somewhat counter-intuitive that the constraint on $f_{\rm dcdm}$ is actually loosened by the inclusion of more data, rather than strengthened. The reason for this shift is, however, easy to understand from figure~\ref{fig:BAO-short}. Given the strong correlation between $f_{\rm dcdm}$ and $\Gamma_{\rm dcdm}$, the shift towards higher values of $\Gamma_{\rm dcdm}$ enforced by the Planck-2018 (and BAO) data removes a large fraction of the low $f_{\rm dcdm}$ parameter space allowed by the Planck-2015 data. This automatically shifts the preferred range of $f_{\rm dcdm}$ upwards and leads to a less restrictive upper bound on this parameter.

\section{Impact on the Hubble and $\sigma_8$ tensions}\label{sec:tensions}

As mentioned in section~\ref{sec:previous}, the idea of decaying dark matter has been able to relieve both the Hubble tension and the $\sigma_8$ tension in matter fluctuations to some extent.

\begin{figure}[t]
    \centering
    \includegraphics{Figures/H0_full_large.pdf}
    \caption{1D-posteriors of $H_0$ using all the data set combinations in both regimes. A posterior for the $\Lambda$CDM model using the newest Planck-2018 data has been included for reference. The shaded areas represent the latest confidence level of measurements in the local Universe of $73.2\pm1.3\,{\rm km\,s^{-1}\,Mpc^{-1}}$ from Ref.~\cite{Riess_2020}.}
    \label{fig:H0}
\end{figure}

From the 1D-posteriors of figure~\ref{fig:H0} we can clearly see that our model is indeed capable of relieving the Hubble tension to some extent using the various data sets. The best result in this regard is obtained using both Planck-2018 data and BAO data. This seems to be the case in both regimes while the short-lived regime is more successful in relieving the tension, which is also apparent from table~\ref{tab:constraints}. Compared to standard $\Lambda$CDM, our best case scenario is, however, only able to relieve the tension with $\sim\!\!1\sigma$, so the simple DCDM model does not offer the solution to this discrepancy. Perhaps one might have expected that the short-lived case would offer a better fit with high values of $H_0$. It is well known that using only CMB temperature data, there is a very strong positive correlation between $N_{\rm eff}$ and $H_0$, and since the short-lived case can be almost exactly mapped to a model with increased $N_{\rm eff}$ (see section~\ref{sec:Neff_mapping}), it is perhaps natural to expect that a high $H_0$ can be accommodated. However, the high $(H_0,N_{\rm eff})$ region is no longer allowed when polarisation data is added, thus disallowing this possibility, and in the end the preferred range of $H_0$ is only shifted marginally towards higher values in this case.

\begin{figure}[t]
    \centering
    \includegraphics{Figures/sigma8_omegaM_full.pdf}
    \caption{2D-posteriors of $\Omega_M$ and $\sigma_8$ using all the data set combinations in both regimes. A posterior from the $\Lambda$CDM model using the newest Planck-2018 data has been included for reference (dotted contours) as well as a posterior from DES Year 1 (purple contour) combining galaxy-galaxy clustering, galaxy-galaxy lensing and cosmic shear (Ref.~\cite{DES_Y1}).}
    \label{fig:sigma8}
\end{figure}

The $\sigma_8$ tension is also not solved with this model, which is apparent from the contours of figure~\ref{fig:sigma8}. This figure includes the contours of the combined DES Year 1 data from local measurements of the matter fluctuations using galaxy clustering and lensing. We can see that the contours of the long-lived regime move towards that of DES, but the tension is only relieved slightly with the 2$\sigma$ contours of the long-lived results being within the 1$\sigma$ contour of DES and vice versa. In the short-lived regime, all contours remain outside the 1$\sigma$ contour of DES, as in the case of standard $\Lambda$CDM, and they only move tangentially along the DES contour in comparison to the $\Lambda$CDM contour, thus yielding no better (nor worse) result.

Since the $\sigma_8$ tension can be relieved (slightly) only in the long-lived regime and since the Hubble tension can be relieved best in the short-lived, we recover a problem stated in Ref.~\cite{Pandey_2020} saying that attempts to relieve one of the tensions often worsens the other. We, however, do not get a worse result for the other when attempting to relieve one of the tension, but since the two tensions are relieved in different decay regimes, the simple DCDM model cannot be a solution to both tensions at once.

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

Using the MCMC sampler \textsc{MontePython} and the most recent CMB data from Planck-2018, we have improved on the previous constraints on the fractional amount of DCDM, $f_{\rm dcdm}$, as well as the parameter $\Gamma f_{\rm dcdm}$ in the two regimes (long-lived and short-lived)  along with the very long-lived sub-regime. In the very long-lived sub-regime, we find that the fractional amount of DCDM cannot be constrained by Planck-2018 data alone due to the degenerate nature of the parameter in this regime. We can however get a constraint for $\Gamma f_{\rm dcdm}$ for which we find an upper bound at 2$\sigma$ of $\Gamma f_{\rm dcdm}<4.01\times10^{-3}\,{\rm Gyr^{-1}}$. In the long-lived regime, we find that the fractional amount of DCDM is much better constrained by Planck-2018 data (as opposed to Planck-2015 data) leading to an upper bound at 2$\sigma$ of $f_{\rm dcdm}<2.44\%$, and the same is true for $\Gamma f_{\rm dcdm}$ for which we find an upper bound at 2$\sigma$ of $\Gamma f_{\rm dcdm}<5.18\times10^{-2}\,{\rm Gyr^{-1}}$. The constraints in the long-lived and very long-lived regimes are thus all tighter than the constraints inferred by Planck-2015 data by a factor of $\sim\!\!2$. For the short-lived regime, the Planck-2018 data leads to more loose constraints than Planck-2015 data does, and we thus find 2$\sigma$ upper bounds of $f_{\rm dcdm}<13.1\%$ and $\Gamma f_{\rm dcdm}<3.01\times10^{4}\,{\rm Gyr^{-1}}$, which are both higher than those of Planck-2015 by a factor of $\sim\!\!1.3$.

To further constrain the decay parameters, we have included BAO data from BOSS DR-12 as well. This leads to an even tighter constraint in the very long-lived sub-regime for $\Gamma f_{\rm dcdm}$ with a 2$\sigma$ upper bound of $\Gamma f_{\rm dcdm}<3.72\times10^{-3}\,{\rm Gyr^{-1}}$, while it is now also possible to constrain $f_{\rm dcdm}$ due to the BAO data lifting the degeneracy. We then get a 1$\sigma$ upper bound of $f_{\rm dcdm}<8.05\%$. In both the long-lived and the short-lived regimes we, however, get looser constraints. In the long-lived regime we get 2$\sigma$ upper bounds of $f_{\rm dcdm}<2.62\%$ and $\Gamma f_{\rm dcdm}<5.84\times10^{-2}\,{\rm Gyr^{-1}}$, while we in the short-lived regime get $f_{\rm dcdm}<14.9\%$ and $\Gamma f_{\rm dcdm}<3.78\times10^{4}\,{\rm Gyr^{-1}}$.

The impact on the Hubble and $\sigma_8$ tensions has also been investigated by comparing our posterior distributions from the MCMC runs with the data from measurements in the local Universe, i.e. the latest value of the Hubble parameter inferred by the distance ladder in Ref.~\cite{Riess_2020} and the combined data of DES Y1 (Ref.~\cite{DES_Y1}) measuring the matter fluctuations using galaxy clustering and lensing. We get a smaller discrepancy in the Hubble parameter in the short-lived regime, reducing the tension by $\sim\!\! 1\sigma$, but in this same regime we get no improvement of the $\sigma_8$ tension. The case is opposite for the long-lived regime, where we are able to relieve the $\sigma_8$ tension slightly, but not as much in the Hubble tension. We must therefore conclude that the simple DCDM model cannot accommodate both tensions at once.

We have in addition to this investigated how the short-lived DCDM (decaying much earlier than matter-radiation equality) is analogous to a universe without DCDM but with a larger initial amount of non-EM radiation (e.g. massless neutrinos). This has been analytically mapped to a correction to the effective number of massless neutrinos species, $N_{\rm eff}$, in eq.~\eqref{eq:neff2}. Using the Boltzmann code \textsc{class}\, we also calculated the actual correction to $N_{\rm eff}$ and saw that this scales in the exact same way as our analytical expression. The analytical and numerical results are also shown to agree with increasing precision as $\Gamma_{\rm dcdm}\to \infty$, where we recover the standard $\Lambda$CDM model.
\\

\noindent\textbf{Reproducibility.} The modified version of \textsc{class}\, used to obtain the results in this paper is available at \url{https://github.com/AarhusCosmology/CLASSpp_public/} as branch 2011.01632 with SHA 767fcdec52f8135dd8cebfcba7e1b2b3cdc7bc6a. The version of \textsc{MontePython} used as well as parameter files and scripts are available at \url{https://github.com/AarhusCosmology/montepython_public/} as branch 2011.01632 with SHA e2a8af41725ce31d 63718eafe7ec614801b291f6.

\section*{Acknowledgements}
The numerical results presented in this work were obtained at the Centre for Scientific Computing, Aarhus \url{http://phys.au.dk/forskning/cscaa/}. A.N. and T.T. was supported by a research grant (29337) from VILLUM FONDEN.

\bibliographystyle{utcaps}