\section{Introduction}

Dark matter (DM) makes up more than $80 \%$ of the matter content of our universe \cite{Planck:2015xua}. Despite the overwhelming quantitative evidence for cold dark matter (CDM) from its gravitational interactions, and its inclusion in the standard $\Lambda$CDM cosmological model, we only have very limited information about possible non-gravitational interactions of DM.
Existing upper bounds on interactions between the DM and the Standard Model (SM) and self-interactions of the DM leave a broad range of possibilities for models of DM. Viable models may include a whole dark sector consisting of different species of particles which interact only with themselves and with DM~(see e.g. \cite{Feng:2008mu,Feng:2009mn,Franca:2013zxa}). Furthermore, there is some tension
between the best fit parameters in $\Lambda$CDM as inferred from the cosmic microwave background (CMB)~\cite{Planck:2015xua}, and direct measurements of $H_0$ using standard candles and $\sigma_8$  using galaxy clusters~\cite{Riess:2011yx,Beutler:2014yhv,Battye:2014qga}). While it is premature to draw definite conclusions about these discrepancies, they motivate the study of interacting dark sectors which go beyond the ordinary {cold dark matter} of $\Lambda$CDM.

In this work, we propose a novel scenario for an interacting dark sector which predicts several important modifications to standard CDM phenomenology. Our dark sector includes a weakly interacting massive particle (WIMP) as the dark matter. As usual, weak interactions between the DM and the SM keep the two sectors in equilibrium at high temperatures and set the abundance of DM through thermal freeze-out. The non-trivial new ingredient of our scenario is that the DM particles are charged under a new non-Abelian ``dark" gauge group.\footnote{%
New visible particles charged under a new non-Abelian gauge group with macroscopic sized confinement scale have been proposed in the past under the name of theta-particles~\cite{Okun:1980mu,Okun:1980kw,Khlopov:1980ar} and more recently as quirks~\cite{Kang:2008ea}.
} Thus our dark sector consists of a degenerate multiplet of massive DM particles coupled to ``dark radiation" in the form of massless non-Abelian gauge bosons. Observations of large scale structure tightly constrain long-range interactions of DM with dark radiation and require very small gauge couplings\footnote{
Such couplings are too small to have appreciable effects on the Bullet Cluster, the shape of galactic dark matter halos or DM protohalos. They are also too small to have any bearing on the ``core vs cusp" or ``missing satellites" problems.},
$g_d < 10^{-3}$, so that the ``dark" confinement temperature is well below the current  CMB temperature.\footnote{
For recent discussions of non-Abelian dark sectors in a different context see~\cite{Feng:2011ik,Jeong:2013eza,Baek:2013dwa,Boddy:2014yra,Blinov:2014nla,Yamanaka:2014pva}.}
Thus confinement in the dark sector is irrelevant, and the ``dark gluons" are a thermal bath of weakly coupled massless particles which are well described as a perfect fluid. Note that this is distinct from neutrinos which free-stream and are therefore not a perfect fluid. The two types of dark radiation can be distinguished through their imprints on the CMB~\cite{Friedland:2007vv,Diamanti:2012tg}.

The general features of our dark sector lead to several distinct observable effects which can be broadly classified as due to ({\it i.}) multiplicity of the dark matter, ({\it ii.}) dark radiation in the form of the dark gluon fluid, and ({\it iii.}) interactions between DM and dark radiation. We briefly summarize the salient features of each of the three effects here. Details will be provided in subsequent sections.

\begin{itemize}

\item {\it DM multiplicity:} Being charged under an unbroken non-Abelian gauge group, DM particles come in degenerate multiplets of $N$ dark colors. This multiplicity leads to several easy-to-understand but important differences from conventional WIMPs. First off, DM annihilation cross sections (into SM particles) are suppressed by $1/N$ from averaging over initial states. Assuming thermal freeze-out, we predict that our dark matter mass is smaller by a factor of $1/\sqrt{N}$ than in models without dark matter multiplicity. Similarly, indirect detection bounds from dark matter annihilation into photons near the galactic center are relaxed because of the $1/N$ in the annihilation cross section.  Second, DM pair production cross sections at colliders are enhanced by a factor of $N$ from the multiplicity of possible final states. Third, direct detection bounds for DM are unchanged because the scattering cross section for DM particles off nuclei does not contain multiplicity factors. We discuss these effects, which would also follow from a non-Abelian global symmetry, in Section 3.

\item {\it Dark radiation:} Through their interactions with the DM, dark gluons come into thermal equilibrium with the SM in the early universe. After freeze-out of the DM the dark gluons decouple from the SM. They maintain a thermal distribution but end up cooling relative to the CMB  because the photons are heated by absorbing the entropy contained in massive SM particles. The dark gluons are observable through their contribution to the radiation density in the universe, an effect which is conventionally expressed in terms of an effective number of neutrino species~(see \cite{Brust:2013ova} for a recent discussion of bounds on new relativistic degrees of freedom from the CMB). We find $\triangle N_{eff} \sim 0.07 (N^2-1)$ for the $N^2-1$ gluons of an \ensuremath{SU(N)_d \ } gauge group. Bounds on $N_{eff}$ from Planck and from nucleosynthesis give $N<4$ at 95\% confidence level. An increase in the radiation density also shifts the best fit values of other cosmological parameters from the CMB. The two most significant shifts are increases in the predicted values of $H_0$ and $\sigma_8$. An increase of $H_0$ would remove the tension between CMB fits and direct measurements of $H_0$. An increase in $\sigma_8$ worsens the tension between CMB data and direct observations of $\sigma_8$ from large scale structure. However, as discussed in the next item, interactions between DM and dark gluons can remove this tension. Finally, the fact that dark gluons are self-interacting allows one to distinguish them from additional neutrino species. Dark gluons are well described by a perfect fluid with no viscosity and therefore lead to less damping of the CMB power spectrum in comparison to extra free-streaming relativistic fluids. Future CMB data will be precise enough to determine the viscosity of any significant new component of dark radiation and can therefore distinguish between additional neutrinos and dark gluons~\cite{Diamanti:2012tg}. We elaborate on the phenomenology of our dark radiation in Section 4.

\item {\it DM - dark radiation interactions:} Scattering of dark radiation (DR) off DM introduces a drag force between  the non-relativistic DM fluid and the relativistic radiation. This drag suppresses gravitational clustering and is therefore observable in the matter power spectrum. Allowing for matter perturbations to grow at least approximately as in $\Lambda$CDM gives a conservative upper bound of $\alpha_d \simlt 10^{-8}$. What we find particularly interesting about our scenario is that the momentum transfer cross section for DM - DR scales with temperature like the Hubble parameter during radiation domination. This means that the effect of the radiation remains equally important throughout radiation domination, leading to a smooth reduction of the power spectrum. This is in contrast with the more frequently studied case where the DM interactions freeze out at a critical temperature, leading to a sharp cutoff in the power spectrum at small scales~\cite{Hofmann:2001bi,Loeb:2005pm,Bertschinger:2006nq,Feng:2009mn}. We show that a very small and smooth reduction of the power spectrum may resolve the tension in the indirect determination of $\sigma_8$ from Planck data and direct measurements of the power spectrum at large scales~(see e.g. \cite{Planck:2015xua,Beutler:2014yhv,Battye:2014qga}).\footnote{%
We thank J. Lesgourgues for alerting us to the discrepancy between determinations of $\sigma_8$ from fits to the CMB and from direct measurements.}
We discuss DM-DR interactions in Section 5.

\end{itemize}

Most of the features associated with non-Abelian dark matter are expected to hold in a broad class of models irrespective of the specific interactions between the DM and the SM. However,  in the interest of concreteness and for clarity of presentation we focus on a specific example model which is described in Section 2. The DM in the model is a ``wino-like" \ensuremath{SU(2)_{weak}\ } triplet Dirac fermion which transforms as a fundamental under a dark \ensuremath{SU(N)_d \ } gauge group.

\section{The Model and DM thermal relic abundance}

For concreteness we focus on a specific realization of a non-Abelian dark matter. We take the dark matter particle to be a Dirac fermion in the $(1,3)_0$ representation of the Standard Model $(SU(3),SU(2))_{U(1)}$ gauge group and in the fundamental representation of dark \ensuremath{SU(N)_d \ }. It has an \ensuremath{SU(2)_{weak}\ } preserving Dirac mass $M_\chi$ and no additional interactions with the Standard Model besides gauge interactions. Electroweak symmetry breaking leads to a loop-induced mass splitting between the charged and neutral components of the \ensuremath{SU(2)_{weak}\ } triplet. In the limit $M_\chi \gg M_W$ the splitting is independent of $M_\chi$ and given by~\cite{Cheng:1998hc,Ibe:2012sx}
\begin{equation}
\delta M_\chi = M_{\chi^{\pm}} - M_{\chi^{0}} \approx 0.16 \, \text{GeV} \, .
\label{eq:mass-splitting}
\end{equation}
The neutral component $\chi^0$ is the lightest particle which carries fundamental \ensuremath{SU(N)_d \ } charge and is therefore stable, it is our DM candidate. From the perspective of SM interactions the dark \ensuremath{SU(N)_d \ } is simply a global symmetry, so that our DM acts like $N$ identical copies of a weakly interacting particle $\chi^0$ with Dirac mass $M_\chi$. In the following, we determine $M_\chi$ by requiring that the correct DM abundance is obtained from thermal freeze-out.

At the time of DM chemical freeze-out $T\sim M_\chi/26 \gg \delta M_\chi$ so that we can ignore the mass splitting between the different components of $\chi^a$ when calculating relic abundances. Once the temperature drops below $\delta M_\chi$, the charged $\chi^\pm$ decay to $\chi^0$ plus standard model particles. Thus the (co-moving) number density of $\chi^0$ particles at low temperatures is simply given by summing over all components of $\chi^a$ at DM freeze-out, i.e.  $\sum_a n_{\chi^a} \rightarrow  n_{\chi^0}$. This explains why we can use the abundance calculation for the whole triplet $\chi^a$ in order to find the DM relic abundance.

We will show in subsequent sections that Cosmology requires the dark gauge coupling to be much smaller than the SM gauge couplings. Therefore the relevant interactions for the thermal relic calculations are of the form $\chi \chi \rightarrow \text{SM SM}$ and independent of the dark gauge coupling. Figure~\ref{fig:annihilation}  shows some of the relevant diagrams for DM annihilating to SM particles. In the limit $M_\chi^2 \gg M_W^2$, the thermally averaged effective annihilation cross-section is
\begin{equation}
\langle \sigma v \rangle = \frac{1}{2 N} \, \frac{37 g_2^4}{192 \pi M_\chi^2} \, ,
\label{eq:annihilation-xsec}
\end{equation}
where $g_2$ is the \ensuremath{SU(2)_{weak}\ } gauge coupling. This cross-section differs from the standard \ensuremath{SU(2)_{weak}\ } triplet ``wino"~\cite{Cirelli:2005uq} by the extra $1/2N$ factor, which comes from the multiplicity associated with the \ensuremath{SU(N)_d \ } representation and the fact that $\chi$ is a 4-component Dirac fermion instead of a 2-component Majorana fermion. The $1/2N$ factor can be easily understood from the fact that any given DM particle carries an dark color charge and can only annihilate if it finds an anti-particle with the corresponding anti-dark color, thus reducing the color averaged annihilation cross-section.

\begin{figure}[ht]
\centering
\includegraphics[width=0.8\textwidth]{dm-annihilation.pdf}
\caption{Annihilation of dark matter into \ensuremath{SU(2)_{weak}\ } gauge bosons or SM fermions. Since the mass splitting between members of \ensuremath{SU(2)_{weak}\ } multiplets is small compared with energy transfer in the annihilation diagrams (i.e. twice the $\chi$ mass) we compute the co-annihilation of full \ensuremath{SU(2)_{weak}\ } multiplets and ignore the mass splittings.}
\label{fig:annihilation}
\end{figure}

If the DM abundance is set by thermal freeze-out, then - to a good approximation - the DM mass density today depends on the DM mass only through its annihilation cross-section~\cite{Kolb:1990vq}. Therefore, holding the cross section in Eq.~\ref{eq:annihilation-xsec} fixed at the Cosmologically preferred value we see that the mass required to get the correct relic abundance decreases as the square root of $N$. In Table~\ref{table:masses} we give the mass of DM for different values of $N$ using the tree level annihilation cross-section. The masses are significantly lower than for the usual ``wino" case where the tree level formula predicts a mass of $2.4$ TeV.~\footnote{In the case of the ``wino", the preferred mass is close to a Sommerfeld resonance of the weak interactions and one must account for Sommerfeld enhancement to obtain $M_{wino} \simeq 3$ TeV ~\cite{Cirelli:2007xd}. In our case, the DM mass is safely below the Sommerfeld resonances and the perturbative annihilation formula \ref{eq:annihilation-xsec} is adequate.}

\begin{table}[ht]
\centering
\begin{tabular}{|c|c|c|c|c|}
\hline
$N=2$ & $N=3$ & $N=4$ & $N=5$ & Generic $N$ \\
\hline \hline
$1.2$ TeV & $1.0$ TeV & $0.9$ TeV & $0.8$ TeV & $\sim {2.4}/{\sqrt{2N}}$ TeV  \\
\hline
\end{tabular}
\caption{Dark Matter masses required to get the correct thermal abundance as a function of $N$.}
\label{table:masses}
\end{table}

\section{Dark Matter multiplicity and experimental searches}

As will be discussed in subsequent Sections, astrophysical and cosmological constraints require very small dark gauge couplings of order $\alpha_d < 10^{-8}$. This is too small to have observable effects for production of DM at colliders, or for direct or indirect detection of DM. However, it is important to take into account that DM particles are in multiplets of the dark \ensuremath{SU(N)_d \ }, i.e. there are multiplicity factors associated with DM processes. In the previous Section we already discussed that the predicted mass for our thermal relic DM is reduced by a multiplicity factor of $1/\sqrt{2N}$, which can lead to large changes in sensitivity of experimental searches. Similar multiplicity/color factors also appear in cross-sections associated with DM detection. Figure~\ref{fig:searches} shows the color factors associated with the different kinds of DM searches. In what follows we briefly describe the effects of \ensuremath{SU(N)_d \ } multiplicity for direct and indirect detection and also for collider searches for DM:

\begin{figure}[h!t]
\centering
\includegraphics[width=0.5\textwidth]{dm-searches.pdf}
\caption{Color factors for the different types of DM search experiments. The different multiplicity factors can be easily understood from the color flow in the figure. For direct detection the color of the incoming dark matter is the same as of the outgoing and so there is no multiplicity factor. For indirect detection it is an annihilation diagram, so just as for the thermal relic calculation there is a $1/2N$ suppression because a DM particle can only annihilate if it finds the anti-particle with the right anti-dark color. For colliders there is an $2N$ enhancement because any of the different N colors can be created and an extra 2 from Dirac vs Majorana.}
\label{fig:searches}
\end{figure}

\subsection{Direct detection}

 For direct detection there is no color factor associated with the multiplicity, thus the only change comes from dark matter being lighter. The spin-independent cross-section for dark matter scattering of the nucleus is approximately $10^{-47} \, \text{cm}^2$, and independent of the DM mass as long as $M_\chi \gg M_W$~\cite{Hill:2013hoa}. This cross-section is an order of magnitude smaller than the projected sensitivity of the next generation direct detection experiments~\cite{Cushman:2013zza}. However, in the mass range of interest (around $1 \, \text{TeV}$) it is above the neutrino background and potentially within reach of future experiments.

\subsection{Indirect detection}

The annihilation cross-section relevant for indirect detection is suppressed by a $1/2N$ factor. Taking into account Sommerfeld enhancement the cross-section is further reduced relative to the standard ``wino" model. This is because there is a resonant Sommerfeld enhancement from weak interactions for dark matter masses in the $2-3 \, \text{TeV}$ range, but the enhancement is much smaller for masses around $1 \, \text{TeV}$. The ``wino" model has been investigated recently~\cite{Cohen:2013ama,Fan:2013faa,Bauer:2014ula,Ovanesyan:2014fwa,Baumgart:2014saa}, and is strongly disfavored by H.E.S.S. data. Our \ensuremath{SU(N)_d \ } model is not yet constrained by either H.E.S.S. or Fermi data for any $N\ge 2$, but the annihilation cross-section is close to H.E.S.S. sensitivity as shown in Figure~\ref{fig:indirect}. It is worth noting that the limits shown in Figure~ \ref{fig:indirect} assume a specific NFW profile and there is a large uncertainty in these limits due to our limited knowledge of the dark matter distribution in the center of the galaxy. Also shown in the figure is the projected reach of CTA~\cite{Funk:2013gxa}, assuming 5 hours of observation time. One can see that CTA should have enough sensitivity to discover or rule out our model up to at least $N=10$, and therefore is the most promising search for discovering dark matter in our model.

\begin{figure}[h!t]
\centering
\includegraphics[width=0.6\textwidth]{dm-indirect-detection.pdf}
\caption{Indirect detection constraints from gamma ray line searches from H.E.S.S. (blue region) and projected sensitivity for gamma ray lines at CTA (gray region) assuming an NFW profile for the dark matter distribution in the galactic center, both taken from~\cite{Ovanesyan:2014fwa}. The red dots are the cross sections for dark matter annihilation into gamma rays in our models for $N=2$ up to $N=10$, which were obtained by appropriately rescaling the NLL cross section with Sommerfeld enhancement from~\cite{Ovanesyan:2014fwa}. For comparison we plotted in black the cross-section for the annihilation cross-section to photons in the standard ``wino'' model as a function of mass, also from~\cite{Ovanesyan:2014fwa}.}
\label{fig:indirect}
\end{figure}

\subsection{LHC and future collider searches}

The multiplicity factor enhances sensitivity of collider searches to our DM in two ways. Most importantly, the predicted DM mass from thermal freeze-out is lowered,  and thus a lower partonic center of mass energy is required to pair produce DM at a collider. For example,  at the 14 TeV LHC and for DM masses near $M_{DM} \sim$ 500 GeV the DM cross section scales as $(1/M_{DM})^6$ because of the strong energy dependence of the parton luminosities. In addition, the cross-section for pair producing dark matter is enhanced by the final state multiplicity factor $2N$.

\begin{figure}[h!t]
\centering
\includegraphics[width=0.5\textwidth]{LHC_reach--a-2_b-0.pdf}
\hskip-5pt
\includegraphics[width=0.5\textwidth]{FHC_reach--a-2_b-0.pdf}
\caption{Expected significance of missing energy (MET) searches for DM at the LHC and a future 100 TeV collider. The solid black lines in each plot correspond to the sensitivity of the collider to ``wino"-like DM, an \ensuremath{SU(2)_{weak}\ } triplet Majorana fermion. The colored dots labeled by different $N$-values correspond to our models in which the DM is a Dirac fermion with multiplicity $N$ and mass chosen to yield the correct abundance from thermal freeze-out.}
\label{fig:LHC}
\end{figure}

In the left panel of Figure~\ref{fig:LHC} we show the expected sensitivity of the high luminosity LHC to the DM in our model. The solid dots correspond to DM with multiplicity $N$ and Dirac masses chosen so that the correct thermal DM abundance is obtained. We see that the 14 TeV LHC is sensitive to DM with $N\ge7$. For comparison, we also show the expected sensitivity to a standard ``wino" \ensuremath{SU(2)_{weak}\ } triplet with a Majorana mass. Note that the ``wino" has the correct thermal abundance only for $M_{wino}\simeq 3$ TeV. Existing monojet searches from ATLAS~\cite{Aad:2015zva} and CMS~\cite{Khachatryan:2014rra} with 8 TeV collisions and a luminosity of 19.5 fb$^{-1}$ already rule out $N\simgt 20$.

In the right panel we show that a 100 TeV future collider can discover our DM above backgrounds even for the smallest non-Abelian multiplet of $N=2$ and perhaps may be able to rule out a Dirac ``wino". The significances for these plots were determined from parton level signal and background events which were computed with MadGraph \cite{Alwall:2011uj}. The main irreducible backgrounds are due to jets plus $Z$ or $W$ with MET from decays to neutrinos. We computed signal and backgrounds to leading order ($\alpha_s \alpha_W$) and assumed that the experiments will be able to limit background systematic uncertainties to 2\%. For the ``wino", our results are consistent with the more sophisticated studies in (\cite{Low:2014cba,Cirelli:2014dsa}). These references also showed that a MET plus ``disappearing track" search for the production of $\chi^\pm$ can improve sensitivity because it is free of irreducible SM backgrounds.

Another observable consequence of this model is a change in the running of the EW gauge coupling~\cite{Alves:2014cda}. The multiplicity of DM leads to an $2N$ enhancement factor in the DM contribution to running of $\alpha_W$ at one loop, which would be observable at the proposed 100 TeV hadron colliders.

\section{Dark gluons as dark radiation}

In this Section we turn our attention to the cosmic evolution of the energy density in dark gluons and its effects on the CMB. An important parameter which determines the effects of dark gluons on the CMB is the ratio of temperatures of the dark gluon plasma $T_d$ over the photon temperature $T$.

Dark gluons are coupled to the thermal bath of SM particles through their couplings to the DM which is in equilibrium with the SM in the early universe. If $\alpha_d$ is not too small the dark gluons also equilibrate with the SM plasma at early times, so that $T_d=T$ before DM freeze-out. To determine the smallest possible coupling $\alpha_d$ for which the dark gluons are in chemical and kinetic equilibrium with the SM we consider the Feynman diagram in Figure~\ref{fig:darkgluon-compton}. At temperatures higher than $M_\chi$, the thermally averaged cross-section for this process times the DM number density is given by
\begin{equation}
n_\chi \langle \sigma v \rangle \sim
T^3\ \frac{\pi \alpha_W \alpha_d}{T^2}\, .
\label{eq:darkgluon-equilibrium}
\end{equation}
Comparing this to the Hubble rate $H$ we find that the dark gluons are in equilibrium with the SM at $T \sim M_\chi$ for
\[
\alpha_d \gtrsim \frac{1}{\alpha_W} \, \frac{M_\chi}{M_\text{Planck}}.
\]

\begin{figure}[ht]
\centering
\includegraphics[width=0.5\textwidth]{dm-darkgluon-compton.pdf}
\caption{The process through which dark gluons maintain equilibrium with the DM and the SM plasma for small $\alpha_d$.}
\label{fig:darkgluon-compton}
\end{figure}

Thus even for $\alpha_d$ as small as $10^{-13}$ the dark gluons come to chemical and thermal equilibrium with the SM at temperatures of order $M_\chi$. However, when the universe cools below $T \sim M_\chi$ dark matter becomes non-relativistic and its number density drops exponentially. Then the rate for the process in Figure~\ref{fig:darkgluon-compton} becomes%
\footnote{Note that for temperatures below the W mass, the W boson in Figure~\ref{fig:darkgluon-compton} should be replaced by a photon and the $\alpha_W$ in Eq.~\ref{eq:gluonfreezeout} by $\alpha_{em}$.}
\begin{equation}
n_\chi \langle \sigma v \rangle \sim
 (M_\chi T)^{3/2} \, e^{-M_\chi/T} \ \ \frac{\pi \alpha_W \alpha_d}{M_\chi^2}\, ,
\label{eq:gluonfreezeout}
\end{equation}
and the dark gluons decouple at temperatures of about $M_\chi$ for $\alpha_d \sim 10^{-13}$ to about $M_\chi/20$ for $\alpha_d \sim 10^{-3}$. Below the decoupling temperature the dark radiation fluid evolves independently with a  temperature $T_d$ which redshifts as $1/a$.

The temperature of the photon fluid also redshifts as $1/a$ for most of the universe's evolution. However, when massive SM particles become non-relativistic they annihilate into the remaining lighter SM particles which effectively heats up the photons compared to the dark gluons (similarly to what happens to photons and neutrinos after neutrino decoupling~\cite{Kolb:1990vq}). The ratio between the photon and dark gluon temperatures can be easily calculated in the instantaneous decoupling approximation by requiring that the entropy per co-moving volume is conserved independently in each fluid,\begin{equation}
\frac{T_d}{T_\gamma} = \left( \frac{g_*^f}{g_*^i} \right)^{1/3},
\label{eq:temperature-ratio}
\end{equation}
where $g_*^i$ is the number of effective degrees of freedom in the SM plasma at the time of dark gluon decoupling and $g_*^f$ the number of effective degrees of freedom at any later time.

The CMB places strong constraints on the energy density in relativistic particles at the time of recombination. This constraint is usually presented in terms of the number of effective neutrino species, $N_{eff}$. The contribution of the dark gluons to $N_{eff}$ is given by
\begin{equation}
\Delta N_{eff} = \frac{8}{7} \, (N^2-1)  \left( T_d/T_\nu \right)^4,
\end{equation}
where the $N^2-1$ is the number of generators of \ensuremath{SU(N)_d \ } and $T_\nu$ is the neutrino temperature. The ratio can be calculated using Eq.~\ref{eq:temperature-ratio} right at neutrino decoupling when neutrino and photon temperatures are still the same. Assuming that the decoupling between dark gluons and the SM happens at temperatures around $50 \, \text{GeV}$ one finds
\[
\Delta N_{eff} = 0.07 \, (N^2-1).
\]

The strongest constrain on $N_{eff}$ comes from the 2015 Planck data~\cite{Planck:2015xua}, which found $N_{eff} = 3. 15 \pm 0.46$ at $95 \% $ confidence to be compared with the SM prediction~\cite{Mangano:2001iu} $N_{eff}^{SM}=3.046$. We see that this rules out $N\geq 4$ and that $N=3$ is within the $2 \sigma$ allowed range.
However, we note that the Planck analysis assumes the $\Lambda$CDM model with one additional parameter, $N_{eff}$. This limit could potentially be relaxed in our scenario where there are important differences as we will now discuss.

The effects of dark radiation can be divided into so-called background effects and perturbation effects. Background effects are due to a change in the average energy density in relativistic degrees of freedom and are not sensitive to any other properties of the dark radiation fluid. The largest background effect of extra radiation from relativistic degrees of freedom is to change the redshift of matter-radiation equality $z_{eq}$. Since $z_{eq}$ is very well measured a fit to the data is forced to maintain the redshift of matter-radiation equality by simultaneously increasing the dark matter density~(see e.g.~\cite{Hu:1998tk,Bowen:2001in,Bashinsky:2003tk}). This change in the matter density in turn requires a change in the Hubble parameter today $H_0$, in order to keep $\Omega_m$ fixed (the ratio between the matter density and the critical density). Thus we see that a fit to the CMB data alone has an approximate flat direction in which an increase in $N_{eff}$ can be compensated for by simultaneous increases in $\rho_{DM}$ and $H_0$.

The perturbation effects are due to perturbations in the dark radiation fluid and thus sensitive to properties of dark radiation. In particular there are two additional parameters which distinguish different types of dark radiation, see e.g.~\cite{Diamanti:2012tg}: the effective sound speed $c_\text{eff}^2$ and the viscosity speed $c_\text{vis}^2$. The dark gluons are a relativistic fluid and have $c_\text{eff}^2 = 1/3$, the same as for neutrinos. However---unlike neutrinos---the dark gluons have self-interactions which come from the non-Abelian gauge kinetic terms. If the rate of dark gluon-gluon scattering is large compared to the Hubble parameter then the dark gluons are well described as an ideal fluid instead of a free-streaming fluid as for neutrinos. For an ideal fluid one has $c_\text{vis}^2 = 0$ instead of $1/3$ as for neutrinos. The interaction rate between dark gluons is approximately given by
\begin{equation}
\tau^{-1} \sim \alpha_d^2 T_d,
\end{equation}
and one sees that as long as $\alpha_d \gtrsim 10^{-13}$ this rate is larger than $H$ during recombination. Thus on the time scale set by Hubble the dark gluons behave as a perfect fluid. Qualitatively, the interactions reduce the damping of overdensities from relativistic particles streaming out of gravitational potential wells. This has the effect that the CMB peaks are not as suppressed as they would be in the case of additional free-streaming dark radiation.

The Planck Collaboration has performed a fit for $c_\text{eff}^2$ and $c_\text{vis}^2$ with $N_{eff}$ fixed to the SM value $3.046$, i.e. no additional dark radiation, and found that the parameters were in perfect agreement with the expected value for neutrinos, $c_\text{eff}^2=c_\text{vis}^2=1/3$. However, Planck has not yet performed a fit for additional radiation $\triangle N_{eff}$ which is allowed to have non-standard values for $c_\text{eff}$ and $c_\text{vis}$. In particular, if future CMB experiments find evidence for a non-zero $\triangle N_{eff}$ then measuring $c_\text{vis}$ of this extra component would allow one to distinguish between dark gluons and free-streaming dark radiation, like dark photons or sterile neutrinos~\cite{Friedland:2007vv,Diamanti:2012tg}.

\section{Dark matter-dark gluon interactions and large scale structure}

In this Section we study how interactions between dark gluons and DM affect the linear evolution of DM overdensities. The interactions have two important effects. One is the transfer of kinetic energy (i.e. temperature) from the DR to the DM.
The other is that DM particles moving through the DR fluid experience a drag force.
This drag slows the growth of large scale structures through gravitational clustering. Both rates can be computed by considering scattering of dark gluons with DM particles.
The process is a generalization of Compton scattering to the non-Abelian case. The most important new feature is that scattering is dominated by the t-channel diagram shown in Figure~\ref{fig:t-channel} which is divergent for small angle scattering and which only exists for non-Abelian gauge bosons. As we will see, this interaction leads to interesting signatures in large scale structure which distinguish our scenario from other models with interactions between dark radiation (DR) and DM (see e.g. \cite{Diamanti:2012tg,Wilkinson:2014ksa}). In particular, for $\alpha_d \simlt 10^{-8.5}$ our  model predicts a smooth suppression of the matter power spectrum at all scales which could resolve the conflict between the indirect measurement of $\sigma_8$ from Planck \cite{Planck:2015xua} and the direct measurement from BOSS \cite{Beutler:2014yhv}.

\begin{figure}[ht]
\centering
\includegraphics[width=0.4\textwidth]{dm-tchannel.pdf}
\caption{t-channel scattering of dark matter with dark gluons.}
\label{fig:t-channel}
\end{figure}

To simplify the calculation of the energy transfer rate we consider the limit in which the DM temperature is negligible so that we can take the DM particles to be at rest. We then compute the rate of energy transfer to the DM from scattering~\cite{McDermott:2010pa,Dvorkin:2013cea} with a thermal bath of dark gluons. As in the well-known case of Coulomb scattering the cross section is dominated by small-angle forward scattering. To significantly impact the energy of a massive DM particle many collisions with the gluons are required. The collisions are uncorrelated so that the resulting momentum of the DM particle performs a random walk with
$$ E=\frac{p^2}{2M_\chi} \simeq \frac{N}{2M_\chi} (\delta p)^2 \simeq \frac1{2M_\chi} \sum (\delta p)^2\,. $$
Here $\delta p$ is the typical momentum transfer in a single collision and $N$ is the number of collisions. After many such random scatters the resulting DM particle distribution is thermal, but not necessarily with the temperature of the gluon bath. The DM temperature depends on the relative size of the Hubble expansion rate and the energy transfer rate.

The rate of energy transfer is calculated by averaging the energy transfer per collision over the initial Bose-Einstein distribution of the dark gluons
\begin{eqnarray}
\dot{E}=a\! \int\! \frac{d^3k}{(2\pi)^3}&& \!\!\!\!2 (N^2\!-\!1) f(k) \times  \\
&&\frac{1}{4 E_p k}\int\!  \frac{d^3k'}{(2\pi)^3 2k'} \frac{d^3p'}{(2\pi)^3 2E'_p}
(2\pi)^4 \delta(p\!+\!k\!-\!p'\!-\!k') |\overline{M}|^2 (E'_p\!-\!E_p) (1\!+\!f(k')) \nonumber \,.
\end{eqnarray}
Here---and for the remainder of this Section---time derivatives are taken with respect to conformal time which is the origin of the scale factor $a$ on the right hand side. Also, $f(k)=1/(\exp(k/T_{dr}-1)$ is the gluon thermal distribution function and it is multiplied by $2(N^2-1)$ for the spin and color of the initial gluon. The $1+f(k')$ final state factor accounts for stimulated emission. Finally, the color and spin summed and averaged matrix element (keeping only the t-channel) is given by
$$|\overline{M}|^2 = \frac12\ 4 g_d^4\, \frac{(s-M_\chi^2)(M_\chi^2-u)}{t^2} \ ,$$
where the $1/2$ is the disappointingly boring color factor of the t-channel diagram and $s$, $t$, and $u$ are the Mandelstam variables.
The integrals are straightforward to evaluate for $p=(M_\chi,\vec0)$ and give
\begin{equation}
\dot{E}= a (N^2\!-\!1)\, \frac{\pi}3\,  \alpha_d^2 \log(\alpha_d^{-1})\ \frac{T_{dr}^3}{M_\chi}=
a \frac5\pi \,  \alpha_d^2 \log(\alpha_d^{-1})\ \frac{\rho_{dr}}{T_{dr} M_\chi}\ .
\label{Eq:energy_rate}
\end{equation}
One subtlety one encounters is a logarithmic divergence at small $t$ which stems from the long range of the interaction mediated by a massless gluon. In the thermal plasma of gluons the interaction range is made finite by screening.  This effect can be parametrized by including a Debye mass for the gluon, $m^2_{Debye}\sim g^2_d T_{dr}^2$ \cite{Arnold:1995bh}. With the Debye mass the logarithmic divergence becomes the $\log(\alpha_d^{-1})$ in Eq.~\ref{Eq:energy_rate}. This is analogous to the ubiquitous ``Coulomb logarithm" in plasma physics.

Using $T=2/3\, E$ we obtain an equation for the evolution of dark matter temperature \cite{Ma:1995ey}
\begin{equation}
\dot T_\chi = -2 \,\frac{\dot{a}}{a}\, T_\chi +  a (N^2\!-\!1) \frac{2\pi}{9}\alpha_d^2 \log(\alpha_d^{-1}) \frac{T_{dr}^2}{M_\chi} (T_{dr}-T_\chi) \ ,
\end{equation}
where the $T_{dr}-T_\chi$ factor generalizes the energy transfer rate in Eq.~\ref{Eq:energy_rate} to the case in which the thermal motion of the DM is not negligible. To understand the possible solutions to this equation one compares the size of the temperature transfer term with the Hubble redshift term $\dot{a}/a\, T_\chi \sim a \, T_\chi T^2/M_{Planck}$ (during radiation domination). For large temperature transfer rates DM and DR are kept in thermal equilibrium by the interactions and $T_\chi=T_{dr}$. Note that both terms scale with the cube of the temperature, thus if DM and DR are in equilibrium at one temperature they will be in equilibrium throughout radiation domination.

For couplings $\alpha_d < 10^{-8}$ and assuming comparable starting temperatures $T_\chi \sim T_{dr}$, the Hubble term initially dominates over the interaction term and the DM temperature plummets relative to the temperature of radiation: $T_\chi \sim a^{-2}$ versus $T_{dr} \sim a^{-1}$. This is the regime where our DM behaves like ordinary cold DM (CDM). However, note that the Hubble term scales proportional to $T_{dr}^2 T_\chi$ whereas the collision term scales like $T_{dr}^3$. Thus for sufficiently small $T_\chi$ the two terms become comparable and the DM temperature switches to scaling proportional to $a^{-1}$, keeping the ratio $T_\chi/T_{dr}$ constant. In Figure~\ref{fig:temps} we show numerical solutions for the DM temperature as a function of scale factor for three representative values of $\alpha_d$.

\begin{figure}[ht]
\centering
\includegraphics[width=0.6\textwidth]{dm-temperature.pdf}
\caption{Dark sector temperatures as a function of scale factor $a$. Shown are the dark gluon temperature (black dashed), and DM temperatures for three representative values of $\alpha_d=10^{-8}$ (upper, blue), $\alpha_d=10^{-9}$ (middle, red), and $\alpha_d=10^{-10}$ (lower, green).}
\label{fig:temps}
\end{figure}

For values of $\alpha_d$ for which DM and DR are in thermal equilibrium DM behaves very differently from ordinary CDM. The pressure from the dark gluons prevents the growth of DM overdensities during radiation domination. This can be seen as the sharp drop in the DM power spectrum in Figure~\ref{fig:power-spectrum}. For smaller values of $\alpha_d$ the dark gluons still influence the evolution of DM overdensities, however the effects are more subtle and we employ the formalism of Ma and Bertschinger \cite{Ma:1995ey} to study them.

Following \cite{Ma:1995ey} we write down the linearized evolution equations for overdensities including the interactions between the DM and DR fluids in conformal Newtonian gauge.
To avoid the complications of solving a full Boltzmann code we work with a simplified scenario. We replace all relativistic energy density in the SM (i.e. neutrinos and photons) with an identical energy density which is made up of only photons. And we replace all matter (dark matter and baryons) with an equivalent energy density of only dark matter. In addition, we approximate by treating the photons as a perfect fluid (zero viscosity). This is only true before recombination but soon after recombination photons contribute only a negligible amount to the energy density and thus to the evolution of DM overdensities. Since our goal is to demonstrate the effect of the interactions on the DM density perturbations we compare our scenario with interactions to the same scenario with $\alpha_d=0$.

In Fourier space the equations for the DM and DR overdensities are
\begin{equation}
\begin{aligned}
& \dot \delta_{DM} = -\theta_{DM} + 3 \dot \psi \\
& \dot \theta_{DM}  = - \frac{\dot a}{a} \theta_{DM} + a \tau^{-1}_c (\theta_{DR} - \theta_{DM} ) + k^2 \psi \\
& \dot \delta_{DR}  = - \frac{4}{3} \theta_{DR} + 4 \dot \psi \\
& \dot \theta_{DR}  = k^2 \frac{\delta_{DR}}{4} + k^2 \psi + \frac{3}{4}\frac{\rho_{DM}}{\rho_{DR}} a \tau^{-1}_c (\theta_{DM} - \theta_{DR}) \\
& \dot \delta_\gamma  = - \frac{4}{3} \theta_\gamma + 4 \dot \psi \\
& \dot \theta_\gamma  = \frac{1}{4} k^2 \delta_\gamma + k^2 \psi \\
& k^2 \psi + 3 \frac{\dot a}{a} \left( \dot \psi + \frac{\dot a}{a} \, \psi \right) =  - \frac{a^2}{2 M_{Pl}^2} \sum_i \rho_i \delta_i \, ,
\end{aligned}
\end{equation}
where the dots represent derivatives with respect to conformal time, $\rho_X$ is the average energy density of fluid $X$ and $\delta_X$ and $\theta_X$ are related to the overdensity and velocity divergence in fluid $X$. We have also set the two metric perturbations $\psi$ and $\phi$ equal because we are treating the photons and dark radiation as ideal fluids (no anisotropic stress) and did not include neutrinos which have sizable anisotropic stresses. The interaction between dark matter and dark radiation is encoded in the momentum transfer rate $\tau^{-1}_c$~\cite{Dvorkin:2013cea}. It is defined as the change in momentum $\dot{\vec{p}}_\chi = - a \tau^{-1}_c\ \vec p_\chi$ which a DM particle with momentum $\vec{P}$ experiences due to friction as it is moving through the dark gluon fluid. Microscopically, the friction arises from collisions between DM particles and dark gluons and to compute it we evaluate
\begin{equation}
\dot{\vec{p}}=a\! \int\! \frac{d^3k}{(2\pi)^3}f(k) \frac{1}{4 E_p k}\int\!  \frac{d^3k'}{(2\pi)^3 2k'} \frac{d^3p'}{(2\pi)^3 2E'_p}
(2\pi)^4 \delta(p\!+\!k\!-\!p'\!-\!k') |M|^2 (\vec{p}{\,'}\!-\!\vec{p}\,) (1\!+\!f(k')) \ ,
\end{equation}
where now the initial DM momentum $\vec{p}$ is non-zero and we expand to first order in $p/M_\chi$. Employing the same approximations as for the energy transfer rate we obtain
\begin{equation}
\tau^{-1}_c = (N^2\!-\!1) \frac{\pi}{9} \alpha_d^2 \log\alpha_d^{-1} \, \frac{T_{dr}^2}{M_\chi}
\label{eq:velocity-rate}
\end{equation}

We integrate the equations for the overdensities from $a=10^{-7}$, when all modes of interest are well outside the horizon, until $a=10^{-3}$. We use initial conditions corresponding to adiabatic perturbations:
\begin{equation}
\begin{aligned}
\delta_\gamma = \delta_{DR} = \frac{4}{3} \delta_{DM} = - 2 \psi = C(k) \, ,
\end{aligned}
\label{eq:initial_conditions}
\end{equation}
where the initial perturbations $C(k) \sim k^{-3/2}$ are determined by the physics of inflation. We define the DM power spectrum equal to the square of the perturbations, $P(k)\equiv \, \delta_{DM}^2$, at scale factor $a=10^{-3}$. To focus only on effects of the coupling between DM and DR we form a ratio where we divide the power spectrum with interactions turned on by the power spectrum with $\alpha_d = 0$. Note that since the equations are linear the initial values for the perturbations, $C(k)$, drop out in the ratio.

\begin{figure}[ht]
\centering
\includegraphics[width=0.65\textwidth]{dm-power-spectrum.pdf}
\caption{Power spectrum including the DM-DR interactions normalized by the power spectrum with interactions turned off. The black dotted curve corresponds to $\alpha_d = 10^{-8}$, the green dashed curve corresponds to $\alpha_d = 10^{-8.5}$ and the red line corresponds to $\alpha_d = 10^{-9}$. The power spectra are defined proportional to $\delta_{DM}^2$ at $a = 10^{-3}$. The vertical yellow band labeled $k_{eq}$ indicates modes which enter the horizon at matter-radiation equality, modes which enter the horizon earlier are to the right (larger $k$). The blue band labeled $\sigma_8$ indicates modes which the observable $\sigma_8$ is most sensitive to.}
\label{fig:power-spectrum}
\end{figure}

The ratios of power spectra for different values of $\alpha_d$ are plotted in Figure~\ref{fig:power-spectrum}.
For the plot we chose the number of colors $N=2$ and $M_\chi = 1.2\,$TeV. One sees that for $\alpha_d = 10^{-8}$ the power spectrum is strongly suppressed for modes which entered the horizon before matter-radiation equality. These are modes with $k > k_{eq} \sim 0.015 \, \text{Mpc}^{-1}$. This should be expected because in this case the DM is in equilibrium with the DR bath throughout radiation domination. For the smaller values of $\alpha_d = 10^{-8.5}$ and $10^{-9}$ the power spectrum is less affected with modes which entered the horizon earlier (larger $k$) suppressed more that those which entered later. Modes which enter the horizon after matter radiation equality are not suppressed for any of the couplings plotted. The light blue vertical band indicates the range of modes which the observable $\sigma_8$ is sensitive to ($\sigma_8$ is a measurement of the matter fluctuations in spheres of radius of $8 h^{-1}$ Mpc).

The smooth suppression of power at all scales that we are finding is special to our scenario and stems from the fact that the momentum transfer rate scales with temperature as $T_{dr}^2$, the same scaling as Hubble. Thus it is possible to arrange for the couplings to have a small effect but act over a large range of scales. This should be contrasted with cases where the interactions scale like a higher power of $T_{dr}$ in which case they are important at high scales and have no effect at low scales. In such scenarios the power spectrum has a sharp cutoff at scales of the size of the horizon at the time when the interactions cease to be important, leaving larger scales unaffected and wiping out the smaller scales that entered the horizon at earlier times.

The smooth suppression of power is exciting because it might help resolve two sources of tension in recent experimental results. A fit to the most recent Planck CMB data is used to fix the parameters of $\Lambda$CDM. Using the model, the Planck collaboration predicts $\sigma_8^{\Lambda CDM}=0.829\pm0.014$ (``TT + lowP", 1$\sigma$ errors). This value is about 2$\sigma$ higher than ``direct" measurements of large scale structure ~\cite{Beutler:2014yhv,Battye:2014qga,Planck:2015xua} (for example, gravitational lensing of the CMB as measured by Planck gives $\sigma_8^{lensing}=0.802\pm0.012$). A reduction of the predicted power spectrum  due to DM-DR interactions as, it occurs in our model for $\alpha_d\sim 10^{-9}$, removes this tension between the Planck fit and LSS data. Interestingly, the tension in $\sigma_8$ is currently also driving the Planck fit for $N_{eff}$ and $H_0$ to lower values, because in $\Lambda$CDM larger values for those parameters would correlate with even larger values for $\sigma_8$. Thus after including DM-DR interactions the Planck fit might prefer $N_{eff}> 3$ which would create more room for dark gluons. Furthermore, larger $N_{eff}$ is correlated in the $\Lambda$CDM fit with a larger value for $H_0$ (to keep the position of the acoustic peaks in the CMB fixed~\cite{Hu:1998tk,Bashinsky:2003tk,Bowen:2001in,Planck:2015xua}). This in turn would allow better agreement between the $H_0$ values from Planck and supernova data~\cite{Planck:2015xua}, another area of mild tension in cosmological data. Clearly, a more quantitative analysis of this issue is desirable and requires including the non-Abelian dark radiation in a full Boltzmann code.

\section{Conclusions}

In this work we have studied a new type of dark sector with massless non-Abelian gauge bosons super-weakly coupled to the DM. There are many different possibilities for the DM coupling to the standard model. As an example, we chose our DM particle to transform as \ensuremath{SU(2)_{weak}\ } triplets and fundamentals under the dark \ensuremath{SU(N)_d \ } gauge group.

Our model has three new parameters, the coupling constant $\alpha_d$, the mass of the DM particles $M_\chi$, and the size of the gauge group $N$. Demanding the correct DM abundance from thermal freeze-out fixes the DM mass in terms of $N$, thus leaving a two dimensional parameter space. Constraints on this parameter space can be derived from several different experiments.

The first set of constraints derives simply from the multiplicity of the dark matter and would even apply if the dark coupling constant were zero. The multiplicity of dark matter affects the usual WIMP searches (direct and indirect detection and collider searches). The effect is simply that dark color multiplicity factors enhance pair production and decrease pair annihilation. Therefore the DM mass required in order to predict the right DM abundance decreases by $\sqrt{2 N}$. It also increases the collider cross-section, placing this type of dark matter within easy reach of the proposed $100$ TeV collider. The decrease in mass and in annihilation cross-section also removes the current tension between thermally produced \ensuremath{SU(2)_{weak}\ } triplet dark matter and H.E.S.S data. Despite the decrease in annihilation cross-section, the \ensuremath{SU(2)_{weak}\ } triplet DM model annihilation cross-section is within the projected reach for CTA.

For $\alpha_d \simgt 10^{-13}$ the dark gluons thermalize with the SM in the early universe and contribute to dark radiation. The limits placed on $N_{eff}$ from Planck constrain $N$ to be at most $3$. However, the self interactions of the dark gluons leave an imprint in the CMB which is distinct from that of free-streaming fluids like neutrinos or dark photons. This can be used to distinguish between the two types of radiation if future experiments establish the need for a non-standard contribution to $N_{eff}$.

Finally we studied the effects of the interactions between DM and DR on the power spectrum. We found that for $\alpha_d \gtrsim 10^{-8}$ the interactions strongly suppress the power spectrum of modes entering the horizon before matter radiation equality and thus such couplings are ruled out. On the other hand, for $\alpha_d \simlt 10^{-8.5}$ the interactions predict a smooth decrease in the power spectrum, which can potentially solve the discrepancy between Planck and large scale structure data and the discrepancy between Planck and Supernova measurements of $H_0$.  We hope to return to this in future work.

\section{Acknowledgments}
We wish to thank Ed Bertschinger, Chris Brust, Mark Hertzberg, David Kaplan, Lloyd Knox, Michele Papucci, Tracy Slatyer, Jesse Thaler, Matt Walters, Kathryn Zurek,  and especially Andy Cohen and Julien Lesgourgues for useful discussions. MS acknowledges support from the Aspen Center for Physics where some of the early work on this paper was performed. This work was supported by the US Department of Energy Office of Science under Award Number DE-SC-0010025. GMT also acknowledges support from a DOE High Energy Physics Graduate Fellowship.

\bibliographystyle{JHEP}


\section{Introduction}
\label{intro}

Ultra-high energy cosmic rays (UHECRs) are the most energetic particles observed in nature, with energies up to a few times $10^{20}$~eV. The first detection of such particles dates back to 1960 and so far, after several decades of observations, the origin of UHECRs still remains an open question.

The most recent and advanced experiments observing UHECRs are the Pierre Auger Observatory \cite{Abraham:2004dt,Aab:2015zoa}, by far the largest experimental set-up devoted to the study of UHECRs, located in the southern hemisphere in Argentina, and the Telescope Array (TA) experiment \cite{AbuZayyad:2012kk,Sagawa:2013pya}, with 1/10 of the Auger statistics at the highest energies, located in the northern hemisphere in the United States. The two experiments adopt a hybrid approach that allows events simultaneously observed with a surface detector and a fluorescence detector to be collected. The method strongly reduces the uncertainties on energy calibration of the large majority of the showers only detected by the surface arrays.

The experimental efforts made in the field of UHECRs produced a number of significant results concerning these particles, such as \cite{Kotera:2011cp,Kampert:2014eja,Blasi:2014roa}: (i) UHECRs are charged particles, with limits on neutral particles around $10^{19}$~eV at the level of few percent for photons and well below for neutrinos \cite{Abraham:2009qb,Abu-Zayyad:2013dii,Abreu:2013zbq}, (ii) the spectra observed on Earth show a slight flattening at energies around $5\times 10^{18}$~eV (called the ankle), with (iii) a steep suppression at the highest energies around $10^{20}$~eV \cite{Abu-Zayyad:2013qwa,ThePierreAuger:2013eja}.

A key piece of information in the physics of UHECRs is the composition of these extremely energetic particles. In this respect experimental observations do not converge to the same conclusions, with different experiments pointing toward different interpretations. The measurements performed by TA are consistent with predictions for a pure proton flux in the whole energy range observed, starting from the lowest energies $10^{18}$~eV up to the highest (see \cite{Tinyakov:2014lla} and references therein), whereas those by Auger signal a richer phenomenology with a flux dominated by protons and light nuclei at low energies but becoming progressively heavier for $E > 3\times10^{18}$~eV \cite{Abreu:2013env,ThePierreAuger:2013eja}, albeit using substantially different analysis
modalities. On the other hand, TA measurements have wider error bars and a recent analysis \cite{Abbasi:2015xga} has shown that they are also compatible with Auger data within their uncertainties.

During their journey from the source to the observer, UHECRs cover cosmological distances and, apart from magnetic fields not considered here, interact with the radiation fields of the astrophysical backgrounds, mainly the cosmic microwave background (CMB) and the extragalactic background light (EBL). The propagated energy spectrum of nucleons\footnote{Hereafter discussing freely propagating UHE nucleons we will always refer only to protons because the decay time of neutrons is much shorter than all other time scales involved \cite{Aloisio:2008pp,Aloisio:2010he,Allard:2011aa}.} is affected almost only by the CMB radiation field and the processes that influence the propagation are: (i) pair production and (ii) photo-pion production \cite{Berezinsky:2002nc}. On the other hand, the propagation of heavier nuclei is affected also by the EBL and the interaction processes relevant are: (i) pair production and (ii) photo-disintegration \cite{Puget:1976nz,Aloisio:2008pp,Aloisio:2010he,Allard:2005ha,Allard:2011aa}.

In this work we will consider both scenarios, proton-dominated and mixed composition, that are adopted in the interpretation of the most recent experimental data. A pure proton composition will be used in the analysis of TA data: under this assumption the observations can be described in the framework of the dip model \cite{Berezinsky:2002nc,Aloisio:2006wv}, which gives an explanation of the ankle observed in the spectrum as an effect of the pair production by protons on the CMB, while the strong flux suppression at the highest energies is the effect of the photo-pion production process, the so-called Greisen-Zatsepin-Kuzmin (GZK) cut-off \cite{Greisen:1966jv,Zatsepin:1966jv}.

The Auger data, on the other hand, indicate that the overall proton contribution to the spectrum is mostly restricted to the lowest energies $\le 3\times 10^{19}$~eV well below the photo-pion production threshold ($\sim 6\times 10^{19}$~eV) \cite{Aloisio:2013hya,Taylor:2013gga}. The ankle feature has to be interpreted as the transition between two different classes of sources, the high-energy one with large metallicity while the flux suppression at the highest energies is generated by  nuclei photo-disintegration and/or the maximum acceleration energy \cite{GR,Aloisio:2009sj,Aloisio:2013hya}.

As was first realised by Berezinsky and Zatsepin \cite{Beresinsky:1969qj}, the proton content of UHECRs, particularly at the highest energies, is a crucial quantity that affects the fluxes of secondary photons and neutrinos produced by UHECR propagation in the intergalactic space \cite{Stecker:1973sy,Berezinsky:1975zz}. In the present paper, we will focus on the connections between UHECR chemical composition and source cosmological evolution and the production of secondary neutrinos. In the framework of different source models, we will give a comprehensive picture of the expected UHECR and neutrino fluxes in comparison with the latest experimental data (previous estimates of these quantities by different authors include \cite{Hill:1983xs,Engel:2001hd,Kalashev:2002kx,Semikoz:2003wv,Allard:2006mv,Takami:2007pp,Gelmini:2011kg,Roulet:2012rv,Stanev:2014asa}). Particularly relevant are the recent observations performed by the IceCube observatory in Antarctica, which detected, for the first time, an extraterrestrial neutrino flux at energies around $10^{15}$~eV \cite{Aartsen:2013jdh,Aartsen:2013bka}, together with the limits on neutrino fluxes fixed by Auger \cite{Aab:2015kma}.

The calculations presented here are based on {\it SimProp}\footnote{{\it SimProp}  is available upon request to \url{SimProp-dev@aquila.infn.it}.} \cite{Aloisio:2012wj,Aloisio:2015sga}, a Monte Carlo (MC) propagation code designed to study UHECR propagation and upgraded here to determine the fluxes of secondary neutrinos, taking into account all possible production channels. The production of secondary neutrinos by the propagation of UHECRs is contributed by the overall universe \cite{Hill:1983xs}, therefore the cosmological evolution of both (i) the sources of UHECRs and (ii) the astrophysical backgrounds play an important role. While the impact on the expected neutrino fluxes of the latter will be one of the main findings of this work, the poor knowledge of the former represents a source of uncertainty that should be quantified. The cosmological evolution of the EBL is inferred from (few) observations at different redshifts implementing {\it ad hoc} models \cite{Stecker:2005qs,Stecker:2006eh,Kneiske:2003tx}. These models show sizeable differences only at high redshift ($z>4$), not actually relevant in the determination of UHE nuclei fluxes but affecting the production of secondary neutrinos \cite{Allard:2011aa}. In the present paper we will use the two recipes proposed in \cite{Stecker:2005qs,Stecker:2006eh} and \cite{Kneiske:2003tx} showing the impact of the considered EBL evolution model on our findings.

The results discussed here, addressing the expected fluxes of secondary neutrinos through a detailed analysis of some common UHECR models, has a twofold interest: on one hand, as stated above, using the latest observations we can already draw interesting conclusions on the sources of UHECRs, on the other hand the study presented should be intended as a benchmark computation to assess the discovery capabilities of the next generation experiments in both fields of UHECR  \cite{Adams:2013hqc} and neutrino observations \cite{Allison:2011wk,Klein:2012bu,Katz:2014aoa}.

The paper is structured as follows: in Section~\ref{sec:prop} we illustrate the physics of the propagation of UHECR through intergalactic space, with a particular emphasis on the processes of secondary particles production and on the related new features implemented in {\it SimProp}; in Section~\ref{sec:scen} we illustrate different scenarios for UHECR sources and their output in terms of secondary neutrinos, focusing on the experimental limits; finally, in Section~\ref{sec:res}, we summarise and discuss our results.

\section{UHECR propagation and secondary neutrinos production}
\label{sec:prop}

Ultra-high energy (UHE) protons or nuclei propagating in the intergalactic space interact with photons of the CMB and EBL; these interactions are particularly important not only as a mechanism of energy losses for the propagating particles but also being the processes responsible for the production of secondary particles such as neutrinos, gamma rays and electron-positron pairs\footnote{In the present paper we will not discuss the case of secondary gamma rays and electron-positron pairs, leaving this subject to a possible forthcoming publication.}. The interaction rate associated to such processes can be written in a very general form as \cite{Stecker:1968uc,Aloisio:2008pp}:
\begin{equation}
\frac{1}{\tau} = \frac{1}{\Gamma^2}\int_{\epsilon'_{\min}}^{+\infty} \epsilon'\sigma(\epsilon')\int_{{\epsilon'}/{2\Gamma}}^{+\infty} \frac{n_\gamma(\epsilon)}{2\epsilon^2} d\epsilon \,d\epsilon' \label{eq:tau}~,
\end{equation}
where $\Gamma$ is the Lorentz factor of the particle, $\sigma(\epsilon')$ is the total cross-section associated to the particle interactions, $\epsilon'$ is the background photon energy in the particle rest frame, $\epsilon'_{\min}$ is the lowest value of~$\epsilon'$ above which the interaction is kinematically possible (threshold), and $n_\gamma(\epsilon)\,d\epsilon$ is the number per unit volume of background photons with energy between $\epsilon$ and $\epsilon + d\epsilon$ in the laboratory reference frame. The photon energy in the particle rest frame is related to that in the laboratory frame by $\epsilon' = \Gamma\epsilon(1-\cos\theta)$, where $\theta$ is the angle between the particle and photon momenta $(0 \le \epsilon' \le 2\Gamma\epsilon)$.

The interaction processes that involve UHECRs with background photons are pair-production and photo-pion production in the case of protons, and also photo-disintegration in the case of heavier nuclei. Given the distribution of background photons and the energies involved, the propagation of protons is substantially affected only by the CMB radiation field, while in the case of nuclei, and only for the photo-disintegration process, also the EBL field plays an important role \cite{Aloisio:2008pp,Aloisio:2010he}. As we will discuss below, the effect of EBL on proton propagation has an important role for the production of secondary neutrinos, but it negligibly affects the expected proton flux.

Let us start from the process of photo-pion production. Nucleons ($N$), whether free or bound in nuclei, with Lorentz factor $\Gamma\gtrsim 10^{10}$ interacting with the CMB photons give rise to the photo-pion production process:
\begin{equation}
N+\gamma \to N + \pi^0   \qquad N+\gamma \to N + \pi^{\pm}.
\end{equation}
At lower energies $\Gamma< 10^{10}$, even if with a lower probability, the same processes can occur on the EBL field. In the case of UHE protons propagating in the CMB, the photo-pion production process involves a sizeable energy loss producing the so-called GZK cut-off \cite{Greisen:1966jv,Zatsepin:1966jv}, a sharp suppression of the flux of protons expected on Earth at $E\simeq 6\times 10^{19}$~eV, which corresponds to the threshold for photo-pion production, that in the nucleon rest frame reads $\epsilon'_{\min} = m_{\pi} + m_{\pi}^2/2m_N \approx 145$ MeV. The photo-pion production cross-section has a complex behavior with a number of peaks corresponding to different hadronic resonances, the largest one being the $\Delta$ resonance placed at $\epsilon' =\epsilon_{\Delta}\approx 340$ MeV \cite{Berezinsky:2002nc}. At energies much larger than $\epsilon_{\Delta}$ the cross-section sits on an approximately constant value \cite{Berezinsky:2002nc}. The photo-pion production process holds also for nucleons bound within UHE nuclei, being the interacting nucleon ejected from the parent nucleus, but this process is subdominant with respect to nucleus photo-disintegration except at extremely high energies \cite{Allard:2011aa}.

UHE nuclei propagating through astrophysical backgrounds can be stripped of one or more nucleons by the interactions with CMB and EBL photons, giving rise to the process of photo-disintegration:
\begin{equation}
(A,Z) + \gamma \to (A-n, Z-n') + nN
\end{equation}
$n$ ($n'$) being the number of stripped nucleons (protons). In the nucleus rest frame the energy involved in such processes is usually much less than the rest mass of the nucleus itself, therefore in the laboratory frame all fragments approximately inherit the same Lorentz factor of the parent nucleus, i.e. we can neglect nucleus recoil. The cross-section is dominated by a smooth peak, the giant dipole resonance, that appears for photon energies close to the threshold ($8~{\rm MeV} \approx \epsilon'_{\min} < \epsilon' \lesssim 30~{\rm MeV}$) \cite{Puget:1976nz}; in this regime nucleons collectively behave as a two fluid system: a proton fluid vibrating against a neutron fluid. The giant dipole resonance corresponds to the extraction of one nucleon and it is the dominating process in UHE nuclei propagation \cite{Puget:1976nz,Allard:2005ha,Aloisio:2008pp,Aloisio:2010he}. At larger energies $\epsilon'>30~{\rm MeV}$ the quasi-deuteron process dominates, in which the photon interacts with one or two nucleons inside the nucleus with the extraction of two or more nucleons. This regime corresponds to an almost constant cross-section and has a small impact on the propagation of UHE nuclei \cite{Puget:1976nz,Allard:2005ha,Aloisio:2008pp,Aloisio:2010he}.

The interaction rate associated to the processes of photo-pion production and photo-disintegration can be written using Eq. (\ref{eq:tau}) specifying the cross-section, the background photon density and all relevant kinematical thresholds of the process.

Protons and nuclei with Lorentz factor $\Gamma\gtrsim 10^{9}$ can undergo the process $p + \gamma \to p + e^{+} + e^{-}$. The mean free path associated to pair production is relatively short compared with all others length scales of UHECR propagation, with a very small amount of energy lost by the propagating particle in each interaction. Therefore, as was shown elsewhere it is possible to neglect the stochasticity of the process assuming a continuous energy loss evolution \cite{Berezinsky:2002nc}.

Particles covering cosmological distances feel the effect of the changes in the background universe due to cosmology. The expansion of the universe produces an adiabatic energy loss mechanism that continuously reduces the energy of the propagating particles. Assuming standard cosmology we can write the energy lost (adiabatically) per unit time by UHECRs as
\begin{equation}
\left(-\frac{1}{\Gamma}\frac{d\Gamma}{dt}\right)_\textrm{ad} = H(z) = H_0\sqrt{(1+z)^3\Omega_\textrm{m} + \Omega_\Lambda}\label{adiabatic}
\end{equation}
where $z$ is the redshift at time~$t$, $H_0 \approx 70~\mathrm{km}/\mathrm{s}/\mathrm{Mpc}$ is the Hubble constant, $\Omega_\textrm{m} \approx 0.3$ is the matter density, and $\Omega_\Lambda \approx 0.7$ is the dark energy density.

Let us now concentrate on the unstable particles (pions, free neutrons and unstable nuclei) produced by the propagation of UHECRs through photo-pion production and photo-disintegration. In most cases the decay length of such particles is much shorter than all other relevant length scales, so these particles decay very soon giving rise to secondary neutrinos.

There are two processes by which neutrinos can be produced in the propagation of UHECRs:
\begin{itemize}
 \item the decay of charged pions produced by photo-pion production, $\pi^{\pm}\to \mu^{\pm} + \nu_\mu(\bar{\nu}_\mu)$, and the subsequent muon decay $\mu^{\pm}\to e^{\pm}+\bar{\nu}_\mu(\nu_\mu)+\nu_e(\bar{\nu}_e)$;
 \item the beta decay of neutrons and nuclei produced by photo-disintegration: $n \to p + e^{-} + \bar{\nu}_e$, $(A,Z) \to (A,Z-1) + e^{+} + \nu_e$, or $(A,Z) \to (A,Z+1) + e^{-} + \bar{\nu}_e$.
\end{itemize}
These processes produce neutrinos in different energy ranges: in the former the energy of each neutrino is around a few percent of that of the parent nucleon, whereas in the latter it is less than one part per thousand (in the case of neutron decay, larger for certain unstable nuclei). This means that in the interactions with CMB photons, which have a threshold around $\Gamma\gtrsim 10^{10}$, neutrinos are produced with energies of the order of $10^{18}$~eV and $10^{16}$~eV respectively. Interactions with EBL photons contribute with a much lower probability than CMB photons, affecting a small fraction of the propagating protons and nuclei. Neutrinos produced through interactions with EBL, characterised by lower thresholds, have energies of the order of $10^{15}$~eV in the case of photo-pion production and $10^{14}$~eV in the case of neutron decay.

The results presented in this work are obtained using {\it SimProp}~v2r2 \cite{Aloisio:2015sga}, an extension for the study of secondary neutrinos of the {\it SimProp} \cite{Aloisio:2012wj} computation scheme. The most relevant technical details and a more detailed discussion of the algorithm and numerical code implemented can be found in Ref.~\cite{Aloisio:2015sga}. The following new features have been implemented in {\it SimProp}:
\begin{itemize}
\item photo-pion production is taken into account as a stochastic process, for protons as well as nuclei, using both the EBL and CMB radiation fields;
\item all allowed isobars for each $A$, stable or not, can be produced, though unstable ones are assumed to decay immediately so only stable isobars for each $A$ are propagated;
\item the decay products of pions, neutrons and unstable nuclei are tracked;
\item several technical improvements were implemented, resulting in a faster execution time;
\end{itemize}

A {\it SimProp} run consists of $N$~events (for a user-supplied value of $N$), each consisting in the generation of a primary particle and its propagation from the source to the observer along with that of any secondary particles produced. All particles are assumed to travel rectilinearly (without taking into account the possible presence of intergalactic magnetic fields) at the speed of light, so only one coordinate, the redshift $z$, is used to keep track of particle positions.

In order to test our algorithm, particularly the new module of the code developed to compute the flux of secondary neutrinos, we have compared our results with those of a simulation performed by Engel, Seckel and Stanev \cite{Engel:2001hd}. In figure  \ref{fig:ESS} we plot the fluxes of secondary neutrinos and anti-neutrinos $\nu_e, \bar{\nu}_e, \nu_\mu, \bar{\nu}_\mu$, where results of \cite{Engel:2001hd} are plotted as smooth lines while our results through histograms.
\begin{figure}
\includegraphics[width=0.5\textwidth]{neutrinos1.pdf}
\includegraphics[width=0.5\textwidth]{neutrinos2.pdf}
\caption{ Fluxes of neutrinos at Earth expected in the Engel--Seckel--Stanev~\cite{Engel:2001hd} scenario. The histograms were obtained via {\it SimProp} and the smooth curves refer to the data from~Ref.~\cite{Engel:2001hd}. [Left panel] $\nu_e$ flux (blue solid) and $\bar{\nu}_e$ flux (red dashed). [Right panel] $\nu_\mu$ flux (blue solid) and $\bar{\nu}_\mu$ flux (red dashed).
\label{fig:ESS} }
\end{figure}
The two simulations show results that are in quite good agreement. The only substantial difference is that our electron antineutrino spectrum lacks the second peak at around $10^{18}$~eV. This follows from our choice of neglecting the sub-dominant interaction channels in which several pions and/or heavier mesons are produced. This simplification is observationally irrelevant, because the phenomenon of neutrino oscillations (not shown here) will mix the flavours so that at Earth all neutrino flavours will show equal fluxes and the total antineutrino production at around $10^{18}$~eV is dominated by $\bar{\nu}_\mu$. Neutrino fluxes presented in the forthcoming discussion are always the sum of neutrino and anti-neutrino fluxes over all flavours.

\section{UHECR scenarios and cosmogenic neutrino fluxes}
\label{sec:scen}

Theoretical models aiming at the explanation of UHECR observations can be distinguished in two different scenarios depending on the chemical composition: on one hand the dip model, which assumes a pure proton composition, on the other hand models based on a mixed composition with the low energy part of the spectrum dominated by protons and the highest energy part dominated by nuclei. As we will show below, the flux of secondary neutrinos, being extremely sensitive to the chemical composition, can be a powerful tool to solve the puzzle.

Neutrinos produced by the interaction of UHECRs, because of their extremely low interaction rate, arrive on Earth unmodified, except for redshift energy losses and flavor oscillations, with the overall universe contributing to their flux. This is an important point that makes neutrinos a viable probe not only of the chemical composition of UHECRs but also of the cosmological evolution of sources that, as we will show below,  can be also constrained by the neutrino flux observed on Earth.
We will consider the case of sources with no cosmological evolution, with the same cosmological evolution as that active galactic nuclei (AGN), an astrophysical object supposed to play a role in particle acceleration to the highest energies \cite{Berezinsky:2002vt}, and with the cosmological evolution of the star formation rate (SFR).

The production of UHECRs at the source will be modelled as a power law in the Lorentz factor of the injected particles $\Gamma$. We will also account for the maximum acceleration energy that sources can provide, $\Gamma_{max}$. This last quantity will be taken to be rigidity dependent, namely given the maximum acceleration energy for protons $E_{max}^p$ the same quantity for any nuclei species will be $E_{max}^A=ZE^p_{max}$,  $Z$ being the atomic number (charge) of the nucleus. The spectrum of cosmic rays effectively injected in the intergalactic medium by a collection of sources can be affected by the convolution of the spectrum of individual sources and their luminosity and/or maximum energy achieved. For instance, under realistic assumptions \cite{Kachelriess:2005xh}, the effective spectrum can have a broken power-law dependence. For soft injection (spectral index $\gamma_g>2$), which is mainly of interest in the case of sources injecting only protons and light nuclei, we will assume as in  \cite{Aloisio:2013hya} an effective injection spectrum of the type
\begin{equation}
Q_{inj}(\Gamma,z) \propto S(z) e^{-\Gamma/\Gamma_{max}} \left\{
\begin{array}{cc}
1/\Gamma^2 & ~~~~~~ \Gamma < \Gamma_0\\
\frac{1}{\Gamma_0^2}\left (\frac{\Gamma}{\Gamma_0} \right
)^{-\gamma_g}&~~~~~~~\Gamma \ge \Gamma_0.
\end{array}\right.
\label{eq:Qinj}
\end{equation}
The value of $\Gamma_0$ depends on the specific adopted model; in the case $\gamma_g>2$ for phenomenological purposes typical values are $\Gamma_0\simeq 10^{6}\div 10^{8}$, as this makes the requirements in terms of the source luminosity less severe \cite{Berezinsky:2002nc}. In this work we use $\Gamma_0 = 10^8$. Following the approach of \cite{Aloisio:2013hya}, in the case of sources also injecting heavier nuclei, we are especially interested in hard injection, where the slope of the generation spectrum is $\gamma_g < 2$; in this work we do not use any spectral break for such sources.

The term $S(z)$ in equation (\ref{eq:Qinj}) accounts for the cosmological evolution of sources. In the forthcoming discussion we will consider separately three different cases: the case of no evolution with $S(z)=1$, the case of the cosmological evolution associated to the SFR \cite{Yuksel:2008cu,Wang:2011qc,Gelmini:2011kg}:
\begin{equation}
S_{\rm SFR}(z)= \left \{
\begin{array}{lll}
(1+z)^{3.4} & \quad & z<1 \\
2^{3.7}(1+z)^{-0.3} & \quad & 1<z<4\\
2^{3.7}5^{3.2}(1+z)^{-3.5} & \quad  & z>4
\end{array} \right.
\label{S_SFR}
\end{equation}
and the case of the evolution typical of the AGN \cite{Hasinger:2005sb,Ahlers:2009rf,Gelmini:2011kg}
\begin{equation}
S_{\rm AGN}(z)= \left \{
\begin{array}{lll}
(1+z)^{5.0} & \quad & z<1.7 \\
(1+1.7)^{5.0} & \quad & 1.7<z<2.7\\
(1+1.7)^{5.0} 10^{(2.7-z)} & \quad & z>2.7
\end{array} \right.
\label{S_AGN}
\end{equation}

All computations presented here are performed under the assumption of a homogenous distribution of sources. This assumption does not affect the expected neutrino spectra because in the case of neutrinos the overall universe, up to the maximum redshift, contributes to the fluxes and possible flux variations due to a local inhomogeneity in source distribution gives a negligible contribution to the total flux. We also fix a maximum redshift of the sources $z_{max}=10$, which is the typical redshift of the first stars (pop III) \cite{Berezinsky:2011bb}. In any case the expected fluxes of primary and secondary particles are almost independent of $z_{max}$ if $z_{max}>3$.

Once produced at cosmological distances neutrinos travel toward the observer almost freely, the opacity of the universe to neutrinos being relevant only at the redshifts $z>10$ \cite{Gondolo:1991rn,Weiler:1982qy}. Therefore, given the assumptions discussed above, in our computations we have neglected any effect due to neutrino propagation apart from the adiabatic energy losses due to the expansion of the universe as discussed in the previous section.

In the following subsections we will consider the two cases of UHECRs consisting of a pure proton composition or a mixed composition, assessing their consequences in terms of secondary neutrino production.

\begin{figure}
\includegraphics[width=0.5\textwidth]{E3J_dip.pdf}
\includegraphics[width=0.5\textwidth]{E2Jnu_dip.pdf}
\caption{[Left Panel] Fluxes of protons expected at Earth in proton-only scenarios with various models for the cosmological evolution of sources (solid red: no evolution; dashed green: SFR evolution; dot-dashed blue: AGN evolution \cite{Gelmini:2011kg}) normalized to TA data (purple squares);  Auger data (olive disks) and KASCADE-Grande data \cite{Apel:2011mi} (blue triangles) are also shown for comparison. [Right Panel] Fluxes of neutrinos in the same scenarios (same color code; from bottom to top: no evolution, SFR, and AGN), with colored bands showing the difference between the Stecker (solid) \cite{Stecker:2005qs,Stecker:2006eh} and Kneiske (dashed) \cite{Kneiske:2003tx} EBL models; thin solid lines are neutrino fluxes obtained taking into account only the interaction with CMB photons.}
\label{fig:TAp}
\end{figure}

\subsection{The dip model}

The dip model is based on the assumption that UHECRs are composed of only protons, therefore in this case we will restrict our analysis to the TA data. In order to describe the overall observed spectrum in terms of a single proton component we need to assume a rather steep injection spectrum with $\gamma_g>2.3$ \cite{Berezinsky:2002nc,Aloisio:2006wv}. The actual value of the injection spectrum depends on the assumptions about the cosmological evolution of sources: increasing the co-moving density of sources at large redshift the required value of $\gamma_g$ decreases \cite{Berezinsky:2002nc,Aloisio:2006wv}.

In figure \ref{fig:TAp} we consider the three different scenarios for cosmological evolution of sources mentioned above, namely: no evolution, SFR and AGN evolution, assuming that contributing sources are restricted to redshifts $z<z_{max}=10$. In order to reproduce the TA observations shown in the left panel in figure \ref{fig:TAp}, the injection power law index $\gamma_g$ needs to be changed according to the cosmological evolution: we used $\gamma_g = 2.6$, $2.5$, and $2.4$ respectively in the three cases of cosmological evolution considered. The maximum acceleration energy used to compute fluxes in figure \ref{fig:TAp} is $E_{max}^p=10^{22}$~eV. The source emissivities, i.e. the energy of UHECRs emitted per unit volume and time, required to reproduce TA data are ${\cal L}_0 \approx 1.5\times10^{46}$, $6\times10^{45}$ and $3.5\times10^{45}$~erg/Mpc$^3$/yr at $z=0$ respectively in the case of $\gamma_g=2.6, 2.5$, and 2.4, assuming $\Gamma_0=10^8$ (see equation (\ref{eq:Qinj})).

Results presented in figure \ref{fig:TAp} are quite interesting, showing that models with cosmological evolution of sources stronger than the case of AGN predict a total neutrino flux (right panel of figure \ref{fig:TAp}) in excess of the experimental limits fixed by Auger and the observations of IceCube. Moreover, as predicted in \cite{Berezinsky:2002vt}, the AGN hypothesis seems to reproduce quite well TA observations in terms of the UHE flux suppression, where our calculations result in a neutrino flux at a level which would have been detected by both IceCube and Auger.

\subsection{Mixed composition}

Let us consider here the case of a mixed composition, which is more suitable to describe the observations of Auger; in this analysis we will refer only to Auger data on both spectrum and chemical composition.

As shown by different authors \cite{Taylor:2013gga,Aloisio:2013hya,Fang:2013cba} scenarios with mixed composition have peculiar characteristics that, in order to describe Auger observations, imply different source families (see \cite{Globus:2015xga,Unger:2015laa} for models with a single type of source). At the highest energies ($E>5\times 10^{18}$~eV) Auger observes predominantly heavy elements, with $A \gtrsim 12$, while at the lowest energies, as a matter of fact, all UHECR experiments observe a light composition dominated by protons. To reproduce such observations two classes of sources are needed \cite{Taylor:2013gga,Aloisio:2013hya,Fang:2013cba}: one that injects only light elements with low maximum energy and steep injection and one that injects also heavy elements with higher maximum energy and a
flatter injection.

In this context, while the extragalactic nature of sources that contribute also heavy elements cannot be questioned, sources providing only light elements can be both galactic or extragalactic, though the galactic hypothesis poses severe problems with both the limits on observed anisotropy and the (galactic) acceleration mechanism \cite{Aloisio:2013hya}. In this case the production of secondary neutrinos will be connected only to the propagation of UHECRs injected by extragalactic sources and will result in a particularly low flux, far below the detection threshold of any experiment. For a different approach with galactic sources providing the low energy light component and secondary EeV neutrinos produced directly at the source the reader can refer to \cite{Fang:2013vla,Murase:2014foa}.

In the present paper, in order to study the production of cosmogenic neutrinos in the framework of the Auger observations, we will consider the case discussed in \cite{Aloisio:2013hya}, but with slightly retuned spectral parameters in order to account  for the different approximations used in {\it SimProp} and those used in \cite{Aloisio:2013hya} and the differences between the 2009 Auger data used in \cite{Aloisio:2013hya} and the 2013 data used here. The model is composed of two different classes of extragalactic cosmic ray sources distributed uniformly in a comoving volume. The first class injects ($75\%$) protons and ($25\%$) helium nuclei with $\gamma_g = 2.6, 2.5$, and 2.4, depending on the cosmological evolution as in the previous section, with a maximum energy $E_{max} = Z\times 2\times10^{18}$~eV and a source emissivity  ${\cal L}_0 \approx 1.5\times 10^{46}$, $5.5\times 10^{45}$ and $2\times10^{45}$~erg/Mpc$^3$/yr at $z=0$ respectively in the three cases $\gamma_g=2.6$, 2.5 and 2.4, assuming $\Gamma_0=10^8$ (see equation (\ref{eq:Qinj})). The second class also injects heavier nuclei, with $35\%$ protons, $30\%$ helium, $25\%$ CNO and $10\%$ MgAlSi (these fractions being defined at fixed energy per nucleus equal to $10^{18}$~eV), with $\gamma = 1.0$, $E_{max} = Z\times 6 \times10^{18}$~eV and a source emissivity ${\cal L}_0\approx5.5\times10^{44}$~erg/Mpc$^3$/yr (above $10^7$~GeV/n).

It is interesting to note here that in the case of nuclei the cosmological evolution of sources has a larger impact on the expected fluxes than in the case of a pure proton composition. This can be understood taking into account that cosmological evolution implies an increased source density at high redshift, through terms of the type $(1+z)^m$ (see above), that corresponds to an increased effect of the photo-disintegration interaction. This fact produces an increased number of secondary light cosmic rays, particularly at intermediate energies $E \sim 10^{18.5}$~eV, spoiling the agreement with both flux and chemical composition observed by Auger. Therefore sources providing heavy elements cannot show any evolution with redshift. This result is actually quite interesting, showing the capability of the Auger data on flux and chemical composition in constraining also the cosmological characteristics of the possible UHECR sources.

Using two classes of sources, with the parameters discussed above, and assuming no cosmological evolution of those providing also heavy elements, we obtain the fluxes at Earth of UHECRs and neutrinos shown in figure \ref{fig:2comp1}, left and right panel respectively. The corresponding chemical composition is shown in figure \ref{fig:2comp1xmax}.
\begin{figure}
\includegraphics[width=0.5\textwidth]{E3J_two.pdf}
\includegraphics[width=0.5\textwidth]{E2Jnu_two.pdf} \\
\caption{Fluxes of UHECRs (left) and neutrinos (right) expected at Earth in the two-class model. UHECR fluxes are grouped according to their mass number at Earth as labeled. Olive disks in the left panel are Auger data.  TA data (purple squares) and KASCADE-Grande data \cite{Apel:2011mi} (blue triangles) are also shown for comparison. Neutrino fluxes are shown (from bottom to top) for the homogeneous (red), SFR (green) and AGN (blue) evolution of the steep (light) component, assuming the Stecker (solid) and Kneiske (dashed) EBL evolution models.  Fluxes of protons and He are summed over the two different classes of sources. For the injection chemical composition and source emissivities see text.}
\label{fig:2comp1}
\end{figure}

\begin{figure}
\includegraphics[width=0.5\textwidth]{Xmax_two.pdf}
\includegraphics[width=0.5\textwidth]{RMS_two.pdf} \\
\caption{Mean (left) and standard deviation (right) of~$X_{\max}$ predicted by EPOS-LHC (solid), QGSJet II-04 (dashed) and SIBYLL 2.1 (dot-dash) for pure protons (top, blue), pure iron (bottom, red), and our two-class model (in the no-evolution scenario), along with Auger data (lines: statistical errors; brackets: systematic errors).}
\label{fig:2comp1xmax}
\end{figure}

The results presented in left panel of figure \ref{fig:2comp1} and in figure \ref{fig:2comp1xmax} confirm the capability of models based on two classes of sources of reproducing the Auger observations, as previously found in \cite{Aloisio:2013hya,Taylor:2013gga,Fang:2013cba}.

The corresponding flux of secondary neutrinos is extremely low at the highest energies as follows from the right panel of figure \ref{fig:2comp1}. This is a direct consequence of the chemical composition that at the highest energies is dominated by heavy nuclei. In this case the production of secondary neutrinos, being mainly connected with the photo-pion production process, involves only nucleons inside nuclei with energy $E/A$. Therefore, assuming a nucleus maximum energy at the level of $E_{max}=Z\times 5\times 10^{18}$~eV, as in figure \ref{fig:2comp1}, implies that only the tails of the cosmic rays spectra contribute to the production of secondary neutrinos with a substantial reduction of their expected flux. Only the light component of UHECRs, being proton dominated, will provide a sizeable amount of secondary neutrinos. Therefore, given the constraints introduced by the Auger data, we can conclude that only the PeV component of cosmogenic neutrinos can be detected, with an expected flux close to the IceCube sensitivity assuming a moderate (SFR) cosmological evolution for the light component of UHECRs and in excess of it for a strong (AGN) evolution, whereas EeV neutrino fluxes are way below the Auger sensitivity.

\section{Conclusions}
\label{sec:res}

Cosmogenic neutrinos mainly originate from the decay of pions (and muons) produced in photo-hadronic interactions of protons and heavy nuclei with CMB and EBL backgrounds. These processes are efficient only in the case of protons, while in the case of heavier nuclei photo-hadronic interactions are significantly suppressed. Therefore the bulk of secondary neutrinos is produced by the lighter components of UHECRs.

There is a solid consensus about the light composition of UHECRs in the low energy part of the observed spectrum. This assures a flux of cosmogenic neutrinos in the PeV energy region, produced by the proton's photo-pion production process on the EBL photons. Our remaining ignorance in predicting this neutrino flux is due to the cosmological evolution of cosmic ray sources and, to a lesser extent, to the uncertainties in the EBL evolution. Models with a strong cosmological evolution, as in the case of AGN, produce a flux of cosmogenic neutrinos almost at the level of the IceCube observations in the PeV region. This fact enables a partial constraining of cosmological evolution of sources, disfavouring models with a too strong evolution: $S(z) \propto (1+z)^m$ with $m \gtrsim 3.5$.

At the highest (EeV) energies the flux of cosmogenic neutrinos mainly originates from the highest energy tail of the UHECR spectrum. If the most energetic cosmic rays, at around $10^{20}$~eV, are mainly protons, then a sizeable flux of cosmogenic EeV neutrinos will be produced. In this case, as before, the experimental limits on EeV neutrinos can be used to constrain the cosmological evolution of sources. For instance, Auger limits on neutrino fluxes already disfavour too strong evolution models, i.e. models with $m\gtrsim 5$, and moderate cosmological evolution, i.e. models with $m<5$, can be tested by future UHE neutrino detectors. Conversely, if the highest energy tail of UHECRs is composed mainly of heavy nuclei, as in models reproducing Auger data on spectrum and chemical composition, then the flux of cosmogenic EeV neutrinos is far below the detection threshold of any running or planned detector.

\section{Acknowledgements}
RA thanks V. Berezinsky and P. Blasi for joint work in the field of UHECRs. The authors thank the Gran Sasso Science Institute where part of this work was developed. The research of DB is supported by SdC Progetto Speciale Multiasse ``La Societ\`a della Conoscenza in Abruzzo'' PO FSE Abruzzo 2007 - 2013.

\bibliographystyle{JHEP}