# Generalizing GradCAM for Embedding Networks

Mudit Bachhawat
<EMAIL_ADDRESS>
###### Abstract

Visualizing CNN is an important part in building trust and explaining model’s
prediction. Methods like CAM[22] and GradCAM[14] have been really successful
in localizing area of the image responsible for the output but are only
limited to classification models. In this paper, we present a new method
EmbeddingCAM, which generalizes the GradCAM for embedding networks. We show
that for classification networks, EmbeddingCAM reduces to GradCAM. We show the
effectiveness of our method on CUB-200-2011 dataset and also present
quantitative and qualitative analysis on the dataset.

Figure 1: Sample results generated using mean proxy method

##  1 Introduction

Convolutional neural networks are widely used for image classification, image
captioning, face recognition, representation learning, open set
classification, visual recommender system. Generally neural networks are
considered as black box and the output generated by the model is always in
question. Thus, visualization tools are very important to build confidence and
trust in the neural networks. Neural networks generally trained end-to-end
makes it hard to debug and explain the output generated. For example: a
product recommender system must know which parts of the image the model is
looking at to generate the recommendations. A bird classification model must
look at the bird area to generate the predictions instead to any background
noise to generate the prediction. This kinds of explanations helps in building
trust in the neural networks.

Interpobaility methods like [22, 14, 2, 20] are widely used for classification
and tasks with classification branches. But these methods fail to generalize
on embedding network. Embedding networks outputs embeddings instead of
discrete class scores. Embeddings, generally in context of neural networks,
are learned vector of continuous variables, which generally represents a
discrete class. Embeddings are useful as they captures low dimensional rich
representation of otherwise high dimentional data. For example the task of
capturing of the body pose can be represented by a embedding [10] rather than
representing it in discrete classes (as body pose is a continuous measure and
the scale of unique poses can be huge).

Few methods have tried to visualize embedding networks. Methods like [21] have
tried to generate visualization heatmaps on embedding network, but have
trained a classifier on top of embeddings. This type of classification head is
not generally available in embedding networks. Another method by [15]
successfully generates heatmap to focus areas in images, which are responsible
for similarity between pair of images. This method requires a pair of images,
which is not a requirement of a embedding network. Chen et al. [3] provides a
method to generate localized heatmap for single image, but requires sampling
of multiple triplets for loss. Additionally, the images needs to be selected
properly as their embeddings should violate the maximum condition for triplet
loss to flow. Our method doesn’t require any such sampling and works with
single image.

Our method EmbeddingCAM generalizes GradCAM and is analogous to GradCAM for
classification type networks. For classification, GradCAM generates class
discriminative map, focusing regions of images which are responsible in
generating the score for that class. The generalization of GradCAM is hard for
embedding network as embeddings cannot be directly used for back-propagating
without the class scores. We solve this problem by using proxies as class
representative and following a similar method as GradCAM. Under special
conditions, i.e. when embeddings are substituted with one-hot encoded vectors,
we show that EmbeddingCAM reduces to GradCAM. In a way, embedding is a
generalization of classification vector, just trained on different loss. Our
method is well represented in figure 2.

In summary our work has following major contribution

  * •

We provide a method to generate a GradCAM type heatmap to highlight regions of
images, which the network is focusing at, to generate the embedding

  * •

We provide a generic method of EmbeddingCAM, which reduces to GradCAM when
adjusted for classification task

  * •

We conduct experiments open set classification method, to show the working of
our method, by generating heatmaps on CUB-200-2011 dataset

  * •

We provide qualitative and quantitative metrics of our method on CUB-200-2011
dataset. We also show the comparison with previous methods

##  2 Related Works

In this section, we explore related works on visual embedding networks and
interopability of neural networks.

Visual Embedding Networks: Visual embedding networks outputs a fixed sized
embedding for a given input image. The network can be then trained using loss
which uses relative distance between these embedding. Additionally, to add
boundedness into the network output, $L_{2}$ normalization is generally added
and cosine similarity is widely used. Such models are trained using Triplet
Loss or Siamese Loss. Finding new losses and efficiently training embedding
networks is an active area of research [11][17][8]. Once the network is
trained, our method can work on any network, trained on any loss.

Most of the methods to train these networks is either a pair-based method or a
proxy-based method. Pair-based methods minimizes the relative distance of the
images in the embedding metric space. For example, triplet loss [18, 13]
samples three images, two from same class and one from different class, and
minimises the relative distance of samples from same class to samples from
different classes. Another example, contrastive loss [6, 4], which aims to
minimize the distance between a pair of data if their class labels are
identical and to separate them otherwise.

Pair-based methods requires sampling of multiple images during training time.
This increases the training complexity and thus creates a need for efficient
mining of hard and semi hard negatives. Proxy-based methods [8] [11] [17] uses
proxy embedding to optimize the sampling complexity and improve the training
process. These methods use proxies to compute the loss between proxy and
embedding, thus preventing use of any other embeddings.

Common methods to train a these models is using triplet loss [18, 13, 9].
Generally these methods use relative distance between images to supervise the
model. Common use cases for these systems includes open set classification,
visual search engine, etc.

Explaining Deep Networks: CAM[22] and Grad-CAM[14] are widely used methods of
visualizing classification networks. These methods highlights important areas
in image which are responsible for generating the class scores for a specific
class. Limitation: These methods are only limited to classification type
networks, and not apply for embedding type networks.

GradCAM [14] and related methods [2] are widely used for visualizing
classification type networks. GradCAM backpropagates the gradient of class
score ($y_{c}$) to the last convolutional layer and reduces those gradient
weights with last convolutional layer output to generate final heatmap. For
simple networks (i.e. Conv -¿ Avg Pool -¿ FC), [14] has shown GradCAM reduces
to CAM [22], which uses fully connected layer weights to generate heatmap.

Visualizing Embedding Networks: Chen et al. [3] proposed a method to generate
Grad-CAM type heatmap for embedding network. Their method requires multiple
training examples to first compute gradient weights for training samples.
Then, during testing phase, they find the nearest training sample to the
testing sample, and use gradient weight of the training sample to generate the
heatmap. This method require multiple images and requires to build an index to
generate heatmap for any arbitrary image.

Figure 2: Diagram showing our method for generating heatmap from embedding
models

##  3 Method

In this section, we go through the method of generating a Grad-CAM style
heatmap to highlight important parts in image. A high level visualization of
our method can be found in figure 2. We use proxy as class embedding to
compute a custom loss[9]. Then, similar to GradCAM we use back-propagated
weights to generate final heatmap. Since the loss depends the proxy used, we
explore couple of method to generate the proxy. The proxy can be imagined as a
ground truth vector, i.e. a vector for the concept we are trying to visualize.
More on this in section 3.2.1.

###  3.1 Preliminary

This section covers the preliminary topics on which we build our method.

####  3.1.1 Visual Metric Learning

Visual metric learning is the task of learning a function $f:I\rightarrow S$
which maps an image $I$ to a semantic space $S=\mathbb{R}^{d}$. This kind of
formulation allows to capture rich representation in the data. The embedding
of semantically similar images should be closer and embeddings of semantically
dissimilar images should be farther in the metric space $\mathbb{R}^{d}$.
Generally the function $f$ is replaced with a convolutional neural network
which outputs an embedding in $\mathbb{R}^{d}$. Euclidean distance is used on
the $L_{2}$ normalized embedding to measure distance between points.

$y=\mathrm{CNN}(I)$ (1)

$\hat{y}=\frac{y}{\left\lVert y\right\rVert}$ (2)

where $\left\lVert.\right\rVert$ is $L_{2}$ norm and
$\hat{y}\in\mathbb{R}^{d}$

####  3.1.2 GradCAM

GradCAM [14] generates a class discriminative map (heatmap) for classification
type network. It uses gradient of the class score, to calculate the area of
image responsible the output score. It has been used to generate class local
heatmaps for various tasks like image classification [14, 7, 5], visual
question answering [12] and image captioning [16].

Mathematically, Let $Y_{c}$ be the class score (pre softmax) of the
discriminative class. We then compute the gradients of $Y_{c}$ with respect to
last convolutional layer of the network $A^{k}_{i,j}$ (generally before global
average pooling)

$g_{c}\left(A^{k}_{i,j}\right)=\frac{\partial y^{c}}{\partial A^{k}_{i,j}}$
(3)

where $(i,j)$ are the spatial indexes of the tensor and Z is the total number
of spatial pixels in each channel. We then reduce the gradient tensor across
spatial dimensions, to generate channel weights.

$\alpha_{k}^{c}=\frac{1}{Z}\sum_{i}\sum_{j}\frac{\partial y^{c}}{\partial
A_{i,j}^{k}}$ (4)

We finally reduce the last convolutional layer output using the channel
weights obtained from equation 4.

$H_{\mathrm{Grad}-\mathrm{CAM}}^{c}=\operatorname{ReLU}\left(\sum_{k}\alpha_{k}^{c}A^{k}\right)$
(5)

Computations for simple networks: For simple networks, i.e. final
convolutional layer followed by global average pooling, followed by dense
layer, followed by softmax layer for classification. The equation 4 reduces to

$\alpha_{k}^{c}=\frac{1}{Z}w_{k,c}$ (6)

Where $w\in\mathbb{R}^{K\times C}$ is the kernel of dense layer, with number
of channels in last convolutional layer K and number of classes C. Ignoring
the $1/Z$ factor, which will be normalized during the visualization, the
$\alpha_{k}^{c}$ is identical to $w_{k,c}$. The proof of equation 6 can be
found in [14]. The channel weights can be arranged in vector as

$\displaystyle\alpha^{c}$
$\displaystyle\propto(w_{1,c},w_{2,c},...,w_{K,c})^{T}$ (7)
$\displaystyle\propto w\cdot o^{c}$ (8)

Where $o^{c}$ is one-hot encoded vector for class $c$. In the following
chapters, we will show how $o^{c}$ can be visualized as proxy vector of
EmbeddingCAM.

###  3.2 Generalizing GradCAM for Metric Learning

As mentioned in [3], applying GradCAM to metric learning have two major
challenges. First the is no class labels ($Y_{c}$). In GradCAM, the gradients
for convolution output is calculated with respect to the class score. Since
most embedding networks are trained using relative distances between
embeddings, there is no notion of class score. Secondly, because of there is
no notion of class, it becomes extremely hard to compute the gradient with one
sample. You will need to take multiple samples to compute the gradient.

Our method solves both the problem by using proxies as class representative.
Formally, let $p_{c}\in\mathbb{R}^{d}$ be the proxy of class c with dimension
d. We define the process of calculating proxies in section 3.2.1. We further
define following as the loss function, for calculating the gradients.

$\mathcal{L}_{c}=y\cdot p_{c}$ (9)

where ’$\cdot$’ represents dot product and $y\in\mathbb{R}^{d}$ is
unnormailzed CNN output. This is analogous to using pre-softmax class score
$Y_{c}$ in GradCAM. Finally, we use $\mathcal{L}_{c}$ to calculate gradients
similar to 3.

$g_{c}\left(A^{k}_{i,j}\right)=\frac{\partial\mathcal{L}_{c}}{\partial
A^{k}_{i,j}}$ (10)

We use equation 4 and 5 to calculate the final heatmap.

####  3.2.1 Calculating Proxies

We explore two schemes for calculating proxy of class C.

Normalized Mean Proxy: For any class c with images ${i_{1},i_{2},...i_{n}}$,
we calculate average proxy for class c as

$p_{c}=\frac{\frac{1}{n}\sum_{j=i}^{n}\hat{y}_{j}}{\left\lVert\frac{1}{n}\sum_{j=i}^{n}\hat{y}_{j}\right\rVert}$
(11)

where

$\hat{y_{j}}=\frac{CNN(i_{j})}{\left\lVert CNN(i_{j})\right\rVert}$ (12)

Single Point Proxy: Single Point Proxy is a special case of normalized mean
proxy where we use only one point i.e. the image embedding itself as class
representative. This method solves the problem of sampling multiple images to
generate a heatmap.

####  3.2.2 Intuition

In GradCAM, we calculate the gradient of class score $y_{c}$ w.r.t the final
convolutional output. The gradients in GradCAM highlights the area, when
magnified, increases the class score. Thus $\alpha_{k}^{c}$ represents
direction of each channel in final convolutional output, which can lead to
increment in the class score $y_{c}$. This direction can be also visualised as
weights for each channel in the generation of final heatmap. In this process
the gradients create class discriminative activation map for the class.

Similar to this, the equation 9 represents the agreement between the model
output and the class proxy. If we think of class proxy as the ground truth
vector of the concept for which we are trying generate visualization, the
equation 9 is analogous to class score of the Grad-CAM.

The equation 9 can be written as $|y|cos\theta$, where $\theta$ is the angle
between model output vector and proxy vector. Thus the region highlighted by
the gradient of this loss is responsible for increasing the model output
confidence ($|y|$) and agreement between proxy and model output vector.
Additionally, if y is an output of classification network with number of
classes $C$ and $p_{c}$ is an one hot encoded vector for class $c$, the
equation 9 reduces to class score of class $c$. This reduces EmbeddingCAM to
GradCAM.

###  3.3 Reduction for Simple Networks

For simple image classification type networks (networks ending with
convolutional layer, followed by average pooling, fully connected and L2
normalization layer to generate final embedding), we simplify the network
calculations and draw similarity between GradCAM and EmbeddingCAM. Let
$w\in\mathbb{R}^{K\times C}$ be the kernel of last fully connected layer.
Continuing from equation 9:

$\frac{\partial L_{c}}{\partial y}=p_{c}$ (13)

Back-propagating the loss through fully connected layer and average pooling
layer, we get following as gradients:

$\frac{\partial\mathcal{L}_{c}}{\partial
A^{k}_{i,j}}=\frac{1}{Z}\sum_{x=1}^{C}w_{k,x}*p_{c,x}$ (14)

Since the above gradient is independent of $i$ and $j$, $\alpha_{k}^{c}$ can
be reduced to 16

$\displaystyle\alpha_{k}^{c}$
$\displaystyle=\frac{1}{Z}\sum_{i}\sum_{j}\frac{\partial\mathcal{L}_{c}}{\partial
A_{i,j}^{k}}$ (15) $\displaystyle=\frac{1}{Z}\sum_{x=1}^{C}w_{k,x}*p_{c,x}$
(16)

Finally, we vectorise the equation 16 to generate 17.

$\alpha^{c}=\frac{1}{Z}w\cdot p_{c}$ (17)

$\alpha^{c}\propto w\cdot p_{c}$ (18)

We again ignore the $1/Z$ factor in equation 17. This factor will get
normalized during heatmap generation. Please note that equation 18 is similar
to 8.

###  3.4 Comparison with other methods

Method [3] requires selective sampling of multiple triplets for a single image
to generate the heatmap. Additionally, this methods uses training data weights
during inference phase. This, the method is not generic and cannot work with
just pretrained model. Our method can work with just pretrained model and can
work with single image (when using single point proxy scheme).

Visual Attention Accuracy
---
| Mean Heatmap Ratio | Accuracy (IoU @ 0.5)
|  | Uniform | All Channels | Top-50 Channels
Chen et al. [3] | BBox | 0.543 | 0.643 $\pm$ 0.265 | 0.760 $\pm$ 0.152 | 0.797
Mask | 0.275 | 0.416 $\pm$ 0.231 | 0.534 $\pm$ 0.146
Our Method (Single Point Proxy) | BBox | 0.522 | 0.776 $\pm$ 0.004 | 0.757 $\pm$ 0.003 | 0.776
Mask | 0.274 | 0.543 $\pm$ 0.004 | 0.520 $\pm$ 0.004
Our Method (Average Proxy) | BBox | 0.522 | 0.789 $\pm$ 0.005 | 0.764 $\pm$ 0.004 | 0.768
Mask | 0.274 | 0.551 $\pm$ 0.006 | 0.524 $\pm$ 0.005
Table 1: Comparison of our method to [3]. The accuracy reported is weakly-
supervised localization accuracy with<EMAIL_ADDRESS>We train 8 models with
different initializations to provide standard deviation of the score.

##  4 Experiments

###  4.1 Dataset: CUB200-2000

We use CUB200-2011 dataset for the evaluation of our method. The dataset
contains 11.8K images of birds across 200 classes. This dataset is commonly
used for fine grained classification, open set classification, and visual
explanations. We trained a ResNet-50 embedding model following the method
details by [19], which outputs embedding of size 128. We use the output of
last convolutional layer (before average pooling layer) of size $7\times
7\times 2048$, for back-propagating the loss and generating the heatmap. The
model is trained on first 100 classes and evaluated on next 100 classes.

Evaluation of the method: Our method takes input a pretrained model which is
used to generates heatmaps on the complete dataset. The heatmaps generated are
evaluated using the metrics defined below. Please note that bounding boxes and
segmentation maps are only used during the evaluation phase, and are not used
during the training of the model.

Metrics: We use the mean heatmap ratio metric to evaluate our method. Mean
heatmap ratio is the ratio of the heatmap score inside the bounding box (or
segmentation map). We calculate this by calculating ratio of heatmap inside
the bounding box (or segmentation map) and then taking mean over all images.
We also present a uniform baseline, which calculates this metric on a constant
heatmap (heatmap that have constant non zero activation in its spatial
dimensions). For a heatmap which has activation on bird area, the heatmap
ratio will be one for both bounding box and segmentation map.

We also report weakly supervised localization accuracy as reported in [3].
This metric requires computation of bounding box from the heatmap and
calculates intersection over union (IoU @ 0.5) accuracy with the ground truth
bounding boxes. This metric measures the bounding box prediction with weakly
supervised training of the model. This metrics is widely used for visual
interpretation tasks [3][1].

To generate the bounding box from heatmap, we first binarize the heatmap with
a threshold $t$. We then find the maximum connected component in binarized
heatmap. Finally we take the smallest bounding box that completely encloses
the largest connected component.

####  4.1.1 Results

We compare our method with [3], which also reports mean heatmap ratio and
weakly supervised localization accuracy. We have attempted our best to keep
the environment same while generating the results. Please note that bounding
boxes and segmentation maps are only used during evaluation phase and not
during training of the model.

Mean Heatmap Ratio: We calculate mean heatmap ratio for all classes using
bounding box and segmentation map. The report our metrics in table 1. Our
method both single point proxy and average point proxy), both performs
competitive to [3] without the need for sampling multiple triplets and test-
time indexing. Our method does not require sampling and thus mean heatmap
ratios are stable (doesn’t have variance) compared to [3]. We additionally
train 8 models with different random initialisation to test the stability of
our model. We report the standard deviation of mean heatmap ratio in the
table. Our single point proxy achieves mean heatmap ratio of 0.776 and 0.543
for bounding box and segmentation map respectively. We observe the standard
deviation of our methods to be less than 1%, demonstrating stability of our
approach.

Weakly Supervised Localization Accuracy: We report weakly supervised
localization accuracy in table 1 using the threshold $t=0.2$. We have found
the accuracy to be stable ($>0.6$) for all values of $t\in[0.05,0.4]$. We
calculate all our metrics

Qualitative Analysis: Figure 1 presents sample images with their heatmaps
generated by our method. We observe that most of our heatmaps are generated on
relevant parts of birds which uniquely identifies the species.

##  5 Conclusion

In this paper, we present a novel approach to explain the prediction of any
metric learning model. A wide variety [22][14][2] of methods have been
developed to explain the prediction of a classification system. But, as the
application of metric learning are increasing there are no simple and
effective method to explain the prediction of the system.

In this paper, we present a novel method to generate GradCAM style heatmap
from just a pretrained metric learning model. We do this without the use of
multiple images and selective sampling of multiple images. We are able to
achieve this by using proxies as class representative and and defining a
custom loss (which can be used for backpropagation, similar to GradCAM).
Finally, we show the working of our method in the results section. We show
that we have achieved comparable mean heatmap ratio and weakly supervised
localization accuracy in comparison to [3] on CUB200-2011 dataset. We also
present some qualitative analysis of the results generated.

The recent advancements in the field of AI has greatly enhanced the accuracy
of the tasks. But, deep learning method primarily being a black box system has
caused havoc when deployed in production. To avoid such incidents, it is
extremely important for ML models to be a white box rather than black box and
developing an explainable AI system plays a huge role in this. We encourage
researchers to apply this method on other embedding networks share the
effectiveness of the approach.

## References

  * [1] S. Belharbi, A. Sarraf, M. Pedersoli, I. Ben Ayed, L. McCaffrey, and E. Granger.  F-cam: Full resolution class activation maps via guided parametric upscaling.  In WACV, 2022.
  * [2] Aditya Chattopadhyay, Anirban Sarkar, Prantik Howlader, and Vineeth N Balasubramanian.  Grad-cam++: Generalized gradient-based visual explanations for deep convolutional networks.  arXiv preprint arXiv:1710.11063, 2017.
  * [3] Lei Chen, Jianhui Chen, Hossein Hajimirsadeghi, and Greg Mori.  Adapting grad-cam for embedding networks.  In Proceedings of the IEEE/CVF Winter Conference on Applications of Computer Vision (WACV), March 2020.
  * [4] Sumit Chopra, Raia Hadsell, and Yann LeCun.  Learning a similarity metric discriminatively, with application to face verification.  In Computer Vision and Pattern Recognition, 2005. CVPR 2005. IEEE Computer Society Conference on, volume 1, pages 539–546. IEEE, 2005.
  * [5] Elias Ennadifi, Sohaib Laraba, Damien Vincke, Benoît Mercatoris, and Bernard Gosselin.  Wheat diseases classification and localization using convolutional neural networks and gradcam visualization.  In 2020 International Conference on Intelligent Systems and Computer Vision (ISCV), pages 1–5, 2020.
  * [6] R. Hadsell, S. Chopra, and Y. LeCun.  Dimensionality reduction by learning an invariant mapping.  In 2006 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR’06), volume 2, pages 1735–1742, 2006.
  * [7] Jeremy Irvin, Pranav Rajpurkar, Michael Ko, Yifan Yu, Silviana Ciurea-Ilcus, Chris Chute, Henrik Marklund, Behzad Haghgoo, Robyn L. Ball, Katie S. Shpanskaya, Jayne Seekins, David A. Mong, Safwan S. Halabi, Jesse K. Sandberg, Ricky Jones, David B. Larson, Curtis P. Langlotz, Bhavik N. Patel, Matthew P. Lungren, and Andrew Y. Ng.  Chexpert: A large chest radiograph dataset with uncertainty labels and expert comparison.  CoRR, abs/1901.07031, 2019.
  * [8] Sungyeon Kim, Dongwon Kim, Minsu Cho, and Suha Kwak.  Proxy anchor loss for deep metric learning.  In IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), June 2020.
  * [9] Gregory Koch, Richard Zemel, and Ruslan Salakhutdinov.  Siamese neural networks for one-shot image recognition.  2015\.
  * [10] Greg Mori, Caroline Pantofaru, Nisarg Kothari, Thomas Leung, George Toderici, Alexander Toshev, and Weilong Yang.  Pose embeddings: A deep architecture for learning to match human poses.  CoRR, abs/1507.00302, 2015.
  * [11] Yair Movshovitz-Attias, Alexander Toshev, Thomas Leung, Sergey Ioffe, and Saurabh Singh.  No fuss distance metric learning using proxies.  pages 360–368, 10 2017.
  * [12] Badri N. Patro, Mayank Lunayach, Shivansh Patel, and Vinay P. Namboodiri.  U-CAM: visual explanation using uncertainty based class activation maps.  CoRR, abs/1908.06306, 2019.
  * [13] Florian Schroff, Dmitry Kalenichenko, and James Philbin.  Facenet: A unified embedding for face recognition and clustering.  CoRR, abs/1503.03832, 2015.
  * [14] R. R. Selvaraju, A. Das, R. Vedantam, M. Cogswell, D. Parikh, and D. Batra.  Grad-CAM: Why did you say that? Visual Explanations from Deep Networks via Gradient-based Localization.  CoRR, abs/1610.02391, 2016.
  * [15] Abby Stylianou, Richard Souvenir, and Robert Pless.  Visualizing deep similarity networks.  CoRR, abs/1901.00536, 2019.
  * [16] Jiamei Sun, Sebastian Lapuschkin, Wojciech Samek, and Alexander Binder.  Understanding image captioning models beyond visualizing attention.  CoRR, abs/2001.01037, 2020.
  * [17] Eu Wern Teh, Terrance DeVries, and Graham W Taylor.  Proxynca++: Revisiting and revitalizing proxy neighborhood component analysis.  arXiv preprint arXiv:2004.01113, 2020.
  * [18] Chao-Yuan Wu, R. Manmatha, Alexander J. Smola, and Philipp Krähenbühl.  Sampling matters in deep embedding learning.  CoRR, abs/1706.07567, 2017.
  * [19] Chao-Yuan Wu, R Manmatha, Alexander J Smola, and Philipp Krähenbühl.  Sampling matters in deep embedding learning.  In ICCV, 2017.
  * [20] Matthew D. Zeiler and Rob Fergus.  Visualizing and understanding convolutional networks.  CoRR, abs/1311.2901, 2013.
  * [21] Meng Zheng, Srikrishna Karanam, Ziyan Wu, and Richard J. Radke.  Re-identification with consistent attentive siamese networks.  CoRR, abs/1811.07487, 2018.
  * [22] Bolei Zhou, Aditya Khosla, Àgata Lapedriza, Aude Oliva, and Antonio Torralba.  Learning deep features for discriminative localization.  CoRR, abs/1512.04150, 2015.


# Hints for decaying dark matter from $S_{8}$ measurements

Guillermo F. Abellán<EMAIL_ADDRESS>Riccardo Murgia
<EMAIL_ADDRESS>Vivian Poulin<EMAIL_ADDRESS>Julien Lavalle<EMAIL_ADDRESS>Laboratoire Univers & Particules de
Montpellier (LUPM), Université de Montpellier (UMR-5299)
Place Eugène Bataillon, F-34095 Montpellier Cedex 05, France

###### Abstract

Recent weak lensing surveys have revealed that the direct measurement of the
parameter combination $S_{8}\equiv\sigma_{8}(\Omega_{m}/0.3)^{0.5}$ –
measuring the amplitude of matter fluctuations on 8 $h^{-1}$Mpc scales – is
$\sim 3\sigma$ discrepant with the value reconstructed from cosmic microwave
background (CMB) data assuming the $\Lambda$CDM model. In this Letter, we show
that it is possible to resolve the tension if dark matter (DM) decays with a
lifetime of $\text{log}_{10}(\Gamma^{-1}/\text{Gyr})=1.75_{-0.95}^{+1.4}$ into
one massless and one massive product, and transfers a fraction
$\varepsilon\simeq 0.7^{+2.7}_{-0.6}\%$ of its rest mass energy to the
massless component. The velocity-kick received by the massive daughter leads
to a suppression of gravitational clustering below its free-streaming length,
thereby reducing the $\sigma_{8}$ value as compared to that inferred from the
standard $\Lambda$CDM model, in a similar fashion to massive neutrino and
standard warm DM. Contrarily to the latter scenarios, the time-dependence of
the power suppression and the free-streaming scale allows the 2-body decaying
DM scenario to accommodate CMB, baryon acoustic oscillation, growth factor and
uncalibrated supernova Ia data. We briefly discuss implications for DM model
building, galactic small-scale structure problems and the recent Xenon-1T
excess. Future experiments measuring the growth factor to high accuracy at
$0\lesssim z\lesssim 1$ can further test this scenario.

††preprint: LUPM:20-026

Introduction$-$ The standard $\Lambda$-cold dark matter ($\Lambda$CDM)
cosmological model provides a remarkable fit to a wide variety of observables,
such as big bang nucleosynthesis (BBN), the cosmic microwave background (CMB),
large-scale structures (LSS), baryonic acoustic oscillations (BAO), and
uncalibrated supernovae of type Ia (SNIa) (see e.g. Scott (2018); Akrami et
al. (2018) for reviews). Nevertheless, tremendous experimental developments
have revealed curious discrepancies between different probes. At the heart of
this study is the growing tension between the cosmological and local
determination of the amplitude of the matter fluctuations on 8 $h^{-1}$/Mpc
scales, typically described through the parameter combination
$S_{8}\equiv\sigma_{8}(\Omega_{m}/0.3)^{0.5}$. Within $\Lambda$CDM, the latest
$S_{8}$ value inferred from a fit to CMB data Aghanim et al. (2018) is $\sim
2-3\sigma$ higher than that measured by a host of weak lensing surveys such as
CFHTLenS Heymans et al. (2013), HSC Hikage et al. (2019), DES Abbott et al.
(2018) and KiDS+Viking Hildebrandt et al. (2020). In particular, the recent
joint analysis of KIDS1000+BOSS+2dfLenS has yielded
$S_{8}=0.766^{+0.020}_{-0.014}$ Heymans et al. (2020), in $\sim 3\sigma$
discrepancy with $\Lambda$CDM from Planck. While an unknown systematic effect
at the origin of this discrepancy is not excluded, the existence of several
independent observations disfavoring the $\Lambda$CDM predictions strengthen
the case for new physics.

In this Letter, we show that the $S_{8}$ tension can be resolved if DM
experiences 2-body decays where the decay products are one massive warm DM
(WDM) particle and one (massless) dark radiation (DR) component. We will refer
to the full model as $\Lambda$DDM. We find that it requires DM to have a
lifetime of $\text{log}_{10}(\Gamma^{-1}/\text{Gyr})=1.75_{-0.95}^{+1.4}$ and
transfers a fraction $\epsilon\simeq 0.7^{+2.7}_{-0.6}\%$ of its rest mass
energy into the DR species. Interestingly, depending on the velocity-kick
received by the massive daughter, this scenario could help resolving some of
the sub-galactic scales issues in $\Lambda$CDM (e.g. Refs. Lin et al. (2001);
Sánchez-Salcedo (2003); Cembranos et al. (2005); Kaplinghat (2005); Strigari
et al. (2007); Borzumati et al. (2008); Peter et al. (2010); Peter and Benson
(2010); Wang et al. (2014)) or could explain the recent Xenon-1T excess The
Xenon-1t Collaboration et al. (2020); Kannike et al. (2020).

Many authors have attempted to explain the $S_{8}$ tension through new
properties of DM (see, e.g., Lesgourgues et al. (2016); Chacko et al. (2016);
Murgia et al. (2016); Kobayashi et al. (2017); Di Valentino et al. (2018);
Buen-Abad et al. (2018); Raveri et al. (2017); Poulin et al. (2019); Lin et
al. (2019); D’Eramo et al. (2018); Dutta et al. (2020); Archidiacono et al.
(2019); Vattis et al. (2019); Bohr et al. (2020); Heimersheim et al. (2020)).
Scenarios where the DM decays only into DR have been discussed in this
context, but have been shown to be at odds with the latest Planck CMB lensing
and BAO data Aoyama et al. (2014); Enqvist et al. (2015); Chudaykin et al.
(2016); Poulin et al. (2016); Bringmann et al. (2018). The extension of these
studies to the case of a massive daughter was recently performed in Ref.
Vattis et al. (2019), where it was suggested that this scenario could resolve
the ‘Hubble tension’ – The $\sim 5\sigma$ discrepancy between the value of the
current expansion rate of the universe inferred from Planck data Aghanim et
al. (2018) under $\Lambda$CDM, and that measured using the cosmic distance
ladder Riess (2019); Verde et al. (2019); Aylor et al. (2019); Wong et al.
(2019); Freedman et al. (2019). However, a recent series of analysis has shown
that a combination of BAO, uncalibrated SNIa data-sets Haridasu and Viel
(2020) and Planck data Clark et al. (2020) excludes this model. Yet, authors
of Refs. Vattis et al. (2019); Haridasu and Viel (2020); Clark et al. (2020)
have limited their analyses to the $\Lambda$DDM background evolution, or used
crude approximations to deal with linear perturbations of the massive daughter
particles.

In this work, we perform the first thorough analysis of the $\Lambda$DDM
including a realistic treatment of linear cosmological perturbations. Our
careful treatment of the warm daughter perturbations allows us to pin down the
space of parameters resolving the $S_{8}$-tension. Indeed, the warm component
produced by decay leads to a suppression of the matter power spectrum at late
times, similar to that of massive neutrinos or standard WDM. However,
contrarily to the latter scenarios, the specific time-dependence of the power
suppression imprinted by the decay allows to accommodate Planck, BAO,
uncalibrated SN1a and S8 measurements, though we confirm that it cannot
simultaneously resolve the Hubble tension. Finally, we briefly discuss
implications of these results for DM model building, Xenon-1T and the ‘small-
scale crisis of CDM’.

Cosmology of 2-body DCDM$-$ Our framework is characterized by two additional
free parameters with respect to $\Lambda$CDM: the DCDM lifetime,
$\Gamma^{-1}$, and the fraction of DCDM rest mass energy converted into DR,
defined as follows Blackadder and Koushiappas (2014):

$\varepsilon=\frac{1}{2}\left(1-\frac{m^{2}_{\rm{wdm}}}{m^{2}_{\rm{dcdm}}}\right),$
(1)

where $0\leq\varepsilon\leq 1/2$. The lower limit corresponds to the standard
CDM case, whereas $\varepsilon=1/2$ corresponds to DM decaying solely into DR.

The general form of the background and linear evolution equations for the
DCDM, WDM and DR components can be found in Ref. Wang and Zentner (2012);
Aoyama et al. (2014). The dynamics of linear perturbations in the DR component
is greatly simplified by integrating its phase-space distribution (PSD) over
all the momentum degrees of freedom Poulin et al. (2016) (after expanding it
over Legendre polynomials). To compute the WDM dynamics one cannot apply the
same strategy, since the integral over the phase space is not analytic. One
should follow the evolution of the full time- and momentum-dependent PSD,
which would require solving $\mathcal{O}(10^{8})$ linear differential
equations for realistic cosmological analyses. We tackle this issue by
devising a new fluid approximation for the WDM species (we demonstrate its
accuracy in Ref. F. Abellán et al. (in prep.)), which allows to integrate out
the dependency on momenta and reduces the hierarchy of equations to the first
three multipoles. The novel approximation scheme, based on the treatment of
massive neutrinos as a viscous fluid by Refs. Lesgourgues and Tram (2011);
Blas et al. (2011), is only valid on sub-Hubble scales, where high- and
low-$\ell$ modes are effectively decoupled. In the synchronous gauge comoving
with the DCDM species (see Ref. Ma and Bertschinger (1995) for all relevant
definitions), it yields the following expressions for the WDM continuity,
Euler and shear equations Wang and Zentner (2012); F. Abellán et al. (in
prep.):

$\displaystyle\dot{\delta}_{\rm{wdm}}$
$\displaystyle=-3\mathcal{H}(c^{2}_{g}-w)\delta_{\rm{wdm}}-(1+w)\left(\theta_{\rm{wdm}}+\frac{\dot{h}}{2}\right)$
(2)
$\displaystyle+(1-\varepsilon)a\Gamma\frac{\bar{\rho}_{\rm{dcdm}}}{\bar{\rho}_{\rm{wdm}}}(\delta_{\rm{dcdm}}-\delta_{\rm{wdm}}),$
$\displaystyle\dot{\theta}_{\rm{wdm}}$
$\displaystyle=-\mathcal{H}(1-3c_{g}^{2})\theta_{\rm{wdm}}+\frac{c^{2}_{g}}{1+w}k^{2}\delta_{\rm{wdm}}-k^{2}\sigma_{\rm{wdm}}$
(3)
$\displaystyle-(1-\varepsilon)a\Gamma\frac{1+c_{g}^{2}}{1+w}\frac{\bar{\rho}_{\rm{dcdm}}}{\bar{\rho}_{\rm{wdm}}}\theta_{\rm{wdm}}.$
$\displaystyle\dot{\sigma}_{\rm{wdm}}$
$\displaystyle=-3\left(\frac{1}{\tau}+\mathcal{H}\xi+(1-\varepsilon)a\Gamma\frac{\bar{\rho}_{\rm{dcdm}}}{\bar{\rho}_{\rm{wdm}}}\frac{1+c_{g}^{2}}{3(1+w)}\right)\sigma_{\rm{wdm}}$
(4)
$\displaystyle+\frac{8wc_{g}^{2}}{1+w}\left[\theta_{\rm{wdm}}+\frac{\dot{h}}{2}\right]-\frac{2}{3}\frac{\varepsilon^{2}a\Gamma}{(1-\varepsilon)}\frac{\bar{\rho}_{\rm{dcdm}}}{\bar{\rho}_{\rm{wdm}}}\frac{\delta_{\rm{dcdm}}}{(1+w)},$

where $\delta_{\rm{wdm}}$, $\theta_{\rm{wdm}}$ and $\sigma_{\rm{wdm}}$
represent the WDM density, velocity divergence and shear perturbations
respectively, $w\equiv\bar{P}_{\rm{wdm}}/\bar{\rho}_{\rm{wdm}}$ is the WDM
equation of state, $\xi\equiv
2/3-c_{g}^{2}-\mathfrak{p}_{\rm{wdm}}/3\bar{P}_{\rm{wdm}}$ and the ‘dot’
refers to derivative with respect to conformal time. The shear equation is
valid in the relativistic limit $\varepsilon\simeq 0.5$, and for
$\varepsilon\leq 0.49$ we find that simply setting $\dot{\sigma}_{\rm wdm}\to
0$ is accurate at the ${\cal O}(1\%)$ level in the matter power spectrum. As
in Ref. Lesgourgues and Tram (2011), these equations rely on the assumption
that the sound speed in the synchronous gauge $c_{\rm syn}$ is nearly equal to
the adiabatic sound speed $c_{g}$, and the latter can be computed as follows

$\displaystyle c_{g}^{2}$
$\displaystyle=w\left(5-\frac{\mathfrak{p}_{\rm{wdm}}}{\bar{P}_{\rm{wdm}}}-\frac{\bar{\rho}_{\rm{dcdm}}}{\bar{\rho}_{\rm{wdm}}}\frac{a\Gamma}{3w\mathcal{H}}\frac{\varepsilon^{2}}{1-\varepsilon}\right)$
(5)
$\displaystyle\times\left[3(1+w)-\frac{\bar{\rho}_{\rm{dcdm}}}{\bar{\rho}_{\rm{wdm}}}\frac{a\Gamma}{\mathcal{H}}(1-\varepsilon)\right]^{-1},$

where the pseudo-pressure $\mathfrak{p}_{\rm{wdm}}$ is a higher momenta
integral of the background WDM PSD, equivalent to the standard pressure in the
relativistic limit Lesgourgues and Tram (2011).

Figure 1:  Reconstructed 2D posterior distribution in the $\Lambda$DDM and in
the $\nu\Lambda$CDM models.

Resolving the $S_{8}$ tension with DCDM$-$ The 2-body $\Lambda$DDM scenario
under study is fully described by the following set of free parameters:

$\left\\{\Omega_{b}h^{2},\text{ln}\left(10^{10}A_{s}\right),n_{s},\tau_{\rm
reio},\Omega_{\rm{dcdm}}^{\rm{ini}},H_{0},{\mathrm{log}_{10}}~{}\Gamma,{\mathrm{log}_{10}}~{}\varepsilon\right\\}.$
(6)

We implement the DDM equations in the publicly available numerical code CLASS
Lesgourgues and Tram (2011); Blas et al. (2011). We use a shooting method that
guarantees that the budget equation is satisfied. We solve the Boltzmann
hierarchy presented in Ref. Aoyama et al. (2014) using 300 momentum bins. At
sub-Hubble scales ($k\tau>25$), we switch to the fluid equations given by Eqs.
(2)-(4). We make use of the code MONTEPYTHON-v3 Audren et al. (2013);
Brinckmann and Lesgourgues (2019) to perform a Monte Carlo Markov chain (MCMC)
analysis with a Metropolis-Hasting algorithm, testing the $\Lambda$DDM model
against the high-$\ell$ CMB TT, TE, EE ‘lite’+low-$\ell$ TT,EE+lensing data
from Planck Aghanim et al. (2018), BAO data from 6dF Beutler et al. (2011),
SDSS DR7 Ross et al. (2015), BOSS DR12 Alam et al. (2017), eBOSS DR14
Ly-$\alpha$ de Sainte Agathe et al. (2019); Blomqvist et al. (2019) and the
Pantheon SNIa catalog Scolnic et al. (2018). We model the
KIDS1000+BOSS+2dfLenS with a split-normal likelihood111The full likelihood is
not yet available. Making use of $S_{8}=0.755^{+0.019}_{-0.021}$ from
KiDS450+DESY1 Joudaki et al. (2020) yields similar results. on
$S_{8}=0.766^{+0.02}_{-0.014}$ Heymans et al. (2020). We adopt the
uninformative priors
$-4\leq{\mathrm{log}_{10}}(\varepsilon)\leq{\mathrm{log}_{10}}(0.5)$,
$-4\leq{\mathrm{log}_{10}}(\Gamma/[{\rm Gyr^{-1}}])\leq 1$ and
$0\leq\Omega_{\rm{dcdm}}^{\rm{ini}}\leq 1$. To gauge the importance of the
late-time decay in the success of the solution, we compare the $\Lambda$DDM
model with another cosmological scenario that features a power suppression at
small scales, namely massive neutrinos222Mixed C+WDM scenarios induce a
similar effect, but the WDM fraction (present at all times) is constrained by
observations on sub-galactic scales Baur et al. (2017); Diamanti et al.
(2017); Gariazzo et al. (2017); Murgia et al. (2018) and they cannot provide a
lower $S_{8}$. ($\nu\Lambda$CDM). We model these as three degenerate states
and vary the total mass $M_{\nu}$, on top of the standard $\Lambda$CDM
parameters.

Our results are reported in Tab. 1 and summarized in Fig. 1: In the
$\Lambda$DDM scenario (blue contours) we find that the best-fit (when
including the $S_{8}$ prior) has $\varepsilon\simeq 0.7\%$ and
$\Gamma^{-1}\simeq 56$ Gyrs, yielding $S_{8}\simeq 0.77$ and $\Omega_{m}\simeq
0.31$, in excellent agreement with the KiDS1000+BOSS+2dfLenS measurement.
Moreover, the decrease in $S_{8}$ is driven by a smaller $\sigma_{8}$, while
$\Omega_{m}$ is not affected, which is also what is favored by the data. We
find a strong negative correlation between $\varepsilon$ and $\Gamma$, which
approximately scales like $\Gamma^{-1}\simeq 56~{}(\varepsilon/0.007)^{1.4}$
Gyrs. On the other hand, the $\nu\Lambda$CDM model can only achieve
$S_{8}\simeq 0.81$, with $M_{\nu}<0.1614$ (95% C.L.). Remarkably, we find a
$\Delta\chi^{2}_{\rm min}\equiv\chi^{2}_{\rm min}(\nu\Lambda{\rm
CDM})-\chi^{2}_{\rm min}(\Lambda{\rm DDM})\simeq-5.4$ in favor of the
$\Lambda{\rm DDM}$ model. The negative $\Delta\chi^{2}$ is driven by the low
$S_{8}$ value. The fit to other data set is barely affected by the inclusion
of a $S_{8}$ prior in the $\Lambda$DDM model (but degrades in the
$\nu\Lambda$CDM case), such that without the $S_{8}$ prior, there is no
preference for DDM and $S_{8}$ seems unchanged. Looking at the reconstructed
$S_{8}=0.821^{+0.016}_{-0.011}$ in the $\Lambda$DDM model without the prior
information, one might naively expect $S_{8}\simeq 0.77$ to be largely
excluded. In fact, the combined $\chi^{2}$ with the prior on $S_{8}$ only
increases by $\sim+1.5$ (as opposed to $\sim+7$ in the $\nu\Lambda$CDM). This
is because $S_{8}$ has a non-Gaussian posterior with a tail extending to low
values due the degeneracy between $\Gamma$ and $\varepsilon$. This degeneracy
becomes clear when incorporating the $S_{8}$ prior, or keeping either $\Gamma$
or $\varepsilon$ fixed. These issues are discussed in a forthcoming
publication F. Abellán et al. (in prep.).

Figure 2:  Top panel $-$ Residuals in the linear matter power spectrum $P(k)$
at redshifts $z=0,3$, with respect to our baseline $\nu\Lambda$CDM model, for
the best-fit $\Lambda$DDM model (red lines) and a $\nu\Lambda$CDM scenario
yielding the same $S_{8}$ (blue lines). The gray band indicates the
approximate range of comoving wavenumbers contributing to $\sigma_{8}$. Bottom
panel $-$ Same as above, for the (lensed) CMB TT, EE and lensing power
spectrum. In this case, the gray bands show Planck 1$\sigma$ errors.

Best-fit cosmology and impact on CMB and LSS$-$ A C+WDM scenario like this one
is expected to produce a suppression in the linear matter power spectrum on
intermediate and small scales, with a non-trivial shape Murgia et al. (2017,
2018); Miller et al. (2019); Bohr et al. (2020). The cut-off scale, as in the
case of massive neutrinos, is determined by the WDM free-streaming scale
$k_{\rm fs}$. At wavenumbers $k>k_{\rm fs}$, pressure becomes important and
structure formation is inhibited. To better understand the $\Lambda$DDM
success in resolving the $S_{8}$ tension compared to the case of massive
neutrinos, we now illustrate the impact of both models on the relevant
cosmological observables for our study.

In the top panel of Fig. 2 we compare the residual differences in linear
matter power spectrum $P(k)$ with respect to our baseline $\nu\Lambda$CDM
model (first column of Tab. 1), for both the best-fit $\Lambda$DDM scenario
(fourth column of Tab. 1) and a $\nu\Lambda$CDM model with three degenerate
massive neutrinos of total mass $M_{\nu}=0.23$ eV (we adjust $\omega_{\rm
cdm}=0.116$ whereas all other parameters are fixed to the baseline
$\nu\Lambda$CDM model), which yields $\sigma_{8}=0.76$ and $\Omega_{\rm
m}=0.31$, in agreement with weak lensing data Heymans et al. (2020). These
scenarios feature two key differences: i) a distinct redshift evolution for
the power suppression. In the $\Lambda$DDM scenario, it is less significant at
higher redshifts, since the abundance of the WDM daughter is smaller; ii) a
time-evolving cut-off scale; in the $\Lambda$DDM model, $k_{\rm
fs}=\sqrt{3/2}\mathcal{H}(a)/c_{g}(a)$, while in the $\nu\Lambda$CDM it is
obtained by evaluating $k_{\rm fs}$ at the redshift at which neutrinos become
non-relativistic Lesgourgues et al. (2016). As a consequence, the CMB power
spectra, well constrained by Planck, are vastly different. This is illustrated
in Fig. 2 bottom panel, for both the best-fit $\Lambda$DDM scenario and the
$\nu\Lambda$CDM model which yield the same $S_{8}$ value. The $\nu\Lambda$CDM
predicts different early-integrated Sachs-Wolfe effects, as well as different
amount of lensing, because of a significant power suppression at $z\sim 2-3$,
where the CMB lensing kernel peaks Abazajian et al. (2016). On the other hand,
the differences between $\Lambda$CDM and $\Lambda$DDM until $z\sim 2$ are very
small, explaining why Planck cannot disentangle between both scenarios.
Detecting the DDM through its impact on CMB power spectra will be challenging,
although CMB lensing measurements accurate at the $\sim 1\%$ level could help
(e.g. with CMB-S4 Abazajian et al. (2016)). Fortunately, the differences
between the growth rate $f\sigma_{8}$ in $\Lambda$DDM scenario and
$\Lambda$CDM (shown in Fig. 3) at $0\lesssim z\lesssim 1$, while below the
sensitivity of current experiments measuring, could be measured by upcoming
surveys such as Euclid Amendola et al. (2018), LSST Alonso et al. (2018), and
DESI Aghamousa et al. (2016).

Figure 3: Growth rate of matter fluctuations for our baseline $\nu\Lambda$CDM
model (solid line), compared to the best-fit $\Lambda$DDM model (dashed line)
and to the $\nu\Lambda$CDM scenario yielding the same $\sigma_{8}$ and
$\Omega_{\rm m}$ (dotted line). The observational constraints are taken from
Ref. Aghanim et al. (2018) and references therein.

Some implications of the DCDM$-$ Concrete realizations of $\Lambda$DDM
scenarios as the ones considered in this work may arise (for instance in the
context of the ‘super weakly interacting massive particle’ (superWIMP) class
of exotic particle physics models Covi et al. (1999); Feng et al. (2003,
2003), whose super weak couplings make them evade many observational
constraints Allahverdi et al. (2015); Feng (2010). To connect with parameters
of these models, we can estimate typical orders of magnitude for the decay
rate, $\Gamma\approx(50\,{\rm Gyr})^{-1}(g_{\rm eff}/(10^{-18}\,{\rm
or}\,10^{-16}))^{2}(m/(0.01\,{\rm or}\,0.02\,{\rm GeV}))$
$\times(\varepsilon/0.007)^{2\,{\rm or}\,4}$, where $g_{\rm eff}$ is the
coupling constant for the decay, $m=m_{\rm{dcdm}}$ is the mass of the decaying
particle, and the model-dependent power index is $2$ for particle of spins
$\leq 1$, and $4$ for e.g. gravitino production. The decaying particle must
have properties similar to CDM candidates, which sets a lower mass bound of
$m\gtrsim 5$ keV from structure formation if it is thermally produced in the
early universe Viel et al. (2005); Boyarsky et al. (2009a); Iršič et al.
(2017); Yeche et al. (2017); Schneider (2018), raising to the MeV mass scale
depending on couplings to standard model particles Bœhm et al. (2013)
(irrelevant for thermal production in hidden sectors Feng et al. (2008)). If
produced non-thermally, the Tremaine-Gunn limit Tremaine and Gunn (1979) of
$m\gtrsim 1$ keV Boyarsky et al. (2009b) applies only to fermions.

Interestingly, late decays of CDM to WDM are among the possible proposed cures
to some observational discrepancies with CDM on small (sub-galactic) scales
after structure formation (e.g. Lin et al. (2001); Sánchez-Salcedo (2003);
Cembranos et al. (2005); Kaplinghat (2005); Strigari et al. (2007); Borzumati
et al. (2008); Peter et al. (2010); Peter and Benson (2010), and e.g. de Blok
(2010); Bullock and Boylan-Kolchin (2017) for reviews on small-scale
issues)—see also e.g. Pontzen and Governato (2014); Read and Erkal (2019);
Ostriker et al. (2019); Zavala and Frenk (2019) for different views. A
detailed inspection of the effective parameter space was performed in Ref.
Peter (2010), which was later supplemented by dedicated cosmological
simulations Wang et al. (2014), showing that this scenario mostly affects
galaxy satellite/subhalo properties. In particular, a daughter particle speed
$3\times 10^{-5}\lesssim v\simeq\varepsilon\lesssim 2\times 10^{-3}$ may
reduce the abundance of subhalos and their concentrations for a large range of
lifetimes, up to $\sim 100$ Gyr. Greater speeds are disfavored from the
existence of dwarf galaxies (see also Wang et al. (2013)), unless
$(\Gamma)^{-1}\gtrsim 100$ Gyr, for which $\Lambda$DDM does not depart from
$\Lambda$CDM as far as structure formation is concerned (this also holds for
very low speeds irrespective of $\Gamma$). Further constraining our favored
range for $\Lambda$DDM models along those lines would require to predict the
non-linear matter power spectrum, which goes beyond the scope of this paper.
However, our best-fitting parameters do fall in ranges that could partly
address small-scale CDM challenges.

From different perspectives, remarking that the higher tail of the daughter
particle speed contour reaches $v\simeq\varepsilon\sim 0.05$ may lead to the
intriguing possibility that this scenario be also connected with the excess
events in the electronic recoils recently reported by the Xenon-1t
Collaboration The Xenon-1t Collaboration et al. (2020). If one resorts to
elastic scattering of our daughter particle off electrons, measured recoils in
the keV energy range translate into a lower bound on $m_{\rm wdm}\sim m\gtrsim
1$ MeV Kannike et al. (2020), and require a local interaction rate with
electrons of $n_{\rm wdm}\sigma_{e}\sim 3.5\times 10^{-44}\,{\rm cm^{-1}}$.
The former constraint is consistent with both thermal and non-thermal
production of the parent CDM particle, while it is not clear whether the
latter can be achieved while preventing the production of WDM particles in the
early universe. Another possibility is that the daughter particle is a dark
photon that kinetically mixes with photons, which could therefore transmit its
entire mass energy to electrons from photoelectric processes Pospelov et al.
(2008). This instead completely sets $m_{\rm wdm}\sim m\sim 2$ keV and no
longer relies on $\varepsilon$, which can therefore only be achieved in non-
thermal production scenarios. Interesting realizations of the latter are e.g.
found in models with a broken $U(1)_{\rm B-L}$ gauge symmetry, also proposed,
incidentally, to solve small-scale CDM issues (e.g. Choi et al. (2020a, b,
c)).

Model |  $\nu\Lambda$CDM |  $\Lambda$DDM
---|---|---
Parameter | w/o $S_{8}$ | w/ $S_{8}$ | w/o $S_{8}$ | w/ $S_{8}$
$100\ \omega_{b}$ | $2.245(2.242)\pm 0.013$ | $2.251(2.253)\pm 0.013$ | $2.243(2.242)_{-0.014}^{+0.013}$ | $2.246(2.240)_{-0.014}^{+0.013}$
$\omega_{\rm cdm}\ \text{or}\ \omega_{\rm dcdm}^{\rm ini}$ | $0.1193(0.1194)\pm 0.0009$ | $0.1182(0.1184)_{-0.0008}^{+0.0009}$ | $0.1195(0.1195)\pm 0.0009$ | $0.1190(0.1194)\pm 0.0009$
$H_{0}/[{\rm km/s/Mpc}]$ | $67.55(67.76)_{-0.44}^{+0.46}$ | $67.85(68.08)_{-0.44}^{+0.47}$ | $67.72(67.68)\pm 0.41$ | $67.94(67.70)_{-0.41}^{+0.42}$
$\text{ln}(10^{10}A_{s})$ | $3.052(3.045)_{-0.016}^{+0.014}$ | $3.047(3.043)_{-0.015}^{+0.014}$ | $3.051(3.046)_{-0.016}^{+0.014}$ | $3.049(3.051)_{-0.016}^{+0.014}$
$n_{s}$ | $0.9676(0.9663)\pm 0.0037$ | $0.9697(0.9683)_{-0.0036}^{+0.0037}$ | $0.9674(0.9673)_{-0.0037}^{+0.0036}$ | $0.9683(0.9673)_{-0.0038}^{+0.0035}$
$\tau_{\rm reio}$ | $0.058(0.055)_{-0.008}^{+0.007}$ | $0.0569(0.0549)_{-0.008}^{+0.007}$ | $0.0576(0.0552)_{-0.0079}^{+0.007}$ | $0.0573(0.0582)_{-0.0081}^{+0.0072}$
$M_{\nu}/\text{eV}$ | $<0.1395$ | $<0.1611$ | $-$ | $-$
$\text{log}_{10}(\varepsilon)$ | $-$ | $-$ | $<-0.79(-2.90)$ | $-2.33(-2.16)_{-0.73}^{+0.86}$
$\text{log}_{10}(\Gamma/[{\rm Gyr}^{-1}])$ | $-$ | $-$ | unconstrained (-3.86) | $-1.88(-1.75)_{-1.4}^{+0.95}$
$\Omega_{\rm m}$ | $0.3127(0.3104)_{-0.0061}^{+0.0057}$ | $0.3083(0.3061)_{-0.006}^{+0.0056}$ | $0.3101(0.3112)_{-0.0056}^{+0.0054}$ | $0.3068(0.3098)_{-0.0055}^{+0.0053}$
$S_{8}$ | $0.824(0.824)\pm 0.011$ | $0.81(0.816)\pm 0.01$ | $0.821(0.826)_{-0.011}^{+0.016}$ | $0.795(0.773)_{-0.015}^{+0.025}$
$\chi^{2}_{\rm min}$ | 2053.4 | 2060.5 | 2053.4 | 2055.1
Table 1: The mean (best-fit) $\pm 1\sigma$ errors of the cosmological
parameters from our combined analysis, with and without a split-normal
likelihood on $S_{8}$ from Ref. Heymans et al. (2020). For each model and
data-set, we also report the best-fit $\chi^{2}$.

Conclusions$-$ In this Letter we have studied the cosmological phenomenology
of a scenario in which CDM decays into one WDM and one DR species. Through a
comprehensive MCMC analysis including up-to-date data, we demonstrated that
the $\Lambda$DDM scenario can resolve the $S_{8}$ tension. The inferred
$S_{8}$ value is in excellent agreement with its direct measurement from WL,
due to the suppression in the gravitational clustering induced by WDM free-
streaming in a similar fashion to massive neutrino or standard WDM
cosmologies. However contrarily to the latter, the specific time dependence of
the power suppression imprinted by the decay allows to accommodate
simultaneously all the data considered in this work. On the other hand, we
confirm that this model cannot resolve the Hubble tension. Finally, we briefly
outlined implications and possible future developments of our study, in view
of DM model building, the $\Lambda$CDM small-scale crisis and the recent
Xenon-1T excess. Future experiments measuring the growth rate of fluctuations
at $0\lesssim z\lesssim 1$ or the CMB lensing power spectrum to $\sim 1\%$
accuracy will further test this scenario.

Acknowledgements$-$ The authors warmly thank Rodrigo Calderón, Zackaria
Chacko, Abhish Dev, Peizhi Du, Yuhsin Tsai, Julien Lesgourgues and Pasquale
Serpico for many useful comments and discussions. This work has been partly
supported by the CNRS-IN2P3 grant Dark21, the ANR project ANR-18-CE31-0006,
the national CNRS-INSU programs PNHE and PNCG, and the European Union’s
Horizon 2020 research and innovation program under the Marie Skłodowska-Curie
grant agreements N∘ 690575 and N∘ 674896\. The authors acknowledge the use of
computational resources from the Ulysses SISSA/ICTP super-computer in Trieste,
the CNRS/IN2P3 Computing Centre (CC-IN2P3) in Lyon, the IN2P3/CNRS and the
Dark Energy computing Center funded by the OCEVU Labex (ANR-11-LABX-0060) and
the Excellence Initiative of Aix-Marseille University - A*MIDEX, part of the
French “Investissements d’Avenir” programme.

## References

  * Scott (2018) D. Scott (2018), eprint 1804.01318.
  * Akrami et al. (2018) Y. Akrami et al. (Planck) (2018), eprint 1807.06205.
  * Aghanim et al. (2018) N. Aghanim et al. (Planck) (2018), eprint 1807.06209.
  * Heymans et al. (2013) C. Heymans et al., Mon. Not. Roy. Astron. Soc. 432, 2433 (2013), eprint 1303.1808.
  * Hikage et al. (2019) C. Hikage, M. Oguri, T. Hamana, S. More, R. Mandelbaum, M. Takada, F. Köhlinger, H. Miyatake, A. J. Nishizawa, H. Aihara, et al., Publications of the Astronomical Society of Japan 71 (2019), ISSN 2053-051X, URL http://dx.doi.org/10.1093/pasj/psz010.
  * Abbott et al. (2018) T. Abbott et al. (DES), Phys. Rev. D 98, 043526 (2018), eprint 1708.01530.
  * Hildebrandt et al. (2020) H. Hildebrandt et al., Astron. Astrophys. 633, A69 (2020), eprint 1812.06076.
  * Heymans et al. (2020) C. Heymans et al. (2020), eprint 2007.15632.
  * Lin et al. (2001) W. B. Lin, D. H. Huang, X. Zhang, and R. Brandenberger, PRL 86, 954 (2001), eprint astro-ph/0009003, URL https://ui.adsabs.harvard.edu/abs/2001PhRvL..86..954L.
  * Sánchez-Salcedo (2003) F. J. Sánchez-Salcedo, Apjl 591, L107 (2003), eprint astro-ph/0305496, URL https://ui.adsabs.harvard.edu/abs/2003ApJ...591L.107S.
  * Cembranos et al. (2005) J. A. R. Cembranos, J. L. Feng, A. Rajaraman, and F. Takayama, PRL 95, 181301 (2005), eprint hep-ph/0507150, URL https://ui.adsabs.harvard.edu/abs/2005PhRvL..95r1301C.
  * Kaplinghat (2005) M. Kaplinghat, Phys. Rev. D 72, 063510 (2005), eprint astro-ph/0507300, URL https://ui.adsabs.harvard.edu/abs/2005PhRvD..72f3510K.
  * Strigari et al. (2007) L. E. Strigari, M. Kaplinghat, and J. S. Bullock, Phys. Rev. D 75, 061303 (2007), eprint astro-ph/0606281, URL https://ui.adsabs.harvard.edu/abs/2007PhRvD..75f1303S.
  * Borzumati et al. (2008) F. Borzumati, T. Bringmann, and P. Ullio, Phys. Rev. D 77, 063514 (2008), eprint hep-ph/0701007, URL https://ui.adsabs.harvard.edu/abs/2008PhRvD..77f3514B.
  * Peter et al. (2010) A. H. G. Peter, C. E. Moody, and M. Kamionkowski, Phys. Rev. D 81, 103501 (2010), eprint 1003.0419, URL https://ui.adsabs.harvard.edu/abs/2010PhRvD..81j3501P.
  * Peter and Benson (2010) A. H. G. Peter and A. J. Benson, Phys. Rev. D 82, 123521 (2010), eprint 1009.1912, URL https://ui.adsabs.harvard.edu/abs/2010PhRvD..82l3521P.
  * Wang et al. (2014) M.-Y. Wang, A. H. G. Peter, L. E. Strigari, A. R. Zentner, B. Arant, S. Garrison-Kimmel, and M. Rocha, MNRAS 445, 614 (2014), eprint 1406.0527, URL https://ui.adsabs.harvard.edu/abs/2014MNRAS.445..614W.
  * The Xenon-1t Collaboration et al. (2020) The Xenon-1t Collaboration, E. Aprile, J. Aalbers, and F. a. Agostini, arXiv e-prints arXiv:2006.09721 (2020), eprint 2006.09721, URL https://ui.adsabs.harvard.edu/abs/2020arXiv200609721A.
  * Kannike et al. (2020) K. Kannike, M. Raidal, H. Veermäe, A. Strumia, and D. Teresi, arXiv e-prints arXiv:2006.10735 (2020), eprint 2006.10735, URL https://ui.adsabs.harvard.edu/abs/2020arXiv200610735K.
  * Lesgourgues et al. (2016) J. Lesgourgues, G. Marques-Tavares, and M. Schmaltz, JCAP 02, 037 (2016), eprint 1507.04351.
  * Chacko et al. (2016) Z. Chacko, Y. Cui, S. Hong, T. Okui, and Y. Tsai, JHEP 12, 108 (2016), eprint 1609.03569.
  * Murgia et al. (2016) R. Murgia, S. Gariazzo, and N. Fornengo, JCAP 04, 014 (2016).
  * Kobayashi et al. (2017) T. Kobayashi, R. Murgia, A. De Simone, V. Irˇsič, and M. Viel, Phys. Rev. D 96, 123514 (2017), eprint 1708.00015.
  * Di Valentino et al. (2018) E. Di Valentino, E. V. Linder, and A. Melchiorri, Phys. Rev. D97, 043528 (2018), eprint 1710.02153.
  * Buen-Abad et al. (2018) M. A. Buen-Abad, M. Schmaltz, J. Lesgourgues, and T. Brinckmann, JCAP 01, 008 (2018), eprint 1708.09406.
  * Raveri et al. (2017) M. Raveri, W. Hu, T. Hoffman, and L.-T. Wang, Phys. Rev. D 96, 103501 (2017), eprint 1709.04877.
  * Poulin et al. (2019) V. Poulin, T. L. Smith, T. Karwal, and M. Kamionkowski, Phys. Rev. Lett. 122, 221301 (2019), eprint 1811.04083.
  * Lin et al. (2019) M.-X. Lin, M. Raveri, and W. Hu, Phys. Rev. D99, 043514 (2019), eprint 1810.02333.
  * D’Eramo et al. (2018) F. D’Eramo, R. Z. Ferreira, A. Notari, and J. L. Bernal, JCAP 1811, 014 (2018), eprint 1808.07430.
  * Dutta et al. (2020) K. Dutta, Ruchika, A. Roy, A. A. Sen, and M. M. Sheikh-Jabbari, Gen. Rel. Grav. 52, 15 (2020), eprint 1808.06623.
  * Archidiacono et al. (2019) M. Archidiacono, D. C. Hooper, R. Murgia, S. Bohr, J. Lesgourgues, and M. Viel, JCAP 10, 055 (2019), eprint 1907.01496.
  * Vattis et al. (2019) K. Vattis, S. M. Koushiappas, and A. Loeb, Phys. Rev. D99, 121302 (2019), eprint 1903.06220.
  * Bohr et al. (2020) S. Bohr, J. Zavala, F.-Y. Cyr-Racine, M. Vogelsberger, T. Bringmann, and C. Pfrommer (2020), eprint 2006.01842.
  * Heimersheim et al. (2020) S. Heimersheim, N. Schöneberg, D. C. Hooper, and J. Lesgourgues (2020), eprint 2008.08486.
  * Aoyama et al. (2014) S. Aoyama, T. Sekiguchi, K. Ichiki, and N. Sugiyama, JCAP 1407, 021 (2014), eprint 1402.2972.
  * Enqvist et al. (2015) K. Enqvist, S. Nadathur, T. Sekiguchi, and T. Takahashi, JCAP 09, 067 (2015), eprint 1505.05511.
  * Chudaykin et al. (2016) A. Chudaykin, D. Gorbunov, and I. Tkachev, EPJ Web Conf. 125, 03004 (2016).
  * Poulin et al. (2016) V. Poulin, P. D. Serpico, and J. Lesgourgues, JCAP 1608, 036 (2016), eprint 1606.02073.
  * Bringmann et al. (2018) T. Bringmann, F. Kahlhoefer, K. Schmidt-Hoberg, and P. Walia, Phys. Rev. D98, 023543 (2018), eprint 1803.03644.
  * Riess (2019) A. G. Riess, Nature Reviews Physics 2, 10 (2019), eprint 2001.03624.
  * Verde et al. (2019) L. Verde, T. Treu, and A. G. Riess, in _Nature Astronomy 2019_ (2019), eprint 1907.10625.
  * Aylor et al. (2019) K. Aylor, M. Joy, L. Knox, M. Millea, S. Raghunathan, and W. K. Wu, Astrophys. J. 874, 4 (2019), eprint 1811.00537.
  * Wong et al. (2019) K. C. Wong et al. (2019), eprint 1907.04869.
  * Freedman et al. (2019) W. L. Freedman et al. (2019), eprint 1907.05922.
  * Haridasu and Viel (2020) B. S. Haridasu and M. Viel (2020), eprint 2004.07709.
  * Clark et al. (2020) S. J. Clark, K. Vattis, and S. M. Koushiappas (2020), eprint 2006.03678.
  * Blackadder and Koushiappas (2014) G. Blackadder and S. M. Koushiappas, Phys. Rev. D90, 103527 (2014), eprint 1410.0683.
  * Wang and Zentner (2012) M.-Y. Wang and A. R. Zentner, Phys. Rev. D 85, 043514 (2012), eprint 1201.2426.
  * F. Abellán et al. (in prep.) G. F. Abellán, V. Poulin, and R. Murgia (in prep.).
  * Lesgourgues and Tram (2011) J. Lesgourgues and T. Tram, JCAP 09, 032 (2011), eprint 1104.2935.
  * Blas et al. (2011) D. Blas, J. Lesgourgues, and T. Tram, JCAP 1107, 034 (2011), eprint 1104.2933.
  * Ma and Bertschinger (1995) C.-P. Ma and E. Bertschinger, Astrophys. J. 455, 7 (1995), eprint astro-ph/9506072.
  * Audren et al. (2013) B. Audren, J. Lesgourgues, K. Benabed, and S. Prunet, JCAP 1302, 001 (2013), eprint 1210.7183.
  * Brinckmann and Lesgourgues (2019) T. Brinckmann and J. Lesgourgues, Phys. Dark Univ. 24, 100260 (2019), eprint 1804.07261.
  * Beutler et al. (2011) F. Beutler, C. Blake, M. Colless, D. Jones, L. Staveley-Smith, L. Campbell, Q. Parker, W. Saunders, and F. Watson, Mon. Not. Roy. Astron. Soc. 416, 3017 (2011), eprint 1106.3366.
  * Ross et al. (2015) A. J. Ross, L. Samushia, C. Howlett, W. J. Percival, A. Burden, and M. Manera, Mon. Not. Roy. Astron. Soc. 449, 835 (2015), eprint 1409.3242.
  * Alam et al. (2017) S. Alam et al. (BOSS), Mon. Not. Roy. Astron. Soc. 470, 2617 (2017), eprint 1607.03155.
  * de Sainte Agathe et al. (2019) V. de Sainte Agathe et al., Astron. Astrophys. 629, A85 (2019), eprint 1904.03400.
  * Blomqvist et al. (2019) M. Blomqvist et al., Astron. Astrophys. 629, A86 (2019), eprint 1904.03430.
  * Scolnic et al. (2018) D. M. Scolnic et al., Astrophys. J. 859, 101 (2018), eprint 1710.00845.
  * Joudaki et al. (2020) S. Joudaki et al., Astron. Astrophys. 638, L1 (2020), eprint 1906.09262.
  * Baur et al. (2017) J. Baur, N. Palanque-Delabrouille, C. Yeche, A. Boyarsky, O. Ruchayskiy, . Armengaud, and J. Lesgourgues, JCAP 12, 013 (2017), eprint 1706.03118.
  * Diamanti et al. (2017) R. Diamanti, S. Ando, S. Gariazzo, O. Mena, and C. Weniger, JCAP 06, 008 (2017), eprint 1701.03128.
  * Gariazzo et al. (2017) S. Gariazzo, M. Escudero, R. Diamanti, and O. Mena, Phys. Rev. D 96, 043501 (2017), eprint 1704.02991.
  * Murgia et al. (2018) R. Murgia, V. Irsic, and M. Viel, Phys. Rev. D 98, 083540 (2018), eprint 1806.08371.
  * Murgia et al. (2017) R. Murgia, A. Merle, M. Viel, M. Totzauer, and A. Schneider, JCAP 11, 046 (2017), eprint 1704.07838.
  * Miller et al. (2019) C. Miller, A. L. Erickcek, and R. Murgia, Phys. Rev. D 100, 123520 (2019), eprint 1908.10369.
  * Abazajian et al. (2016) K. N. Abazajian et al. (CMB-S4) (2016), eprint 1610.02743.
  * Amendola et al. (2018) L. Amendola et al., Living Rev. Rel. 21, 2 (2018), eprint 1606.00180.
  * Alonso et al. (2018) D. Alonso et al. (LSST Dark Energy Science) (2018), eprint 1809.01669.
  * Aghamousa et al. (2016) A. Aghamousa et al. (DESI) (2016), eprint 1611.00036.
  * Covi et al. (1999) L. Covi, J. E. Kim, and L. Roszkowski, PRL 82, 4180 (1999), eprint hep-ph/9905212, URL https://ui.adsabs.harvard.edu/abs/1999PhRvL..82.4180C.
  * Feng et al. (2003) J. L. Feng, A. Rajaraman, and F. Takayama, PRL 91, 011302 (2003), eprint hep-ph/0302215, URL https://ui.adsabs.harvard.edu/abs/2003PhRvL..91a1302F.
  * Feng et al. (2003) J. L. Feng, A. Rajaraman, and F. Takayama, Phys. Rev. D 68, 063504 (2003), eprint hep-ph/0306024, URL https://ui.adsabs.harvard.edu/abs/2003PhRvD..68f3504F.
  * Allahverdi et al. (2015) R. Allahverdi, B. Dutta, F. S. Queiroz, L. E. Strigari, and M.-Y. Wang, Phys. Rev. D 91, 055033 (2015), eprint 1412.4391, URL https://ui.adsabs.harvard.edu/abs/2015PhRvD..91e5033A.
  * Feng (2010) J. L. Feng, ARAA 48, 495 (2010), eprint 1003.0904, URL https://ui.adsabs.harvard.edu/abs/2010ARA%26A..48..495F.
  * Viel et al. (2005) M. Viel, J. Lesgourgues, M. G. Haehnelt, S. Matarrese, and A. Riotto, Phys. Rev. D 71, 063534 (2005), eprint astro-ph/0501562, URL https://ui.adsabs.harvard.edu/abs/2005PhRvD..71f3534V.
  * Boyarsky et al. (2009a) A. Boyarsky, J. Lesgourgues, O. Ruchayskiy, and M. Viel, JCAP 5, 012 (2009a), eprint 0812.0010, URL https://ui.adsabs.harvard.edu/abs/2009JCAP...05..012B.
  * Iršič et al. (2017) V. Iršič, M. Viel, M. G. Haehnelt, J. S. Bolton, S. Cristiani, G. D. Becker, V. D’Odorico, G. Cupani, T.-S. Kim, T. A. M. Berg, et al., Phys. Rev. D 96, 023522 (2017), eprint 1702.01764, URL https://ui.adsabs.harvard.edu/abs/2017PhRvD..96b3522I.
  * Yeche et al. (2017) C. Yeche, N. Palanque-Delabrouille, J. . Baur, and H. du Mas des BourBoux, ArXiv e-prints (2017), eprint 1702.03314, URL https://ui.adsabs.harvard.edu/abs/2017arXiv170203314Y.
  * Schneider (2018) A. Schneider, Phys. Rev. D 98, 063021 (2018), eprint 1805.00021, URL https://ui.adsabs.harvard.edu/abs/2018PhRvD..98f3021S.
  * Bœhm et al. (2013) C. Bœhm, M. J. Dolan, and C. McCabe, JCAP 8, 041 (2013), eprint 1303.6270, URL https://ui.adsabs.harvard.edu/abs/2013JCAP...08..041B.
  * Feng et al. (2008) J. L. Feng, H. Tu, and H.-B. Yu, JCAP 10, 043 (2008), eprint 0808.2318, URL https://ui.adsabs.harvard.edu/abs/2008JCAP...10..043F.
  * Tremaine and Gunn (1979) S. Tremaine and J. E. Gunn, PRL 42, 407 (1979), URL https://ui.adsabs.harvard.edu/abs/1979PhRvL..42..407T.
  * Boyarsky et al. (2009b) A. Boyarsky, O. Ruchayskiy, and D. Iakubovskyi, JCAP 3, 005 (2009b), eprint 0808.3902, URL https://ui.adsabs.harvard.edu/abs/2009JCAP...03..005B.
  * de Blok (2010) W. J. G. de Blok, Advances in Astronomy 2010, 789293 (2010), eprint 0910.3538, URL https://ui.adsabs.harvard.edu/abs/2010AdAst2010E...5D.
  * Bullock and Boylan-Kolchin (2017) J. S. Bullock and M. Boylan-Kolchin, ARAA 55 (2017), eprint 1707.04256, URL https://ui.adsabs.harvard.edu/abs/2017arXiv170704256B.
  * Pontzen and Governato (2014) A. Pontzen and F. Governato, Nature (London) 506, 171 (2014), eprint 1402.1764, URL https://ui.adsabs.harvard.edu/abs/2014Natur.506..171P.
  * Read and Erkal (2019) J. I. Read and D. Erkal, MNRAS 487, 5799 (2019), eprint 1807.07093, URL https://ui.adsabs.harvard.edu/abs/2019MNRAS.487.5799R.
  * Ostriker et al. (2019) J. P. Ostriker, E. Choi, A. Chow, and K. Guha, Astrophys. J. 885, 97 (2019), eprint 1904.10471, URL https://ui.adsabs.harvard.edu/abs/2019ApJ...885...97O.
  * Zavala and Frenk (2019) J. Zavala and C. S. Frenk, Galaxies 7, 81 (2019), eprint 1907.11775, URL https://ui.adsabs.harvard.edu/abs/2019Galax...7...81Z.
  * Peter (2010) A. H. G. Peter, Phys. Rev. D 81, 083511 (2010), eprint 1001.3870, URL https://ui.adsabs.harvard.edu/abs/2010PhRvD..81h3511P.
  * Wang et al. (2013) M.-Y. Wang, R. A. C. Croft, A. H. G. Peter, A. R. Zentner, and C. W. Purcell, Phys. Rev. D 88, 123515 (2013), eprint 1309.7354, URL https://ui.adsabs.harvard.edu/abs/2013PhRvD..88l3515W.
  * Pospelov et al. (2008) M. Pospelov, A. Ritz, and M. Voloshin, Phys. Rev. D 78, 115012 (2008), eprint 0807.3279, URL https://ui.adsabs.harvard.edu/abs/2008PhRvD..78k5012P.
  * Choi et al. (2020a) G. Choi, M. Suzuki, and T. T. Yanagida, arXiv e-prints arXiv:2006.12348 (2020a), eprint 2006.12348, URL https://ui.adsabs.harvard.edu/abs/2020arXiv200612348C.
  * Choi et al. (2020b) G. Choi, M. Suzuki, and T. T. Yanagida, Phys. Rev. D 101, 075031 (2020b), eprint 2002.00036, URL https://ui.adsabs.harvard.edu/abs/2020PhRvD.101g5031C.
  * Choi et al. (2020c) G. Choi, M. Suzuki, and T. T. Yanagida, arXiv e-prints arXiv:2004.07863 (2020c), eprint 2004.07863, URL https://ui.adsabs.harvard.edu/abs/2020arXiv200407863C.