#  ProtGNN: Towards Self-Explaining Graph Neural Networks

Zaixi Zhang1, Qi Liu1∗, Hao Wang1, Chengqiang Lu1, Cheekong Lee2 Qi Liu and
Chee-Kong Lee are corresponding authors.

# ProtGNN: Towards Self-Explaining Graph Neural Networks

Zaixi Zhang1, Qi Liu1∗, Hao Wang1, Chengqiang Lu1, Cheekong Lee2 Qi Liu and
Chee-Kong Lee are corresponding authors.

###### Abstract

Despite the recent progress in Graph Neural Networks (GNNs), it remains
challenging to explain the predictions made by GNNs. Existing explanation
methods mainly focus on _post-hoc_ explanations where another explanatory
model is employed to provide explanations for a trained GNN. The fact that
_post-hoc_ methods fail to reveal the original reasoning process of GNNs
raises the need of building GNNs with _built-in_ interpretability. In this
work, we propose Prototype Graph Neural Network (ProtGNN), which combines
prototype learning with GNNs and provides a new perspective on the
explanations of GNNs. In ProtGNN, the explanations are naturally derived from
the case-based reasoning process and are actually used during classification.
The prediction of ProtGNN is obtained by comparing the inputs to a few learned
prototypes in the latent space. Furthermore, for better interpretability and
higher efficiency, a novel conditional subgraph sampling module is
incorporated to indicate which part of the input graph is most similar to each
prototype in ProtGNN+. Finally, we evaluate our method on a wide range of
datasets and perform concrete case studies. Extensive results show that
ProtGNN and ProtGNN+ can provide inherent interpretability while achieving
accuracy on par with the non-interpretable counterparts.

## Introduction

Graph Neural Networks (GNNs) have become increasingly popular since many real-
world relational data can be represented as graphs, such as social networks
(Bian et al. 2020), molecules (Gilmer et al. 2017) and financial data (Yang et
al. 2020). Following a message passing paradigm to learn node representations,
GNNs have achieved state-of-the-art performance in node classification, graph
classification, and link prediction (Kipf and Welling 2017; Veličković et al.
2017; Xu et al. 2019). Despite the remarkable effectiveness of GNNs,
explaining predictions made by GNNs remains a challenging open problem.
Without understanding the rationales behind the predictions, these black-box
models cannot be fully trusted and widely applied in critical areas such as
medical diagnosis. In addition, model explanations can facilitate model
debugging and error analysis. These indicate the necessity of investigating
the explainability of GNNs.

Recently, extensive efforts have been made to study explanation techniques for
GNNs (Yuan et al. 2020b). These methods can explain the predictions of node or
graph classifications of trained GNNs with different strategies. For example,
GNNExplainer (Ying et al. 2019) and PGExplainer (Luo et al. 2020) are proposed
to select a compact subgraph structure that maximizes the mutual information
with the GNN’s predictions as the explanation. PGM-Explainer (Vu and Thai
2020) firstly obtains a local dataset by random node feature perturbation.
Then it employs an interpretable Bayesian network to fit the local dataset and
to explain the predictions of the original GNN model. In addition, XGNN (Yuan
et al. 2020a) generates graph patterns to maximize the predicted probability
for a certain class and provides model-level explanation. Despite the
tremendous developments in the interpretation of GNNs, most existing
approaches can be classified as _post-hoc_ explanations where another
explanatory model is used to provide explanations for a trained GNN. _Post-
hoc_ explanations can be inaccurate or incomplete in revealing the actual
reasoning process of the original model (Rudin 2018). Therefore, it is more
desirable to build models with inherent interpretability where the
explanations are generated by the model themselves.

We leverage the concept of prototype learning to construct GNNs with _built-
in_ interpretability (i.e. self-explaining GNNs). In contrast to _post-hoc_
explanation methods, the explanations generated by self-explaining GNNs are
actually used during classification and are not generated _post-hoc_.
Prototype learning is a form of case-based reasoning (Kolodner 1992; Schmidt
et al. 2001), which makes the predictions for new instances by comparing them
with several learned exemplar cases (i.e. prototypes). It is a natural
practice in solving problems with graph-structured data. For example, chemists
identify potential drug candidates based on known functional groups (i.e. key
subgraphs) in molecular graphs (He et al. 2010; Zhang et al. 2021). Prototype
learning provides better interpretability by imitating such a human problem-
solving process. Recently the concept of the prototype has been incorporated
in convolutional neural networks to build interpretable image classifiers
(Chen et al. 2018; Rymarczyk et al. 2021). However, so far prototype learning
is not yet explored for explaining GNNs.

Building self-explaining GNNs based on prototype learning brings unique
challenges. First, the discreteness of the edges makes the projection and
visualization of the graph prototypes difficult. Second, the combinatorial
nature of graph structure makes it hard to build self-explaining models with
both efficiency and high accuracy for graph modeling.

In this paper, we tackle the aforementioned challenges and propose Prototype
Graph Neural Network (ProtGNN), which provides a new perspective on the
explanations of GNNs. Specifically, various popular GNN architectures can be
employed as the graph encoder in ProtGNN. Prediction on a new input graph is
performed based on its similarity to the prototypes in the prototype layer.
Furthermore, we propose to employ the Monte Carlo tree search algorithm
(Silver et al. 2017) to efficiently explore different subgraphs for prototype
projection and visualization. In addition, in ProtGNN+, we design a
conditional subgraph sampling module to identify which part of the input graph
is most similar to each prototype for better interpretability and efficiency.
Finally, extensive experiments on several real-world datasets show that
ProtGNN/ProtGNN+ provides _built-in_ interpretability while achieving
comparable performance with the non-interpretable counterparts.

## Related Work

### Graph Neural Networks

Graph neural networks have demonstrated their effectiveness on various graph
tasks. Let $G=(V,E)$ denotes a graph with node attributes $X_{v}$ for $v\in V$
and a set of edges $E$. GNNs leverage the graph connectivity as well as node
and edge features to learn a representation vector (i.e., embedding) $h_{v}$
for each node $v\in V$ or a vector $h_{G}$ for the entire graph $G$.
Generally, GNNs follows a message passing paradigm, in which the
representation of node $v$ is iteratively updated by aggregating the
representations of $v$’s neighboring nodes $\mathcal{N}(v)$. Here we use Graph
Convolutional Network (GCN) (Kipf and Welling 2017) as an example to
illustrate such message passing procedures:

$h_{v}^{k+1}=\sigma\bigg{(}\mathop{\sum}_{u\in\mathcal{N}(v)}\big{(}W^{k}h_{u}^{k}\tilde{A}_{uv}\big{)}\bigg{)},$
(1)

where $h_{u}^{k}$ is the representation vector of node $u$ at the $k$-th layer
and $\tilde{A}=\hat{D}^{-\frac{1}{2}}\hat{A}\hat{D}^{-\frac{1}{2}}$ is the
normalized adjacency matrix. $\hat{A}=A+I$ is the adjacency matrix of the
graph $G$ with self connections added and $\hat{D}$ is a diagonal matrix with
$\hat{D}_{ii}=\sum_{j}\hat{A}_{ij}$. $\sigma(\cdot)$ in Eq. (1) is the ReLU
function and $W^{k}$ is the trainable weight matrix of the $k$-th layer.

### Explainability in Graph Neural Networks

As the application of GNNs grows, understanding why GNNs make such predictions
becomes increasingly critical. Recently, the study of the explainability in
GNNs is experiencing rapid developments. As Suggested by a recent survey (Yuan
et al. 2020b), existing methods for explaining GNNs can be categorized into
several classes: gradients/features-based methods (Baldassarre and Azizpour
2019; Pope et al. 2019), perturbation-based methods (Ying et al. 2019; Luo et
al. 2020; Yuan et al. 2021; Schlichtkrull, De Cao, and Titov 2020),
decomposition methods (Schwarzenberg et al. 2019; Schnake et al. 2020), and
surrogate methods (Vu and Thai 2020; Huang et al. 2020).

Specifically, the gradients/features-based methods employ the gradients or the
feature values to indicate the importance of different input features. These
methods simply adapt existing explanation techniques in the image domain to
the graph domain without incorporating the properties of graph data.
Perturbation-based methods monitor the changes in the predictions by
perturbing different input features and identifies the most influential
features. Decomposition methods explain GNNs by decomposing the original model
predictions into several terms and associating these terms with graph nodes or
edges. Given an input example, surrogate methods firstly sample a dataset from
the neighborhood of the given example and then fit a simple and interpretable
model, e.g., a decision tree to the sampled dataset. The surrogate models are
usually easier to interpret, shedding light into the inner-workings of more
complex models.

However, all the above methods are _post-hoc_ explanation methods. Compared
with _post-hoc_ explanation methods, _built-in_ interpretability (Chen et al.
2018; Ming et al. 2019) is more desirable since _post-hoc_ explanations
usually do not fit the original model precisely (Rudin 2018). Therefore, it is
necessary to build models with inherent interpretability and high accuracy.

## The Proposed ProtGNN

In this section, We introduce the architecture of ProtGNN/ProtGNN+, formulate
the learning objective and describe the training procedures.

### ProtGNN Architecture

Figure 1: The architecture of our proposed ProtGNN/ProtGNN+. The model mainly
consists of three parts: GNN encoder $f$, prototype layer $g_{P}$, and the
fully connected layer $c$ appended by softmax to output probabilities in
multi-class classification tasks. ProtGNN calculates the similarity score
($sim(p_{k},\cdot)$ in the illustration) between the graph embedding and the
learned prototypes in the prototype layer. For further interpretability, the
conditional subgraph sampling module (in the dashed bounding box) is
incorporated in ProtGNN+ to output subgraphs most similar to each learned
prototype.

We let $\\{x_{i},y_{i}\\}_{i=1}^{n}$ be a labeled training dataset, where
$x_{i}$ is the input attributed graph and $y_{i}\in\\{1,...,C\\}$ is the label
of the graph. We aim to learn representative prototypical graph patterns that
can be used for classification references and analogical explanations. For a
new input graph, its similarities with each prototype are measured in the
latent space. Then, the prediction of the new instance can be derived and
explained by its similar prototype graph patterns.

In Figure 1, we show the overview of the architecture of our proposed ProtGNN.
The network consists of three key components: a graph encoder $f$, a prototype
layer $g_{P}$, and a fully connected layer $c$ appended by softmax to output
the probabilities in multi-class classification tasks.

For a given input graph $x_{i}$, the graph encoder $f$ maps the whole graph
into a single graph embedding $h$ with a fixed length. The encoder could be
any backbone GNN e.g., GCN, GAT or GIN. The graph embedding $h$ could be
obtained by taking a sum or max pooling of the last GNN layer.

In the prototype layer, we allocate a pre-determined number of prototypes $m$
for each class. In the final trained ProtGNN, each class can be represented by
a set of learned prototypes. The prototypes should capture the most relevant
graph patterns for identifying graphs of each class. For each input graph
$x_{i}$ and its embedding $h$, the prototype layer computes the similarity
scores:

$sim(p_{k},h)=log(\frac{\|p_{k}-h\|_{2}^{2}+1}{\|p_{k}-h\|_{2}^{2}+\epsilon})$
(2)

where $p_{k}$ is the $k$-th prototype with the same dimension as the graph
embedding $h$. The similarity function is designed to be monotonically
decreasing to $\|p_{k}-h\|_{2}$ and always greater than zero. In experiments,
$\epsilon$ is set to a small value e.g., 1e-4. Finally, with the similarity
scores, the fully connected layer with softmax computes the output
probabilities for each class.

### Learning Objective

Our goal is to learn a ProtGNN with both accuracy and inherent
interpretability. For accuracy, we minimize the cross-entropy loss on the
training dataset: $\frac{1}{n}\sum_{i=1}^{n}{\rm CrsEnt}(c\circ g_{p}\circ
f(x_{i}),y_{i})$. For better interpretability, we consider several constraints
in constructing prototypes for the explanation. Firstly, the cluster cost
(Clst) encourages that each graph embedding should at least be close to one
prototype of its own class. Secondly, the separation cost (Sep) encourages
that each graph embedding should stay far away from prototypes not of its
class. Finally, we found in experiments that some learned prototypes are very
close to each other in the latent space. We encourage the diversity of the
learned prototypes by adding the diversity loss (Div) which penalizes
prototypes too close to each other.

To sum up, the objective function we aim to minimize is

$\frac{1}{n}\sum_{i=1}^{n}{\rm CrsEnt}(c\circ g_{p}\circ
f(x_{i}),y_{i})+\lambda_{1}{\rm Clst}+\lambda_{2}{\rm Sep}+\lambda_{3}{\rm
Div},$ (3) ${\rm Clst}=\frac{1}{n}\sum_{i=1}^{n}\mathop{\rm
min}\limits_{j:p_{j}\in P_{y_{i}}}\|f(x_{i})-p_{j}\|_{2}^{2}$ (4) ${\rm
Sep}=-\frac{1}{n}\sum_{i=1}^{n}\mathop{\rm min}\limits_{j:p_{j}\notin
P_{y_{i}}}\|f(x_{i})-p_{j}\|_{2}^{2}$ (5) ${\rm
Div}=\sum_{k=1}^{C}\sum_{\mathop{i\neq j}\limits_{p_{i},p_{j}\in P_{k}}}{\rm
max}(0,cos(p_{i},p_{j})-s_{max})$ (6)

where $\lambda_{1}$, $\lambda_{2}$, and $\lambda_{3}$ are hyper-parameters
controlling the weights of the losses. $P_{y_{i}}$ is the set of prototypes
belonging to class $y_{i}$. $s_{max}$ is the threshold of the cosine
similarity measured by $cos(\cdot,\cdot)$ in the diversity loss.

### Prototype Projection

The learned prototypes are embedding vectors that are not directly
interpretable. For better interpretation and visualization, we design a
projection procedure performed in the training stage. Specifically, we project
each prototype $p_{j}$ ($p_{j}\in P_{k}$) onto the nearest latent training
subgraph from the same class as that of $p_{j}$ (see Eq. (7)). In this way, we
can conceptually equate each prototype with a subgraph, which is more
intuitive and human-intelligible. To reduce the computational cost, the
projection step is only performed every few training epochs:

$\displaystyle p_{j}$ $\displaystyle\leftarrow{\rm arg}~{}\mathop{\rm
min}\limits_{\widetilde{h}\in\mathcal{H}_{j}}\|\widetilde{h}-p_{j}\|_{2},$ (7)
$\displaystyle\mathcal{H}_{j}=\\{\widetilde{h}:f(\widetilde{x}),\widetilde{x}$
$\displaystyle\in{\rm Subgraph}(x_{i})~{}\forall i~{}s.t.~{}y_{i}=k\\}.$

Unlike grid-like data such as images, the combinatorial characteristic of
graph makes it unrealistic to find the nearest subgraph by enumeration (Chen
et al. 2018). In graph prototype projection, we employ the Monte Carlo tree
search algorithm (MCTS) (Silver et al. 2017) as the search algorithm to guide
our subgraph explorations (see Figure 2). We build a search tree in which the
root is associated with the input graph and each of other nodes corresponds to
an explored subgraph. Formally, we define each node in the search tree
$\mathcal{T}$ as $\mathcal{N}_{i}$ and $\mathcal{N}_{0}$ denotes the root
node. The edges in the search tree represent the pruning actions. In the
search tree, the graph associated with a child node can be obtained by
performing node-pruning from the graph corresponding to its parent node. To
limit the search space, we have added two additional constraints:
$\mathcal{N}_{i}$ has to be a connected subgraph and the size of the projected
subgraph should be small.

During the search process, the MCTS algorithm records the statistics of
visiting counts and rewards to guide the exploration and reduce the search
space. Specifically, for the node and pruning action pair ($\mathcal{N}_{i}$,
$a_{j}$), we assume that the subgraph $\mathcal{N}_{j}$ is obtained by action
$a_{j}$ from $\mathcal{N}_{i}$. The MCTS algorithm records four variables for
($\mathcal{N}_{i}$, $a_{j}$):

  * •

C($\mathcal{N}_{i}$, $a_{j}$) denotes the number of counts for selecting
action $a_{j}$ for node $\mathcal{N}_{i}$.

  * •

W($\mathcal{N}_{i}$, $a_{j}$) is the total reward for all ($\mathcal{N}_{i}$,
$a_{j}$) visits.

  * •

Q($\mathcal{N}_{i}$, $a_{j}$) is the averaged reward for multiple visits.

  * •

R($\mathcal{N}_{i}$, $a_{j}$) is the immediate reward for selecting $a_{j}$ on
$\mathcal{N}_{i}$, which is measured by the similarity between the prototype
and the subgraph embedding in this paper. The subgraph embedding is obtained
by encoding the subgraph with the GNN encoder $f$.

Figure 2: An illustration of graph prototype projection with Monte Carlo Tree
Search. The bottom shows one selected path from the root to leaves in the
search tree, which corresponds to one iteration of MCTS. Nodes that are not
selected are ignored for simplicity. For each node, its subgraph is evaluated
by computing the similarity score via GNN Encoder and the similarity function.
In this figure, we show the computation of similarity score for the node
$\mathcal{N}_{1}$ (shown in blue dashed box). In the backward pass, the model
updates the statistics of each node.

Guided by these statistics, MCTS searches for the nearest subgraphs in
multiple iterations. Each iteration consists of two phases. In the forward
pass, MCTS selects a path starting from the root $\mathcal{N}_{0}$ to a leaf
node $\mathcal{N}_{l}$. To keep subgraphs connected, we select to prune
peripheral nodes with minimum degrees. The leaf node can be defined based on
the numbers of nodes in subgraphs such that $|\mathcal{N}_{l}|\leq N_{\rm
min}$. The action selection criteria at node $\mathcal{N}_{i}$ is:

$a^{*}=\mathop{\rm
argmax}\limits_{a_{j}}Q(\mathcal{N}_{i},a_{j})+U(\mathcal{N}_{i},a_{j})$ (8)
$U(\mathcal{N}_{i},a_{j})=\lambda
R(\mathcal{N}_{i},a_{j})\frac{\sqrt{\sum_{k}C(\mathcal{N}_{i},a_{k})}}{1+C(\mathcal{N}_{i},a_{j})},$
(9)

where $\lambda$ is a hyper-parameter to control the trade-off between
exploration and exploitation. The strategy initially prefers to select child
nodes with low visit counts to explore different pruning actions, but
asympotically prefers actions leading to higher similarity scores.

In the backward pass, the statistics of all node and action pairs selected in
this path are updated:

$C(\mathcal{N}_{i},a_{j})=C(\mathcal{N}_{i},a_{j})+1$ (10)
$W(\mathcal{N}_{i},a_{j})=W(\mathcal{N}_{i},a_{j})+sim(p_{k},h_{\mathcal{N}_{l}}),$
(11)

where $h_{\mathcal{N}_{l}}$ is the embedding of the subgraph associated to the
leaf node $\mathcal{N}_{l}$. In the end, we select the subgraph with the
highest similarity score from all the expanded nodes as the new projected
prototype.

### Conditional Subgraph Sampling module

We further propose ProtGNN+ with a novel conditional subgraph sampling module
to provide better interpretation. In ProtGNN+, we not only show the similarity
scores to prototypes, but also identify which part of the input graph is most
similar to each prototype as part of the reasoning process. In Figure 1, the
conditional subgraph sampling module outputs different subgraph embeddings for
each prototype. While this task can also be accomplished by MCTS, the
exponentially-growing time complexity to the graph size and the difficulty of
parallelization and generalization make MCTS algorithm an undesirable choice.
Instead, we adopt a parameterized method for efficient similar subgraph
selection conditioned on given prototypes.

Formally, we let $e_{ij}\in\\{0,1\\}$ be the binary variable indicating
whether the edge between node $i$ and $j$ is selected. The matrix of $e_{ij}$
is denoted as $\mathcal{E}$. The optimization objective of the conditional
subgraph sampling module is:

$\mathop{\rm
max}\limits_{\mathcal{E}}sim(p_{k},f(\mathcal{G}_{s}))~{}s.t.~{}|\mathcal{G}_{s}|\leq
B,$ (12)

where $\mathcal{G}_{s}$ is the selected subgraph whose adjacency matrix is
$\mathcal{E}$. $B$ is the maximum size of the subgraph.

The combinatorial and discrete nature of graph makes the direct optimization
of the above objective function intractable. We first consider a relaxation by
assuming that the explanatory graph is a Gilbert random graph (Gilbert 1959)
where the state of each edge is independent to each other. Furthermore, for
ease of gradient computation and update, we relax
$\mathcal{E}\in\\{0,1\\}^{N\times N}$ into convex space
$\mathcal{E}\in[0,1]^{N\times N}$. $N$ is the number of nodes in the input
graph. For efficiency and generalizability, we adopt deep neural networks to
learn $e_{ij}$:

$e_{ij}=\sigma({\rm MLP}_{\theta}([z_{i};z_{j};p_{k}])),$ (13)

where $\sigma(\cdot)$ here is the Sigmoid function. MLP is a multi-layer
neural network parameterized with $\theta$ and [·; ·; ·] is the concatenation
operation. $z_{i}$ and $z_{j}$ are node embedding obtained from the GNN
Encoder, which encodes the feature and structure information of the nodes’
neighborhood. Then the objective in Eq. (12) becomes

$\displaystyle\mathop{\rm
max}\limits_{\theta}sim(p_{k},f(\mathcal{G}_{s}))-\lambda_{b}R_{b}$ (14)
$\displaystyle R_{b}={\rm ReLU}(\sum_{e_{ij}\in\mathcal{E}}e_{ij}-B),$

where $\lambda_{b}$ is the weight for the budget regularization $R_{b}$. In
our experiments, we adopt stochastic gradient descent to optimize the
objective function.

Comparison with MCTS: Our designed conditional subgraph sampling module is
much more efficient than MCTS and easier for parallel computation. The
parameters of our conditional subgraph sampling module are fixed and
independent of the graph size. To sample from a graph with $|\mathcal{E}|$
edges, the time complexity of our method is $\mathcal{O(|E|)}$. One limitation
of the conditional subgraph sampling module is that it requires additional
training. Therefore, MCTS is still used in the prototype projection step of
ProtGNN+ for the stability of optimization.

### Theorem on Subgraph Sampling

To provide more understandable visualization, ProtGNN+ prunes the input graph
to find the subgraphs most similar to prototypes and then calculates the
similarity scores. Compared with ProtGNN, the subgraph sampling procedure may
affect the classification accuracy. The following theorem provides some
theoretical understanding of how input graph sampling affects classification
accuracy.
Theorem 1: Let $c\circ g_{p}\circ f$ be a ProtoGNN. The embedding of the input
graph is $h$. We assume that the number of prototypes is the same for each
class, and is denoted as $m$. For each class $k$, the weight connection in the
last layer $c$ between a class $k$ prototype and the class $k$ logit is 1, and
that between a non-class $k$ prototype and the class $k$ logit is 0. We denote
$p_{l}^{k}$ as the $l$-th prototype for class $k$ and $h_{l}^{k}$ the
embedding of the pruned subgraph. ProtGNN and ProtGNN+ has the same graph
encoder $f$. We make the following assumptions: there exists some $\delta$
with $0<\delta<1$,

  * •

for the correct class, we have
$\|h-h_{l}^{k}\|_{2}\leq(\sqrt{1+\delta}-1)\|h-p_{l}^{k}\|_{2}$ and
$\|h-p_{l}^{k}\|_{2}\leq\sqrt{1-\delta}$;

  * •

for the incorrect classes,
$\|h-h_{l}^{k}\|_{2}\leq\theta\|h-p_{l}^{k}\|_{2}-\sqrt{\epsilon}$,
$\theta={\rm min}(\sqrt{1+\delta}-1,1-\frac{1}{\sqrt{2-\delta}})$.

For one correctly classified input graph in ProtGNN, if the output logits
between the top-2 classes are at least $2mlog((1+\delta)(2-\delta))$, then
ProtGNN+ can classify the input graph correctly as well.

The intuition behind Theorem 1 is that if the subgraph sampling does not
change the graph embedding too much, ProtGNN+ will have the same correct
predictions as ProtGNN. The proof is included in the appendix.

### Training Procedures

In Algorithm 1, we show the training procedure of ProtGNN/ProtGNN+. Before
training starts, we randomly initialize the model parameters. We let $w_{c}$
be the weight matrix of the fully connected layer $c$ and $w_{c}^{(k,j)}$ be
the weight connection between the output of the $j$-th prototype and the logit
of class $k$. In particular, for a class k, we set $w_{c}^{(k,j)}=1$ for all
$j$ with $p_{j}\in P_{k}$ and $w_{c}^{(k,j)}=0$ for all $j$ with $p_{j}\notin
P_{k}$. Intuitively, such initialization of $w_{h}$ encourages prototypes
belonging to class $k$ to learn semantic concepts that are characteristic to
class $k$. After training begins, we employ gradient descents to optimize the
objective function in Eq. (3). If the training epoch is larger than the
projection epoch $T_{p}$, we perform the prototype projection step every few
training epochs. Furthermore, if we train ProtGNN+, the conditional subgraph
sampling module and ProtGNN are iteratively optimized after the warm-up epoch
$T_{w}$ when the optimization of GNN encoder and prototypes are stabilized.

Algorithm 1 Overview of ProtGNN/ProtGNN+ Training

Input: Training dataset $\\{x_{i},y_{i}\\}_{i=1}^{n}$
Parameter: Training epochs $T$, Warm-up epoch $T_{w}$, Projection epoch
$T_{p}$, Prototype projection period $\tau$, ProtGNN+

1: Initialize model parameters.

2: for training epochs $t=1,2,\cdots,T$ do

3: Optimizing objective function in Eq. (3)

4: if $t\textgreater T_{p}$ and $t\%\tau=0$ then

5: Performing prototype projection with MCTS

6: end if

7: if ProtGNN+ enabled and $t\textgreater T_{w}$ then

8: Optimizing the objective function in Eq. (14).

9: end if

10: end forOutput: Trained model, prototype visualization

### ProtGNN for Generic Graph Tasks

In the above sections and illustrations, we have described ProtGNN/ProtGNN+
using graph classification as an example. It is worth noting that
ProtGNN/ProtGNN+ can be easily generalized to other graph tasks, such as node
classification and link prediction. For example, in the node classification
task, the explanation target is to provide the reasoning process behind the
prediction of node $v_{i}$. Assuming the GNN encoder has $L$ layers, the
prediction of node $v_{i}$ only relies on its $L$-hop computation graph.
Therefore, prototype projection and conditional subgraph sampling are all
performed in the $L$-hop computation graph.

## Experimental Evaluation

Figure 3: The reasoning process of ProtGNN+ in deciding whether the molecular
graph is mutagenic (a) or the sentiment of input text graph is positive (b).
The predictions are based on the similarity between the latent input
representations against the prototypes. The network tries to find evidence by
looking at which subgraph was mostly similar to the prototypes. The selected
subgraphs are highlighted. Due to space constraint, we only show several
prototypes with the largest weights for each class.

### Datasets and Experimental Settings

Table 1: The classification accuracies and standard deviations ($\%$) of ProtGNN, ProtGNN+, and the original GNNs.  Datasets | GCN | GIN | GAT
---|---|---|---
Original | ProtGNN | ProtGNN+ | Original | ProtGNN | ProtGNN+ | Original | ProtGNN | ProtGNN+
MUTAG | 73.3$\pm 5.8$ | 76.7$\mathbf{\pm 6.4}$ | 73.3$\pm 2.9$ | 93.3$\mathbf{\pm 2.9}$ | 90.7$\pm 3.2$ | 91.7$\pm 2.9$ | 75.0$\pm 5.0$ | 78.3$\pm 4.2$ | 81.7$\mathbf{\pm 2.9}$
BBBP | 84.6$\pm 3.4$ | 89.4$\mathbf{\pm 4.1}$ | 88.0$\pm 4.6$ | 86.2$\pm 1.1$ | 86.5$\mathbf{\pm 1.6}$ | 85.9$\pm 4.0$ | 83.0$\pm 2.6$ | 85.9$\mathbf{\pm 2.5}$ | 85.5$\pm 0.8$
Graph-SST2 | 89.7$\pm 0.5$ | 89.9$\mathbf{\pm 2.4}$ | 89.0$\pm 3.0$ | 92.2$\pm 0.3$ | 92.0$\pm 0.2$ | 92.3$\mathbf{\pm 0.4}$ | 88.1$\pm 0.8$ | 89.1$\mathbf{\pm 1.2}$ | 88.7$\pm 0.9$
Graph-Twitter | 67.5$\pm 1.9$ | 68.9$\mathbf{\pm 5.9}$ | 66.1$\pm 6.5$ | 66.2$\pm 1.3$ | 75.2$\pm 2.8$ |  76.5 $\mathbf{\pm 3.4}$ | 69.6$\mathbf{\pm 6.5}$ | 64.8$\pm 4.0$ | 66.4$\pm 3.3$
BA-Shape | 91.9$\pm 1.7$ | 95.7$\mathbf{\pm 1.4}$ | 94.3$\pm 3.7$ | 92.9$\pm 0.5$ | 95.2$\pm 1.3$ |  95.5 $\mathbf{\pm 2.4}$ | 92.9$\pm 1.2$ | 93.4$\mathbf{\pm 3.4}$ | 93.2$\pm 2.0$

Datasets: We conduct extensive experiments on different datasets and GNN
models to demonstrate the effectiveness of our proposed model. These datasets
are listed as below:

  * •

MUTAG (Debnath et al. 1991) and BBBP (Wu et al. 2018) are molecule datasets
for graph classification. In these datasets, nodes represent atoms and edges
denote chemical bonds. The labels of molecular graphs are determined by the
molecular compositions.

  * •

Graph-SST2 (Socher et al. 2013) and Graph-Twitter (Dong et al. 2014) are
sentiment graph datasets for graph classification. They convert sentences to
graphs with Biaffine parser (Gardner et al. 2018) that nodes denote words and
edges represent the relationships between words. The node embeddings are
initialized with Bert word embeddings (Devlin et al. 2018). The labels are
determined by the sentiment of text sentences.

  * •

BA-Shape is a synthetic node classification dataset. Each graph contains a
base graph obtained from the Barabási-Albert (BA) mode (Albert and Barabási
2002) and a house-like five-node motif attached to the base graph. Each node
is labeled based on whether it belongs to the base graph or the different
spatial locations of the motif.

Experimental Settings: In our evaluation, we use three variants of GNNs, i.e.
GCN, GAT, and GIN. The split for train/validation/test sets is
$80\%:10\%:10\%$. All models are trained for 500 epochs with an early stopping
strategy based on accuracy on the validation set. We adopt the ADAM optimizer
with a learning rate of 0.005. In Eq.(3), the hyper-parameters $\lambda_{1}$,
$\lambda_{2}$, and $\lambda_{3}$ are set to 0.10, 0.05, and 0.01 respectively.
$s_{max}$ is set to 0.3 in Eq. (6). The number of prototypes per class $m$ is
set to 5. In MCTS for prototype projection, we set $\lambda$ in Eq. (9) to 5
and the number of iterations to 20. Each node in the Monte Carlo Tree can
expand up to 10 child nodes and $N_{\rm min}$ is set to 5. The prototype
projection period $\tau$ is set to 50 and the projection epoch $T_{p}$ is set
to 100. In the training of ProtGNN+, the warm-up epoch $T_{w}$ is set to 200.
We employ a three-layer neural network to learn edge weights. In Eq. (14),
$\lambda_{b}$ is set to 0.01 and $B$ is set to 10. We select hyper-parameters
based on related works or grid search, an analysis on hyper-parameters is
included in the appendix. All our experiments are conducted with one Tesla
V100 GPU.

### Evaluations on ProtGNN/ProtGNN+

#### Comparison with Baselines

In Table 1, we compare the classification accuracy of ProtGNN/ProtGNN+ with
the original GNNs. We apply 3 independent runs on random data splitting and
report the means and standard deviations. In the following sections, we use
GCN as the default backbone model. As we can see, ProtGNN and ProtGNN+ achieve
comparable classification performance with the corresponding original GNN
models, which also empirically verifies Theorem 1.

Figure 4: Visualization on BBBP dataset of the graph (dots) and prototype
(stars) embeddings using the t-SNE method. Different colors indicate different
classes.

#### Reasoning Process of Our Network

In Figure 3, we perform case studies on MUTAG and Graph-SST2 to qualitatively
evaluate the performance of our proposed method. We visualize the prototypes
and show the reasoning process of our ProtGNN+ in reaching a classification
decision on input graphs. In particular, given an input graph $x$, the network
finds the likelihood to be in each class by comparing it with prototypes from
each class. The conditional subgraph sampling module finds the most similar
subgraphs in $x$. These similarity scores are calculated, weighted, and summed
together to give a final score for $x$ belonging to each class. For example,
Figure 3(a) shows the reasoning process of ProtGNN+ in deciding whether the
input molecular graph is mutagenic. Based on chemical domain knowledge
(Debnath et al. 1991), carbon rings and $NO_{2}$ groups tend to be mutagenic.
In the Prototype column of the mutagenic class, we can observe that the
prototypes can capture the structures of $NO_{2}$ and carbon rings well.
Moreover, in the column of Similar Subgraphs, the conditional subgraph
sampling module can effectively identify the most similar subgraphs. For
instance, in the first row of the mutagenic class, the $NO_{2}$ group and part
of the carbon ring can be identified, which is quite similar to the prototype
on the right.

Compared with biochemistry datasets, examples on text data could be more
understandable since no special domain knowledge is required. In Figure 3(b),
the input graph “can take the grandkids or the grandparents and never worry
about anyone being bored” is positive. Our method can effectively capture the
key phrase/subgraph leading to positiveness, “never worry about bored”.
Furthermore, we can observe that the similarity score between the input graph
with the positive prototypes e.g., “kind of entertainment love to have” is
much larger than negative prototypes e.g., “embarrassed by invention”.

Overall, our method provides interpretable evidence to support
classifications. Such explanations participate in the actual model computation
and is always faithful to the classification decisions. More examples and case
studies are reported in appendix.

#### t-SNE Visualization of Prototypes

In Figure 4 we show the visualization on BBBP dataset of the graph and
prototype embeddings using t-SNE method. We can observe that the prototypes
can occupy the centers of graph embeddings, which verifies the effectiveness
of prototype learning.

Table 2: Efficiency studies of different methods on BBBP Methods | GCN | ProtGNN | ProtGNN+ | ProtGNN+*
---|---|---|---|---
Time | 177.9 s | 506.3 s | 632.7 s |  $\textgreater$ 2 hrs

#### Efficiency Studies

Finally, we study the efficiency of our proposed methods. In Table 2, we show
the time required to finish training for each model. Here ProtGNN+* denotes
using MCTS for subgraph sampling in the training of ProtGNN+. The time
complexity of ProtGNN+* is extremely high due to the complexity of MCTS. The
proposed conditional subgraph sampling module can effectively reduce the time
cost of ProtGNN+. Although ProtGNN and ProtGNN+ have a larger time cost
compared to GCN (largely due to prototype projection with MCTS), the time cost
is still acceptable considering the provided built-in interpretability.

## Conclusion

While extensive efforts have been made to explain GNNs from different angles,
none of existing methods can provide _built-in_ explanations for GNNs. In this
paper, we propose ProtGNN/ProtGNN+ which provides a new perspective on the
explanations of GNNs. The prediction of ProtGNN is obtained by comparing the
inputs to a few learned prototypes in the prototype layer. For better
interpretability and higher efficiency, a novel conditional subgraph sampling
module is proposed to indicate the subgraphs most similar to prototypes.
Extensive experimental results show that our method can provide a human-
intelligible reasoning process with acceptable classification accuracy and
time-complexity.

## References

  * Albert and Barabási (2002) Albert, R.; and Barabási, A.-L. 2002.  Statistical mechanics of complex networks.  _Reviews of modern physics_ , 74(1): 47.
  * Baldassarre and Azizpour (2019) Baldassarre, F.; and Azizpour, H. 2019.  Explainability techniques for graph convolutional networks.  _ICML workshop_.
  * Bian et al. (2020) Bian, T.; Xiao, X.; Xu, T.; Zhao, P.; Huang, W.; Rong, Y.; and Huang, J. 2020.  Rumor detection on social media with bi-directional graph convolutional networks.  In _AAAI_ , volume 34, 549–556.
  * Chen et al. (2018) Chen, C.; Li, O.; Tao, C.; Barnett, A. J.; Su, J.; and Rudin, C. 2018.  This looks like that: deep learning for interpretable image recognition.  _NeurIPS_.
  * Debnath et al. (1991) Debnath, A. K.; Lopez de Compadre, R. L.; Debnath, G.; Shusterman, A. J.; and Hansch, C. 1991.  Structure-activity relationship of mutagenic aromatic and heteroaromatic nitro compounds. correlation with molecular orbital energies and hydrophobicity.  _Journal of medicinal chemistry_ , 34(2): 786–797.
  * Devlin et al. (2018) Devlin, J.; Chang, M.-W.; Lee, K.; and Toutanova, K. 2018.  Bert: Pre-training of deep bidirectional transformers for language understanding.  _arXiv preprint arXiv:1810.04805_.
  * Dong et al. (2014) Dong, L.; Wei, F.; Tan, C.; Tang, D.; Zhou, M.; and Xu, K. 2014.  Adaptive recursive neural network for target-dependent twitter sentiment classification.  In _ACL_ , 49–54.
  * Gardner et al. (2018) Gardner, M.; Grus, J.; Neumann, M.; Tafjord, O.; Dasigi, P.; Liu, N.; Peters, M.; Schmitz, M.; and Zettlemoyer, L. 2018.  Allennlp: A deep semantic natural language processing platform.  _arXiv preprint arXiv:1803.07640_.
  * Gilbert (1959) Gilbert, E. N. 1959.  Random graphs.  _The Annals of Mathematical Statistics_ , 30(4): 1141–1144.
  * Gilmer et al. (2017) Gilmer, J.; Schoenholz, S. S.; Riley, P. F.; Vinyals, O.; and Dahl, G. E. 2017.  Neural message passing for quantum chemistry.  In _ICML_ , 1263–1272. PMLR.
  * He et al. (2010) He, Z.; Zhang, J.; Shi, X.-H.; Hu, L.-L.; Kong, X.; Cai, Y.-D.; and Chou, K.-C. 2010\.  Predicting drug-target interaction networks based on functional groups and biological features.  _PloS one_ , 5(3): e9603.
  * Huang et al. (2020) Huang, Q.; Yamada, M.; Tian, Y.; Singh, D.; Yin, D.; and Chang, Y. 2020.  Graphlime: Local interpretable model explanations for graph neural networks.  _arXiv preprint arXiv:2001.06216_.
  * Kipf and Welling (2017) Kipf, T. N.; and Welling, M. 2017.  Semi-supervised classification with graph convolutional networks.  _ICLR_.
  * Kolodner (1992) Kolodner, J. L. 1992.  An introduction to case-based reasoning.  _Artificial intelligence review_ , 6(1): 3–34.
  * Luo et al. (2020) Luo, D.; Cheng, W.; Xu, D.; Yu, W.; Zong, B.; Chen, H.; and Zhang, X. 2020.  Parameterized explainer for graph neural network.  _NeurIPS_.
  * Ming et al. (2019) Ming, Y.; Xu, P.; Qu, H.; and Ren, L. 2019.  Interpretable and steerable sequence learning via prototypes.  In _SIGKDD_.
  * Pope et al. (2019) Pope, P. E.; Kolouri, S.; Rostami, M.; Martin, C. E.; and Hoffmann, H. 2019.  Explainability methods for graph convolutional neural networks.  In _CVPR_ , 10772–10781.
  * Rudin (2018) Rudin, C. 2018.  Please stop explaining black box models for high stakes decisions.  _stat_ , 1050: 26.
  * Rymarczyk et al. (2021) Rymarczyk, D.; Struski, Ł.; Tabor, J.; and Zieliński, B. 2021.  ProtoPShare: Prototype Sharing for Interpretable Image Classification and Similarity Discovery.  _SIGKDD_.
  * Schlichtkrull, De Cao, and Titov (2020) Schlichtkrull, M. S.; De Cao, N.; and Titov, I. 2020.  Interpreting graph neural networks for nlp with differentiable edge masking.  _arXiv preprint arXiv:2010.00577_.
  * Schmidt et al. (2001) Schmidt, R.; Montani, S.; Bellazzi, R.; Portinale, L.; and Gierl, L. 2001.  Cased-based reasoning for medical knowledge-based systems.  _International Journal of Medical Informatics_ , 64(2-3): 355–367.
  * Schnake et al. (2020) Schnake, T.; Eberle, O.; Lederer, J.; Nakajima, S.; Schütt, K. T.; Müller, K.-R.; and Montavon, G. 2020.  XAI for graphs: explaining graph neural network predictions by identifying relevant walks.  _arXiv e-prints_ , arXiv–2006.
  * Schwarzenberg et al. (2019) Schwarzenberg, R.; Hübner, M.; Harbecke, D.; Alt, C.; and Hennig, L. 2019.  Layerwise relevance visualization in convolutional text graph classifiers.  _arXiv preprint arXiv:1909.10911_.
  * Silver et al. (2017) Silver, D.; Schrittwieser, J.; Simonyan, K.; Antonoglou, I.; Huang, A.; Guez, A.; Hubert, T.; Baker, L.; Lai, M.; Bolton, A.; et al. 2017.  Mastering the game of go without human knowledge.  _nature_ , 550(7676): 354–359.
  * Socher et al. (2013) Socher, R.; Perelygin, A.; Wu, J.; Chuang, J.; Manning, C. D.; Ng, A. Y.; and Potts, C. 2013.  Recursive deep models for semantic compositionality over a sentiment treebank.  In _EMNLP_ , 1631–1642.
  * Veličković et al. (2017) Veličković, P.; Cucurull, G.; Casanova, A.; Romero, A.; Lio, P.; and Bengio, Y. 2017.  Graph attention networks.  _arXiv preprint arXiv:1710.10903_.
  * Vu and Thai (2020) Vu, M.; and Thai, M. T. 2020.  PGM-Explainer: Probabilistic Graphical Model Explanations for Graph Neural Networks.  In Larochelle, H.; Ranzato, M.; Hadsell, R.; Balcan, M. F.; and Lin, H., eds., _NeurIPS_ , volume 33.
  * Wu et al. (2018) Wu, Z.; Ramsundar, B.; Feinberg, E. N.; Gomes, J.; Geniesse, C.; Pappu, A. S.; Leswing, K.; and Pande, V. 2018.  MoleculeNet: a benchmark for molecular machine learning.  _Chemical science_ , 9(2): 513–530.
  * Xu et al. (2019) Xu, K.; Hu, W.; Leskovec, J.; and Jegelka, S. 2019.  How Powerful are Graph Neural Networks?  In _ICLR_.
  * Yang et al. (2020) Yang, S.; Zhang, Z.; Zhou, J.; Wang, Y.; Sun, W.; Zhong, X.; Fang, Y.; Yu, Q.; and Qi, Y. 2020.  Financial Risk Analysis for SMEs with Graph-based Supply Chain Mining.  In _IJCAI_ , 4661–4667.
  * Ying et al. (2019) Ying, R.; Bourgeois, D.; You, J.; Zitnik, M.; and Leskovec, J. 2019.  Gnnexplainer: Generating explanations for graph neural networks.  _NeurIPS_ , 32: 9240.
  * Yuan et al. (2020a) Yuan, H.; Tang, J.; Hu, X.; and Ji, S. 2020a.  Xgnn: Towards model-level explanations of graph neural networks.  In _SIGKDD_ , 430–438.
  * Yuan et al. (2020b) Yuan, H.; Yu, H.; Gui, S.; and Ji, S. 2020b.  Explainability in graph neural networks: A taxonomic survey.  _arXiv preprint arXiv:2012.15445_.
  * Yuan et al. (2021) Yuan, H.; Yu, H.; Wang, J.; Li, K.; and Ji, S. 2021.  On explainability of graph neural networks via subgraph explorations.  _ICML_.
  * Zhang et al. (2021) Zhang, Z.; Liu, Q.; Wang, H.; Lu, C.; and Lee, C.-K. 2021.  Motif-based Graph Self-Supervised Learning for Molecular Property Prediction.  _NeurIPS_ , 34.

Table 3: Statistics of five datasets | Datasets
---|---
MUTAG | BBBP | Graph-SST2 | Graph-Twitter | BA-Shape
$\\#$ of Edges (avg) | 19.79 | 25.95 | 9.20 | 20.10 | 2055
$\\#$ of Nodes (avg) | 17.93 | 24.06 | 10.19 | 21.10 | 700
$\\#$ of Graphs | 188 | 2039 | 70042 | 6940 | 1
$\\#$ of Classes | 2 | 2 | 2 | 3 | 4

##  Appendix A Dataset Statistics

In Table 3, we show the detailed statistics of five datasets. These datasets
include biological data, text data, and synthetic data. The first four
datasets are used for graph classification tasks while BA-Shape is used for
node classfication.

Figure 5: The reasoning process of ProtGNN+ in deciding whether the input
molecular graph is penetrative.

Figure 6: The reasoning process of ProtGNN+ in deciding whether the sentiment
of the text graph is positive, neutral, or negative.

##  Appendix B Proof of Theorem 1

In this section, we provide a proof for Theorem 1 in the main paper.
Theorem 1: Let $c\circ g_{p}\circ f$ be a ProtoGNN. The embedding of the input
graph is $h$. We assume that the number of prototypes is the same for each
class, which is denoted as $m$. For each class $k$, the weight connection in
the last layer $c$ between a class $k$ prototype and the class $k$ logit is 1,
and that between a non-class $k$ prototype and the class $k$ logit is 0. Let
$p_{l}^{k}$ denote the $l$-th prototype for class $k$ and $h_{l}^{k}$ the
embedding of the pruned subgraph. ProtGNN and ProtGNN+ has the same graph
encoder $f$.
We make the following assumptions: there exists some $\delta$ with
$0<\delta<1$,

  * •

for the correct class, we have
$\|h-h_{l}^{k}\|_{2}\leq(\sqrt{1+\delta}-1)\|h-p_{l}^{k}\|_{2}$ and
$\|h-p_{l}^{k}\|_{2}\leq\sqrt{1-\delta}$;

  * •

for the incorrect classes,
$\|h-h_{l}^{k}\|_{2}\leq\theta\|h-p_{l}^{k}\|_{2}-\sqrt{\epsilon}$,
$\theta={\rm min}(\sqrt{1+\delta}-1,1-\frac{1}{\sqrt{2-\delta}})$.

For one correctly classified input graph in ProtGNN, if the output logits
between the top-2 classes are at least $2mlog((1+\delta)(2-\delta))$, then
ProtGNN+ can classify the input graph correctly as well.

Proof: For any class $k$, let $L_{k}(x,\\{p_{l}^{k}\\}_{l=1}^{m})$ denotes the
summed contributed scores for graph $x$ belonging to class $k$ in ProtGNN.
According to Eq. (2) and the assumption:

$L_{k}(x,\\{p_{l}^{k}\\}_{l=1}^{m})=\sum_{l=1}^{m}log(\frac{\|h-p_{l}^{k}\|_{2}^{2}+1}{\|h-p_{l}^{k}\|_{2}^{2}+\epsilon}).$
(15)

Let $L^{\prime}_{k}(x,\\{p_{l}^{k}\\}_{l=1}^{m})$ denotes the summed
contributed scores in ProtGNN+:

$L^{\prime}_{k}(x,\\{p_{l}^{k}\\}_{l=1}^{m})=\sum_{l=1}^{m}log(\frac{\|h_{l}^{k}-p_{l}^{k}\|_{2}^{2}+1}{\|h_{l}^{k}-p_{l}^{k}\|_{2}^{2}+\epsilon}).$
(16)

Then, the gap between the summend contributed scores denoted by $\Delta_{k}$
is:

$\displaystyle\Delta_{k}$
$\displaystyle=L^{\prime}_{k}(x,\\{p_{l}^{k}\\}_{l=1}^{m})-L_{k}(x,\\{p_{l}^{k}\\}_{l=1}^{m})$
(17)
$\displaystyle=\sum_{l=1}^{m}log(\frac{\|h_{l}^{k}-p_{l}^{k}\|_{2}^{2}+1}{\|h-p_{l}^{k}\|_{2}^{2}+1}\cdot\frac{\|h-p_{l}^{k}\|_{2}^{2}+\epsilon}{\|h_{l}^{k}-p_{l}^{k}\|_{2}^{2}+\epsilon}).$

Correct class: We first derive the lower bound of $\Delta_{k}$ for the correct
class. Firstly, we have

$\frac{\|h_{l}^{k}-p_{l}^{k}\|_{2}^{2}+1}{\|h-p_{l}^{k}\|_{2}^{2}+1}\geq\frac{1}{\|h-p_{l}^{k}\|_{2}^{2}+1}\geq\frac{1}{2-\delta}$
(18)

Then, by the triangle inequality , we have
$\|h_{l}^{k}-p_{l}^{k}\|\leq\|h-p_{l}^{k}\|+\|h-h_{l}^{k}\|$. As a result, we
have:

$\displaystyle\frac{\|h-p_{l}^{k}\|_{2}^{2}+\epsilon}{\|h_{l}^{k}-p_{l}^{k}\|_{2}^{2}+\epsilon}$
$\displaystyle\geq\frac{\|h-p_{l}^{k}\|_{2}^{2}+\epsilon}{(\|h-p_{l}^{k}\|_{2}+\|h-h_{l}^{k}\|_{2})^{2}+\epsilon}$
(19)
$\displaystyle\geq\frac{\|h-p_{l}^{k}\|_{2}^{2}+\epsilon}{(1+\delta)\|h-p_{l}^{k}\|_{2}^{2}+\epsilon}$
$\displaystyle\geq\frac{1}{1+\delta}$

Combining the above two inequalities, $\Delta_{k}$ for the correct class is
$m{\rm log}(\frac{1}{(1+\delta)(2-\delta)})=-m{\rm
log}((1+\delta)(2-\delta))$.
Wrong class: Now we begin to derive an upper bound of $\Delta_{k}$ for
incorrect classes. First,

$\displaystyle\frac{\|h_{l}^{k}-p_{l}^{k}\|_{2}^{2}+1}{\|h-p_{l}^{k}\|_{2}^{2}+1}$
$\displaystyle\leq\frac{(\|h-p_{l}^{k}\|_{2}+\|h-h_{l}^{k}\|)^{2}+1}{\|h-p_{l}^{k}\|_{2}^{2}+1}$
(20)
$\displaystyle\leq\frac{(\|h-p_{l}^{k}\|_{2}+(\sqrt{1+\delta}-1)\|h-p_{l}^{k}\|_{2})^{2}+1}{\|h-p_{l}^{k}\|_{2}^{2}+1}$
$\displaystyle\leq 1+\delta$

Then,

$\displaystyle\frac{\|h-p_{l}^{k}\|_{2}^{2}+\epsilon}{\|h_{l}^{k}-p_{l}^{k}\|_{2}^{2}+\epsilon}$
$\displaystyle\leq\frac{\|h-p_{l}^{k}\|_{2}^{2}+\epsilon}{(\|h-p_{l}^{k}\|_{2}-\|h-h_{l}^{k}\|_{2})^{2}+\epsilon}$
(21)
$\displaystyle\leq(\frac{\|h-p_{l}^{k}\|_{2}+\sqrt{\epsilon}}{\|h-p_{l}^{k}\|_{2}-\|h-h_{l}^{k}\|_{2}})^{2}$

According to the assumption for incorrect classes, we have:

$\|h-h_{l}^{k}\|_{2}\leq(1-\frac{1}{\sqrt{2-\delta}})\|h-p_{l}^{k}\|_{2}-\sqrt{\epsilon}$
(22)
$\frac{1}{\sqrt{2-\delta}}\|h-p_{l}^{k}\|_{2}+\sqrt{\epsilon}\leq\|h-p_{l}^{k}\|_{2}-\|h-h_{l}^{k}\|_{2}$
(23)

Combining Eq. (21) and Eq. (23), we have:

$\frac{\|h-p_{l}^{k}\|_{2}^{2}+\epsilon}{\|h_{l}^{k}-p_{l}^{k}\|_{2}^{2}+\epsilon}\leq(\sqrt{2-\delta})^{2}=2-\delta$
(24)

For incorrect classes, $\Delta_{k}\geq m{\rm log}((1+\delta)(2-\delta))$.

Finally, suppose the summed contributed scores of the correct class is at
least $2m{\rm log}((1+\delta)(2-\delta))$ larger than any other classes in
ProtGNN, the input graph will still be correctly classified by ProtGNN+.
$\hfill\qed$

##  Appendix C Architecture of the Conditional Subgraph Sampling Module

Table 4: Architecture of the Conditional Subgraph Sampling Module Layer | Size
---|---
Input | 128 + 128 + 128
Fully Connected + ReLU | 64
Fully Connected + ReLU | 8
Sigmoid | 1

In the conditional subgraph sampling module, we adopt deep neural networks to
learn $e_{ij}$:

$e_{ij}=\sigma({\rm MLP}_{\theta}([z_{i};z_{j};p_{k}])).$ (25)

In Table 4, we show the details of architecture. In our experiments, the node
embedding size and prototype size are set to 128. To make sure the selected
adjacency matrix is symmetric, we set $\mathcal{E}$ as
$\frac{\mathcal{E}+\mathcal{E}^{T}}{2}$ in experiments.

##  Appendix D More Case Studies

(a) BBBP (b) Graph-Twitter

Figure 7: The influence of the number of prototypes per class $m$ on
performance. The accuracy of ProtGNN firstly increase dramatically as $m$
increases. Then the increasing slope flattens after $m$ exceeds 5.

(a) w/ Div Class 0 (b) w/ Div Class 1 (c) w/o Div Class 0 (d) w/o Div Class 1

Figure 8: The influence of the diversity regularization in Eq. (6). In the
first row, we show the cosine similarities among prototypes learned with the
diversity loss. In the Second row, the diversity loss is removed from the
learning objective. The cosine similarities among prototypes without diversity
regularization are much larger than those with diversity loss. In some extreme
cases, the similarities are close to 1, which means the learned prototypes are
nearly the same. Figure 9: Effects of the cluster and separation losses on
BBBP. When $\lambda_{1}=0.10$ and $\lambda_{2}=0.05$, ProtGNN achieves the
best performance on BBBP dataset.

In Figure 5 and Figure 6, we show more case studies on BBBP and Graph-Twitter.
Note that Graph-Twitter is a 3-class dataset and we show the prototypes for
negativeness, neutrality, and positiveness. The input graph in Figure 6 is
positive. Our method can effectively capture the key phrase/subgraph leading
to positiveness, “amazing lady gaga I love”.

##  Appendix E Hyper-parameters Analysis

In this section, we provide some analysis on hyper-parameters in
ProtGNN/ProtGNN+.

### Choosing the Number of Prototypes per Class

We first investigate how would the number of prototypes per class $m$
influence the performance of ProtGNN using BBBP and Graph-Twitter. With the
default setting of hyper-parameters, we train ProtGNN with varying $m$. In
Figure 7, we observe that the accuracy of ProtGNN firstly increase
dramatically as $m$ increases. Then the increasing slope flattens after $m$
exceeds 5.

Actually, there is one trade-off between accuracy and interpretability when
choosing $m$. The accuracy increases when $m$ increases. However, a large
number of prototypes makes the model difficult to train and comprehend. In
experiments, we choose $m=5$ since the increasing $m$ only brings marginal
improvement to the performance.

### Influence of Diversity Loss

We further show the effectiveness of the diversity loss (Eq. 6). In Figure 8,
we plot the cosine similarity matrices of the learned prototypes on BBBP. The
first row are similarity matrices with the diversity loss while the second row
without the diversity loss. We can observe that the cosine similarities among
prototypes without diversity regularization are much larger than those with
diversity loss. In some extreme cases, the similarities are close to 1, which
means the learned prototypes are nearly the same. Therefore, the diversity
loss can help ProtGNN learn more diverse and evenly distributed prototypes.

### Influence of the Cluster and Separation Loss

Here we want to show the influence of $\lambda_{1}$ and $\lambda_{2}$ which
controls the weights of the cluster loss and separation loss respectively. In
Figure 9, we can observe that the cluster and separation constraints play
important roles in ProtGNN. When $\lambda_{1}=0.10$ and $\lambda_{2}=0.05$,
ProtGNN achieves the best performance on BBBP dataset.


# WIMP Dark Matter Hidden behind its Companion

Jun Guo<EMAIL_ADDRESS>College of Physics and Communication Electronics,
Jiangxi Normal University, Nanchang 330022, China  Zhaofeng Kang
<EMAIL_ADDRESS>School of physics, Huazhong University of Science and
Technology, Wuhan 430074, China  Peng Zhang<EMAIL_ADDRESS>School
of physics, Huazhong University of Science and Technology, Wuhan 430074, China

###### Abstract

The WIMP dark matter (DM) hypothesis now is in an awkward position, owing to
the stronger and stronger exclusion from DM direct detection. In this article
we design a mechanism to evade this constraint. The idea is simple. DM has a
companion, and they are both charged under the DM protecting symmetry $G$;
they admit the trilinear coupling DM-DM-companion, so the latter provides a
portal to the standard model (SM) via, for instance, the coupling to Higgs
doublet. Then, DM semi-annihilates into the companion to arrive correct relic
density, without leaving DM-nucleon scattering signal. The idea can be
realized for $Z_{N}$ symmetric models with $N>2$. We stress that this
mechanism has the characteristics of co-annihilation, and as a matter of fact
its effect becomes necessary near or above the TeV region. This means that it
may be difficult to detect our dark matter directly or indirectly.

###### pacs:

12.60.Jv, 14.70.Pw, 95.35.+d

##  I Introduction

The existence of dark matter (DM) is one of the widely accepted evidences
beyond the particle standard model (SM) Bertone:2004pz . But so far, the only
accurate quantitative information we have about DM is its relic density,
$\Omega h^{2}\simeq 0.12$ Planck:2018vyg . This illustrative value inspires
the well-known weakly interacting massive particle (WIMP) picture for DM,
which is supposed to be a neutral particle around the weak scale, interacting
with the SM members with a strength $\sim{\mathcal{O}}(g_{2})$, thus
demonstrating the WIMP miracle. In such a paradigm, barring parameters fine-
tuning and peculiar structure for interactions, $\sigma_{\rm SI}$, the cross
section of DM-nucleon spin-independent scattering should be large enough,
leaving detectable signals which would have already been hunted by the
underground experiments such as XENON Aprile:2018dbl  and PandaX Meng:2021mui
, etc. Nevertheless, we just see the experimental exclusion lines on
$\sigma_{\rm SI}$ approaching the neutrino floor, below which DM signals will
be buried in the neutrino backgrounds and become nearly undetectable.

The decreasing upper bound on $\sigma_{\rm SI}$ means that we may have to face
such a fact: the WIMP DM paradigm probably is not correct, at least does not
assemble the typical one. This may be a good news to the feebly interacting
(with the members of the SM) massive particle scenario McDonald:2001vt ;
Hall:2009bx ; Bernal:2017kxu ; Kang:2014cia , which gives up the freeze-out
dynamics to determine DM relic density and assumes that DM has never reached
thermal equilibrium with the SM bath, but it is still capable of producing
correct relic density via the IR or UV-freeze-in mechanism Chen:2017kvz . Even
within the WIMP framework, we still have many ways to maintain the WIMP
miracle while avoiding the very stringent direct detection bounds, for
instance, the isospin-violating DM Feng:2011vu ; Kang:2010mh ; Gao:2011ka ,
and the pseudo-Goldstone DM Gross:2017dan ; Okada:2020zxo ; Jiang:2019soj .

This article is devoted to exploring a new mechanism. The idea is inspired by
the structure of the Lagrangian for semi-annihilating DM models DEramo:2010keq
; Belanger:2012vp ; Ko:2014nha , to which a generalization admits the
separation between DM annihilation cross section and $\sigma_{\rm SI}$. For
illustration, let us consider the $Z_{3}$ symmetric model, where we introduce
two identical scalars $S_{1,2}$, rather than merely one in the minimal
version; both of them are transformed as $e^{ik2\pi/3}S_{1,2}$. Let $S_{1}$ be
the DM, and its companion $S_{2}$ is properly heavier than $S_{1}$. These
ingredients open a possibility to hide $S_{1}$ DM: $S_{1}$ has a negligible
Higgs portal coupling, so its scattering with nuclei is highly suppressed; the
companion $S_{2}$ serves as the portal (dubbed as the companion portal)
between $S_{1}$ and the SM Higgs sector, so the correct relic density of
$S_{1}$ can be achieved by the semi-annihilation mode $S_{1}S_{1}\rightarrow
S_{2}(\rightarrow S_{1}+\bar{f}f)h$. This structure can be easily generalized
to $Z_{N}$ with $N>3$. In particular, if $S_{1}$ instead is a spin-1/2
fermion, its direct coupling with the SM Higgs sector is forbidden,
naturalizing the absence of $\sigma_{\rm SI}$, which is unlike the previous
version that requires one to squeeze the Higgs portal by hand.

We must point out a related work Saez:2021oxl , which appeared on arxiv at the
final stage of our article. Authors studied the assisted freeze-out mechanism
Belanger:2011ww , which actually is a two components DM model with the
dominant one $X_{1}$ feebly interacting with the SM while the subdominant one
$X_{2}$ opening the annihilation channel for $X_{1}$ freeze-out. Their
hypothesis has a similar picture to ours in treating the freeze-out dynamics
of WIMP DM. However, the difference in the model building is obvious: we start
from a single symmetry whereas they should impose more, and in our approach
the companion for assisted freeze-out also transforms under this symmetry. In
fact, we find that the assisted freeze-out mechanism can be easily integrated
into our mechanism, for $N>3$.

If direct detection is no longer feasible, we must rely on other detection
methods, such as looking for relevant signals in cosmic rays or on colliders.
However, this requires professional technology, which we do not fully have at
present, so the detailed analysis is postponed to the follow-up work. For the
indirect detection, the main difficulty is to obtain the precise spectra of DM
annihilation, a pair of Higgs boson with one (on or off shell) from the
companion cascade decay; such spectra is not regular and needs simulation.
Note that our mechanism may have a strong co-annihilation effect, in
particular for DM near the TeV scale, then the resulting WIMP DM is also
hidden in the sky.

The paper is organized as follows: In Section II we present the mechanism and
models. In Section III we discuss the closely related phenomenology, focusing
on the ways to obtain correct relic density, stressing the role of
coannihilation. Section IV contains the conclusions and discussions. We add
the details of calculations for the fermionic DM annihilation cross section in
the appendix, which may involve some subtleties.

##  II The mechanism and the models

Since our mechanism is motivated by the semi-annihilating DM models, we will
start with the models with $Z_{3}$ and then generalize it to $Z_{N}$.

###  II.1 The $Z_{3}$ models

####  II.1.1 Model-A: Scalar DM with suppressed Higgs-portal by hand

Let us demonstrate the mechanism based on the well-known $Z_{3}$ symmetric
model, whose minimal form contains just a complex scalar $S_{1}$. But we will
see one more such a scalar $S_{2}$ makes a real difference. For concreteness,
both fields transform under $Z_{3}$ as $S_{i}\rightarrow
e^{ik_{i}2\pi/3}S_{i}$ with $k_{i}=1$ or $2$, then this model (model-A
hereafter) has the (incomplete) Lagrangian incorporating the following terms

$\displaystyle-{\cal L}_{Z_{3}}\supset$ $\displaystyle
m_{1}^{2}S_{1}S_{1}^{*}+m_{2}^{2}S_{2}S_{2}^{*}+\lambda_{1h}|S_{1}|^{2}|H|^{2}+\lambda_{2h}|S_{2}|^{2}|H|^{2}$
(1)
$\displaystyle+\left(\frac{A_{1}S_{1}^{3}}{3}+\frac{A_{2}S_{2}^{3}}{3}+\frac{1}{2}A_{12}S_{1}^{2}S_{2}+\frac{1}{2}A_{21}S_{1}S_{2}^{2}+c.c\right),$

In the present work, we impose the CP symmetry, and therefore all couplings
are real. Later, we will comment on the consequence of relaxing this
assumption. $S_{1}$ and $S_{2}$ receive mass terms both from the bare mass
term $m_{i}^{2}S_{i}S_{i}^{*}$ and the Higgs portal terms
$\lambda_{ih}|S_{i}|^{2}|H|^{2}$, so their squared masses are given by

$\displaystyle m_{S_{i}}^{2}=m_{i}^{2}+\lambda_{ih}v^{2},$ (2)

where we have written $H=v+h/\sqrt{2}$ with $v=174$GeV. $S_{1}$ is assumed to
be the lighter component thus the dark matter candidate, while $S_{2}$ can
decay into $S_{1}+h(\rightarrow\bar{f}f)$ after including a term
$\lambda_{12}S_{1}S_{2}^{*}|H|^{2}$, with the coupling strength $\lambda_{12}$
set to be sufficiently small, irrelevant to other dynamics 111As a matter of
fact, $S_{1}$ and $S_{2}$ are indistinguishable and therefore in general they
are expected to have a mixing via the term $m_{12}^{2}S_{1}S_{2}^{*}$ and
$\lambda_{12}S_{1}S_{2}^{*}|H|^{2}$; the mixing must be small, otherwise
$S_{1}$ will have a sizeable coupling to the Higgs doublet, which should be
suppressed in our mechanism. So, in writing the Lagrangian Eq. (II.1.1), we
are working in the mass basis, leaving the term $\lambda_{12}S_{1}S_{2}^{*}h$
for $S_{2}$ to decay..

The mechanism of hiding the WIMP dark matter candidate $S_{1}$ is simple. Let
us assume that $S_{1}$ is properly heavier than its companion $S_{2}$, namely
having mass within the window

$\displaystyle m_{S_{2}}+m_{h}<2m_{S_{1}}<2m_{S_{2}}.$ (3)

The window is wide, and it opens for the mass gap $\Delta=m_{S_{2}}-m_{S_{1}}$
satisfying $0<\Delta<m_{S_{1}}-m_{h}$; DM mass must be heavier than the Higgs
boson mass $m_{h}\approx 125$GeV. Eq. (3) ensures that the annihilation
channel $S_{1}S_{1}\rightarrow S_{2}^{*}h$ opens, with a cross section
proportional to $\propto A_{12}^{2}\lambda_{2h}^{2}$ rather than $\propto
A_{12}^{2}\lambda_{1h}^{2}$. Thus, we can turn off the Higgs portal coupling
of the dark matter candidate $S_{1}$ without turning off its annihilation
channels. In this way, we disconnect the interactions for DM annihilating from
scattering with quarks, to avoid the extremely tight constraints from the
latter. Note that the WIMP scenario is maintained, both $S_{1}$ and $S_{2}$
can be around the weak scale, with weak couplings strength.

As a matter of fact, $S_{2}$ plays the role of mediator between the dark
matter and the Higgs sector in the SM. But different than the usual models
where the DM purely annihilate into a pair of mediator (may or may be not the
Higgs boson pair), our model is characterized by the semi-mediator final
states, due to the fact that the mediator is not neutral under the DM
protecting symmetry.

####  II.1.2 Model-B: Fermionic DM without Higgs portal

One may be dissatisfied with the assumption that the Higgs portal coupling of
$S_{1}$ should be very small. Is it possible to forbid it in a natural and
simple way? To that end, we instead consider that $S_{1}$ is not a spin-0 but
a spin-1/2 fermion, denoted as $\Psi_{L}$. Its $Z_{3}$ charge is not changed.
As a consequence, many unwanted terms go away due to the Lorentz symmetry.
Owing to the $Z_{3}$ symmetry, $\Psi_{L}$ cannot be a Majorana fermion. So, to
make it massive, the chiral partner $\Psi_{R}$ is needed. Then, the general
model takes the form of

$\displaystyle-{\cal L}_{Z_{3}}\supset$
$\displaystyle+m_{S}^{2}|S|^{2}+M_{\Psi}\bar{\Psi}\Psi+\lambda_{sh}|S|^{2}|H|^{2}$
(4)
$\displaystyle+\left(\frac{A_{s}}{3}S^{3}+\lambda_{L}\overline{\Psi^{C}}P_{L}\Psi{S}+\lambda_{R}\overline{\Psi^{C}}P_{R}\Psi{S}+c.c.\right),$

where $\Psi=(\Psi_{L},\Psi_{R})$ is the Dirac fermion and the companion field
$S_{2}$ is renamed as $S$; the projector $P_{L,R}=\frac{1\mp\gamma^{5}}{2}$
and $\Psi^{C}=C\Psi^{T}$ with $C=i\gamma^{0}\gamma^{2}$ is the charge
conjugate field of $\Psi$. As expected, the direct connection between $\Psi$
and SM is absent, by virtue of the symmetry rather than by hand.

This model (model-B hereafter) respects an accidental parity, namely the
fermion number $\Psi\rightarrow-\Psi$. Therefore, both $S$ and $\psi$ are
stable without introducing other interactions. But we have no interest in
multi-component DM, and want to avoid such trouble. It is not a big problem,
and in the proper extensions to the above minimal model, this accidental
parity is supposed to be violated, admitting decay like
$S\rightarrow\Psi+\nu_{L}$. As proof of existence, let us consider the SM
extended by the type-I seesaw mechanism with right handed neutrinos $\nu_{R}$,
which are $Z_{3}$ neutral. Then couplings $S\bar{\Psi}_{L}\nu_{R}$ leads to
the decay $S\rightarrow\Psi_{L}+\nu_{R}$ if $\nu_{R}$ is sufficiently light,
otherwise $S\rightarrow\Psi_{L}+\nu_{L}$ via the mixing between $\nu_{R}$ and
$\nu_{L}$ 222It is of interest to consider a quite different realization of
the mechanism within the type-I seesaw extension to SM. Now, the term like
$\bar{\Psi}P_{R}\nu S$ provide a portal to DM and opens the annihilation
channel to a pair of sterile neutrinos. Such a way to avoid the direct
detection bound heavily relies on the presence of $\nu_{R}$, and moreover even
the $Z_{2}$ symmetry works. .

Another way to break the accidental symmetry is introducing an additional
Higgs doublet $\eta$ with quantum number like $H$ but carrying one unit of
$Z_{3}$ charge, and thus in general it mixes with $S$ due to the term
$\eta^{\dagger}HS$; further, we couple it to the SM left-handed lepton doublet
via $\eta\bar{\Psi}P_{L}L_{i}$, which then leads to the desired decays of $S$
and $\eta$ to $\Psi$. Although for different purposes, the Lagrangian Eq. (4),
along with those terms to break the accidental symmetry for $\Psi$, is a
subset of the model with radiative neutrino masses at two loop level Ma:2007gq
. The impact of semi-annihilation is studied for the first time in Ref.
Aoki:2014cja , which essentially covers our scenario of freeze-out for $\Psi$.
But their studies involve much more degrees of freedom and parameters, so, in
this sense, our study is a simplification of their analysis. Another obvious
difference is that, due to $\eta$ in the loop, the Dirac DM in their model
gains effective vertex $\bar{\Psi}\sigma^{\mu\nu}\Psi F_{\mu\nu}$, which
contributes to DM-nucleon scattering.

###  II.2 Generalization to $Z_{N}$ or beyond

From the previous model building, it is seen that the cubic term between DM
and its companion fields, such as $S_{1}^{2}S_{2}$ and
$\bar{\Psi}^{c}P_{L}\Psi{S}$, is the crucial structure for our mechanism to
work, while the portal term between the companion field and SM fields
like$|S_{2}|^{2}|H|^{2}$ is trivial viewing from symmetry, and it can be
replaced by other portals. Since those cubic terms are compatible with
$Z_{3}$, it explains why we start from this symmetry. But the cubic term is
allowed by many other symmetries, and this subsection is devoted to exploring
simple generalizations of $Z_{3}$.

It is illustrative to analyze the model from the pattern of symmetry, which is
helpful to generalize the mechanism from $Z_{3}$ to other possible symmetry
$G$. Here, the annihilation of DM is characterized by non-$G$ invariant: The
initiate state is the DM-DM pair, carrying a non-vanishing $G-$charge, and the
final states, containing the companion field, are not $G$ invariant neither,
carrying the same $G-$charge. Such a situation is in contrast to the familiar
case like $Z_{2}$, where the initial DM pair and the final states must be
$Z_{2}$ invariant. Although the minimal, $Z_{3}$ is definitely not the unique
symmetry to realize our scenario of DM annihilation. For simplicity, let us
focus on the case that $G$ is an Abelian symmetry $Z_{N}={e^{i2\pi
k/N}}(k=0,1,...N-1)$with $N>3$. In the following, we will derive the condition
for $Z_{N}$ to work. Suppose the DM and its companion field carry $Z_{N}$
charges $k_{1}$ and $k_{2}$, respectively. Therefore, charge conservation
yields the following equation

$\displaystyle 2k_{1}+k_{2}=0~{}{\rm module}~{}N,~{}~{}(k_{i}=1,...N-1).$ (5)

In the language of Lagrangian, one readily sees that the crucial cubic term
$S_{1}^{2}S_{2}$ along with the trivial Higgs portal term $|S_{2}|^{2}|H|^{2}$
are always present (Obviously, it is true also for the fermionic model.). For
$N>3$, the equation admits many solutions, for instance,

$\displaystyle N=4:~{}~{}k_{1}=1,~{}k_{2}=2;\quad k_{1}=3,~{}k_{2}=2,$ (6)
$\displaystyle N=5:~{}~{}k_{1}=1,~{}k_{2}=3;\quad k_{1}=2,~{}k_{2}=1;\quad
k_{1}=3,~{}k_{2}=4;\quad k_{1}=4,~{}k_{2}=2;$ (7) $\displaystyle
N=6:~{}~{}k_{1}=1,~{}k_{2}=4;\quad k_{1}=2,~{}k_{2}=2;\quad
k_{1}=4,~{}k_{2}=4;\quad k_{1}=5,~{}k_{2}=2.$ (8)

By the way, for $N=3$, the solutions are $k_{1}=k_{2}=1$ or $2$, thus in
either case the two fields carrying the same charge 333One can prove that only
the $Z_{N}$ which has the $Z_{3}$ subgroup admits such solutions: Let the
solution be $1\leq k_{1}=k_{2}=k=nN/3\leq N-1$, with $n=1,2$. So, to make $k$
be an integer, one must require $N=3m$. Then, in $Z_{3m}$, the two solutions
$k=m$ and $2m$ exactly correspond to the two solutions in $Z_{3}$..

But generically speaking, for $N>3$, these solutions require two fields with
different charges. As a consequence, they can be distinguished and cannot be
mixed with each other even if the spins are the same (both are scalar fields).
Then it is likely that both particles will be stable, leading to multi-
component dark matter; a related study is Ref. Yaguna:2019cvp . For
concreteness, let us show a model for $N=5$, with two scalars $S_{1}$ and
$S_{2}$ which respectively carry charges $k_{1}=1$ and $k_{2}=3$. Then, the
interactions of the most general (also typical for the larger $N$)
renormalizable Lagrangian is

$\displaystyle{\cal
L}_{Z^{5}}=\lambda_{1}|S_{1}|^{2}|H|^{2}+\lambda_{2}|S_{2}|^{2}|H|^{2}+\left(A_{12}S_{2}^{2}S_{1}+c.c.\right),$
(9)

$S_{2}$ can decay only via the term $S_{1}^{2}S_{2}$, however, it is
kinematically forbidden in the mass window Eq. (3) we are considering. In this
case, we must consider the rather strong coupling $\lambda_{2}$, so that the
annihilation cross section of the companion $S_{2}$ is large enough to make
its final abundance sufficiently small. Interestingly, this picture is nothing
but just the one mentioned in the introduction Saez:2021oxl ; Belanger:2011ww
.

Note for the higher $N$ case, some irrelevant terms such as $S_{i}^{3}$ may be
forbidden by symmetry, but other irrelevant terms like $S^{4}_{i}$ may arise.
Anyway, in addition to the crucial terms, the other terms appearing in the
Lagrangian depends on the concrete $N$ and as well charge assignments.

Any $Z_{N}$ discrete symmetry can be a residual symmetry of $U(1)$; see an
example Guo:2015lxa . Therefore, $G$ can be $U(1)$, either global or local. In
particular, if $U(1)$ is not spontaneously broken, then Eq. (5) should be
understood as the charge neutrality condition.

##  III Phenomenology of the DM-companion dark sector

In our mechanism, the main concern on DM is its relic density. In the present
study, we focus on the $Z_{3}$-like case, where the companion can decay. Short
comments on the possible probes to the model are made.

###  III.1 Freeze-out via semi-annihilation, including co-annihilation

####  III.1.1 Not standard semi-annihilation

In our models, the annihilation modes for DM are very limited, and in either
model, they are dominated by a pair of DM particles annihilating into the SM
Higgs boson plus the companion particle. Since the companion will eventually
contribute to the final DM relic density, essentially the DM freezes out by
means of semi-annihilation, which halves the efficiency of DM annihilation.

In the model-A, the Feymann diagrams for DM semi-annihilation are shown in
Fig. 3. There are two modes. We did not show the ordinary semi-annihilation
mode, $S_{1}S_{1}\rightarrow hS_{1}^{*}$, which is supposed to be negligible
because of the suppressed Higgs portal for DM. But in the numerical study we
will also include this contribution, to see if the state-of-the-art DM direct
detection still allows a pure Higgs portal DM scenario. The left panel shows
the dominant mode, $S_{1}S_{1}\rightarrow hS_{2}^{*}$ in the $s$-channel
exchange of the companion. Its cross section $\sigma_{s,1}$ times the DM
relative velocity is give by

$\displaystyle\sigma_{s,1}v_{r}=\frac{1}{64\pi}\frac{|\vec{p}_{h}|}{m_{S_{1}}^{3}}\frac{|A_{12}|^{2}(\lambda_{2h}v)^{2}}{(4m_{S_{1}}^{2}-m_{S_{2}}^{2})^{2}},$
(10)

where ${\vec{p}}_{h}$ is the momentum of the Higgs boson, in the center-of-
mass (CM) frame, taking the form of

$\displaystyle|\vec{p}_{h}|\approx
m_{S_{1}}\lambda(1,m_{S_{2}}^{2}/4m_{S_{1}}^{2},m_{h}^{2}/4m_{S_{1}}^{2}),$
(11)

with $\lambda(1,x,y)=\sqrt{(1-x-y)^{2}-4xy}$, which usually results in a
moderate phase space suppression. One can rewrite Eq. (10) as

$\displaystyle\sigma_{s,1}v_{r}=\frac{1}{64\pi}\frac{|A_{12}|^{2}(\lambda_{2h}v)^{2}}{16m_{S_{1}}^{6}}f(1,m_{S_{2}}^{2}/4m_{S_{1}}^{2},m_{h}^{2}/4m_{S_{1}}^{2}),$
(12)

where $f(x,y)=\lambda(1,x,y)/(1-x)^{2}$ is an order one number except for
$x\rightarrow 1$, hitting the resonant pole. Since we are considering CP
conservation in the dark sector, the CP conjugate process
$S^{*}_{1}S^{*}_{1}\rightarrow hS_{2}$ has an equal annihilation rate. At the
end of this article, we will briefly discuss the consequence if this is not
true.

Figure 1: Left: The Feymann diagram for DM semi-annihilation
$S_{1}S_{1}\rightarrow S_{2}^{*}h$ in the $s$-channel and DM-companion
coannihilation $S_{1}S_{2}^{*}\rightarrow S_{1}^{*}h$ in the $t$-channel;
Right: The diagram for companion annihilation into the SM memebers such as the
vector bosons via the Higgs-portal.

In the model-B, the DM is instead a Dirac fermion $\Psi$, its only
annihilation mode is $\Psi\Psi\rightarrow hS^{*}$, in the $s$-channel exchange
of the companion. In the appendix. A, we give the details of calculating the
annihilation cross section $\sigma_{f,1}$ times the DM relative velocity
$v_{r}$:

$\displaystyle\sigma_{f,1}v_{r}=\frac{1}{64\pi}\frac{(\lambda_{L}-\lambda_{R})^{2}(\lambda_{sh}v)^{2}}{16M_{\Psi}^{4}}f(1,m_{S}^{2}/4M_{\Psi}^{2},m_{h}^{2}/4M_{\Psi}^{2}).$
(13)

This cross section is almost identical to the one given in Eq. (10), except
that here the parameter $(\lambda_{L}-\lambda_{R})M_{\Psi}$ plays the role of
$A_{12}$. In addition, $\sigma_{f,1}v_{r}$ anishes as
$\lambda_{L}\rightarrow\lambda_{R}$, due to the interference effect. One can
understand this limit by combining the two Yukawa couplings in the Lagrangian
Eq. (4), to yield a Majorana-like Yukawa coupling
$\lambda\overline{\Psi^{C}}\Psi{S}+c.c.$; such a coupling is well-known to
give an annihilation cross section without $s$-wave.

####  III.1.2 Impacts of DM-companion co-annihilation

Although semi-annihilation has been widely studied, in our mechanism, it is
additionally characterized by the possibility of co-annihilation Griest:1990kh
due to the presence of companion. Despite of the more or less tuning of
parameters, the physical consequence of semi & co-annihilation DM is still
very significant. The co-annihilation modes involving the companion in the
initial state are inactive today, and therefore, if co-annihilation is the
dominant mechanism to determine DM relic density, the DM detectors, not only
direct detectors installed in the ground, but also indirect detectors flying
in the sky, turn a blind eye to DM.

In the following, we analyze the underlying reasons for the presence of co-
annihilation, and as well its impacts to the ordinary semi-annihilation. Our
example is again model-A, but the discussions completely apply to model-B.

  * •

First, the mass of the companion is supposed to be close to that of the DM, in
particular for DM mass not far above the Higgs boson mass, which is clearly
seen from the allowed window: $\Delta/m_{S_{1}}<1-m_{h}/m_{S_{1}}$. The
coannihilation effect becomes nonnegligible when the ratio
$\Delta/m_{S_{1}}\lesssim 1/x_{f}\sim 5\%-10\%$, with
$x_{f}=m_{S_{1}}/T_{f}\sim 20$ denoting the decoupling temperature of DM
freeze-out; the freeze-out process of two almost degenerate thermal particles
cannot be simply separated.

  * •

Next, the DM-companion co-annihilation mode is simply from the Feymann diagram
that leads to DM-DM annihilation; see the left panel of Fig. 3. The co-
annihilation $S_{1}S_{2}\rightarrow S^{*}_{1}h$ is just the $t$-channel
exchange of the companion. We denote its cross section as $\sigma_{s,2}$,
which is related to $\sigma_{s,1}$ by the cross symmetry, and therefor they
have almost equal area. This means that its impact is not very dramatic,
considering the exponential suppressing factor shown in Eq. (15).

  * •

Moreover, $\sigma_{3}$, the cross section of the companion-companion
annihilation $S_{2}S_{2}^{*}\rightarrow hh/VV/f\bar{f}$ with $V=Z$ or $W$ (see
the right panel of Fig. 3), can be enhanced by its relatively strong Higgs
portal coupling and the trilinear term coupling $A_{2}$ 444Such a strong Higgs
portal is beneficial to make the electroweak phase transition strong first
order, a sufficient condition for the electroweak baryogenesis mechanism. As a
matter of fact, in the $Z_{3}$ symmetric model, this possibility has been
studied by several groups Kang:2017mkl ; Kannike:2019mzk ; Chiang:2019oms . In
addition to that, if the companion is also the component of DM relic, which
probably is the situation in the models with $Z_{N}$ with $N>3$, the large
coupling is required to sufficiently reduce the fraction of companion in the
total DM relic density., and consequently it may even dominate the effective
cross section. Note that the companion Higgs portal thus $\sigma_{3}$ is
common to model-A and -B, so we do not need a further subscript to distinguish
the two models.

The precise relic density of DM is determined by the Boltzmann equations which
describe the evolution of DM number density. In the coannihilation scenario,
the relevant number density is the sum of DM and companion,
$n=n_{S_{1}}+n_{S_{2}}$, and its evolution resembles the ordinary one, with
the cross section replaced by the effective cross section,

$\displaystyle\frac{\mathrm{d}n}{\mathrm{d}t}$
$\displaystyle=-3Hn-\langle\sigma_{eff}v\rangle(n^{2}-n_{eq}^{2}),$ (14)
$\displaystyle\sigma_{eff}$
$\displaystyle\simeq\frac{1}{2}\sigma_{s,1}+\sigma_{s,2}(1+\delta)^{\frac{3}{2}}e^{-x_{f}\delta}+\sigma_{3}(1+\delta)e^{-2x_{f}\delta},$
(15)

where $H$ is the Hubble parameter. The cross sections $\sigma_{s,1/2}$ and
$\sigma_{3}$ are introduced previously, and in $\sigma_{eff}$ the
contributions from $\sigma_{s,1/2}$ have been halved, because each semi-
annihilation just reduces one dark particle (DM or the companion). The ratio
$\delta\equiv\Delta/m_{S_{1}}\ll 1$ measures the DM-companion mass degeneracy.
If their mass difference is large, the co-annihilation effect is not important
and then the Boltzmann equation above will reduce to the standard form.

###  III.2 Comments on probes to the dark sector

Although in the model-B the fermionic DM candidate $\Psi$ does not scatter
with the atom, in the model-A, the scalar DM candidate $S_{1}$ does, due to
the non-vanishing Higgs portal coupling $\lambda_{1h}|S_{1}|^{2}|H|^{2}$.
Then, the Higgs boson mediates $S_{1}$-nucleon spin-independent scattering,
with a cross section estimated to be Jungman:1995df

$\displaystyle\sigma^{(n)}_{\rm
S_{1}}=\frac{\lambda_{1h}^{2}}{4\pi}\frac{\mu_{n}^{2}m_{n}^{2}}{m_{S_{1}}^{2}m_{h}^{4}}\left(\sum_{q=u,d,s}f^{(n)}_{T_{q}}+\frac{2}{9}f^{(n)}_{T_{G}}\right)^{2},$
(16)

where $\mu_{n}\simeq m_{n}$ is the DM-nucleon reduced mass, and
$f_{T_{q}}^{(n)}$ the nucleon parameters. In the case that the companion is
also a tiny component of dark matter, the companion-nucleon scattering has a
cross section similar to the above expression. But the corresponding
contribution to the event rate should be multiplied by the fraction of this
component.

In general, our dark matter candidate leaves hints in the cosmic ray via the
following signature:

$\displaystyle{\rm DM}+{\rm DM}\rightarrow{\rm companion}(\rightarrow{\rm
DM}+h~{}~{}{\rm or}~{}~{}h^{*})+h,$ (17)

where the Higgs boson from the companion decay may be on-shell or off shell,
depending on the mass splitting $\Delta$. In either case, the resulting cosmic
ray spectra, such as of anti-proton and photon, need detailed simulations,
which is beyond the current project and will be present elsewhere. As a fairly
rough approximation, we may use the spectra from ${\rm DM}+{\rm DM}\rightarrow
h+h$, and there are available constraints in the literature Abazajian:2020tww
. Later, we will show this schematic constraint. Note that the indirect
detection signal of DM may be greatly weakened in the co-annihilation
mentioned region.

If the companion has a relatively strong interaction with the Higgs sector,
the dark sector may leave imprints at the collider such as LHC. As noted in
footnote 4, such a companion can play the role of trigger for strongly first
order of electroweak phase transition. Then, one can probe its effects on the
Higgs physics by measuring the Higgs self-coupling, which has been widely
discussed. But this approach can not distinguish different bosonic triggers.
To identify the companion, it is better to produce it directly at the
collider. They can be pairly produced via gluon gluon fusion at the LHC,
giving rise to signals of missing energy plus a pair of (on shell or off
shell) Higgs bosons. Note that similar signals are present in supersymmetry,
where the next-to-lightest sparticle (LSP) may decay into the LSP plus Higgs
boson Martin:1997ns . However, if the companion itself behaves as missing
energy at the detector, one should turn to different strategies.

###  III.3 Parameter Scan and Result

We display the numerical results of the two benchmark models, using the
program micrOMEGAs Belanger:2020gnr . The model-A given in Eq. (II.1.1)
contains 8 parameters, but only four of them are relevant, which are chosen to
be: 1) masses of DM and the companion,
$m_{i}=\sqrt{M_{i}^{2}+\lambda_{ih}v^{2}}$; 2) the companion-Higgs portal
coupling $\lambda_{2h}\sim{\cal O}(1)$, whereas the DM-Higgs portal coupling
is suppressed $\lambda_{1h}\ll 1$; 3) the DM-companion coupling $A_{12}$ and
the trilinear term coupling for the companion, $A_{2}$, which is relevant in
the coannihilation scenario to enhance the annihilation cross section for the
companion. Similar analysis on the parameter space apply to model-B, and one
can find that it also involves four key parameters. To learn how well the
parameter space of our model fits the experimental observations, we use
Metropolis-Hastings algorithm mickay  to scan the parameter space. In Table.
III.3 and Table. III.3 we show the range and the random walk size for each of
parameters.

model-A
---
parameters | range | step size
$A_{1}$ | $10\rm\ GeV$ | $0\rm\ GeV$
$A_{2}$ | $[0,1000]\rm\ GeV$ | $200\rm\ GeV$
$A_{12}$ | $[0,1000]\ \rm GeV$ | $200\rm\ GeV$
$A_{21}$ | $10\rm\ GeV$ | $0\rm\ GeV$
$\log_{10}\lambda_{1h}$ | $[-3,0]$ | $1$
$\lambda_{2h}$ | $[0.1,2.5]$ | $0.6$
$m_{S_{1}}$ | $[m_{h},2000]\rm\ GeV$ | $400\rm\ GeV$
$\log_{10}\delta$ | $[-2,0]$ | $1$
Table 1: Information of the scanned parameters with
$\delta=\frac{m_{S_{2}}}{m_{S_{1}}}-1$ model-B
---
parameters | range | step size
$A_{S}$ | $[0,1000]\rm\ GeV$ | $200\rm\ GeV$
$\lambda_{sh}$ | $[0.1,2.5]$ | $0.6$
$\lambda_{L}$ | $[0.1,2.5]$ | $0.6$
$\lambda_{R}$ | $0.1$ | $0$
$M_{\Psi}$ | $[m_{h},2000]\rm\ GeV$ | $400\rm\ GeV$
$\log_{10}\delta$ | $[-2,0]$ | $1$
|  |
|  |
Table 2: Information of the scanned parameters with
$\delta=\frac{m_{S}}{m_{\Psi}}-1$

The experimental observations we consider are relic density and two most
stringent DM direct detection experimental results, XENON. Aprile:2018dbl  and
PandaX. Meng:2021mui . We set the observed relic density in log scale
$\log_{10}\widehat{\Omega h^{2}}=-1$; to get a larger allowed parameter space,
the uncertainty of $\log_{10}\Omega h^{2}$ is set to be $\sigma=0.5$. Then,
the posterior log-likelihood for a certain parameter point could be written
as:

$\displaystyle\ln
P(Obs|x)=\ln\frac{1}{\sqrt{2\pi}\sigma}\exp\left(\frac{(\log_{10}\widehat{\Omega
h^{2}}-\log_{10}\Omega h^{2})^{2}}{2\sigma^{2}}\right)+\ln(p_{val}),$ (18)

where $p_{val}$ is the probability of acceptance by all recent DM direct
detection experiments, and it is calculated by micrOMEGAs Belanger:2020gnr
using the maxsimum gap method gap method . After applying MCMC process for
$1\times 10^{6}$ times, we get the posterior distribution of parameter space,
and we select the last $1\times 10^{5}$ points for analysis.

Let us end up the main text with some figures, to show the preferred parameter
space and its features. For the scanning parameters given in Table. III.3 and
Table. III.3, the correct relic density can be achieved in a wide parameter
space, where the companion-Higgs portal coupling and the DM-companion coupling
should be sufficiently large; see the left panels in Fig. 2 and Fig. 4. The
viable DM mass mostly lies above about 200 GeV, and can extend above the TeV
region. However, for the upper bound of parameters that we impose, assistant
from coannihilation is required, as is explicit in the example of model-B
(similar in model-A although not present) shown in the right panel of Fig. 4:
For $M_{\Psi}\gtrsim 1$TeV, the required degeneracy $\delta\lesssim 4\%$.
Therefore, the natural parameter space predicts a sub-TeV DM, if we do not
want so large couplings that the theory hits the Landau pole at a quite low
scale.

Figure 2: Correct DM relic density in model-A. Left: the viable parameter
space on the $m_{S_{1}}-\lambda_{2h}$ plane, with colors encoding the value of
DM-companion coupling $A_{12}$. Right: for the low value of $\lambda_{1h}$,
our model escape the strong constrain of Xenon 1T

To show the indirect detection prospects, in the left panels of Fig. 3 and
Fig. 5 we plot the annihilation cross sections of ${\rm DM+DM}\rightarrow{\rm
companion}+h$ time DM velocity today, whose horizontal axis is the degeneracy
parameter $\delta$, to highlight the coannihilation effect. As expected, this
effect has a significant impact on the indirect detection signals, and from
these panels it is seen that $\langle\sigma v\rangle$ may obviously lie below
the typical thermal cross section$\sim 1$pb. As explained before, we
approximate the spectra of DM annihilation as that of the DM directly
annihilating into a pair of Higgs boson. Then, in the left panels of Fig. 3
and Fig. 5, we schematically constrain the parameter space. Of interest, the
available data is already able to rule out the lighter DM space without
coannihilation effect. This situation makes further detailed investigation on
the indirect signals of the model necessary, which needs other tools to
simulate the spectra.

By the way, in the numerical study on model-A, we did not limit to a very
small $\lambda_{1h}$, so we show the tight bounds from the DM direct detection
in the right panel of Fig. 2. Even for $m_{S_{1}}$ in the TeV region, where
the bound is relatively loose, still the allowed $\lambda_{1h}\lesssim 0.5$ .
As a consequence, the latest PandaX-4T data Meng:2021mui  imposing a bound up
to 10 TeV (compared to the XENON1T 2018, its exclusion ability is improved by
an order of magnitude), rules out the conventional Higgs-portal scenario,
which otherwise survived in the high mass region.

Figure 3: Left: In model-A, plot of DM-DM annihilation cross section times DM
relative velocity today, depending on the degeneracy $\delta$; the color bar
is for the companion-Higgs portal coupling. Right: The schematic constraints
on DM via indirect detection, approximating the DM annihilating spectra as
DM+DM$\rightarrow hh$.

Figure 4: Correct DM relic density in model-B. Left: the viable parameter
space on the $M_{\Psi}-\lambda_{sh}$ plane, with colors encoding the value of
$\lambda_{L}$. Right: Plots on the $M_{\Psi}-\delta$ plane with colors
encoding the value of $\lambda_{sh}$.

Figure 5: Left: In model-B, plot of DM-DM annihilation cross section times DM
relative velocity today, depending on the degeneracy $\delta$; the color bar
is for the companion-Higgs portal coupling. Right: The schematic constraints
on DM via indirect detection, approximating the DM annihilating spectra as
DM+DM$\rightarrow hh$.

##  IV Conclusions and discussions

Because the correct relic abundance can be obtained naturally, WIMP has been
expected to be the leading candidate for dark matter. Nevertheless, facing the
more and more stringent constraints from many DM direct detection experiments,
the typical WIMP dark matter candidate becomes less and less possible. Hence,
in order to save the WIMP DM scenario, we should consider a more complicated
structure of the WIMP dark sector.

In this article we propose a novel approach to overcome the strong DM direct
detection bound on the WIMP DM, inspired by the semi-annihilating DM models.
The crucial assumption is that the DM has a companion which is also charged
under the DM protecting symmetry, and moreover the symmetry allows the
companion (usually a spin-0 particle) to play the role of bridge between DM
and SM, via the trilinear coupling DM-DM-companion. Then, DM semi-annihilates
into the companion to reduce its relic density, having no DM-neulceon
scattering signals. Such an idea can be easily realized in the $Z_{N}$
symmetric models with $N>2$. If DM is a scalar field, one may have to suppress
its Higgs portal coupling by hand, whereas the fermionic DM naturally has
vanishing coupling to the Higgs. We stress that this mechanism has the
characteristics of co-annihilation, and as a matter of fact its effect becomes
necessary near or above the TeV region. This means that maybe it is difficult
to detect our dark matter directly and indirectly.

There are some open and pending issues. First, in the scenario of semi-
annihilation, DM is a complex field, consisting of DM and anti-DM. However, DM
freeze-out is characterized by DM-DM instead of DM and anit-DM annihilation,
which is traced back to DM number violation from the cubic terms. Therefore,
the model provides all three elements to generate an asymmetry between DM and
anti-DM if CP conservation is abandoned. Such a very interesting possibility
deserves a further study elsewhere. Next, the companion may be also a DM
component, which brings a remarkable difference to DM phenomenologies. Third,
from the side of model building, it is of interest to explore if the idea can
be realized in the case that DM protecting symmetry is a non-Abelian group.
Last but not least, looking for imprints of our model in the cosmic ray and
accelerator such as LHC are also on the agenda. We commented both aspects in
the main text, and moreover used a direct di-Higgs boson to approximate the
cosmic ray spectra from DM annihilation,. However, a lot of work needs to be
done to obtain reliable spectral lines

##  Appendix A DM-DM semi-annihilation cross section in model-B

In this appendix we give the details of the cross section of DM annihilation
in model-B with a fermionic DM. The Lagrangian is given by Eq. (4). We list
below the relevant interactions and the corresponding Feynman rules:

$\displaystyle\lambda_{L}\bar{\Psi}^{c}P_{L}\Psi{S}:\lambda_{L}(1-\gamma 5),$
(19) $\displaystyle\lambda_{R}\bar{\Psi}^{c}P_{R}\Psi{S}:\lambda_{R}(1+\gamma
5),$ (20)
$\displaystyle\lambda_{L}^{*}\bar{\Psi}P_{L}\Psi^{c}{S}:\lambda_{L}^{*}(1+\gamma
5),$ (21)
$\displaystyle\lambda_{R}^{*}\bar{\Psi}P_{R}\Psi^{c}{S}:\lambda_{R}^{*}(1-\gamma
5).$ (22)

With them we can calculate the usual spin averaged amplitude squared for DM-DM
annihilation:

$\displaystyle\frac{1}{4}{\sum_{spins}|M|^{2}}$
$\displaystyle=2\frac{\lambda^{2}_{sh}(\lambda_{L}^{2}+\lambda_{R}^{2}-2\lambda_{L}\lambda_{R})v^{2}m^{2}_{\chi}}{(q^{2}-m^{2}_{s})^{2}},$
(23)

and then the cross section times the DM relative velocity is

$\displaystyle\sigma_{0}v$
$\displaystyle=\frac{|\vec{p}_{h}|(\lambda_{L}^{2}+\lambda_{R}^{2}-2\lambda_{L}\lambda_{R})\lambda_{sh}^{2}v^{2}m_{\chi}^{2}}{2\pi\sqrt{s}^{3}(s-m_{s}^{2})^{2}},$
(24)

where $s$ is the square of the energy of the center of mass. Noe that the
above expression includes both (putatively equal) cross sections for DM-DM and
${\rm DM}^{*}-{\rm DM}^{*}$ annihilation.

But for the fermonic DM, to obtain the cross section that appears in the
formula for calculating relic density,

$\displaystyle\Omega
h^{2}=1.69\times\frac{x_{f}}{20}\sqrt{100/g_{*}}\left(\frac{10^{-10}\rm
GeV^{-2}}{\langle\sigma v\rangle}\right)$ (25)

with $g_{*}$ the total degree of freedom, we should handle the spin degrees of
freedom properly. In the density expression of the initial state, a Dirac
fermion has 4 degrees of freedom, including 2 spin states for each of the
fermion $\psi$ and an antifermion $\bar{\psi}$. According to the equipartition
theorem, we know $n_{i}=\frac{n}{4}$ for every spin type particle, with $n$
the total particle number density. Moreover, the number density of each spin
degree evolves as
$\frac{\mathrm{d}n_{i}}{\mathrm{d}t}=-\sum_{j=1}^{4}\langle\sigma_{ij}v\rangle
n_{i}n_{j}$. Then, summing over $i$, we get the following equation
Srednicki:1988ce ,

$\displaystyle\frac{\mathrm{d}n}{\mathrm{d}t}$
$\displaystyle=-\frac{1}{16}\sum_{i,j=1}^{4}\langle\sigma_{ij}v\rangle
n^{2}=-\langle\sigma v\rangle n^{2},$ (26)

where we have defined the cross section

$\displaystyle\sigma\equiv\frac{1}{16}\sum_{i,j=1}^{4}\sigma_{ij}=\frac{1}{16}\sum_{spins}\left(\sigma_{\bar{\Psi}\Psi}+\sigma_{{\Psi}\bar{\Psi}}+\sigma_{{\Psi}\Psi}+\sigma_{\bar{\Psi}\bar{\Psi}}\right).$
(27)

Specific to our model, there are only two types of reactions in our model,
namely $\sigma_{\bar{\Psi}\Psi}=\sigma_{\Psi\bar{\Psi}}=0$, and moreover we
assumed $\sigma_{{\Psi}\Psi}=\sigma_{\bar{\Psi}\bar{\Psi}}$. Therefore, we can
establish the relationship between $\sigma$ and $\sigma_{0}$ as follows:

$\displaystyle\sigma$
$\displaystyle=\frac{1}{16}\sum_{spins}(\sigma_{\bar{\Psi}\bar{\Psi}}+\sigma_{\Psi\Psi})=\frac{1}{4}\times\frac{1}{4}\sum_{spins}(\sigma_{\bar{\Psi}\bar{\Psi}}+\sigma_{\Psi\Psi})=\frac{1}{4}\sigma_{0}.$
(28)

Eventually, we obtain the cross section adopt in the relic density formula Eq.
(25),

$\displaystyle\sigma v$
$\displaystyle=\frac{1}{4}\sigma_{0}v=\frac{|\vec{p}_{h}|(\lambda_{L}^{2}+\lambda_{R}^{2}-2\lambda_{L}\lambda_{R})\lambda_{sh}^{2}v^{2}m_{\chi}^{2}}{8\pi\sqrt{s}^{3}(s-m_{s}^{2})^{2}}.$
(29)

Acknowledgements

This work is supported in part by the National Science Foundation of China
(11775086).

## References

  * (1) G. Bertone, D. Hooper and J. Silk, Phys. Rept. 405, 279-390 (2005).
  * (2) N. Aghanim et al. [Planck], Astron. Astrophys. 641, A6 (2020) [erratum: Astron. Astrophys. 652, C4 (2021)].
  * (3) E. Aprile et al. [XENON], Phys. Rev. Lett. 121, no.11, 111302 (2018).
  * (4) Y. Meng, Z. Wang, Y. Tao, A. Abdukerim, Z. Bo, W. Chen, X. Chen, Y. Chen, C. Cheng and Y. Cheng, et al. [arXiv:2107.13438 [hep-ex]].
  * (5) J. McDonald, Phys. Rev. Lett. 88, 091304 (2002) [hep-ph/0106249].
  * (6) L. J. Hall, K. Jedamzik, J. March-Russell and S. M. West, JHEP 1003, 080 (2010).
  * (7) N. Bernal, M. Heikinheimo, T. Tenkanen, K. Tuominen and V. Vaskonen, Int. J. Mod. Phys. A 32, no.27, 1730023 (2017).
  * (8) F. D’Eramo and J. Thaler, JHEP 06, 109 (2010).
  * (9) C. Gross, O. Lebedev and T. Toma, Phys. Rev. Lett. 119, no.19, 191801 (2017).
  * (10) N. Okada, D. Raut and Q. Shafi, Phys. Rev. D 103, no.5, 055024 (2021).
  * (11) X. M. Jiang, C. Cai, Z. H. Yu, Y. P. Zeng and H. H. Zhang, Phys. Rev. D 100, no.7, 075011 (2019).
  * (12) G. Belanger, K. Kannike, A. Pukhov and M. Raidal, JCAP 04, 010 (2012).
  * (13) G. Bélanger, K. Kannike, A. Pukhov and M. Raidal, JCAP 06, 021 (2014).
  * (14) P. Ko and Y. Tang, JCAP 05, 047 (2014).
  * (15) S. L. Chen and Z. Kang, JCAP 05, 036 (2018).
  * (16) Z. Kang, Eur. Phys. J. C 75, no.10, 471 (2015).
  * (17) J. L. Feng, J. Kumar, D. Marfatia and D. Sanford, Phys. Lett. B 703, 124-127 (2011).
  * (18) Z. Kang, T. Li, T. Liu, C. Tong and J. M. Yang, JCAP 01, 028 (2011).
  * (19) X. Gao, Z. Kang and T. Li, JCAP 01, 021 (2013).
  * (20) B. Díaz Sáez, K. Möhling and D. Stöckinger, [arXiv:2103.17064 [hep-ph]].
  * (21) G. Belanger and J. C. Park, JCAP 03, 038 (2012).
  * (22) C. E. Yaguna and Ó. Zapata, JHEP 03, 109 (2020).
  * (23) E. Ma, Phys. Lett. B 662, 49-52 (2008).
  * (24) M. Aoki and T. Toma, JCAP 09, 016 (2014).
  * (25) J. Guo, Z. Kang, P. Ko and Y. Orikasa, Phys. Rev. D 91, no.11, 115017 (2015).
  * (26) K. Griest and D. Seckel, Phys. Rev. D 43, 3191-3203 (1991) doi:10.1103/PhysRevD.43.3191
  * (27) Z. Kang, P. Ko and T. Matsui, JHEP 02, 115 (2018).
  * (28) K. Kannike, K. Loos and M. Raidal, Phys. Rev. D 101, no.3, 035001 (2020).
  * (29) C. W. Chiang and B. Q. Lu, JHEP 07, 082 (2020).
  * (30) G. Jungman, M. Kamionkowski and K. Griest, Phys. Rept. 267, 195-373 (1996).
  * (31) David J. C. MacKay, Information Theory, Inference, and Learning Algorithms
  * (32) G. Belanger, A. Mjallal and A. Pukhov, Eur. Phys. J. C 81, no.3, 239 (2021)
  * (33) S.Yellin, Phys. Rev. D66, 032005 (2002).
  * (34) K. N. Abazajian, S. Horiuchi, M. Kaplinghat, R. E. Keeley and O. Macias, Phys. Rev. D 102 (2020) no.4, 043012
  * (35) S. P. Martin, Adv. Ser. Direct. High Energy Phys. 18, 1-98 (1998).
  * (36) M. Srednicki, R. Watkins and K. A. Olive, Nucl. Phys. B 310, 693 (1988).