text:
\section{Appendix A: Construction of likelihood function}
Here, we provide the individual contributions to the full likelihood function defined in Eq. (3)
The likelihood that a given event $i$ is astrophysical is the following
\begin{eqnarray}
P_{{\rm{astro}},i} = \frac{\sum_{f = e,\mu,\tau} P_{\rm{t}}(f|t_i)\int dE  A_{\rm{eff}}^f(E,\vec{x}) R(E,E_{\rm{dep}}E_{\rm{dep}})  \Phi_{\rm{astro}}}{\int dE_{\rm{dep}} \int d^2\vec{x}\sum_{f = e,\mu,\tau} P_{\rm{t}}(f|t_i)\int dE  A_{\rm{eff}}^f(E,\vec{x}) R(E,E_{\rm{dep}})  \Phi_{\rm{astro}}},
\label{eq:pastro}
\end{eqnarray}
where $P_{\rm{t}}$ is the probability that the observed topology $t =$ (track, shower) originated from incoming neutrino flavor $f$. We introduce the detector energy resolution $R_E$ as normal distributions with uncertainties given in~\cite{Aartsen:2015rwa}. The neutrino effective areas for each flavor, $A_{\rm{eff}}^f(E,\vec{x})$, parametrize the detector geometry as well as the neutrino-nucleus (and neutrino-electron) cross sections and the effects of attenuation from propagation within the earth, atmosphere and ice sheet; they thus depend on energy and zenith angle~\cite{IceCubeScience}. $\Phi_{\rm{astro}}$ is the solution to the cascade equation (1) and represents the incoming neutrino flux at Earth, after DM-neutrino interactions have taken place. It depends on neutrino energy and arrival direction $\vec x = (b,l)$. The trivial solution (where $\sigma_{\nu-DM} =0$, or $\tau = 0$) is an isotropic power law in energy. The denominator of Eq. \eqref{eq:pastro} ensures the likelihood is normalized to one. The integral over $E_{\rm{dep}}$ goes from 10 TeV to 10 PeV; this means that the non-observation of events in energy ranges where a $P_{\rm{astro}}$ is high is implicitly taken into account in the likelihood.

The atmospheric likelihood per event is the following
\begin{eqnarray}
P_{{\rm{atm}},i} = \frac{\sum_{f = e,\mu} P_{\rm{veto}} P_{\rm{t}}(f|t_i)\int dE  A_{\rm{eff}}^f(E,\vec{x}) R_E(E,E_{\rm{dep}})  \Phi_{\rm{atm}}}{\int dE_{\rm{dep}} \int d^2\vec{x}\sum_{f = e,\mu,\tau} P_{\rm{veto}} P_{\rm{t}}(f|t_i)\int dE  A_{\rm{eff}}^f(E,\vec{x}) R(E,E_{\rm{dep}})  \Phi_{\rm{atm}}}.
\label{eq:patmo}
\end{eqnarray}
The energy distribution and flavor ratio of the atmospheric flux are specified by the Honda model~\cite{Honda:2006qj,Sanuki:2006yd} with the Gaisser cosmic ray knee correction~\cite{Gaisser:2013bla}. The probability of passing the atmospheric self-veto $P_{\rm{veto}}$ was calculated using~\cite{Gaisser:2014bja}.

Finally, the energy dependence of the atmospheric muon likelihood uses the empirical parametrization given by Vincent et al.~\cite{Vincent:2016nut} and the angular distribution from Aartsen et al.~\cite{Aartsen:2015zva}:
\begin{equation}
P_{\mu,i}(t_i,E_i,\vec x_i) = \frac{\int dE_{\mu} R(E,E_{\rm{dep}}) P_{\mu}(E_\mu)P_{\mu}(\vec x)\delta_{t,{\rm{tr}}}}{\int dE_{\rm{dep}} \int d^2\vec{x} \int dE_{\mu} P_{\mu}(E_\mu)P_{\mu}(\vec x)\delta_{t,{\rm{tr}}}},
\end{equation}
where the Kroenecker delta $\delta_{t,{\rm{tr}}}$ ensures only tracks can be counted as atmospheric muons.

To account for the finite possibility that track events have in fact been misclassified as showers, we perform the following substitution for each observed shower event:
\begin{equation}
P_{i,t = {\rm{sh}}} \rightarrow 0.9 P_{i,t={\rm{sh}}} + 0.1 P_{i,t={\rm{tr}}}.
\end{equation}

\section{Appendix B: Changes in incoming astrophysical power law\label{sec:gamma_var}}
In this appendix we show the results of Fig. 3
, fixing the spectral index of the extragalactic isotropic flux to three different values: $E^{-2}$, corresponding to the prediction from Fermi acceleration, $E^{-2.6}$, the best fit from \cite{Aartsen:2015zva}, and $E^{-2.9}$ in line with latest fits such as~\cite{Vincent:2016nut,ClaudioTalk}.
\begin{figure}
\begin{tabular}{c c}
\includegraphics[width=0.5\textwidth]{gconstraints_vec.pdf}& \includegraphics[width=0.5\textwidth]{gconstraints_SF.pdf} \\
\includegraphics[width=0.5\textwidth]{gconstraints_vec26.pdf}& \includegraphics[width=0.5\textwidth]{gconstraints_SF26.pdf} \\
\includegraphics[width=0.5\textwidth]{gconstraints_vec29.pdf}& \includegraphics[width=0.5\textwidth]{gconstraints_SF29.pdf}
\end{tabular}
\caption{Same as Fig.
, but for different fixed values of the spectral index $\gamma$ of the isotropic extragalactic flux $\phi_0 \propto E^{-\gamma}$. Top: $\gamma = 2$; middle: $\gamma = 2.6$; bottom: $\gamma = 2.9$. Left: fermionic DM with a vector mediator; right: scalar DM with a fermionic mediator.}
\end{figure}

\end{document}