11institutetext: Jerzy Proficz 22institutetext: Gdansk University of
Technology
Centre of Informatics - Tricity Academic Supercomputer & networK (CI TASK)
11/12 Gabriela Narutowicza Street, 80-233 Gdańsk, Poland
22email<EMAIL_ADDRESS>Tel.: +48 58 3486343

# Improving all-reduce collective operations for imbalanced process arrival
patterns

Jerzy Proficz

###### Abstract

Two new algorithms for the all-reduce operation, optimized for imbalanced
process arrival patterns (PAPs) are presented: (i) sorted linear tree (SLT),
(ii) pre-reduced ring (PRR) as well as a new way of on-line PAP detection,
including process arrival time (PAT) estimations and their distribution
between cooperating processes was introduced. The idea, pseudo-code,
implementation details, benchmark for performance evaluation and a real case
example for machine learning are provided. The results of the experiments were
described and analyzed, showing that the proposed solution has high
scalability and improved performance in comparison with the usually used ring
and Rabenseifner algorithms.

###### Keywords:

All-reduce, Pre-reduced ring, Sorted linear tree, Process arrival pattern, MPI

##  1 Introduction

Collective communication MPICollective  is frequently used by the programmers
and designers of parallel programs, especially in high performance computing
(HPC) applications related to scientific simulations and data analysis,
including machine learning calculations. Usually, collective operations, e.g.
implemented in MPI MPI , are based on algorithms optimized for the
simultaneous entering of all participants into the operation, i.e. they do not
take into consideration possible differences in process arrival times (PATs),
thus, in real environment, where such imbalances are ubiquitous, they can have
significant performance issues. It is worth to note that well performing
algorithms for the balanced times, work poorly in the opposite case Faraj2008
.

Fig. 1 presents an example of a typical execution of a distributed program,
where after the computation phase all processes exchange data with each other
using some kind of collective communication operation, e.g. all-reduce. We can
observe that even for the same computation volume, different processes arrive
at the communication phase in different time, in the example processes $1$ is
the slowest. Sometimes such differences can be observed in the communication,
where the data exchange time is shorter for process 1 than for the other
processes.

Figure 1: Example of a process arrival pattern with the process #1 delayed,
where $P$ is the number of cooperating processes, $a_{i}$, $f_{i}$ and $e_{i}$
are respectively arrival, finish and elapsed times of process $i$ for the
performed collective communication operation

As a contribution of this paper, we present two new algorithms for the all-
reduce operation, optimized for imbalanced process arrival patterns (PAPs):
sorted linear tree (SLT) and pre-reduced ring (PRR). We described their idea,
pseudo-code, implementation details, benchmark for their evaluation as well as
a real case example related to machine learning. Additionally we introduced a
new way of on-line PAP detection, including PAT estimations and their
distribution among cooperating processes.

The following section presents the related works in the subject, the next one
describes the used computation and communication model, section 4 presents the
proposed algorithms, section 5 provides the evaluation of the algorithms using
a benchmark, section 6 shows a real case example of the algorithms’
utilization, and the last section presents the final remarks.

##  2 Related works

We grouped the related works into three areas: the all-reduce operation in
general, i.e. the review of the currently used algorithms in different
implementations of MPI MPI , then we describe the current state of the art in
process arrival patterns (PAPs), and finally we present the works related to
process arrival times (PATs) on-line monitoring and estimation.

###  2.1 All-reduce operation

All-reduce operation is one of the most common collective operations used in
HPC software Faraj2008 . We can define it as a reduction in a vector of
numbers using a defined operation, e.g. sum, which needs to be cumulative, and
then distribution of the result into all participating processes, in short:
all-reduce = reduce + broadcast. There are plenty of all-reduce algorithms,
Table 1 summarizes the ones used in two currently most popular, open source
MPI implementations: OpenMPI OpenMPI  and MPICH MPICH .

Table 1: All-reduce algorithms implemented in OpenMPI OpenMPI  and MPICH MPICH  Selection criteria |  OpenMPI |  MPICH
---|---|---
Short messages |  Recursive doubling Thakur2005  |  Recursive doubling Thakur2005
Long messages |  Ring Thakur2005  or Segmented ring OpenMPI  |  Rabenseifner Rabenseifner2004
Non-commutive reduce operation |  Nonoverlapping OpenMPI  |  Recursive doubling Thakur2005
Others |  Basic linear OpenMPI , Segmented linear tree OpenMPI  |

The _basic linear_ algorithm performs linear reduce (flat tree: the root
process gathers and reduces all data from the other processes) followed by the
broadcast without any message segmentation OpenMPI . _Segmented linear tree_
creates pipeline between the participating processes, where the data are split
into segments and sent from process $0$ to $1$ to $2$ …to $P-1$. The
_nonoverlapping_ algorithm uses default reduce followed by the broadcast,
these operations are not overlapping: the broadcast is performed sequentially
after the reduce, even if both use segmentation and some segments are ready
after the reduce OpenMPI . _Recursive doubling_ (a.k.a. _butterfly_ and
_binary split_) is performed in rounds, in every round each process exchanges
and reduces its data with another, corresponding process, whose rank changes
round by round, after $log_{2}P$ rounds ($P$ is the total number of
participating processes) all data are reduced and distributed to the processes
Thakur2005 . _Ring_ is also performed in rounds, the data are divided into $P$
segments and each process sends one segment per round (the segment index
depending on the round number and the process rank) to the next process (the
processes are usually ordered according to their ranks). At the beginning,
there are performed $P-1$ rounds, where a segment delivery is followed by its
data reduction, thus after these rounds, each process has exactly one segment
fully reduced, and needs to forward it to other processes. This is performed
during next $P-1$ rounds, where the data are gathered, i.e. the processes
transfer the reduced segments, thus the final result is delivered to all
processes Thakur2005 . Fig. 2 presents an example of a ring execution.
_Segmented ring_ is similar to the ring, but after the segmentation related to
process number, an additional one is performed for pipeline effect between
corresponding processes OpenMPI . In the _Rabenseifner_ algorithm the
reduction is performed in two phases, first the scatter and reduce of data are
executed (using recursive doubling on divided data) and then the data
gathering takes place (using recursive halving on gathered data), where the
reduced data come back to all processes Rabenseifner2004 .

[width=8cm]Fig2

Figure 2: Example of a ring algorithm execution, showing data distribution
between the processes (one square corresponds to one data segment), a solid
arrow represents send and reduce a data segment, an empty arrow represents
send and override a data segment, and $\tau$ is the time of transfer and
reduction/overriding of a data segment

###  2.2 Imbalanced process arrival times

Not much work has been done for imbalanced process arrival patterns analysis.
To the best of our knowledge, the following four papers cover the research
performed in this area.

In Faraj2008  Faraj et al. performed advanced analysis of _process arrival
patterns_ (PAPs) observed during execution of typical MPI programs. They
executed a set of tests proving that the differences between _process arrival
times_ (PATs) in operations of the collective communication are significantly
high and they influence the performance of the underlying computations. The
authors defined a PAP to be imbalanced for a given collective operation with a
specific message length when its imbalance factor (a ratio between the highest
difference between the arrival times of the processes and time of the simple
(point-to-point) message delivery between each other) is larger than 1. The
authors provided examples of typical HPC benchmarks, e.g. NAS, LAMMPS or
NBODY, where imbalance factor, during their execution in a typical cluster
environment, equals 100 or even more. They observed that, such behavior
usually cannot be controlled directly by the programmer, and the imbalances
are going to occur in any typical HPC environment. The authors proposed a
mini-benchmark for testing various collective operations and found out the
conclusion that the algorithms which perform better with balanced PAPs tend to
behave worse when dealing with imbalanced ones. Finally, they proposed
solution: their self-tuning framework – STAR-MPI Faraj2006 , which includes a
large set of various implementations of collective operation and can be used
for different PAPs, with automatic detection of the most suitable algorithm.
The framework efficiency was proved by an example of tuned all-to-all
operations, where the performance of the set of MPI benchmarks was
significantly increased.

As a continuation of the above work Patarasuk et al. proposed a new solution
for broadcast operation used in MPI application concerning imbalanced PAPs of
the cooperating processes Patarasuk2008 . The authors proposed a new metric
for the algorithm performance: competitive ratio – $PERF(r)$, which describes
the influence of the imbalanced PATs for the algorithm execution, regarding
the behavior for its worst case PAT. They evaluated well known broadcast
algorithms, using the above metric, and presented two new algorithms, which
have constant (limited) value of the metric. The algorithms are meant for
large messages, and use the subsets of cooperating processes to accelerate the
overall process: the data are sent to the earliest processes first. One of the
algorithms is dedicated for non-blocking (_arrival_nb_) and other for blocking
(_arrival_b_) message passing systems. The authors proposed a benchmark for
algorithms’ evaluation, which introduces random PAPs and measures their impact
on the algorithm performance. The experiments were performed using two
different 16-node compute clusters (one with Infiniband and other with
Ethernet interconnecting network), and 5 broadcast algorithms, i.e.
_arrival_b_ , flat, linear, binomial trees and the one native to the machine.
The results of the experiments showed the advantage of the _arrival_b_
algorithm for large messages and imbalanced PAPs.

In Marendic2012  Marendic et al. focused on an analysis of reduce algorithms
working with imbalanced PATs. They assumed atomicity of reduced data (the data
cannot be split into segments and reduced piece by piece), as well as the
Hockney Hockney1994  model of message passing (time of message transmission
depends on the link bandwidth: $\beta$ and constant latency: $\alpha$, with an
additional computation speed parameter: $\gamma$) and presented related works
for typical reduction algorithms. They proposed a new static load balancing
optimized reduction algorithm requiring a priori information about current
PATs of all cooperating processes. The authors performed a theoretical
analysis proving the algorithm is nearly optimal for the assumed model. They
showed that the algorithm gives the minimal completion time under the
assumption that the corresponding point-to-point operations start exactly at
the same time for any two communicating processes. However if the model
introduces a delay of the receive operation in comparison with the send one,
which seems to be the case in real systems, the algorithm does not utilize
this additional time in receiving process, although, in some cases, it could
slightly improve the performance of the overall reduce operation. The other
proposed algorithm, presented by the authors: a dynamic load balancing can
operate under the limited knowledge about PATs, being able to atomically
reconfigure the message passing tree structure while performing reduce
operation using auxiliary short messages for signaling the PATs between the
cooperating processes. The overhead is minimal in comparison with the gains of
the PAP optimization. Finally a mini-benchmark was presented and some typical
PAPs were examined, the results showed the advantage of the proposed dynamic
load balancing algorithm versus other algorithms: binary tree and all-to-all
reduce.

In Marendic2016  Marendic et al. continued the work with optimization of the
MPI reduction operations dealing with the imbalanced PAPs. The main
contribution is a new algorithm, called Clairvoyant, scheduling the exchange
of data segments (fixed parts of reduced data) between reducing processes,
without assumption of data atomicity, and taking into account PATs, thus
causing as many as possible segments to be reduced by the early arriving
processes. The idea of the algorithm bases on the assumption that the PAP is
known during process scheduling. The paper provided a theoretical background
for the PAPs, with its own definition of the time imbalances, including a PAT
vector, absolute imbalance, absorption time as well as their normalized
versions, followed by the analysis of the proposed algorithm, and its
comparison to other typically used reduction algorithms. Its pseudo-code was
described and the implementation details were roughly provided with two
examples of its execution for balanced and imbalanced PAPs. Afterwards the
performed experiments were described, including details about used mini-
benchmark and the results of practical comparison with other solutions
(typical algorithms with no support for imbalanced PAPs) were provided.
Finally, the results of the experiments showing advantage of the proposed
algorithm were presented and discussed.

###  2.3 Process arrival time estimation

The PAP collective communication algorithms require some knowledge about the
PATs for their execution. Table 2 presents the summary of the approaches used
in the works described above.

Table 2: Approaches for PAP detection Operation |  PAP detector
---|---
All-to-all Faraj2008  |  Uses STAR-MPI overall efficiency indicators
Broadcast Patarasuk2008  |  Uses its own messaging for PAT signaling
Local redirect Marendic2012  |  Uses its own messaging for PAT signaling
Clairvoyant Marendic2016  |  None, suggested usage of the static analysis or SMA

In Faraj2008  (dedicated for an all-to-all collective operation), there is
assumption about the call site (a place in the code where the MPI collective
operation is called) paired with the message size that they have a similar
PAPs for the whole program execution, or at least their behavior changes
infrequently. The proposed STAR-MPI Faraj2006  system periodically assesses
the call site performance (exchanging measured times between processes) and
adapts a proper algorithm, trying one after another. The authors claim that it
requires 50–200 calls to provide desired optimization. This is a general
approach and it can be used even for other performance issues, e.g. network
structure adaptation.

In Patarasuk2008  (dedicated for a broadcast operation), the algorithm uses
additional, short messages sent to root process signaling process readiness
for the operations. In case the some processes are ready, the root performs
sub-group broadcast, thus the a priori PATs are not necessary for this
approach. Similar idea is used in Marendic2012  (dedicated for a reduce
operation) where the additional messages are used not only to indicate
readiness, but also to redirect delayed processes.

In Marendic2016  (dedicated for a reduce operation), the algorithm itself does
not include any solution for the PAT estimation, that is why it is called
Clairvoyant, but the authors assume recurring PAPs and give the suggestion
that there can be used simple moving averages (SMA) approximation. This
solution requires the additional communication to exchange the PAT values,
what is performed every $k$ iterations, thus introducing the additional
communication time. The authors claim that the speedup introduced by the usage
of the algorithm overcomes this cost, and provide some experimental results
showing the total time reduction in the overall computations.

##  3 Computation and communication model

We assume usage of the message-passing paradigm, within a homogeneous
environment, where each communicating process is placed on a separated compute
node. The nodes are connected by a homogeneous communication network. Every
process can handle one or more threads of control communicating and
synchronizing with each other using shared memory mechanisms. However there is
no shared memory accessible simultaneously by different processes.

As a _process arrival pattern_ (PAP) we understand the timing of different
processes arrivals for a concrete collective operation, e.g. all-reduce in an
MPI program. We can evaluate a given PAP by measuring the _process arrival
time_ (PAT) for each process. Formally, a PAP is defined as the tuple
$(a_{0},a_{1},\ldots a_{P-1})$, where $a_{i}$ is a measured PAT for process
$i$, while $P$ is the number of processes participating in the collective
operation. Similarly _process exit pattern_ (PEP) is defined as the tuple
$(f_{0},f_{1},\ldots f_{P-1})$, where $f_{i}$ is the time when process $i$
finishes the operation Faraj2008 . Fig. 1 presents an example of arrival and
exit patterns.

For each process participating in a particular operation Faraj et al. define
the _elapsed time_ as $e_{i}=f_{i}-a_{i}$, and the _average elapsed time_ for
the whole PAP: $\bar{e}=\frac{1}{P}\sum_{i=0}^{P-1}e_{i}$ Faraj2008 . This is
a mean value of time spent for communication by each process, the rest of the
time is used for computations. Thus minimizing the elapsed times of the
participating processes decreases the total time of program execution and, in
our case, is the goal of the optimization.

For an all-reduce operation, assuming $\delta$ to be the time of sending the
reduced data between any two processes and only one arbitrary chosen process
$k$ is delayed (others have the same arrival time, see an example in Fig. 1,
where $k=1$), we can estimate the lower bound of $\bar{e}$ as
$\bar{e}_{lo}=a_{k}-a_{o}+\delta$, where $a_{o}$ is the time of arrival of all
processes except $k$. On the other hand, assuming $\Delta$ to be time of an
all-reduce operation for perfectly balanced PAT ($a_{i}=a_{j}$, for all
$i,j\in\langle 0,k-1\rangle$), we can estimate the upper bound of $\bar{e}$ as
$\bar{e}_{up}=a_{k}-a_{o}+\Delta$. Thus using PAP optimized all-reduce
algorithm can decrease the elapsed time by $\Delta-\delta$ or less. E.g. for a
typical ring all-reduce algorithm, working on 16 processes, 4 MB data size and
1 Gbps Ethernet network we can measure $\Delta=$45.7 ms and $\delta=$18.2 ms,
thus using PAP optimized algorithm can save at most 27.5 ms of average elapsed
time, no matter how slow the delayed process is.

Furthermore, we assume a typical iterative processing model with two phases:
the computation phase where every process performs independent calculations,
and the communication phase where the processes exchange the results, in our
case, using all-reduce collective operation. These two phases are combined
into an iteration, which is repeated sequentially during the processing. We
assume that the whole program execution consists of $N$ iterations.

Normally, during a computation phase, the message communication between
processes/nodes is suspended. Nevertheless, each process can contain many
threads carrying on the parallel computations exploiting shared memory for
synchronization and data exchange. Thus, during this phase the communication
network connecting nodes is unused and can be utilized for exchange of
additional messages, containing information about a progress of computations
or other useful events (e.g. a failure detection).

Thus, we introduce an additional thread which is responsible for inter-node
communication during computation phase. It monitors the progress of the phase
on its node, estimates the remaining computation time and exchanges this
knowledge with other processes on the cooperating nodes. Gathering this
information, every thread can approximate a process arrival pattern (PAP) for
itself and other processes (see Fig. 3).

Figure 3: Iterative processing model (solid lines), extended by progress
monitoring (dashed lines)

For monitoring purposes the computation threads need to pass the status of
current iteration processing calling a special function: edge(). We assume
that for all processes this call is made after a defined part of performed
computations e.g. 50%, while the exact value is passed as the function
parameter. In our implementation the edge() function is executed in some kind
of callback function, reporting a status of the iteration progress.

Beside the PAP monitoring and estimation functions the additional thread can
be used for other purposes. In our case, the proposed all-reduce algorithms,
described in section 4, are working much better (performing faster message
exchange) when connections between the communicating processes/nodes are
already established, thus we introduced an additional warm-up procedure where
a series of messages are transferred in the foreseen directions of
communication. This operation is performed by the thread after the exchange of
the progress data.

Thus the analysis of the above requirements implemented in the additional
thread, shows that the thread should react accordingly to the following events
(see Fig. 4):

  * •

the beginning of processing, when the thread is informed about the computation
phase start and when it stores the timestamp for the further time estimation,

  * •

the edge, in the middle of processing, when the thread estimates the ongoing
computation phase time for its process, exchanges this information with other
processes and performs the warm up procedure for establishing connections to
speed up the message transfer in the coming collective operation,

  * •

the finish of processing, when the thread is informed about the end of the
computation phase and when it stores the timestamp for the further time
estimations.

[width=8cm]Fig4

Figure 4: Events and actions performed in the additional thread

##  4 New all-reduce algorithms optimized for PAPs

In this section we introduce two new algorithms for all-reduce operations,
optimized for a PAP observed during the computation phase: (i) sorted linear
tree (SLT), (ii) pre-reduced ring (PRR). Both of them are based on the well
known, and widely used regular all-reduce algorithms: linear tree OpenMPI  and
ring Thakur2005 , respectively, and have similar communication and time
complexity.

###  4.1 Sorted linear tree (SLT)

The algorithm is an extension of the linear tree OpenMPI , which transfers the
data segments sequentially through processes exploiting the pipeline parallel
computation model. The proposed modification causes the processes to be sorted
by their arrival times. While the faster processes start the communication
earlier, the later ones have more time to finish the computations.

Fig. 5 presents pseudo-code of the algorithm. At the beginning new identifiers
are assigned to the processes according to the arrival order, then the data
are split into equal segments, in our case we assume $P$ segments, where $P$
is a number of the cooperating processes. Afterwards the reduce loop is
started (lines: 1-6), where the data segments are transferred and reduced, and
then the override loop is executed (lines: 7-12), where the segments are
distributed back to all processes.

input parameters:
$P$ – number of processes or nodes, we assume one process per node
$a_{i}$ – arrival time of process $i$
$d_{x}$ – $x$-th of $P$ data segments to be reduced, e.g. $d_{0}$ \- first
data segment
$id$ – a new process rank after sorting according to arrival times
variables:
$si$ – segment index
// reduce loop
1\. for $i:=0$ to $P-1$
2\. $si:=i\ mod\ P$
3\. if $id\neq 0$ then // not the first process
4\. $rd$ = receive_segment() // blocking
5\. reduce($d_{si}$, $rd$)
6\. send_segment($d_{si}$, $(id+1)\ mod\ P$) // non-blocking
// override loop
7\. for $i:=0$ to $P-1$
8\. $si:=i\ mod\ P$
9\. if $id\neq P-1$ then // not the last process
10\. $d_{si}:=$ receive_segment() // blocking
11\. if $id<P-2$ then // not the last or penultimate process
12\. send_segment($d_{si}$, $(id+1)\ mod\ P$) // non-blocking

Figure 5: Pseudo-code of the sorted linear tree (SLT) algorithm

Let’s assume $\tau$ is a time period required for transfer and
reduction/overriding of one data segment between any two processes. Fig. 6b
presents an example of SLT execution, where the process $0$ arrival time is
delayed for $4\tau$ in comparison with all other processes. Using the regular
linear tree reduction, the total time of the execution would be $14\tau$ (see
Fig. 6a), while the knowledge of the PAP and procedure of sorting the
processes by their arrival times make the delayed process to be the last one
in the pipeline and cause the whole operation time to be decreased to
$12\tau$.

Figure 6: Example of (a) linear tree and (b) sorted linear tree (SLT)
algorithm executions, showing data distribution between the processes (one
square corresponds to one data segment), where process #0 is delayed by
$4\tau$, a solid arrow represents send and reduce a data segment, an empty
arrow represents send and override a data segment, and $\tau$ is the time of
transfer and reduction/overriding of a data segment

###  4.2 Pre-reduced ring (PRR)

This algorithm is an extended version of ring Thakur2005 , where each data
segment is reduced and then passed to other processes in synchronous manner.
The idea of the algorithm is to perform a number of so-called, reducing pre-
steps, between faster processes (with lower arrival times), and then the
regular processing, like in the typical ring algorithm, is performed.

Fig. 7 presents pseudo-code of the algorithm. First of all, the processes are
sorted by their arrival times and the new ids are assigned, i.e. initially the
processes perform communication in direction from the earliest to the latest
ones. Then data are split into $P$ equal segments and for each such a segment
the number of pre-steps is calculated (lines: 1-6), its value depending
directly on the estimated process arrival times (PATs): $a_{i}$ and the time
of transfer and reduction/overriding of one segment: $\tau$. Knowing the
above, the algorithm resolves where each segment starts and finishes its
processing, thus two process ids are assigned for each segment: $sp_{i}$ and
$rp_{i}$, respectively (lines: 7-13).

input parameters:
$\tau$ – time of transfer and reduction/overriding of one segment
$P$ – number of processes or nodes, we assume one process per node
$a_{i}$ – arrival time of process $i$
$d_{x}$ – $x$-th of $P$ data segments to be reduced, e.g. $d_{0}$ \- first
data segment
$id$ – a new process rank after sorting according to arrival times
variables:
$k_{i}$ – number of pre-steps in process $i$
$si$ – segment index
$sp_{j}$ – id of the first process where the slice $j$ is sent
$rp_{j}$ – process id where the slice j is reduced last time
// calculate a number of pre-steps for all processes
1\. $k_{P-1}:=0$
2\. for $i:=P-2$ to 0
3\. if $a_{P-1}-a_{i+1}\geq(k_{i+1}+1)\times\tau$ then
4\. $k_{i}:=k_{i+1}+1$
5\. else
6\. $k_{i}:=k_{i+1}$
// calculate ids of the processes where a slice id send first time
// and ids of the processes where a slice is reduced last time
7\. $i:=0$
8\. for $j:=0$ to $P-1$
9\. if $i+k_{i}\geq j$ then
10\. $sp_{j}:=i$
11\. else
12\. $i:=i+1$
13\. $rp_{j}:=(sp_{j}+P-1)\ mod\ P$
// calculate what segment start for
14\. $si:=id+k_{id}$
// reduce loop
15\. for $i:=0$ to $P-1$
16\. if $sp_{si}\neq id$ then
17\. $rd:=$ receive_segment() // blocking
18\. reduce($d_{si}$, $rd$)
19\. if $(rpsi+P-1)\ mod\ P\neq is$ then
20\. send_segment($d_{si}$, $(id+1)\ mod\ P$) // non-blocking
21\. $si:=(si+P-1)\ mod\ P$
// override loop
22\. for $i:=0$ to $P-1$
23\. if $rp_{si}\neq id$ then
24\. $d_{si}:=$ receive_segment() // blocking
25\. if $(rp_{si}+P-1)\ mod\ P\neq id$ then
26\. send_segment($d_{si}$, $(id+1)\ mod\ P$) // non-blocking
27\. $si:=(si+P-1)\ mod\ P$

Figure 7: Pseudo-code of the pre-reduced ring (PRR) algorithm

Afterwards, initial value of the segment index: $s_{i}$ is calculated from
which every process starts processing. In the case of the regular ring
algorithm, its value depends on the process $id$ only, however for PRR it also
regards the pre-steps performed by the involved processes (line: 14). Then the
reduce loop is started, the pre-steps and the regular steps are performed in
one block, where the variable $s_{i}$ controls which segments are sent,
received and reduced (lines: 15-21). Similarly, the override loop is
performed, being controlled by the same variable (lines: 22-27).

In the regular ring every process performs $2P-2$ receive and send operations,
thus total number is $P\times(2P-2)$. In case of PRR, this total number is the
same, however the faster processes (the ones which finished computation
earlier) tend to perform more communication while executing the pre-steps. In
the case, when the arrival of the last process is largely delayed to the next
one (i.e. more than $P\times\tau$) it performs only $P-1$ send and receive
operations, while every other process does $2P-1$ ones.

Fig. 8b presents an example of PRR execution, where the process $0$ arrival
time is delayed for $2\tau$ ($2\times$ time of transfer and
reduction/overriding of one data segment) in comparison with all other
processes. Using the regular ring reduction, the total time of the execution
would be $8\tau$ (see Fig. 8b), while the knowledge of the PAP, performing two
additional reduction pre-steps and setting the delayed process to be the last
one in the pipeline, causes the whole operation time to be reduced to $7\tau$.

Figure 8: Example of (a) ring and (b) pre-reduced ring (PRR) algorithm
executions, showing data distribution between the processes (one square
corresponds to one data segment), where process #0 is delayed by $2\tau$, a
solid arrow represents send and reduce a data segment, an empty arrow
represents send and override a data segment, and $\tau$ is the time of
transfer and reduction/overriding of a data segment

##  5 The experiments

The above algorithms were implemented and tested in a real HPC environment,
the following subsections describe a proposed benchmark (including its pseudo-
code and implementation details), the experiment’s setup and provide the
discussion about the observed results.

###  5.1 The benchmark

Fig. 9 presents pseudo-code of the proposed benchmark evaluating the
performance of the proposed algorithms in the real HPC environment using MPI
MPI  for communication purposes. For every execution there is a sequence of
repeated iterations consisting of the following actions:

line 2:

data generation, where the data are randomly assigned,

line 3:

calculation of the emulated computation time, ensuring that the progress
monitoring communication (the additional 100 ms) is completely covered by the
computation phase,

lines 4–8:

the delay mode is applied, in ‘one-late’ mode only one process (id: 1) is
delayed for $maxDelay$, while in ‘rand-late’ mode all processes are delayed
for random time up to $maxDelay$,

lines 9–10:

two MPI_Barrier() calls, making sure all the processes are synchronized,

lines 11–13:

the emulation of computation phase by using sleep() (usleep() in the
implementation) function, including call to edge() function (see section 3),
providing the progress status to the underlying monitoring thread,

lines 14–16:

call to the all-reduce algorithm implementation, including the commands for
the time measurements, for the benchmark purposes we assumed sum as the reduce
operator,

line 17:

checking the correctness of the performed all-reduce operation, using regular
MPI_Allreduce() function,

lines 18–19:

calculation of the elapsed time for the current and for all processes using
MPI_Allgather() function,

line 20:

saving the average elapsed time into the $results$ vector.

input parameters:
$size$ – number of elements (floats) in reduced data
$N$ – number of iterations
$mode$ – mode of delay: one-late/rand-late
$maxDelay$ – maximal delay of the process(es)
$algorithm$ – tested algorithm, one of ring, Rabenseifner, PRR, SLT
$P$ – number of processes
$id$ – process id – MPI rank: $0$…$P-1$
output:
$results$ – vector of measured average elapsed times for the given algorithm
variables:
$halfTime$ – 50% of the emulated computation time
$startTime$ – start time of measurement
$endTime$ – end time of measurement
$myElapsedTime$ – elapsed time measured in the current process
$allElapsedTimes$ – vector of the measured times of all processes
$data$ – vector of data to be reduced
1\. for $i:=1$ to $N$
2\. $data$ := generateRandomData($size$)
3\. $halfTime=(maxDelay+100\,ms)/2$
4\. if $mode=$ ”one-late” then
5\. if $id=1$ then
6\. $halfTime:=halfTime+maxDelay$
7\. else // $mode=$ ”rand-late”
8\. $halfTime:=halfTime+$ random($0$…$maxDelay$)$/2$
9\. MPI_Barrier()
10\. MPI_Barrier()
11\. sleep($halfTime$)
12\. edge($0.5$)
13\. sleep($halfTime$)
14\. $startTime:=$ MPI_Wtime()
15\. makeAllReduce($algorithm$, $data$)
16\. $endTime:=$ MPI_Wtime()
17\. checkCorrectness($data$)
18\. $myElapsedTime:=endTime-startTime$
19\. MPI_Allgather($allElapsedTimes$, $myElapsedTime$, …)
20\. results[i] := average($allElapsedTimes$)

Figure 9: Pseudo-code of the performance benchmark

For the comparison purposes we used two typical all-reduce algorithms: ring
Thakur2005  and Rabenseifner Rabenseifner2004 , they are implemented in the
two most popular open source MPI implementations: OpenMPI OpenMPI  and MPICH
MPICH  respectively, and used for large input data size. To the best knowledge
of the author there are no PAP optimized all-reduce algorithms described in
the literature.

The benchmark was implemented in C language v. C99, compiled using GCC v.
7.1.0, with the maximal code optimization (-O3). The program uses OpenMPI v.
2.1.1 for processes/nodes message exchange and POSIX Threads v. 2.12 for
intra-node communication and synchronization, moreover for managing of dynamic
data structures GLibc v. 2.0 library was used. In this implementation the
reduce operation is based on a sum (equivalent of using MPI_SUM for $op$
parameter in MPI_Allreduce()).

###  5.2 Environment and test setup

The benchmark was executed in a real HPC environment using cluster
supercomputer Tryton, placed in Centre of Informatics - Tricity Academic
Supercomputer and NetworK (CI TASK) at Gdansk University of Technology in
Poland Krawczyk2015 . The cluster consists of homogeneous nodes, where each
node contains 2 processors (Intel Xeon Processor E5 v3, 2.3GHz, Haswell), with
12 physical cores (24 logical ones, due to Hyperthreading technology) and 128
GB RAM memory. In total the supercomputer consists of 40 racks with 1,600
servers (nodes), 3,200 processors, 38,400 compute cores, 48 accelerators and
202 TB RAM memory. It uses fast FDR 56 Gbps InfniBand in fat tree topology and
1 Gbps Ethernet for networking. Total computing power is 1.48 PFLOPS. The
cluster weighs over 20 metric tons.

The experiments were performed using a subset of the nodes, with HT switched
off, grouped in one rack and connected to each other trough a 1 Gbps Ethernet
switch (HP J9728A 2920-48G). The benchmark input parameters were set up to the
following values:

$algorithm$:

ring, Rabenseifner, pre-reduced ring (PRR) and sorted linear tree (SLT);

$size$:

(of data vector) 128 K, 512 K, 1 M, 2 M, 4 M, 8 M of floats (4 bytes long);

$mode$:

(of process delay) one-late (where only one process is delayed by $maxDelay$)
and rand-late (where all processes are delayed randomly up to $maxDelay$);

$maxDelay$:

(of processes arrival times) 0, 1, 5, 10, 50, 100, 500, 1000 ms;

$P$:

(number of processes/nodes) 4, 6, 8, 10, 12, 16, 20, 24, 28, 32, 36, 40, 44,
48;

$N$:

(number of iterations) 64–256, depending on $maxDelay$ (more for lower delay);

###  5.3 Benchmark results

Table 3 presents the results of the benchmark execution for 1 M of floats of
reduced data, 1 Gbps Ethernet network, where only one process was delayed on
48 nodes in a cluster environment of Tryton Krawczyk2015  HPC computer. The
results are presented as absolute values of average elapsed time:
$\bar{e}_{alg}$ and speedup: $s_{alg}$, in comparison with ring algorithm
$s_{alg}=\frac{\bar{e}_{ring}}{\bar{e}_{alg}}$, where $alg$ is the evaluated
algorithm.

In this setup the ring algorithm seems to be more efficient than the
Rabenseifner, thus in further analysis we use the former for reference
purposes. For more balanced PAPs, where the delay is below 5 ms the proposed
algorithms perform worse than the ring, it is especially visible for SLT
(about 25% slower), however PRR shows only slight difference (below 5%). On
the other hand, when the PAP is more imbalanced with the delay over 10 ms,
both algorithms perform much better (up to 15% faster than the ring). For the
really high delays, the results stabilize providing about 17 ms of average
elapsed time savings, causing the total speedup to be lower. The comparison of
influence of changes in the (maximum) delay on the algorithms’ performance is
presented in Fig. 10.

Table 3: Benchmark results for 1 M of floats of reduced data, 1 Gbps Ethernet network, only one process delayed and 48 processes/nodes. Maximum delay is measured in ms and each result consists of two values: the average elapsed times in ms and speedup in comparison with the ring algorithm ($\frac{\bar{e}_{ring}}{\bar{e}_{alg}}$). The bold values indicate better performance in comparison with the ring algorithm Max delay$\rightarrow$ |  0 |  1 |  5 |  10 |  50 |  100 |  500 |  1000
---|---|---|---|---|---|---|---|---
Algorithm$\downarrow$ |  |  |  |  |  |  |  |
Ring |  67.8 |  67.6 |  70.1 |  75.6 |  115.9 |  165.5 |  558.6 |  1047.0
|  1.00 |  1.00 |  1.00 |  1.00 |  1.00 |  1.00 |  1.00 |  1.00
Rabenseifner |  76.3 |  75.9 |  82.9 |  85.7 |  125.0 |  181.6 |  572.0 |  1061.8
|  0.89 |  0.89 |  0.85 |  0.88 |  0.93 |  0.91 |  0.98 |  0.99
SLT |  90.3 |  90.3 |  89.7 |  89.7 |  101.0 |  149.8 |  542.1 |  1031.9
|  0.75 |  0.75 |  0.78 |  0.84 |  1.15 |  1.10 |  1.03 |  1.01
PRR |  70.9 |  70.8 |  71.1 |  73.0 |  101.0 |  149.7 |  541.7 |  1031.5
|  0.96 |  0.95 |  0.99 |  1.04 |  1.15 |  1.11 |  1.03 |  1.01

[width=8cm]Fig10

Figure 10: Benchmark results showing influence of the increasing maximum delay
on average elapsed times of the tested algorithms. The experiments were
performed on 48 nodes for 1 Gbps Ethernet network, one process delayed, and
message size: 1024 K of float numbers. The error bars are set to
$\pm{}2\sigma$ (95% of the measurements for the normal distribution)

Table 4 presents the benchmark results related to PRR in comparison with the
ring algorithm, executed on 48 nodes with 1 Gbps Ethernet interconnecting
network. While the absolute average elapsed time savings of PRR algorithm are
higher for longer messages (up to 82 ms for 8 M), the high speedup values
occur for all message sizes providing relative savings up to 15%. We can
observe that the lower delays cause lager speedup losses (up to 22%, but for
absolute time: 1.5 ms only), which is especially visible for the smallest
message size: 128 K. In general, for the lower delays, the PRR is comparable
with the ring (see Fig. 10 for 1 M of float numbers reduced data size).

Table 4: Comparison of ring and PRR algorithms for 1 Gbps Ethernet and 48 processes/nodes. Maximum delay is measured in ms and size in Kfloats (4$\times$KB). Each entry consists of two values: a difference of the average elapsed times in ms ($\bar{e}_{PRR}-\bar{e}_{ring}$) and speedup: a quotient of elapsed times ($\frac{\bar{e}_{ring}}{\bar{e}_{PRR}}$). The bold values indicate better performance in comparison with the ring algorithm Max delay$\rightarrow$ |  0 |  1 |  5 |  10 |  50 |  100 |  500 |  1000
---|---|---|---|---|---|---|---|---
Size$\downarrow$ | Only one process delayed
|  1.0 |  1.5 |  -5.0 |  -4.1 |  -4.3 |  -3.2 |  -0.9 |  -0.8
128 |  0.85 |  0.78 |  1.53 |  1.26 |  1.08 |  1.03 |  1.00 |  1.00
|  -1.0 |  -0.3 |  -10.3 |  -12.6 |  -16.4 |  -11.6 |  -1.5 |  -1.7
512 |  1.04 |  1.01 |  1.37 |  1.41 |  1.26 |  1.10 |  1.00 |  1.00
|  2.9 |  3.7 |  1.2 |  -3.2 |  -14.1 |  -15.2 |  -15.8 |  -16.0
1024 |  0.96 |  0.95 |  0.98 |  1.04 |  1.14 |  1.10 |  1.03 |  1.02
|  7.8 |  6.4 |  3.4 |  -1.7 |  -8.5 |  -23.9 |  -27.1 |  -27.7
2048 |  0.93 |  0.94 |  0.97 |  1.02 |  1.06 |  1.13 |  1.05 |  1.03
|  7.6 |  6.0 |  5.2 |  2.8 |  -20.9 |  -29.6 |  -38.1 |  -46.8
4096 |  0.96 |  0.97 |  0.97 |  0.99 |  1.10 |  1.12 |  1.06 |  1.04
|  5.8 |  17.0 |  13.7 |  11.7 |  -29.0 |  -52.1 |  -75.7 |  -82.1
8192 |  0.98 |  0.95 |  0.96 |  0.97 |  1.08 |  1.14 |  1.10 |  1.07
Size$\downarrow$ | All processes delayed randomly
|  1.2 |  1.0 |  -1.1 |  -1.1 |  -2.7 |  -1.0 |  -0.9 |  -1.3
128 |  0.82 |  0.85 |  1.11 |  1.09 |  1.09 |  1.02 |  1.00 |  1.00
|  2.3 |  -0.3 |  -9.5 |  -9.3 |  -16.3 |  -11.5 |  -2.3 |  -2.9
512 |  0.91 |  1.01 |  1.37 |  1.35 |  1.38 |  1.18 |  1.01 |  1.01
|  2.9 |  3.8 |  1.7 |  0.7 |  -13.0 |  -13.2 |  -14.1 |  -14.5
1024 |  0.96 |  0.95 |  0.98 |  0.99 |  1.17 |  1.13 |  1.05 |  1.03
|  7.1 |  8.3 |  4.8 |  3.0 |  -8.8 |  -21.2 |  -23.9 |  -24.5
2048 |  0.94 |  0.92 |  0.96 |  0.97 |  1.08 |  1.16 |  1.07 |  1.04
|  6.1 |  0.1 |  3.2 |  5.6 |  -21.8 |  -19.0 |  -42.6 |  -39.8
4096 |  0.97 |  1.00 |  0.98 |  0.97 |  1.11 |  1.09 |  1.11 |  1.06
|  13.8 |  16.3 |  6.8 |  7.1 |  0.4 |  -34.4 |  -75.9 |  -75.1
8192 |  0.96 |  0.95 |  0.98 |  0.98 |  1.00 |  1.10 |  1.15 |  1.10

The mode of the introduced PAT delay influences slightly the measured values.
For larger delays and message sizes, when only one process is delayed the PRR
algorithm provides slightly smaller relative savings than in the case when the
delays were introduced for all processes, with the uniform probabilistic
distribution, it is especially visible for 500–1000 ms delays and message
sizes of 2–8 M of floats. However, in general, the PRR works fine for both
modes of PAT delay.

Fig. 11 presents the behavior of the tested algorithms in a context of the
changing data size. While the PRR algorithm performs well for a wide range of
size values, the SLT lags for the larger data, the threshold depends on the
(maximum) delay, e.g. for 50 ms the SLT is worse than the regular ring for 4+
M of float numbers data size.

[width=8cm]Fig11

Figure 11: Benchmark results showing the algorithms’ behavior related to the
reduced data size. The experiments were performed on 48 nodes for 1 Gbps
Ethernet network, one process delayed, and the (maximum) delay equals 50 ms.
The error bars are set to $\pm{}2\sigma$ (95% of the measurements for the
normal distribution)

Fig. 12 shows the tested algorithms performance related to the increase of the
number of the nodes/processes exchanging reduced data. For 32 nodes, we
observed interesting behavior of the Rabenseifner algorithm, which was
designed for power of 2 node/process numbers. Both the PRR and SLT algorithms
present stable speedup over the regular ring algorithm, proving good usability
and high scalability for imbalanced PAPs.

[width=8cm]Fig12

Figure 12: Benchmark scalability results for 1 GB Ethernet network, one
process delayed for 100 ms, 1 M of float numbers data size. The error bars are
set to $\pm{}2\sigma$ (95% of the measurements for the normal distribution)

##  6 All-reduce PAP optimization for training of a deep neural network

In this section we present a practical application of the proposed method for
a deep learning iterative procedure, implemented using tiny-dnn open-source
library tiny-dnn . The example is focused on training of a convolutional
neural network to classify graphical images (photos).

For the experiments we used a training dataset usually utilized for
benchmarking purposes: CIFAR-10 cifar10 , it contains 60,000 32$\times$32
color images grouped in 10 classes. There are 50,000 training images, exactly
5,000 images per class. The other 10,000 test images are used to evaluate
training results. In our example we assess the performance of the distributed
processing, especially collective communication, thus we did not need to use
the test set.

We test the network architecture 8 layers: 3 convolutional, 3 average pooling,
and 2 fully connected ones, the whole model has 145,578 parameters of float
size. Each process trains such a network and after each training iteration it
is averaged over the other processes, a similar procedure was described in
Dean2012 , with the distinction in using a separated parameter server. The
mini-batch size was set to 8 images per node, what gives 390 iterations in
total.

The training program was implemented in C++ language (v. C++11) using tiny-dnn
library tiny-dnn . It was chosen because it is very lightweight: header only,
easy to install, dependency-free and has an open-source license: BSD 3-Clause.
The above features made it easy to introduce the required modification: a
callback function for each neural network layer, called during the distributed
SGD (Stochastic Gradient Descent Dean2012 ) training. The function is used for
progress monitoring, where the egde() function (see section 3) is called just
before the last layer is processed, which takes about 44% of computation time
of the iteration. Additionally, the computation part of each iteration is
performed in parallel, using POSIX threads pthreads  executed on available
(24) cores (provided by two processors with Hyper Threading switched off).

We tested 3 all-reduce algorithms: ring, sorted linear tree (SLT) and pre-
reduced ring (PRR). The PAP framework, including estimation of computation
time and warm-up was initiated only for the latter two. The benchmark was
executed 128 times for each algorithm, and each execution consists of 390
training iterations with all-reduce function calls. The tests were executed in
an HPC cluster environment, consisting of 16 nodes with 1 Gbps Ethernet
interconnecting network (the configuration is described in section 5.2).

The results of the benchmark are presented in Table 5. The PRR has the lowest
time of communication: 31.9 ms (average elapsed time of all-reduce) as well as
the total time of training: 78.6 s, SLT is slightly worse having 33.1 ms and
78.8 s respectively. These algorithms were compared to a typical ring
Thakur2005  implementation, where the average elapsed time equals 35.7 ms and
the total training took 81.9 s. In context of the average elapsed time, the
PRR and SLT algorithms are faster for 12.1% and 7.9%, while for the training
total time 4.2% and 4.0% respectively.

Table 5: Times and speedup of Cifar10 cifar10  benchmark execution. The times are presented in ms, and the speedup is calculated in comparison with a ring algorithm implementation Algorithm |  All-reduce average elapsed time |  All-reduce speedup |  Training total time |  Training speedup
---|---|---|---|---
Ring |  35.7 |  1.000 |  81,900 |  1.000
SLT |  33.1 |  1.079 |  78,768 |  1.040
PRR |  31.9 |  1.121 |  78,604 |  1.042

The above results seem to be just a slight improvement, however in the
current, massive processing systems, e.g. neural networks, which are trained
using thousands of compute nodes, consuming MegaWatts of energy, with budgets
of millions of dollars, introducing 4% computing time reduction, without
additional resource demand can provide great cost savings.

##  7 Final remarks

The proposed algorithms provide optimizations for all-reduce operations
executed in environment of imbalanced PAPs. The experimental results show
improved performance and good scalability of the proposed solution over
currently used algorithms. The real case example (machine learning using
distributed SGD Dean2012 ) shows usability of the prototype implementation
with sum as the reducing operation. The solution can be used in a wide
spectrum of applications using iterative computation model, including many
machine learning algorithms.

The future works cover the following topics:

  * •

evaluation of the method for a wider range of interconnecting network speeds
and larger number of nodes using a simulation tool e.g. Czarnul2017 ;
Proficz2016 ,

  * •

expansion of the method for other collective communication algorithms, e.g.
all-gather,

  * •

a framework for automatic PAP detection and proper algorithm selection, e.g.
providing a regular ring for balanced PAPs and PRR for imbalanced ones,

  * •

introduction of the presented PAT estimation method for other purposes e.g.
asynchronous SDG training Dean2012  or deadlock and race detection in
distributed programs Krawczyk2000 ,

  * •

deployment of the solution in a production environment.

We believe that the ubiquity of the imbalanced PAPs in HPC environments
Faraj2008  will cause a fast development of new solutions related to this
subject.

## References

  * [1] CIFAR-10 and CIFAR-100 datasets.  URL: https://www.cs.toronto.edu/~kriz/cifar.html.  Accessed: 2018-01-04.
  * [2] MPI 3.1 collective communication.  URL: http://mpi-forum.org/docs/mpi-3.1/mpi31-report/node95.htm.  Accessed: 2018-01-26.
  * [3] MPICH high-performance portable MPI.  URL: https://www.mpich.org/.  Accessed: 2017-09-07.
  * [4] Open MPI: Open source high performance computing.  URL: https://www.open-mpi.org/.  2017-08-27.
  * [5] POSIX threads programming.  URL: https://computing.llnl.gov/tutorials/pthreads/.  Accessed: 2018-01-05.
  * [6] The standarization forum for Messsage Passing Interface (MPI).  URL: http://mpi-forum.org/.  Accessed: 2018-01-24.
  * [7] tiny-dnn header only, dependency-free deep learning framework in c++.  URL: https://github.com/tiny-dnn/tiny-dnn.  Accessed: 2018-01-04.
  * [8] P. Czarnul, J. Kuchta, M. Matuszek, J. Proficz, P. Rościszewski, M. Wójcik, and J. Szymański.  MERPSYS: An environment for simulation of parallel application execution on large scale HPC systems.  Simulation Modelling Practice and Theory, 77:124–140, 2017.
  * [9] Jeffrey Dean, Greg Corrado, Rajat Monga, Kai Chen, Matthieu Devin, Quoc V Le, Mark Mao, Marc’Aurelio Ranzato, Andrew Senior, Paul Tucker, Ke Yang, and Andrew Y Ng.  Large Scale Distributed Deep Networks.  In Advances in Neural Information Processing Systems 25: 26th Annual Conference on Neural Information Processing Systems 2012, pages 1223–1231. Curran Associates, Inc., 2012.
  * [10] A. Faraj, X. Yuan, and D. Lowenthal.  STAR-MPI: self tuned adaptive routines for MPI collective operations.  Proceedings of the 20th annual international conference on Supercomputing, pages 199—-208, 2006.
  * [11] Ahmad Faraj, Pitch Patarasuk, and Xin Yuan.  A Study of Process Arrival Patterns for MPI Collective Operations.  International Journal of Parallel Programming, 36(6):543–570, dec 2008.
  * [12] Roger W. Hockney.  The communication challenge for MPP: Intel Paragon and Meiko CS-2.  Parallel Computing, 20(3):389 – 398, 1994.
  * [13] Henryk Krawczyk, Bartosz Krysztop, and Jerzy Proficz.  Suitability of the time controlled environment for race detection in distributed applications.  Future Generation Computer Systems, 16(6):625–635, apr 2000.
  * [14] Henryk Krawczyk, Michał Nykiel, and Jerzy Proficz.  Tryton Supercomputer Capabilities for Analysis of Massive Data Streams.  Polish Maritime Research, 22(3):99–104, jan 2015.
  * [15] P. Marendic, J. Lemeire, D. Vucinic, and P. Schelkens.  A novel MPI reduction algorithm resilient to imbalances in process arrival times.  Journal of Supercomputing, 72:1973––2013, 2016.
  * [16] Petar Marendić, Jan Lemeire, Tom Haber, Dean Vučinić, and Peter Schelkens.  An Investigation into the Performance of Reduction Algorithms under Load Imbalance.  In Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics), volume 7484 LNCS, pages 439–450. 2012.
  * [17] Pitch Patarasuk and Xin Yuan.  Efficient MPI_Bcast across different process arrival patterns.  IPDPS Miami 2008 - Proceedings of the 22nd IEEE International Parallel and Distributed Processing Symposium, Program and CD-ROM, (1), 2008\.
  * [18] Jerzy Proficz and Paweł Czarnul.  Performance and Power-Aware Modeling of MPI Applications for Cluster Computing.  In Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics), volume 9574, pages 199–209. 2016.
  * [19] Rolf Rabenseifner.  Optimization of Collective Reduction Operations.  In Lecture Notes in Computer Science, volume 3036, pages 1–9. 2004\.
  * [20] Rajeev Thakur, Rolf Rabenseifner, and William Gropp.  Optimization of Collective Communication Operations in MPICH.  The International Journal of High Performance Computing Applications, 19(1):49–66, feb 2005.


# Self-interacting dark baryons

James M. Cline  McGill University, Department of Physics, 3600 University St.,
Montréal, QC H3A2T8, Canada  Cédrick Perron  University of Toronto, Department
of Physics, 60 St George St, Toronto, ON M5S1A7, Canada

###### Abstract

Using results from lattice QCD, it is possible to quantitatively design models
of dark baryons leading to velocity-dependent self-interaction cross sections
that match the values needed for solving small-scale structure problems of
standard cold dark matter. However it is not obvious that the main dark matter
component in such models will be nucleons rather than large nuclei, or dark
pions or atoms, whose scattering properties would be different. We first
identify the parameters of a dark SU(3) sector analogous to QCD—the
confinement scale $\Lambda$ and pion mass $m_{\pi}$—needed to reproduce
desired self-interaction cross sections. Then we show that these values can
generically be compatible with the absence of a sufficiently stable deuteron
bound state, and hence leading to no heavier nuclei, thus establishing the
consistency of the scenario for self-interacting dark nucleons. The range of
dark photon masses needed to avoid dominant pion or atomic dark matter is
determined, as well as allowed values for the kinetic mixing parameter. The
dark proton might be detected directly in future searches, by dark photon
exchange.

##  I Introduction

Over two decades ago, it was proposed that dark matter (DM) with strong self-
interactions Spergel:1999mh ; Dave:2000ar  could address a discrepancy between
gravitational $N$-body simulations of galaxy formation, which predict cuspy
central density profiles Navarro:1995iw ; Navarro:1996gj , versus observations
that indicate otherwise 1994ApJ…427L…1F ; Walker_2011 ; 10.1093/mnras/stu474 ;
de_Blok_2001 ; 2002A&A…385..816D ; 10.1111/j.1365-2966.2004.07836.x ;
Newman_2013 . Subsequent to the early observations, several other
discrepancies were identified, known as the missing satellites
10.1093/mnras/264.1.201  and too-big-too-fail 10.1111/j.1745-3933.2011.01074.x
; 10.1093/mnras/stu1477  problems, that could be ameliorated by including DM
self-interactions in the numerical simulations 10.1093/mnras/sts514 ;
10.1093/mnras/sts535 . The missing satellites problem has tended to disappear
over time (even becoming a “too many satellites problem”), as observations and
simulations have improved Kim:2017iwr ; 2019MNRAS.486.4545F ; Kim:2021zzw ,
but the core-cusp problem seems more robust.

The small scale structure issue is complicated by baryonic feedback,
originating from supernova shock waves sweeping material out of the dense
inner regions of galaxies 10.1093/mnras/stz1890 , or active galactic nuclei
10.1111/j.1365-2966.2012.20879.x , which were neglected in early simulations.
These effects are difficult to simulate from first principles, due to the vast
range of distance scales that need to be considered in structure formation.
The results to date leave doubt as to whether baryonic physics by itself can
resolve the small-scale structure puzzles. Dark matter self-interactions thus
continue to be an interesting possible resolution.

The magnitude of the interaction cross section needed is not far below the
upper limit implied by Bullet Cluster observations Randall:2008ppe ;
Markevitch:2003at ; Robertson:2016xjh ; Wittman:2017gxn , $\sigma/m\sim 2$
cm2/g $\cong 4$ b/GeV. These large values are suggestive of the strong
interactions of the standard model (SM), making it natural to consider
composite dark matter candidates similar to SM nucleons. An early study
Cline:2013zca  showed that it is possible to rescale lattice QCD results to
make quantitative predictions for the nucleon scattering cross sections in a
dark SU(3) model with different confinement scale $\Lambda$ and pion mass
$m_{\pi}$ than in the SM. Curves in the $m_{\pi}$-$\Lambda$ parameter space
consistent with the desired value of $\sigma/m$ were identified, assuming a
constant (velocity-independent) cross section.

Since that time, it was shown that the cusp-core problem is less pronounced on
the scales of galactic clusters than in individual galaxies Newman:2012nv ;
Newman:2012nw , suggesting that a unified solution requires velocity-dependent
DM scattering, which would be weaker for the higher DM velocity dispersions
found in clusters Kaplinghat:2015aga . As was discussed in Refs.
Mahbubani:2019pij ; Chu:2019awd , velocity-dependent DM scattering is
naturally accommodated within the effective range formalism for parametrizing
scattering amplitudes. In this work we extend the analysis of Ref.
Cline:2013zca  to include the velocity dependence in a QCD-like dark sector,
which reduces the degeneracy of the favored values of $m_{\pi}$ and $\Lambda$
to discrete regions of parameter space.

These conclusions are only relevant if the dark sector consists primarily of
nucleons and not higher-mass dark nuclei, whose scattering properties would
differ from those being calculated here. For example, Ref. Krnjaic:2014xza
(see also Ref. Redi:2018muu ) showed that dark nucleosynthesis can easily be
dominated by high-mass nuclei in a generic confining dark sector. Moreover if
there are dark electrons, then nucleon-nucleon scattering might be superseded
by atom-atom scattering. One must further ensure that dark pions do not
dominate the DM. Hence a further goal of this work is to identify the other
conditions needed to establish that nucleons constitute the dominant DM
component, while fulfilling their QCD-like nature that allows us to
incorporate results from lattice QCD. We will show that these requirements put
a lower bound on the dark photon mass.

Figure 1: Data points taken from Ref. Kaplinghat:2015aga  for DM self-
scattering cross sections $\langle\sigma v\rangle/m$ versus characteristic
relative velocity of dark matter particles in galaxies or galactic clusters,
and predictions from the dark baryon model (black, present work) and a
previously studied DM model with dark photon exchange (orange, Ref.
Kaplinghat:2015aga ). Diagonal lines show $\langle\sigma v\rangle/m$ assuming
different values of a velocity-independent cross section. Brown line shows
best-fit phenomenological power law.

##  II Dark nucleon framework

We assume a dark sector which, like in the SM, has two flavors of light
quarks111and possibly a third heavier (strange) quark, included in some of the
lattice simulations relevant for this study Beane:2005rj , whose masses are
sufficiently below the confinement scale so that pions can be treated as
pseudo-Goldstone bosons. This is sufficient for utilizing determinations of
nucleon scattering lengths $a_{s}$ from lattice data, presented in Ref.
Chen:2010yt . The general form of the $S$-wave scattering amplitude is
Chu:2019awd

${\cal A}_{s}={4\pi\over m_{N}(-ip-a_{s}^{-1}+{\textstyle{1\over
2}}r_{0,s}p^{2}+O(p^{4}))}\,,$ (1)

where $s$ is the spin state of scattering particles (either singlet $s=0$ or
triplet $s=1$). The parameter $r_{0,s}$ is the effective range, which was
neglected in Ref. Cline:2013zca , but a priori it could be relevant in the
present study, since it can affect the velocity-dependence of the cross
section at low energy, via

$\sigma=\sum_{s=0}^{1}(2s+1)\,{m_{N}^{2}|{\cal A}_{s}|^{2}\over 16\pi}\,.$ (2)

Figure 2: Best-fit regions of parameter space for velocity-dependent dark
baryon scattering. Upper (lower) row: scattering dominated by pole in singlet
(triplet) scattering length, showing 68%, 95% and 99% confidence regions.
Figure 3: Contours of $\log_{10}[(\sigma/m_{N})/(0.6\,{\rm cm}^{2}/{\rm g})]$
in the $v\to 0$ limit. (colored). The (blue) ones labeled “0” represent the
degenerate solutions for the small-scale structure problems, neglecting
velocity dependence of the cross section. Black regions show the more narrow
predictions accounting for the velocity dependence, at 68% confidence.

The parameters $a_{s}$ and $r_{0,s}$ have been fitted to lattice QCD data as a
function of the pion mass in Ref. Chen:2010yt . By dimensional analysis, these
results can be generalized to a QCD-like theory with a different confinement
scale, by assuming that $a_{s}^{-1}$ scales linearly with $\Lambda$, and
taking a specific value for QCD to determine the dimensionless proportionality
constant; we adopt $\Lambda=250\,$MeV.

We find that the scattering lengths, determined in Ref. Chen:2010yt  using
nuclear effective field theory Beane:2008bt ), can be accurately represented
using the simple analytic functions222This corrects formulas in Ref.
Cline:2013zca  that had erroneous digitization of the curves

$a_{0}={0.71\,\Lambda^{-1}\over m_{\pi}/\Lambda-0.58},\quad
a_{1}={0.96\,\Lambda^{-1}\over m_{\pi}/\Lambda-0.42}\,,$ (3)

given our choice of $\Lambda=250\,$MeV for QCD. In the present study we found
no such simple formulas for fitting $r_{0,s}$ ; instead we digitized the
results from Fig. 2 of Ref. Chen:2010yt . However it turns out that
$r_{0,s}\ll a_{s}$ in the regions of parameter space of interest; hence
$r_{0,s}$ can be neglected without appreciably affecting our results for
scattering.

The cross section depends on the relative velocity through $p=m_{N}v/2$ in
(1), where $m_{N}\cong 3.76\,\Lambda$ by normalizing to QCD (in the limit of
light quarks whose mass has a negligible effect on $m_{N}$). To fit the two
free parameters $\Lambda$ and $m_{\pi}$, we used the inferred values of
$\langle\sigma v\rangle/m$ for a selection of galaxies and galactic clusters
from Fig. 1 of Ref. Kaplinghat:2015aga  to construct a $\chi^{2}$ function,
and minimized it. This requires performing the phase space average over the DM
velocity distributions, to compute $\langle\sigma v\rangle$. Assuming a
Maxwellian $e^{-v^{2}/v_{0}^{2}}$ distribution, the averaging can be carried
out analytically, in the approximation of ignoring $r_{0,s}$ in Eq. (1),
giving (in $\hbar=c=1$ units)

$\displaystyle\langle\sigma v\rangle$ $\displaystyle=$
$\displaystyle\sum_{s}(2s+1){a_{s}\over m_{N}}F(b_{s});$ $\displaystyle F(b)$
$\displaystyle=$ $\displaystyle\sqrt{\pi/2}\left(2b+b^{3}e^{b^{2}/2}\,{\rm
Ei}(-b^{2}/2)\right)$ (4)

where $b_{s}\equiv 2/(m_{N}a_{s}v_{0})$, Ei is the exponential integral
function, and $v_{0}=\sqrt{\pi}\langle v\rangle/2$ in terms of the average DM
velocity $\langle v\rangle$.

Fig. 1 shows the data points and one of our best-fit predictions (black solid
line), as well as the prediction from an alternative model involving exchange
of a light dark photon (solid orange) Kaplinghat:2015aga . The data points
were determined in Ref. Kaplinghat:2015aga  using a semianalytic model of halo
profiles calibrated with N-body simulations, assuming a large enough value of
$\langle\sigma v\rangle/m$ to explain observed coring of the DM profile in
correlation with the velocity dispersion of the system.

We find four best-fit regions of parameter space, shown in Fig. 2. These fall
on the previously identified best-fit regions (blue curves) from Ref.
Cline:2013zca , shown in Fig. 3, where the velocity-dependence of $\sigma$ was
neglected. The four-fold degeneracy can be understood from Eq. (3): a large
enough cross section requires being close to one of the poles of the two
scattering lengths, either slightly above or below. $m_{\pi}/\Lambda$ must be
tuned to one part in $\sim 50$ at the 99% confidence level (C.L.). On the
other hand, the value of $\Lambda$ is more weakly constrained, $\Lambda\cong
3.0{+1.7\atop-1.0}\,$GeV near the singlet pole or $\Lambda\cong
2.1{+1.2\atop-0.6}\,$GeV near the triplet pole. We find a minimum $\chi^{2}$
of 22, lower than that of the best-fit dark photon model shown in Ref.
Kaplinghat:2015aga  (which also has two free parameters, the DM and photon
masses, with $\alpha^{\prime}$ held fixed), whose $\chi^{2}$ is 37. The latter
is dominated by contributions from three data points that lie well below the
curve, which cause our best fit curve to be higher.

For comparison, we show the best fit for a hypothetical power law ansatz,
$\langle\sigma v\rangle/m=N(v/100\,{\rm km/s})^{\alpha}$, which also has two
free parameters. It has approximately the same minimum $\chi^{2}$ as the dark
baryon model, with $N=63$ cm2/g$\cdot$ km/s and $\alpha=0.28$. However, we are
not aware of a particle physics model corresponding to this ansatz.

##  III Cosmological consistency constraints

For these results to be relevant, the dark baryons must be the primary DM
components. In the following, we outline the simplest scenario for ensuring
this criterion, and consistency with other constraints. First, the dark pions
(including other pseudoscalar mesons, in the case of three light flavors) and
vector mesons must be unstable or subdominant. Although one could introduce an
analog of weak interactions into the dark sector, a simpler possibility is to
require that the $\pi^{\pm}$ mesons that are charged under a dark U(1)′,
analogous to electromagnetism, annihilate into dark photons $\gamma^{\prime}$
to a low abundance.333We assume that the mechanism that generates the dark
baryon asymmetry does not also create a dark pion asymmetry. The case of
symmetric dark baryons created by freeze-in has been considered in Ref.
Garani:2021zrr .  This implies Steigman:2012nb

$\langle\sigma
v\rangle_{\pi^{+}\pi^{-}\to\gamma^{\prime}\gamma^{\prime}}\sim{g^{4}\over\pi
m_{\pi}^{2}}\gg\langle\sigma v\rangle_{0}\equiv 3\times 10^{-26}{\rm
cm^{3}/s}\,,$ (5)

where $g$ is the U(1)′ coupling. Taking $m_{\pi}\sim 1$ GeV, we find the
modest requirement $g\gg 10^{-4}$. For example if
$\alpha^{\prime}=g^{2}/4\pi=10^{-3}$, dark mesons would constitute only
$10^{-5}$ of the DM. Similar conclusions hold for dark $K^{0}$ and $K^{+}$ in
the three-flavor case (assuming the strange quark mass is not much greater
than $m_{\pi}$).

Even a small residual rate of dark
$\pi^{\pm}\to\gamma^{\prime}\gamma^{\prime}$ annihilations may have observable
effects at late times, since the $\gamma^{\prime}$ decays can distort the
cosmic microwave background (CMB) Slatyer:2015jla , 21-cm signal Liu:2018uzy ,
and big bang nucleosynthesis (BBN) Depta:2019lbe . For GeV-scale DM, the CMB
constrains $\langle\sigma v\rangle\lesssim 10^{-27}$ cm3/s, which is not far
below the fiducial cross section (5). A subdominant DM component like our dark
$\pi^{\pm}$ would thus be unconstrained.

Considering the neutral $\pi^{0}$ and $\eta$, their decays into
$\gamma^{\prime}\gamma^{\prime}$ will be fast, so long as
$m_{\gamma}^{\prime}<m_{\pi}/2\sim 500\,$MeV. Similarly, dark glueballs have a
mass of order $7\Lambda$ in QCD, which is far above the threshold for decay
into two or three mesons in our scenario. They are therefore short-lived and
pose no risk as DM relics. One expects vector mesons to decay to pions, like
in QCD.

If the dark photon is still massless when the dark proton asymmetry is
generated in the early universe, there should be a compensating dark electron
($e^{\prime}$) asymmetry to ensure U(1)′ charge neutrality. There is then the
danger that dark H atoms constitute a large fraction of the DM, and self-
interact more strongly than nucleons. This can be avoided if the dark photon
mass is large enough to inhibit the formation of bound states Petraki:2016cnz
; Cline:2021itd :

$m_{\gamma^{\prime}}\gtrsim\alpha^{\prime}{m_{e^{\prime}}m_{N}\over
m_{e^{\prime}}+m_{N}}\,,$ (6)

where $\alpha^{\prime}=g^{2}/4\pi$. Since this depends on the dark electron
mass, it does not give an independent constraint on $m_{\gamma^{\prime}}$. In
fact, for the preferred regions of $m_{N}$ and the small value
$\alpha^{\prime}\sim 10^{-3}$ adopted below, the constraint (6) is satisfied
for any $m_{e^{\prime}}$.

Next we consider how to prevent dark nucleons from binding significantly into
larger nuclei. It is interesting that the criterion of being close to the pole
of the scattering length implies that there is no bound deuteron $d$ in the
case where $m_{\pi}/\Lambda=0.42+\epsilon$, which provides an obstruction to
nucleosynthesis from proceeding through the reaction $n+p\to
d+\gamma^{\prime}$. More generally, one can block dark BBN (big bang
nucleosynthesis) by taking $m_{\gamma^{\prime}}$ to exceed the binding energy
of the deuteron, or of the spin singlet dinucleon state in the isotriplet
channel $N_{i}+N_{j}\to D_{k}+\gamma^{\prime}$, where $N_{i}$ are the nucleon
states and $D_{k}$ is the dinucleon isotriplet.

The binding energies of $d$ and $D_{k}$ as determined by lattice QCD studies,
combined with effective field theory, are given as a function of $m_{\pi}$ in
Ref. Chen:2010yt . In fact they are given in terms of the effective range
parameters by
$E_{b,s}=(m_{N}r_{0,s}^{2})^{-1}(1-\sqrt{1-2r_{0,s}/a_{s}})^{2}\cong(m_{N}a_{s})^{-2}$
for $|r_{0,s}|\ll|a_{s}|$, and binding can only occur for the ranges of
$m_{\pi}/\Lambda$ where $a_{s}>0$ and $0<2r_{0,s}/a_{s}<1$. (These ranges are
such that there is no overlap between the singlet and triplet: at most one
channel can have a bound state.) We find that the binding energies are of
order $0.01$ MeV. Therefore if $m_{\gamma^{\prime}}\gtrsim 0.01\,$MeV, dark
BBN is generally inhibited. The dark pions themselves are much heavier than
this scale, and therefore cannot serve as a light mediator for carrying away
binding energy to form $d$ or $D_{k}$. Aside from kinematically blocking $d$
or $D_{k}$ production, a sufficiently small coupling $\alpha^{\prime}$ will
impede these processes, as was quantified in Ref. Redi:2018muu .

Beyond dark BBN considerations, the fusion of dark nucleons into dinucleon
bound states plus dark photon is constrained by indirect signals, notably the
effect of energy injection from dark photon decays into the cosmic microwave
background (CMB) Mahbubani:2019pij . By blocking bound state production we
also satisfy these constraints.

A further requirement is that nucleon self-interactions mediated by dark
photon exchange should be subdominant to the strong interactions. This puts a
stronger constraint on $m_{\gamma^{\prime}}$ than does dark BBN. The
scattering can be computed nonrelativistically using the Yukawa potential
$\alpha^{\prime}e^{-m_{\gamma^{\prime}}r}/r$. Depending on the values of the
parameters $R=(m_{N}v/m_{\gamma^{\prime}})^{2}$ and
$Q=\alpha^{\prime}m_{N}/m_{\gamma^{\prime}}$, where $v$ is the relative
velocity, the cross section can be Sommerfeld enhanced. If both $R,Q\ll 1$,
the Born approximation is valid, and the momentum transfer cross section is
given by Tulin:2013teo

$\sigma_{T}={8\pi\alpha^{\prime 2}\over m_{N}^{2}v^{4}}\left[\ln(1+R)-{1\over
1+R^{-1}}\right]\,.$ (7)

For $\alpha^{\prime}=10^{-3}$ and $m_{N}\sim 10\,$GeV, this approximation is
adequate for $m_{\gamma^{\prime}}\gg 10\,$MeV; otherwise a numerical solution
of the Schrödinger equation going beyond perturbation theory is required. For
simplicity, we use Eq. (7) to estimate the lower bound on
$m_{\gamma^{\prime}}$ (taking benchmark value $\alpha^{\prime}=10^{-3}$), and
cross-check our conclusion against the numerical results of Ref. Tulin:2013teo
. The $v$-dependence of $\sigma_{T}v/m$ is plotted in Fig. 4 for several
values of $m_{\gamma^{\prime}}$, and assuming $m=10\,$GeV, compatible with
models in Table 1. Comparison with the data suggests that if
$m_{\gamma^{\prime}}\gtrsim 35\,$MeV, this contribution to dark nucleon
scattering is unimportant. This is corroborated by Fig. 6 of Ref.
Tulin:2013teo , which takes account of nonperturbative effects.

Figure 4: Perturbative predictions of the dark-photon mediated contribution to
scattering of dark photons, for $\alpha^{\prime}=10^{-3}$ and
$m_{\gamma^{\prime}}=10,\,20,\,35,\,50$ MeV, from top to bottom. The DM mass
is taken to be 10 GeV. Experimental data points from Fig. 1 are also plotted.

Lastly, the dark photons must be unstable, to avoid overclosing the Universe
Kaplinghat:2013yxa . This is usually accomplished by introducing a kinetic
mixing Lagrangian $-(\epsilon/2)F^{\mu\nu}F^{\prime}_{\mu\nu}$ with the SM
hypercharge; then $A^{\prime}_{\mu}$ couples to $\epsilon e$ times the SM
currents of charged particles, allowing $\gamma^{\prime}\to e^{+}e^{-}$
decays. For $m_{\gamma^{\prime}}\sim 35\,$MeV, there are windows of allowed
$\epsilon$ in the vicinity of $10^{-4}$ Fabbrichesi:2020wbt  and $\lesssim
10^{-10}$, where the latter is the bound from supernova 1987A Chang:2016ntp ;
Li:2020roy . The former value $\sim 10^{-4}$ is excluded in the present model
beam dump experiments (see Fig. 3.4 of Fabbrichesi:2020wbt ). The cross
section for $pp^{\prime}\to pp^{\prime}$ elastic scattering is given by

$\sigma_{p}\cong{(\epsilon egm_{p})^{2}\over\pi
m_{\gamma^{\prime}}^{4}}\lesssim 10^{-45}\,{\rm cm}^{2}\,,$ (8)

where $m_{p}$ is the ordinary proton mass and the upper limit is from Refs.
XENON:2018voc ; PandaX-II:2018xpz . Taking $\epsilon=10^{-10}$ and
$\alpha^{\prime}=10^{-3}$, we find $\sigma_{p}=2\times 10^{-46}\,$cm2. Even
though this is only a factor of 5 less than the current limit, it is below the
sensitivity of currently planned experiments if the dark nucleon mass is
$m_{N}\lesssim 9\,$GeV Billard:2021uyg , which is compatible with the allowed
regions from our scattering fits. However for somewhat larger masses $\gtrsim
10\,$GeV, also allowed, this cross section would be observable through direct
detection, due to the quickly rising sensitivity as the recoil energy
increases.444 Alternatively, one could turn off the kinetic mixing and allow
$\gamma^{\prime}$ to decay into dark radiation (e.g., massless fermions),
thereby circumventing all of these bounds.

##  IV Conclusions

In summary, we find that the strong interactions of a QCD-like dark nucleon
sector could have the desired velocity-dependent self-scattering cross section
for solving the cusp-core problem of dark matter halos, on scales from dwarf
galaxies to galactic clusters. At the same time, a variety of self-consistency
requirements must be satisfied for this scenario to be realistic,
necessitating the existence of dark photons and electrons in addition to the
baryons. A portal between the dark sector and the SM is needed to deplete dark
pions, through annihilation or decay, which we took to be kinetic mixing of
$\gamma^{\prime}$ with hypercharge. The $\gamma^{\prime}$ mass should be
between $\sim 35$ MeV (to sufficiently deplete dark pions via annihilation)
and $m_{\pi}/2$ (to allow decay of neutral pions into
$\gamma^{\prime}\gamma^{\prime}$). The $\gamma^{\prime}$ must decay into SM
$e^{+}e^{-}$ or dark radiation; in the former case, kinetic mixing at the
level of $\epsilon\sim 10^{-10}$ is needed to avoid supernova and
nucleosynthesis constraints on $\gamma^{\prime}$. The model predicts a level
of dark proton-proton scattering that could be detectable in direct searches,
depending on $m_{N}$. Benchmark values of allowed parameters are summarized in
Table 1.

${m_{\pi}/\Lambda}$ | $\Lambda$ | $m_{N}$ | $m_{\pi}$ | $m_{\gamma^{\prime}}$ | $\alpha^{\prime}$ | $\epsilon$
---|---|---|---|---|---|---
$\\!\\!\\!\\!{}_{\phantom{|_{|}}}{0.58\atop 0.42}^{\phantom{|^{|}}}\\!\\!\\!\\!$ | ${1.5{\rm-}3.3\atop 2.0{\rm-}4.7}$ | ${5.6{\rm-}12\atop 7.5{\rm-}18}$ | ${0.87{\rm-}1.9\atop 0.84{\rm-}2.0}$ | $\gtrsim 0.035$ | $10^{-3}$ | $10^{-10}$
Table 1: Values of model parameters in benchmark models. $\Lambda$, $m_{N}$
$m_{\pi}$ and $m_{\gamma^{\prime}}$ are in GeV/$c^{2}$ units.

One peculiarity of this scenario is that the ratio $m_{\pi}/\Lambda$ is tuned
at the level of 2% to give a large enough scattering length, in either the
singlet or triplet channel, which has the consequence of making the deuteron
or dinucleon being very close being a zero-energy (quasi)bound state.
Curiously, a similar coincidence occurs in the SM QCD sector, where $m_{\pi}$
is within 5% of the boundary for no deuteron bound state. This is the origin
of the famous deuterium bottleneck of BBN.

In our study we chose a particular EFT matching to rather old lattice QCD data
in Ref. Chen:2010yt . However there are other EFT choices that one could make,
reviewed in Ref. Hammer:2019poc , and newer data. It has been argued in Ref.
Nicholson:2021zwi  that there remains a large degree of theoretical
uncertainty in EFT extrapolations to lower pion masses than can currently be
achieved on the lattice. We have not attempted to quantify these
uncertainties, but we do not expect them to affect the qualitative features of
our conclusions. Namely, all such studies predict scattering lengths of the
general form in Eq. (3). Hence the ultimate values of $\Lambda$ and
$m_{\pi}/\Lambda$ needed to fit cosmological data may differ by factors of
order 1 from those we have found, but the predicted shape of the
$\langle\sigma v\rangle/m$ curve shown in Fig. 1 should remain substantially
unchanged. In any case, the best fit values of the dark QCD parameters will
also depend on the light quark mass spectrum, which is an additional model-
building input.

We have not tried to address how the needed asymmetry between dark baryons and
antibaryons could be generated. Ref. Frandsen:2011kt  proposed a model of dark
baryons whose relic density is naturally achieved through dynamical
electroweak symmetry breaking. It was recently shown that a dark sector with
confining SU(2) interactions can be unified with SU(3)color to provide a UV-
complete explanation of simultaneous baryogenesis in both sectors
Murgui:2021eqf . In this framework, the relative closeness of the asymmetry in
the two sectors, as observed in CMB data, is a consequence of a mild hierarchy
$\Lambda/\Lambda_{QCD}\sim 6$ between the confinement scales. Our model has a
similar hierarchy. It might be interesting to investigate whether unification
of SU$(3)_{\rm dark}\times$SU(3)QCD to SU(6) could provide a more complete
framework for self-interacting dark baryons of the kind we have considered.

Acknowledgments. We thank G. Alonso-Álvarez, K. Moorthy, M. Redi, J.-S. Roux
and A. Tesi for helpful discussions. We thank anonymous referees for very
useful suggestions. This research was supported by the Natural Sciences and
Engineering Research Council (NSERC) of Canada.

## References

  * (1) D. N. Spergel and P. J. Steinhardt, “Observational evidence for self-interacting cold dark matter,” Phys. Rev. Lett. 84 (2000) 3760–3763, arXiv:astro-ph/9909386.
  * (2) R. Dave, D. N. Spergel, P. J. Steinhardt, and B. D. Wandelt, “Halo properties in cosmological simulations of selfinteracting cold dark matter,” Astrophys. J. 547 (2001) 574–589, arXiv:astro-ph/0006218.
  * (3) J. F. Navarro, C. S. Frenk, and S. D. M. White, “The Structure of cold dark matter halos,” Astrophys. J. 462 (1996) 563–575, arXiv:astro-ph/9508025.
  * (4) J. F. Navarro, C. S. Frenk, and S. D. M. White, “A Universal density profile from hierarchical clustering,” Astrophys. J. 490 (1997) 493–508, arXiv:astro-ph/9611107.
  * (5) R. A. Flores and J. R. Primack, “Observational and Theoretical Constraints on Singular Dark Matter Halos,” ”Astrophys. J. Lett.” 427 (May, 1994) L1, arXiv:astro-ph/9402004 [astro-ph].
  * (6) M. G. Walker and J. Peñarrubia, “A Method for measuring (slopes of) the mass profiles of dwarf spheroidal galaxies,” Astrophys. J. 742 no. 1, (Nov, 2011) 20. https://doi.org/10.1088/0004-637x/742/1/20.
  * (7) E. J. Tollerud, M. Boylan-Kolchin, and J. S. Bullock, “M31 satellite masses compared to $\Lambda$CDM subhaloes,” MNRAS 440 no. 4, (04, 2014) 3511–3519. https://doi.org/10.1093/mnras/stu474.
  * (8) W. J. G. de Blok, S. S. McGaugh, A. Bosma, and V. C. Rubin, “Mass density profiles of low surface brightness galaxies,” Astrophys. J. 552 no. 1, (May, 2001) L23–L26. https://doi.org/10.1086/320262.
  * (9) W. J. G. de Blok and A. Bosma, “High-resolution rotation curves of low surface brightness galaxies,” Astron. & Astrophys. 385 (Apr., 2002) 816–846, arXiv:astro-ph/0201276 [astro-ph].
  * (10) G. Gentile, P. Salucci, U. Klein, D. Vergani, and P. Kalberla, “The cored distribution of dark matter in spiral galaxies,” MNRAS 351 no. 3, (07, 2004) 903–922. https://doi.org/10.1111/j.1365-2966.2004.07836.x.
  * (11) A. B. Newman, T. Treu, R. S. Ellis, and D. J. Sand, “The density profiles of massive relaxed galaxy clusters. ii. separating luminous and dak matter in cluster cores,” Astrophys. J. 765 no. 1, (Feb, 2013) 25. https://doi.org/10.1088/0004-637x/765/1/25.
  * (12) G. Kauffmann, S. D. M. White, and B. Guiderdoni, “The formation and evolution of galaxies within merging dark matter haloes,” MNRAS 264 no. 1, (09, 1993) 201–218. https://doi.org/10.1093/mnras/264.1.201.
  * (13) M. Boylan-Kolchin, J. S. Bullock, and M. Kaplinghat, “Too big to fail? The puzzling darkness of massive Milky Way subhaloes,” MNRAS Lett. 415 no. 1, (07, 2011) L40–L44.
  * (14) S. Garrison-Kimmel, M. Boylan-Kolchin, J. S. Bullock, and E. N. Kirby, “Too big to fail in the Local Group,” MNRAS 444 no. 1, (08, 2014) 222–236. https://doi.org/10.1093/mnras/stu1477.
  * (15) M. Rocha, A. H. G. Peter, J. S. Bullock, M. Kaplinghat, S. Garrison-Kimmel, J. Oñorbe, and L. A. Moustakas, “Cosmological simulations with self-interacting dark matter – I. Constant-density cores and substructure,” MNRAS 430 no. 1, (01, 2013) 81–104. https://doi.org/10.1093/mnras/sts514.
  * (16) A. H. G. Peter, M. Rocha, J. S. Bullock, and M. Kaplinghat, “Cosmological simulations with self-interacting dark matter – II. Halo shapes versus observations,” MNRAS 430 no. 1, (01, 2013) 105–120. https://doi.org/10.1093/mnras/sts535.
  * (17) S. Y. Kim, A. H. G. Peter, and J. R. Hargis, “Missing Satellites Problem: Completeness Corrections to the Number of Satellite Galaxies in the Milky Way are Consistent with Cold Dark Matter Predictions,” Phys. Rev. Lett. 121 no. 21, (2018) 211302, arXiv:1711.06267 [astro-ph.CO].
  * (18) C. E. Fielder, Y.-Y. Mao, J. A. Newman, A. R. Zentner, and T. C. Licquia, “Predictably missing satellites: subhalo abundances in Milky Way-like haloes,” MNRAS 486 no. 4, (July, 2019) 4545–4568, arXiv:1807.05180 [astro-ph.GA].
  * (19) S. Y. Kim and A. H. G. Peter, “The Milky Way satellite velocity function is a sharp probe of small-scale structure problems,” arXiv:2106.09050 [astro-ph.GA].
  * (20) A. Benítez-Llambay, C. S. Frenk, A. D. Ludlow, and J. F. Navarro, “Baryon-induced dark matter cores in the eagle simulations,” MNRAS 488 no. 2, (07, 2019) 2387–2404. https://doi.org/10.1093/mnras/stz1890.
  * (21) D. Martizzi, R. Teyssier, B. Moore, and T. Wentz, “The effects of baryon physics, black holes and active galactic nucleus feedback on the mass distribution in clusters of galaxies,” MNRAS 422 no. 4, (05, 2012) 3081–3091. https://doi.org/10.1111/j.1365-2966.2012.20879.x.
  * (22) S. W. Randall, M. Markevitch, D. Clowe, A. H. Gonzalez, and M. Bradac, “Constraints on the Self-Interaction Cross-Section of Dark Matter from Numerical Simulations of the Merging Galaxy Cluster 1E 0657-56,” Astrophys. J. 679 (2008) 1173–1180, arXiv:0704.0261 [astro-ph].
  * (23) M. Markevitch, A. H. Gonzalez, D. Clowe, A. Vikhlinin, L. David, W. Forman, C. Jones, S. Murray, and W. Tucker, “Direct constraints on the dark matter self-interaction cross-section from the merging galaxy cluster 1E0657-56,” Astrophys. J. 606 (2004) 819–824, arXiv:astro-ph/0309303.
  * (24) A. Robertson, R. Massey, and V. Eke, “What does the Bullet Cluster tell us about self-interacting dark matter?,” Mon. Not. Roy. Astron. Soc. 465 no. 1, (2017) 569–587, arXiv:1605.04307 [astro-ph.CO].
  * (25) D. Wittman, N. Golovich, and W. A. Dawson, “The Mismeasure of Mergers: Revised Limits on Self-interacting Dark Matter in Merging Galaxy Clusters,” Astrophys. J. 869 no. 2, (2018) 104, arXiv:1701.05877 [astro-ph.CO].
  * (26) J. M. Cline, Z. Liu, G. D. Moore, and W. Xue, “Composite strongly interacting dark matter,” Phys. Rev. D 90 no. 1, (2014) 015023, arXiv:1312.3325 [hep-ph].
  * (27) A. B. Newman, T. Treu, R. S. Ellis, D. J. Sand, C. Nipoti, J. Richard, and E. Jullo, “The Density Profiles of Massive, Relaxed Galaxy Clusters: I. The Total Density Over 3 Decades in Radius,” Astrophys. J. 765 (2013) 24, arXiv:1209.1391 [astro-ph.CO].
  * (28) A. B. Newman, T. Treu, R. S. Ellis, and D. J. Sand, “The Density Profiles of Massive, Relaxed Galaxy Clusters: II. Separating Luminous and Dark Matter in Cluster Cores,” Astrophys. J. 765 (2013) 25, arXiv:1209.1392 [astro-ph.CO].
  * (29) M. Kaplinghat, S. Tulin, and H.-B. Yu, “Dark Matter Halos as Particle Colliders: Unified Solution to Small-Scale Structure Puzzles from Dwarfs to Clusters,” Phys. Rev. Lett. 116 no. 4, (2016) 041302, arXiv:1508.03339 [astro-ph.CO].
  * (30) R. Mahbubani, M. Redi, and A. Tesi, “Indirect detection of composite asymmetric dark matter,” Phys. Rev. D 101 no. 10, (2020) 103037, arXiv:1908.00538 [hep-ph].
  * (31) X. Chu, C. Garcia-Cely, and H. Murayama, “A Practical and Consistent Parametrization of Dark Matter Self-Interactions,” JCAP 06 (2020) 043, arXiv:1908.06067 [hep-ph].
  * (32) G. Krnjaic and K. Sigurdson, “Big Bang Darkleosynthesis,” Phys. Lett. B 751 (2015) 464–468, arXiv:1406.1171 [hep-ph].
  * (33) M. Redi and A. Tesi, “Cosmological Production of Dark Nuclei,” JHEP 04 (2019) 108, arXiv:1812.08784 [hep-ph].
  * (34) NPLQCD Collaboration, S. R. Beane, P. F. Bedaque, K. Orginos, and M. J. Savage, “I = 2 pi-pi scattering from fully-dynamical mixed-action lattice QCD,” Phys. Rev. D 73 (2006) 054503, arXiv:hep-lat/0506013.
  * (35) J.-W. Chen, T.-K. Lee, C. P. Liu, and Y.-S. Liu, “On the Quark Mass Dependence of Two Nucleon Observables,” Phys. Rev. C 86 (2012) 054001, arXiv:1012.0453 [nucl-th].
  * (36) S. R. Beane, D. B. Kaplan, and A. Vuorinen, “Perturbative nuclear physics,” Phys. Rev. C 80 (2009) 011001, arXiv:0812.3938 [nucl-th].
  * (37) R. Garani, M. Redi, and A. Tesi, “Dark QCD matters,” JHEP 12 (2021) 139, arXiv:2105.03429 [hep-ph].
  * (38) G. Steigman, B. Dasgupta, and J. F. Beacom, “Precise Relic WIMP Abundance and its Impact on Searches for Dark Matter Annihilation,” Phys. Rev. D 86 (2012) 023506, arXiv:1204.3622 [hep-ph].
  * (39) T. R. Slatyer, “Indirect dark matter signatures in the cosmic dark ages. I. Generalizing the bound on s-wave dark matter annihilation from Planck results,” Phys. Rev. D 93 no. 2, (2016) 023527, arXiv:1506.03811 [hep-ph].
  * (40) H. Liu and T. R. Slatyer, “Implications of a 21-cm signal for dark matter annihilation and decay,” Phys. Rev. D 98 no. 2, (2018) 023501, arXiv:1803.09739 [astro-ph.CO].
  * (41) P. F. Depta, M. Hufnagel, K. Schmidt-Hoberg, and S. Wild, “BBN constraints on the annihilation of MeV-scale dark matter,” JCAP 04 (2019) 029, arXiv:1901.06944 [hep-ph].
  * (42) K. Petraki, M. Postma, and J. de Vries, “Radiative bound-state-formation cross-sections for dark matter interacting via a Yukawa potential,” JHEP 04 (2017) 077, arXiv:1611.01394 [hep-ph].
  * (43) J. M. Cline, “Dark atoms and composite dark matter,” in Les Houches summer school on Dark Matter.  8, 2021.  arXiv:2108.10314 [hep-ph].
  * (44) S. Tulin, H.-B. Yu, and K. M. Zurek, “Beyond Collisionless Dark Matter: Particle Physics Dynamics for Dark Matter Halo Structure,” Phys. Rev. D 87 no. 11, (2013) 115007, arXiv:1302.3898 [hep-ph].
  * (45) M. Kaplinghat, S. Tulin, and H.-B. Yu, “Direct Detection Portals for Self-interacting Dark Matter,” Phys. Rev. D 89 no. 3, (2014) 035009, arXiv:1310.7945 [hep-ph].
  * (46) M. Fabbrichesi, E. Gabrielli, and G. Lanfranchi, “The Dark Photon,” arXiv:2005.01515 [hep-ph].
  * (47) J. H. Chang, R. Essig, and S. D. McDermott, “Revisiting Supernova 1987A Constraints on Dark Photons,” JHEP 01 (2017) 107, arXiv:1611.03864 [hep-ph].
  * (48) J.-T. Li, G. M. Fuller, and E. Grohs, “Probing dark photons in the early universe with big bang nucleosynthesis,” JCAP 12 (2020) 049, arXiv:2009.14325 [astro-ph.CO].
  * (49) XENON Collaboration, E. Aprile et al., “Dark Matter Search Results from a One Ton-Year Exposure of XENON1T,” Phys. Rev. Lett. 121 no. 11, (2018) 111302, arXiv:1805.12562 [astro-ph.CO].
  * (50) PandaX-II Collaboration, X. Ren et al., “Constraining Dark Matter Models with a Light Mediator at the PandaX-II Experiment,” Phys. Rev. Lett. 121 no. 2, (2018) 021304, arXiv:1802.06912 [hep-ph].
  * (51) J. Billard et al., “Direct Detection of Dark Matter – APPEC Committee Report,” arXiv:2104.07634 [hep-ex].
  * (52) H. W. Hammer, S. König, and U. van Kolck, “Nuclear effective field theory: status and perspectives,” Rev. Mod. Phys. 92 no. 2, (2020) 025004, arXiv:1906.12122 [nucl-th].
  * (53) A. Nicholson et al., “Toward a resolution of the NN controversy,” PoS LATTICE2021 (2022) 098, arXiv:2112.04569 [hep-lat].
  * (54) M. T. Frandsen, S. Sarkar, and K. Schmidt-Hoberg, “Light asymmetric dark matter from new strong dynamics,” Phys. Rev. D 84 (2011) 051703, arXiv:1103.4350 [hep-ph].
  * (55) C. Murgui and K. M. Zurek, “Dark Unification: a UV-complete Theory of Asymmetric Dark Matter,” arXiv:2112.08374 [hep-ph].


* [44] E. Lubetzky and Y. Zhao, _On the variational problem for upper tails in sparse random graphs_ , Random Structures & Algorithms 50 (2017), no. 3, 420–436.
  * [45] R.J. McCann, _A convexity principle for interacting gases_ , Advances in mathematics 128 (1997), no. 1, 153–179.
  * [46] S. Mukherjee and S. Sen, _Variational inference in high-dimensional linear regression_ , arXiv preprint arXiv:2104.12232 (2021).
  * [47] R. Oliveira and G. Reis, _Interacting diffusions on random graphs with diverging average degrees: Hydrodynamics and large deviations_ , Journal of Statistical Physics 176 (2019), no. 5, 1057–1087.
  * [48] F. Otto and C. Villani, _Generalization of an inequality by talagrand and links with the logarithmic sobolev inequality_ , Journal of Functional Analysis 173 (2000), no. 2, 361–400.
  * [49] K. Ray and B. Szabó, _Variational bayes for high-dimensional linear regression with sparse priors_ , Journal of the American Statistical Association (2021), 1–12.
  * [50] A. Saumard and J.A. Wellner, _Log-concavity and strong log-concavity: a review_ , Statistics surveys 8 (2014), 45.
  * [51] C. Villani, _Topics in optimal transportation_ , no. 58, American Mathematical Soc., 2003.
  * [52] M. J. Wainwright and M. I. Jordan, _Graphical models, exponential families, and variational inference_ , Foundations and Trends in Machine Learning 1 (2008), no. 1–2, 1–305.
  * [53] J. Yan, _Nonlinear large deviations: Beyond the hypercube_ , The Annals of Applied Probability 30 (2020), no. 2, 812–846.
  * [54] A. Y. Zhang and H. H. Zhou, _Theoretical and computational guarantees of mean field variational inference for community detection_ , The Annals of Statistics 48 (2020), no. 5, 2575–2598.