# A Bi-level Globalization Strategy for Non-convex Consensus ADMM and ALADIN

Xu Du*, Jingzhe Wang, Xiaohua Zhou and Yijie Mao *Corresponding authorXu Du,
is with the Institute of Mathematics HNAS, Henan Academy of Science,
Zhengzhou, China. (e-mail: duxu@hnas.ac.cn).Jingzhe Wang is with the
Department of Informatics & Networked Systems, School of Computing and
Information, University of Pittsburgh, Pittsburgh, PA, USA. (e-mail:
jiw148@pitt.edu).Xiaohua Zhou and Yijie Mao are with the School of Information
Science and Technology, ShanghaiTech University, China. (e-mail:
{zhouxh3,maoyj}@shanghaitech.edu.cn).

###### Abstract

In this paper, we formally analyze global convergence in the realm of
distributed consensus optimization. Current solutions have explored such
analysis, particularly focusing on consensus alternating direction method of
multipliers (C-ADMM), including convex and non-convex cases. While such
efforts on non-convexity offer elegant theory guaranteeing global convergence,
they entail strong assumptions and complicated proof techniques that are
increasingly pose challenges when adopted to real-world applications. To
resolve such tension, we propose a novel bi-level globalization strategy that
not only guarantees global convergence but also provides succinct proofs, all
while requiring mild assumptions. We begin by adopting such a strategy to
perform global convergence analysis for the non-convex cases in C-ADMM. Then,
we employ our proposed strategy in consensus augmented Lagrangian based
alternating direction inexact Newton method (C-ALADIN), a more recent and
generalization of C-ADMM. Surprisingly, our analysis shows that C-ALADIN
globally converges to local optimizer, complementary to the prior work on
C-ALADIN, which had primarily focused on analyzing local convergence for non-
convex cases.

###### Index Terms:

Distributed Consensus Optimization, Non-convex, Globalization, C-ADMM,
C-ALADIN

##  I Introduction

Alternating direction method of multipliers (ADMM) is the most well-known
algorithm within the realm of distributed optimization. The primary idea of
ADMM is to decompose a complex optimization problem into simpler sub-problems,
each can be solved in parallel by different agents. Since its first proposal
in [3, 4], there has been much effort to develop ADMM.  Such efforts have
mainly focused on the following two lines of research. The first line of work,
including but not limited to [5, 1], has explored potential applications that
could harness the advantages of ADMM, such as machine learning [6], signal
processing [7], wireless communication [8]. To enrich the theoretical ground
of ADMM, the second line of work, represented by [9, 10, 11, 12], has turned
its attention to revisiting the structure of ADMM. In particular, the most
influential work is that of He and Yuan [10], which first analyzed the
convergence rate of ADMM for general convex cases and investigated the sub-
linear convergence result. Later, Shi et al. extended this trajectory by
introducing the assumption of strongly convex [11], leading to a linear
convergence rate of consensus ADMM (C-ADMM). The aforementioned advancements,
however, have been exclusively to deal with convex problems. The application
of for non-convex consensus problems remains limited.

In tackling the convergence challenge [13, 14] caused by non-convexity,
several authors, including Hong and others [9, 15, 16, 17], have suggested
enforcing global convergence of ADMM with sufficiently large parameter of the
augmented term. As a special case, Wang et al. proposed a variant of ADMM,
known as Bregman ADMM [16] (an extension of proximal point method). Later,
Wang et al. [18] explored some non-convex non-smooth special cases that can be
effectively managed by ADMM. All the aforementioned methods are special
instances of the _proximal ADMM_ algorithm proposed in [17]. We refer [19] as
a comprehensive summary of ADMM variations to interested readers. We are aware
of that most of the existing works of non-convex ADMM rely on the following
two different methods: a) showing monotonic decreasing of a critical Lyapunov
function [18, 17], b) showing a corresponding (augmented) Lagrangian function
[9] converging to a limit point. Though applicable to C-ADMM, these two
approaches inevitably entail quite complicated proofs. Specifically, the
former method requires positive definiteness of the Lyapunov function,
constructed hard in nature, while the latter requires the smoothness of the
objective function, thereby narrowing the scope of applicability of C-ADMM.

Fortunately, in pursing concise proofs, Houska et al. proposed the notion of
T-ALADIN (typical augmented Lagrangian based alternating direction inexact
Newton method) [13], a solution that provides theoretical local convergence
guarantees when tackling general distributed non-convex optimization problems,
assuming a good initial point is available. It was not until seven years later
that Du and Wang realized that T-ALADIN might have a more efficient form of
implementation [2] (C-ALADIN, consensus ALADIN) for distributed consensus
problem in the ALADIN framework. While the approaches showcased promising
results, their guarantees were limited to local convergence.

In this paper, depart from the conventional _Gauss-Seidel decomposition_ ,
instead, we delve into the parallelizable implementation of C-ADMM and
C-ALADIN, based on _Jacobian decomposition_.Inspired by the globalization
approach proposed in [13, Section 6], we initiate the study of bi-level
globalization, which ensures that C-ADMM and C-ALADIN globally converge to
local optimizer. In summary, this paper has the following key contributions:

a) A novel _bi-level 111In fact, there are several papers also works on bi-
level optimization structure, see [20, 21, 22, 23, 24], however, that have
different purpose compare with this paper. proximal term_ is proposed to
handle the non-convexity in the classical distributed consensus optimization
optimization problem. Two distributed optimization algorithms, named CADMM-
Prox and CALADIN-Prox are then proposed for combining bi-level globalization
with C-ADMM and C-ALADIN, respectively. To our best knowledge, this is the
first work that applies bi-level globalization to C-ADMM and C-ALADIN. It’s
worth noting the proposed bi-level globalization approach is promising as it
ensures the strong convexity of the lower level problem and is guaranteed to
converge to the (local) optimal solution of the original non-convex problem.
This contrasts with the globalization employed in [13, Section 6], which is
used to tuning the step size only and only local optimal solution is
guaranteed.

b) Surprisingly, we obtain a simple and straightforward proof procedure thanks
to the use of globalization.

The rest of this paper is organized as follows. In Section II, we review the
preliminaries regarding globalization strategy. In Section III, with the
$L_{2}$-norm proximal term and the $L_{1}$-norm merit function, we propose our
globalization strategy for C-ALADIN and C-ADMM. In Section IV, we establish
the global convergence theory. In the end, Section V concludes this paper.

##  II Preliminaries

A distributed consensus optimization problem222In the rest of this paper, we
assume that, for all local minimizers of Problem (1), the second order
sufficient condition (SOSC) and linear independence constraint qualification
(LICQ) are satisfied. can be summarized as follows

$\begin{split}\min_{x,y}\;\;&\mathop{\sum}_{i=1}^{N}f_{i}(x_{i})\\\
\mathrm{s.t.}\;\;\;\;&\;x_{i}=y\;|\;\lambda_{i},\quad\forall
i\in\\{1,2,\dots,N\\},\end{split}$ (1)

where $f_{i}$s are potentially non-smooth, non-convex, and bounded from below.
Here, $x_{i}\in\mathbb{R}^{n}$ denotes the local variable of agent $i$,
$x\in\mathbb{R}^{N\cdot n}$ collects $x_{i}$s as a vector,
$y\in\mathbb{R}^{n}$ represents the global variable, and $\lambda_{i}$s
indicate the corresponding Lagrangian multiplier of agent $i$.

###  II-A An overview of C-ADMM and C-ALADIN

The augmented Lagrangian function of Problem (1) can be formalized as follows

$\begin{split}\mathscr{L}_{\rho}(x,y,\lambda)\overset{\text{def}}{=}&\;\mathop{\sum}_{i=1}^{N}f_{i}(x_{i})\\\
&+\mathop{\sum}_{i=1}^{N}\lambda_{i}^{\top}\left(x_{i}-y\right)+\frac{\rho}{2}\mathop{\sum}_{i=1}^{N}\|x_{i}-y\|^{2}.\end{split}$
(2)

Based on (2), C-ADMM can be then represented as (3)

$\left\\{\begin{split}x_{i}^{+}=&\;\mathop{\arg\min}_{x_{i}}f_{i}(x_{i})+\lambda_{i}^{\top}\left(x_{i}-y\right)+\frac{\rho}{2}\|x_{i}-y\|^{2};\\\
y^{+}=&\;\frac{1}{N}\mathop{\sum}_{i=1}^{N}\left(x_{i}^{+}+\frac{1}{\rho}\lambda_{i}\right);\\\
\lambda_{i}^{+}=&\;\lambda_{i}+\rho\left(x_{i}^{+}-y^{+}\right).\\\
\end{split}\right.$ (3)

If $f_{i}$s are convex, (3) is guaranteed to converge globally within finite
iterations [1].

By combining C-ADMM with sequential quadratic programming (SQP) method333 The
main difference between (3) and (4) is that C-ALADIN updates the dual
$\lambda$ together with the global variable $y$, while C-ADMM updates the all
the variables separately., we get C-ALADIN [2], as per (4). It only guarantees
local convergence for non-convex problems when starting from a good initial
guess.

$\left\\{\begin{split}&x_{i}^{+}=\mathop{\arg\min}_{x_{i}}f_{i}(x_{i})+\lambda_{i}^{\top}\left(x_{i}-y\right)+\frac{\rho}{2}\|x_{i}-y\|^{2};\\\
&\left(y^{+},\lambda^{+}\right)\\\
=&\left\\{\begin{split}\vspace{-1cm}&\mathop{\mathrm{\arg}\mathrm{\min}}_{\Delta\bar{x},y}\;\mathop{\sum}_{i=1}^{N}\left(\frac{1}{2}\Delta\bar{x}_{i}^{\top}B_{i}\Delta\bar{x}_{i}+\partial
f_{i}(x_{i}^{+})^{\top}\Delta\bar{x}_{i}\right)\\\
&\quad\mathrm{s.t.}\qquad\Delta\bar{x}_{i}+x_{i}^{+}=y\;|\lambda_{i}\\\
\end{split}\right\\}.\\\ \end{split}\right.$ (4)

Here $B_{i}\succ 0$ and $\partial f_{i}$ denote the local Hessian
approximation and (sub)gradient of $f_{i}$, respectively.

###### Remark 1

It is worth noting that, in solving large-scale optimization problems, C-ADMM
is a _first-order_ algorithm that is computationally tractable. However, it
lacks theoretical guarantees for non-convex problems. Conversely, C-ALADIN
does provide rigorous convergence analysis, but it only has local convergence
guarantees instead of global ones. Our proposed globalization strategy
(detailed in IV) solved the two challenges.

###  II-B The $L_{2}$-norm Proximal Term $\&$ the $L_{1}$-norm Merit Function

Since $f_{i}$s are in general non-convex, we introduce the $L_{2}$-norm
proximal term [25, 26] as follows

$\frac{\gamma}{2}\|x_{i}-z\|^{2}.$ (5)

By incorporating (5) with an upper global variable $z$ and a sufficiently
large parameter $\gamma>0$, we guarantee the strong convexity of $F_{i}^{{z}}$
with respect to $x_{i}$, formulated in Equation (6)

$F_{i}^{z}(x_{i})\overset{\text{def}}{=}f_{i}(x_{i})+\frac{\gamma}{2}\|x_{i}-z\|^{2}.$
(6)

Instead of solving (1) directly, in the following sections, two optimization
algorithms are proposed for (7) based on bi-level optimization. Specifically,
the update of $x$ and $y$ will be performed in the _lower level_ while the
update of $z$ will be done in the _upper level_. In the lower level, with a
fixed $z$, we focus on searching the optimizer of (7).

$\begin{split}\min_{x,y}\;\;&\mathop{\sum}_{i=1}^{N}F^{z}_{i}(x_{i})\\\
\mathrm{s.t.}\;\;&\;\;x_{i}=y\;|\;\lambda_{i},\quad\forall
i\in\\{1,2,\dots,N\\}.\end{split}$ (7)

Note that, $\gamma$ guarantees the strong convexity of the _lower level
problem_ (7), making it easier to be solved. Then, the update of $z$ – a key
ingredient within the _upper level_ – enforces the solution of Problem (7)
asymptotically reaching towards the (local) optimal solution of Problem (1)
eventually.

With the lower global variable $y$, the $L_{1}$-norm _merit function_ [27,
Chapter 18.3] of the sub-problem can be then expressed as Equation (8)

$\Phi_{i}^{(z,y)}(x_{i})\overset{\text{def}}{=}F^{z}_{i}(x_{i})+\sigma_{i}\|x_{i}-y\|_{1},$
(8)

that aims to measure the performance of the objective and the equality
constraints violation of each sub-problem. Here $\sigma_{i}$ is sufficiently
large, such that, $\sigma_{i}\geq\|\lambda_{i}\|_{\infty}\;\forall i$ [27].
Obviously, with Equation (8), the $L_{1}$-norm merit function of problem (7)
can be described as (9)

$\Phi^{(z,y)}(x)\overset{\text{def}}{=}\;\sum_{i=1}^{N}F^{z}_{i}(x_{i})+\sum_{i=1}^{N}\sigma_{i}\|x_{i}-y\|_{1}.$
(9)

In the following sections, _merit function_ only relates to Equation (9). Note
that,

$\sum_{i=1}^{N}f_{i}(z)=\Phi^{(z,z)}(z)$ (10)

is an important property that helps with establishing the global convergence
theory through this paper. The analysis in the following sections depends on
the update of the triple $(x,y,z)$. Formally, we define

$\left\\{\begin{split}\Delta x_{i}&\overset{\text{def}}{=}x_{i}^{+}-y;\\\
\Delta\tilde{x}_{i}&\overset{\text{def}}{=}y^{+}-x_{i}^{+}.\end{split}\right.$
(11)

Here $(\cdot)^{+}$ denotes the update of corresponding decision variables
while $(\cdot)^{-}$ represents such variables from the latest iteration.

##  III Proposed Optimization Frameworks

In this section, we take advantage of our proposed globalization strategy and
propose two novel algorithms to solve Equation (1), namely CADMM-Prox and
CALADIN-Prox. Specifically, in Subsection (III-A), we first propose the
updating strategy of the local variable $x_{i}$ and global variable $y$. Then,
we detail the entire optimization frameworks in Subsection (III-B).

###  III-A Conservative Update of the Lower Level

The lower level of C-ALADIN and C-ADMM has two main steps: a) local nonlinear
programming (NLP) minimization of $x_{i}$ (Section III-A1), b) consensus step
(Section III-A2).

####  III-A1 Local Primal Variable Update

With the proximal term (5), we suggest that updating the local primal variable
$x_{i}$ through approximated reformulation. Concretely, for C-ALADIN, we
linearize the upper objective $F^{z}_{i}(x_{i})$, such that

$\begin{split}{x_{i}}^{+}=&\mathop{\arg\min}_{x_{i}}F^{z}_{i}(x_{i})+\lambda_{i}^{\top}x_{i}+\frac{1}{2}\|x_{i}-y\|^{2}_{B_{i}}\\\
\approx&\mathop{\arg\min}_{x_{i}}\left(\partial
f_{i}(y)+\gamma(y-z)\right)^{\top}(x_{i}-y)\\\
&+\lambda_{i}^{\top}x_{i}+\frac{1}{2}\|x_{i}-y\|^{2}_{B_{i}}\\\
=&y+B_{i}^{-1}\left(\gamma(z-y)-\lambda_{i}-\partial
f_{i}(y)\right)\end{split}$ (12)

with $B_{i}\succ 0$. Here $\|\cdot\|_{B_{i}}$ denotes the Mahalanobis distance
[16] with symmetric positive definite matrix $B_{i}$.

For C-ADMM, by setting $B_{i}=\rho I$, Equation (12) degrades into (13)

$\begin{split}{x_{i}}^{+}=y+\frac{1}{\rho}\left(\gamma(z-y)-\lambda_{i}-\partial
f_{i}(y)\right).\end{split}$ (13)

Indeed, alternative approaches to updating the local variables may also work.
We detail this part in Appendix A.

####  III-A2 Lower Level Global Variable Update (Consensus Step)

In the C-ALADIN framework, if we have the local Hessian approximation matrices
$B_{i}\approx\nabla^{2}F^{z}(x_{i}^{+})\succ 0$, the previous lower global
variable $y^{-}$, $\beta>0$, and the local (sub)gradient $g_{i}$ formalized as
follows

$\begin{split}g_{i}&\overset{\text{def}}{=}\nabla F^{z}(x_{i}^{+})\\\
&=\partial f_{i}(x_{i}^{+})+\gamma(x_{i}^{+}-z)\\\ &\approx
B_{i}(y-x_{i}^{+})-\lambda_{i},\end{split}$ (14)

we can then derive the update of the lower global variable $y^{\text{ALADIN}}$
in Equation (15)

$\displaystyle\left(y^{\text{ALADIN}}\right)^{+}$ (15) $\displaystyle=$
$\displaystyle\left\\{\begin{split}&\mathop{\mathrm{\arg}\mathrm{\min}}_{\Delta\bar{x},y}\quad\mathop{\sum}_{i=1}^{N}\left(\frac{1}{2}\Delta\bar{x}_{i}^{\top}B_{i}\Delta\bar{x}_{i}+g_{i}^{\top}\Delta\bar{x}_{i}\right)\\\
&\qquad\qquad+\frac{\beta}{2}\|y-y^{-}\|^{2}\\\
&\quad\mathrm{s.t.}\qquad\Delta\bar{x}_{i}+x_{i}^{+}=y\;|\;\lambda_{i}\\\
\end{split}\right\\}.$

Note that, Equation (15) has a close form as follows

$\begin{split}&\left(y^{\text{ALADIN}}\right)^{+}\\\
=&\left(\sum_{i=1}^{N}B_{i}+\beta I\right)^{-1}\left(\beta
y^{-}+\sum_{i=1}^{N}\left(B_{i}x_{i}^{+}-g_{i}\right)\right).\end{split}$ (16)

By setting $B_{i}=\rho I$ and $\lambda_{i}=-g_{i}$, Equation (16) boils down
to (17) that represents the lower global variable update of C-ADMM.

$\begin{split}\left(y^{\text{ADMM}}\right)^{+}=\frac{1}{N\rho+\beta}\left(\beta
y^{-}+\sum_{i=1}^{N}\left(\rho
x_{i}^{+}+\lambda_{i}\right)\right).\end{split}$ (17)

###  III-B Proposed Algorithms

####  III-B1 CALADIN-Prox

With the approximate update of local variable $x_{i}$ (Equation (12)) and
global variable $y$ (Equation (15)), the linearized form of C-ALADIN, equipped
with globalization strategy, can be then expressed as Algorithm 1.

Algorithm 1 CALADIN-Prox

Initialization: Initial guess $\lambda_{i}$, choose upper global variable $z$,
set lower global variable $y=z$, local Hessian approximation $B_{i}\succ 0$.
Set $\gamma>0,\rho>0,\beta\geq 0,\sigma_{i}=\|\lambda_{i}\|_{\infty}$.
Repeat:

  1. 1.

Inexactly updates the local variable $x_{i}$ with Equation (12) or (33) where
$\lambda_{i}=B_{i}(x_{i}^{-}-y)-g_{i}$.

  2. 2.

Evaluate the new sub-gradient

$\begin{split}g_{i}=B_{i}(y-x_{i}^{+})-\lambda_{i}.\end{split}$

from each ${x_{i}}^{+}$ and upload to the master.

  3. 3.

Update the lower global variable $y$:

$\begin{split}y^{+}=&\left(\sum_{i=1}^{N}B_{i}+\beta I\right)^{-1}\left(\beta
y+\sum_{i=1}^{N}\left(B_{i}x_{i}^{+}-g_{i}\right)\right).\end{split}$

  4. 4.

Set $\sigma_{i}^{+}=\|B_{i}(x_{i}^{+}-y^{+})-g_{i}\|_{\infty}$ if

$\sigma_{i}<\|B_{i}(x_{i}^{+}-y^{+})-g_{i}\|_{\infty}.$

  5. 5.

Update the upper global variable $z^{+}=y^{+}$ if

$\Phi^{(z,y^{+})}(y^{+})<\Phi^{(z,y)}(y).$

  6. 6.

Download the updated variables:

$\left\\{\begin{array}[]{l}\begin{split}y&\leftarrow y^{+};\\\
0&\;\text{or}\;1\qquad\text{1 for updated $z$ and 0 for no reaction
{}}.\end{split}\end{array}\right.$

###### Remark 2

The numerical value of $\left(\sum_{i=1}^{N}B_{i}+\beta I\right)^{-1}$ can be
prepared for acceleration. As discussed in [2], when setting $B_{i}=\rho I$,
C-ALADIN is still different from C-ADMM in updating the dual variables.

####  III-B2 CADMM-Prox

Similar to Algorithm 1, C-ADMM with globalization strategy can be summarized
as Algorithm 2, supported by Equation (13) and (17).

###### Remark 3

A proximal term has also been introduced in _Proximal ADMM_ [17] as
$\frac{\gamma}{2}\|x_{i}-y\|^{2}$, however, only a single-level optimization
structure is considered in that work. The advantage of the bi-level structure
in Algorithm 2 is that the upper level global variable $z$ in
$\frac{\gamma}{2}\|x_{i}-z\|^{2}$ is not updated per iteration in the lower
level, which gives the strong convexity of (7). In fact, for single-level
ADMM, the strong convexity of the augmented sub-problem in Equation (3) is not
a sufficient condition for ADMM to converge. See a contrary example [13,
Example 2.1].

##  IV Convergence Analysis

In this section, we establish the global convergence theory of Algorithm 1 and
2, which is the major contribution of this paper. In Subsection (IV-A), we
begin with analyzing the convergence behavior of the lower level of Algorithms
1 and 2. Later, in Subsection (IV-B), we establish the global convergence
theory with our proposed globalization strategy in the sense of upper level.
The combination of the two finally yields the global convergence.

Algorithm 2 CADMM-Prox

Initialization:Initial the upper and lower global parameter $z,y$ and the dual
$\lambda_{i}$. Set $\gamma>0,\rho>0,\beta\geq
0,\sigma_{i}=\|\lambda_{i}\|_{\infty}$.
Repeat:

  1. 1.

Evaluate $\lambda_{i}^{+}=\lambda_{i}+\rho(x_{i}-y)$, then inexactly update
the local variable as $x_{i}^{+}$ with Equation (13) or (35) on the agent side
and upload to the master.

  2. 2.

Update the lower global variable $y$:

$\begin{split}y^{+}&=\frac{1}{N\rho+\beta}\left(\beta
y+\sum_{i=1}^{N}\left(\rho x_{i}^{+}+\lambda_{i}\right)\right).\end{split}$

  3. 3.

Set $\sigma_{i}^{+}=\|\lambda_{i}+\rho(x_{i}^{+}-y^{+})\|_{\infty}$ if

$\sigma_{i}<\|\lambda_{i}+\rho(x_{i}^{+}-y^{+})\|_{\infty}.$

  4. 4.

Update the upper global variable $z=y^{+}$ if

$\Phi^{(z,y^{+})}(y^{+})<\Phi^{(z,y)}(y).$

  5. 5.

Broadcast the updated variables:

$\left\\{\begin{array}[]{l}\begin{split}y&\leftarrow y^{+}\\\
0&\;\text{or}\;1\qquad\text{1 for updated $z$ and 0 for no
reaction.}\end{split}\end{array}\right.$

###  IV-A Lower Level Merit Function Decreasing Guarantees

Given a nonlinear mapping $\Lambda(\xi):\mathbb{R}^{n}\rightarrow\mathbb{R}$,
and the increment of the decision variable $p\in\mathbb{R}^{n}$, we have the
concept of _directional derivative_ , defined as (18)

$\text{D}\Lambda(\xi)[p]\overset{\text{def}}{=}\lim_{t\rightarrow
0^{+}}\frac{\Lambda(\xi+tp)-\Lambda(\xi)}{t}.$ (18)

Before providing the decent condition of the merit function (9), we first give
the following lemma.

###### Lemma 1

With the definition of $\text{D}\Lambda(\xi)[p]$ in Equation (18), assume that
$\Delta\tilde{x}$ is generated by the QP iteration from (15), then the
directional derivative of Equation (9) is represented as

$\begin{split}&\text{D}\Phi^{(z,y^{+})}(x^{+})[\Delta\tilde{x}]\\\
=&\sum_{i=1}^{N}\nabla
F_{i}^{z}(x_{i}^{+})^{\top}\Delta\tilde{x}_{i}-\sum_{i=1}^{N}\sigma_{i}\|\Delta\tilde{x}_{i}\|_{1}+O(t).\end{split}$
(19)

Here, $y^{+}$ can be pre-computed by following Equation (16) in Algorithm 1 or
(17) in Algorithm 2.

###### Proof 1

See Appendix B. $\blacksquare$

Moreover, for the $L_{1}$ norm gap $\|x-y\|_{1}$, the directional derivative
with respect to $x$ at point $y$, in direction $p$, can be denoted as follows

$\text{D}\|x-y\|_{1}(y)[p]=\|p\|_{1}.$ (20)

With Equation (19) and (20), now we can show that the merit function (9) is
decreasing in the lower level of Algorithm 1 and 2 , detailed in Theorem 1 and
2. Note that, the following analysis mainly relates to the increment with
respect to the local variables $x_{i}$s introduced in Equation (11).

####  IV-A1 Merit Function Analysis at the Local NLP Update Step

The following theorem guarantees the decrement of the merit function with the
local primal $x_{i}$ update for CALADIN-Prox.

###### Theorem 1

With the local update (12), the directional derivative of the lower merit
function (9) with respect to the local NLP (12) has the following property:

$\text{D}\Phi^{(z,y)}(y)[\Delta x]<0$ (21)

will hold if

$\begin{split}\sum_{i=1}^{N}\lambda_{i}^{\top}\Delta
x_{i}+\sum_{i=1}^{N}\|\Delta
x_{i}\|^{2}_{B_{i}}>\sum_{i=1}^{N}\sigma_{i}\|\Delta x_{i}\|_{1}.\end{split}$
(22)

###### Proof 2

See Appendix C. $\blacksquare$

Since $x_{i}^{+}=y+\Delta x_{i}$, Theorem 1 implies

$\Phi^{(z,y)}(x_{i}^{+})<\Phi^{(z,y)}(y).$ (23)

By setting $B_{i}=\rho I$ and updating with Equation (13), the above statement
also works for CADMM-Prox.

Note that, the decreasing analysis of the merit function on local NLP update,
described by Equation (32), (33), (34) and (35), are similar to Theorem 1. We
thus do not show detail here.

####  IV-A2 Merit Function Analysis at the Lower Level Global Variable Update
Step

The following theorem provides the descent condition of the merit function
(9). Such a condition can evaluate the results of the consensus step in
Algorithm 1.

###### Theorem 2

With Lemma 1, the directional derivative of the lower merit function (9) with
respect to the consensus QP (15) has the following property:

$\begin{split}&\text{D}\Phi^{(z,y^{+})}(x^{+})[{\Delta\tilde{x}}]<0,\\\
\end{split}$ (24)

if

$\sigma_{i}>\|\lambda_{i}^{+}\|_{\infty}.$ (25)

###### Proof 3

See Appendix D. $\blacksquare$

After updating the lower global variable $y$, the linear consensus coupling of
Equation (15) holds, such that $x_{i}^{+}+\Delta\tilde{x_{i}}=y^{+}$.
Therefore,

$\Phi^{(z,y^{+})}(y^{+})<\Phi^{(z,y)}(x^{+})$ (26)

is guaranteed with the Theorem 2. The above theorem shows that Equation (15)
will give a decent direction of the lower merit function. Similar to the above
discussion, by replacing $B_{i}$ with $\rho I$ and applying Appendix D again,
(26) also works for CADMM-Prox.

By combining Theorem 1 and 2, the merit function is guaranteed to decrease
with respect to the lower level of Algorithm 1 and 2 in only one step, such
that

$\Phi^{(z,y^{+})}(y^{+})<\Phi^{(z,z)}(z).$ (27)

###### Remark 4

If one of (22) and (25) fail to meet, after one step of the lower level, the
update of the upper global variable $z$ of Algorithm 1 and 2 can not be
guaranteed. Fortunately, with the global convergence theory of convex C-ADMM
[11], Equation (27) still holds within finite steps for the lower level of
Algorithm 2.

For convex C-ALADIN consisting of a special case such that $B_{i}=\rho I$, [2]
shows the linear convergence of strongly convex problems. The following
statement will show that C-ALADIN can also globally converge within finite
iteration for convex problems in a more general setting. The strongly convex
distributed scenario relates to the lower level of Algorithm 1 such that (27)
also holds.

###### Theorem 3

Assume $(y^{*},\lambda^{*})$ denotes the KKT (Karush-Kuhn-Tucker) point of
Problem (7), for general local Hessian approximation $B_{i}\succ 0$, the
Lyapunov function (28)

$\mathcal{L}(y,\lambda)=\sum_{i=1}^{N}\|y-y^{*}\|_{B_{i}}^{2}+\sum_{i=1}^{N}\|\lambda_{i}-\lambda_{i}^{*}\|_{B_{i}^{-1}}^{2}$
(28)

is monotonely decreasing with the lower level of Algorithm 1, such that

$\mathcal{L}(y^{+},\lambda^{+})-\mathcal{L}(y,\lambda)\leq 0.$

###### Proof 4

See Appendix E. $\blacksquare$

Since Equation (28) is positive definite, it converges to $0$ in finite
iterations by adopting C-ALADIN. For strongly convex problems, we can further
show the linear rate of C-ALADIN by applying [2, Theorem 7] again.

###  IV-B Enforcing Global Convergence to Local Optimizer

With the discussion of Equation (27) in Subsection (IV-A), now the global
convergence theory of Algorithm 1 and 2 can be established with the following
theorem.

###### Theorem 4

With the lower level of Algorithm 1 and 2, the decrement of
$\Phi^{(\cdot,\cdot)}(\cdot)$ always holds, such that

$\Phi^{(z,z)}(z)-\Phi^{(y^{+},y^{+})}(y^{+})>\frac{\gamma
N}{2}\|y^{+}-z\|^{2}>0.$ (29)

Here, $y^{+}$ is a strict global minimizer of Problem (7) generated by the
lower level of Algorithm 1 or 2 with ${z}$ ($y^{+}\neq{z}$) as an initial
point.

###### Proof 5

See Appendix F. $\blacksquare$

By setting $z^{+}=y^{+}$, Equation (29) implies

$\Phi^{(z^{k},z^{k})}(z^{k})<\Phi^{(z^{0},z^{0})}(z^{0})-\frac{\gamma
N}{2}\sum_{j=0}^{k-1}\|z^{j}-z^{j+1}\|^{2},$ (30)

where $k$ is the iteration index of the upper level. Since
$\sum_{i=1}^{N}f_{i}(x_{i})$ is bounded below, with the proposed globalization
routine, the following limit always exists

$\Phi^{(z^{*},z^{*})}(z^{*})=\lim_{k\rightarrow\infty}\Phi^{(z^{k},z^{k})}(z^{k})>-\infty.$

Here $z^{*}$ donates the limit point of $\Phi^{(\cdot,\cdot)}(\cdot)$. By
taking the limit of Equation (30) on both sides, we have

$\Phi^{(z^{*},z^{*})}(z^{*})<\Phi^{(z^{0},z^{0})}(z^{0})-\frac{\gamma
N}{2}\sum_{j=0}^{\infty}\|z^{j}-z^{j+1}\|^{2}$

which implies

$\begin{split}&\sum_{j=0}^{\infty}\|z^{j}-z^{j+1}\|^{2}\\\ <&\frac{2}{\gamma
N}\left(\Phi^{(z^{0},z^{0})}(z^{0})-\Phi^{(z^{*},z^{*})}(z^{*})\right).\end{split}$
(31)

Notice that, the left hand side of Inequality (31) must be bounded such that

$\lim_{j\rightarrow\infty}\|z^{j}-z^{j+1}\|^{2}=0.$

With the above analysis,

$\|z^{j}-z^{j+1}\|^{2}\leq\epsilon$

holds for sufficiently small $\epsilon>0$ with finite index $j$ for the upper
level of Algorithm 1 and 2. Moreover, with the definition of merit function
(9), the limit point of $\Phi^{(\cdot,\cdot)}(\cdot)$ is, at least, a
_critical point_ of Problem (1). To test whether the critical point is a local
optima is trivial. By applying a perturbation
$\varsigma\in\mathbb{R}^{n}\rightarrow 0$ to $z^{*}$, serving as a new initial
point, Algorithm 1 and 2 restart iterating until they converge to $z^{*}$
again, which means that $z^{*}$ is a _local minimizer_. Otherwise, a new
_critical point_ will be found, which indicates that $z^{*}$ is a saddle
point. Note that the above global convergence analysis is also applicable to
general ADMM and T-ALADIN which is not shown here.

Note that, for convex cases, the proposed globalization strategy is useless.
For non-convex cases, the proposed globalization strategy can not guarantee to
find a global minimizer of the non-convex distributed problem (1) but can, at
least, find a local optimizer.

###### Remark 5

In Algorithm 1 and 2, due to the strongly convexity of Problem (7), there
exists a parameter $\delta(z^{k})>0$ (with fixed $z^{k}$) such that the lower
global variable $y$ converges to the global optimizer of Problem (7) with
linear rate $\left(\frac{1}{\sqrt{1+\delta(z^{k})}}\right)$ [28, 2]. Assume
that there exists a $\bar{\delta}\geq\delta(z^{k})\;\forall z^{k}$, then, with
our proposed bi-level globalization strategy, the upper global variable $z$
will, at least, globally converge to a local optimizer of Problem (1) with
linear convergence rate $\left(\frac{1}{\sqrt{1+\bar{\delta}}}\right)$.

###### Remark 6

If we skip the local update of Algorithm 1 and update the lower global
variable $y$ in a decentralized fashion, then it relates to another algorithm,
called a consensus version of _decentralized bi-level SQP_ [22]. In this
mentioned bi-level structure, our proposed globalization strategy is redundant
for convex cases. However, the proposed globalization strategy still helps to
establish the global convergence theory for non-convex cases.

##  V Conclusion

In this paper, a novel bi-level globalization strategy of C-ADMM and C-ALADIN
for non-convex problems is proposed. Without the help of the Lyapunov function
or the corresponding Lagrangian function, we establish the global convergence
theory for non-convex problems based on the state-of-the-art convex analysis
of C-ADMM and C-ALADIN.Future work will concentrate on non-convex distributed
consensus optimization with constrained sub-problems. Decentralized variation
of the proposed globalization strategy will also be considered.

##  VI Acknowledgments

We are grateful to Boris Houska, Bingsheng He, Shijie Zhu, Alexander
Engelmann, Kai Wang, Yuning Jiang and Xiaojun Yuan for their helpful
discussion.

##  Appendix A Comments on Local Variable Update in Section III-A

Other approximation techniques can be also implemented here. For C-ALADIN, a
linearization of the lower objective $f_{i}(x_{i})$ is shown as Equation (32).

$\begin{split}{x_{i}}^{+}\approx&\mathop{\arg\min}_{x_{i}}\partial
f_{i}(y)^{\top}(x_{i}-y)+\frac{\gamma}{2}\|x_{i}-z\|^{2}\\\
&+\lambda_{i}^{\top}x_{i}+\frac{1}{2}\|x_{i}-y\|^{2}_{B_{i}}\\\
=&\left(B_{i}+\gamma I\right)^{-1}\left(\gamma z+B_{i}y-\lambda_{i}-\partial
f_{i}(y)\right).\end{split}$ (32)

The local primal variable update can be also designed as a multiple iteration
version, as shown in Equation (33) where $x_{i}$ can be updated within finite
iterations

$\left\\{\begin{split}x_{i}^{-}&=x_{i}^{+}\quad\text{(take $x_{i}^{-}=y$ as an
initial)};\\\ x_{i}^{+}&=\left(B_{i}+\gamma I\right)^{-1}\left(\gamma
z+B_{i}y-\lambda_{i}-\partial f_{i}(x_{i}^{-})\right).\end{split}\right.$ (33)

For C-ADMM, by setting $B_{i}=\rho I$ and with $\rho>0$, Equation (32) boils
down to (34) and has a similar expression as (35)

$\begin{split}{x_{i}}^{+}=\frac{1}{\gamma+\rho}\left(\gamma z+\rho
y-\lambda_{i}-\partial f_{i}(y)\right).\end{split}$ (34)

Similar as (33), (34) can be also replaced by (35)

$\left\\{\begin{split}x_{i}^{-}&=x_{i}^{+}\quad\text{(take $x_{i}^{-}=y$ as an
initial)};\\\ x_{i}^{+}&=\frac{1}{\gamma+\rho}\left(\gamma z+\rho
y-\lambda_{i}-\partial f_{i}(x_{i}^{-})\right)\end{split}\right.$ (35)

where $x_{i}$ is updated with finite iterations. A similar idea has been
proposed in [28] and shows a linear convergence rate of strongly convex ADMM.
Moreover, a quadratic approximation method can be also considered [29].

##  Appendix B Proof of Lemma 1

The following prove inherits from [27, Theorem 18.2]. Assume that the merit
function with a increment has the following linear representation

$\begin{split}&\Phi^{(z,y^{+})}(x^{+}+t\Delta\tilde{x})\\\
=&\sum_{i=1}^{N}F_{i}^{z}(x_{i}^{+}+t\Delta\tilde{x}_{i})+\sum_{i=1}^{N}\sigma_{i}\|x_{i}^{+}-y^{+}+t\Delta\tilde{x}_{i}\|_{1}\\\
=&\sum_{i=1}^{N}F^{z}_{i}(x_{i}^{+})+t\sum_{i=1}^{N}\nabla
F^{z}_{i}(x_{i}^{+})^{\top}\Delta\tilde{x}_{i}\\\
&+\sum_{i=1}^{N}\sigma_{i}\|x_{i}^{+}-y^{+}-t(x_{i}^{+}-y^{+})\|_{1}+O(t^{2}).\end{split}$
(36)

By isolating the first order block with respect to $t$, Equation (36) can be
reformulated as follows

$\begin{split}&\Phi^{(z,y^{+})}(x^{+}+t\Delta\tilde{x})\\\
=&t\left\\{\sum_{i=1}^{N}\nabla
F^{z}_{i}(x_{i}^{+})^{\top}\Delta\tilde{x}_{i}-\sum_{i=1}^{N}\sigma_{i}\|\Delta\tilde{x}_{i}\|_{1}\right\\}\\\
&+\underbrace{\left\\{\sum_{i=1}^{N}F^{z}_{i}(x_{i}^{+})+\sum_{i=1}^{N}\sigma_{i}\|x_{i}^{+}-y^{+}\|_{1}\right\\}}_{\Phi^{(z,y^{+})}(x^{+})}+O(t^{2}).\\\
\end{split}$ (37)

Then, from the definition of directional derivative (18), we have the
resultant derivation as

$\begin{split}&\text{D}\Phi^{(z,y^{+})}(x^{+})[\Delta\tilde{x}]\\\
=&\lim_{t\rightarrow
0^{+}}\frac{\Phi^{(z,y^{+})}(x^{+}+t\Delta\tilde{x})-\Phi^{(z,y^{+})}(x^{+})}{t}\\\
=&\sum_{i=1}^{N}\nabla
F^{z}_{i}(x_{i}^{+})^{\top}\Delta\tilde{x}_{i}-\sum_{i=1}^{N}\sigma_{i}\|\Delta\tilde{x}_{i}\|_{1}+O(t).\end{split}$
(38)

This completes the proof.

##  Appendix C Proof of Theorem 1

First, we have a first order Taylor expansion of Equation (12),

$\begin{split}&\mathop{\min}_{x_{i}}F^{z}_{i}(x_{i})+\lambda_{i}^{\top}x_{i}+\frac{1}{2}\|x_{i}-y\|^{2}_{B_{i}}\\\
\approx&\nabla F^{z}_{i}(y)^{\top}\Delta x_{i}+\lambda_{i}^{\top}(y+\Delta
x_{i})+\frac{1}{2}\|\Delta x_{i}\|^{2}_{B_{i}}.\end{split}$ (39)

Clearly, the first order optimality condition of Equation (39) shows as

$\begin{split}&\nabla F^{z}_{i}(y)+\lambda_{i}+B_{i}\Delta
x_{i}=0.\end{split}$ (40)

By multiplying $\Delta x_{i}$ to the left hand side of (40) and then
performing summation over the agent index, we get the following equality

$\begin{split}&\sum_{i=1}^{N}\nabla F^{z}_{i}(y)^{\top}\Delta
x_{i}=\sum_{i=1}^{N}\left(-\lambda_{i}^{\top}\Delta x_{i}-\|\Delta
x_{i}\|^{2}_{B_{i}}\right).\\\ \end{split}$ (41)

With Equation (20) and (41), if (22) holds, then (21) is guaranteed, which
completes our proof.

##  Appendix D Proof of Theorem 2

From the first order optimal condition of Equation (15) with respect to
$\Delta\tilde{x}_{i}$, the follow equality always holds.

$B_{i}\Delta\tilde{x}_{i}+\lambda_{i}^{+}+\nabla F^{z}_{i}(x_{i}^{+})=0.$ (42)

We then multiply $\Delta\tilde{x}_{i}$ to the left hand side and get

$\Delta\tilde{x}_{i}^{\top}B_{i}\Delta\tilde{x}_{i}+(\lambda_{i}^{+})^{\top}\Delta\tilde{x}_{i}+\nabla
F^{z}_{i}(x_{i}^{+})^{\top}\Delta\tilde{x}_{i}=0.$ (43)

This implies

$\begin{split}&\sum_{i=1}^{N}\nabla
F^{z}_{i}(x_{i}^{+})^{\top}\Delta\tilde{x}_{i}\\\
=&-\sum_{i=1}^{N}\left(\Delta\tilde{x}_{i}^{\top}B_{i}\Delta\tilde{x}_{i}+(\lambda_{i}^{+})^{\top}\Delta\tilde{x}_{i}\right)\\\
\leq&\sum_{i=1}^{N}\|\lambda_{i}^{+}\|_{\infty}\|\Delta\tilde{x}_{i}\|_{1}-\sum_{i=1}^{N}\|\Delta\tilde{x}_{i}\|_{B_{i}}^{2}.\end{split}$
(44)

With Equation (25) and $B_{i}\succ 0$, once we plug (44) into (19), Equation
(24) holds. This completes the proof.

##  Appendix E Proof of Theorem 3

The following proof is similar to that of [2, Appendix C].

We then show that, by proving a non-positive equation (45), Equation (28) is
monotonically decreasing. Specifically,

$4\sum_{i=1}^{N}\left(\lambda_{i}-\lambda_{i}^{*}+B_{i}(x_{i}^{+}-y)\right)^{\top}\left(x_{i}^{+}-y^{*}\right)$
(45)

is equivalent to

$\mathcal{L}(y^{+},\lambda^{+})-\mathcal{L}(y,\lambda).$

Equation (45) can also be expressed as

$\begin{split}&4\sum_{i=1}^{N}\left(\lambda_{i}-\lambda_{i}^{*}+B_{i}(x_{i}^{+}-y)\right)^{\top}\left(x_{i}^{+}-y^{*}\right)\\\
=&4\sum_{i=1}^{N}\left(\lambda_{i}-\lambda_{i}^{*}\right)^{\top}\left(x_{i}^{+}-y^{*}\right)\\\
&+4\sum_{i=1}^{N}\left(B_{i}(x_{i}^{+}-y)\right)^{\top}\left(x_{i}^{+}-y^{*}\right).\end{split}$
(46)

Notice that from the KKT optimality condition of Equation (7), we have

$\sum_{i=1}^{N}\lambda_{i}=0.$ (47)

By applying (47), Equation (46) becomes as

$\begin{split}&4\sum_{i=1}^{N}\left(\lambda_{i}-\lambda_{i}^{*}\right)^{\top}\left(x_{i}^{+}-y\right)\\\
&+4\sum_{i=1}^{N}\left(B_{i}(x_{i}^{+}-y)\right)^{\top}\left(x_{i}^{+}-y^{*}\right)\\\
=&4\sum_{i=1}^{N}\left(\lambda_{i}-\lambda_{i}^{*}+B_{i}\left(x_{i}^{+}-y^{*}\right)\right)^{\top}\left(x_{i}^{+}-y\right).\end{split}$
(48)

In the end, by introducing an important formula appeared in [2]

$x_{i}^{+}=\frac{B_{i}^{-1}}{2}\left(\lambda_{i}^{+}-\lambda_{i}\right)+\frac{1}{2}\left(y+y^{+}\right),$

we simplify Equation (48) as

$\begin{split}&4\sum_{i=1}^{N}\left(\lambda_{i}-\lambda_{i}^{*}+B_{i}\left(x_{i}^{+}-y^{*}\right)\right)^{\top}\left(x_{i}^{+}-y\right)\\\
=&\sum_{i=1}^{N}\|\lambda_{i}^{+}-\lambda_{i}^{*}\|^{2}_{B_{i}^{-1}}-\sum_{i=1}^{N}\|\lambda_{i}-\lambda_{i}^{*}\|^{2}_{B_{i}^{-1}}\\\
&+\sum_{i=1}^{N}\|y^{+}-y^{*}\|^{2}_{B_{i}}-\sum_{i=1}^{N}\|y-y^{*}\|^{2}_{B_{i}}\\\
=&\mathcal{L}(y^{+},\lambda^{+})-\mathcal{L}(y,\lambda).\end{split}$ (49)

As we mentioned, Equation (45) is non-positive, therefore

$\mathcal{L}(y^{+},\lambda^{+})-\mathcal{L}(y,\lambda)\leq 0$

is guaranteed. This completes our proof.

##  Appendix F Proof of Theorem 4

In the lower level, Problem (7) of both algorithms is strongly convex.
Algorithm 1 and 2 will guarantee that, with finite iteration of the lower
level, a better minimizer $y^{+}$ can be found such that

$\begin{split}&\sum_{i=1}^{N}f_{i}(z)\\\
=&\underbrace{\sum_{i=1}^{N}\left(f_{i}(z)+\frac{\gamma}{2}\|z-z\|^{2}\right)+\sum_{i=1}^{N}\sigma_{i}\|z-z\|_{1}}_{\overset{\eqref{eq:
Merit}}{=}\Phi^{(z,{z})}(z)}\\\ \overset{\eqref{eq:
lower}}{>}&\underbrace{\sum_{i=1}^{N}\left(f_{i}(y^{+})+\frac{\gamma}{2}\|y^{+}-z\|^{2}\right)+\sum_{i=1}^{N}\sigma_{i}\|y^{+}-y^{+}\|_{1}}_{\overset{\eqref{eq:
Merit}}{=}\Phi^{(z,y^{+})}(y^{+})}.\end{split}$ (50)

By adding $\frac{\gamma N}{2}\|y^{+}-y^{+}\|^{2}$ on the right hand side of
(50), we can easily show

$\Phi^{({z},{z})}(z)-\Phi^{(y^{+},y^{+})}(y^{+})>\frac{\gamma
N}{2}\|y^{+}-z\|^{2}.$

This completes the proof.

## References

  * [1] Stephen Boyd, Neal Parikh, and Eric Chu.  Distributed optimization and statistical learning via the alternating direction method of multipliers.  Now Publishers Inc, 2011.
  * [2] Xu Du and Jingzhe Wang.  Consensus aladin: A framework for distributed optimization and its application in federated learning, 2023.
  * [3] Sur l’approximation, par éléments finis d’ordre un, et la résolution, par pénalisation-dualité d’une classe de problèmes de dirichlet non linéaires.  Revue française d’automatique, informatique, recherche opérationnelle. Analyse numérique, 9(R2):41–76, 1975.
  * [4] Daniel Gabay and Bertrand Mercier.  A dual algorithm for the solution of nonlinear variational problems via finite element approximation.  Computers & Mathematics with Applications, 2(1):17 – 40, 1976\.
  * [5] Zaiwen Wen, Donald Goldfarb, and Wotao Yin.  Alternating direction augmented lagrangian methods for semidefinite programming.  Mathematical Programming Computation, 2(3-4):203–230, 2010.
  * [6] Shenglong Zhou and Geoffrey Ye Li.  Federated learning via inexact admm.  arXiv preprint arXiv:2204.10607, 2022.
  * [7] Stephen Boyd, Neal Parikh, Eric Chu, Borja Peleato, and Jonathan Eckstein.  Distributed optimization and statistical learning via the alternating direction method of multipliers.  Found. Trends Mach. Learn., 3(1):1–122, 2011.
  * [8] Chengcheng Xu, Bruno Clerckx, Shiwa Chen, Yijie Mao, and Jianyun Zhang.  Rate-splitting multiple access for multi-antenna joint radar and communications.  IEEE Journal of Selected Topics in Signal Processing, 15(6):1332–1347, 2021.
  * [9] Mingyi Hong, Zhi-Quan Luo, and Meisam Razaviyayn.  Convergence analysis of alternating direction method of multipliers for a family of nonconvex problems.  SIAM Journal on Optimization, 26(1):337–364, 2016.
  * [10] Bingsheng He and Xiaoming Yuan.  On the o(1/n) convergence rate of the douglas–rachford alternating direction method.  SIAM Journal on Numerical Analysis, 50(2):700–709, 2012.
  * [11] W. Shi, Q. Ling, K. Yuan, G. Wu, and W. Yin.  On the linear convergence of the admm in decentralized consensus optimization.  IEEE Transactions on Signal Processing, 62(7):1750–1761, 2014.
  * [12] Qing Ling, Wei Shi, Gang Wu, and Alejandro Ribeiro.  Dlm: Decentralized linearized alternating direction method of multipliers.  IEEE Transactions on Signal Processing, 63(15):4051–4064, 2015\.
  * [13] B. Houska, J. Frasch, and M. Diehl.  An augmented Lagrangian based algorithm for distributed nonconvex optimization.  SIAM Journal on Optimization, 26(2):1101–1127, 2016.
  * [14] M. Hong, Z.-Q. Luo, and Meisam Razaviyayn.  Convergence analysis of alternating direction method of multipliers for a family of nonconvex problems.  SIAM Journal on Optimization, 26(1):337–364, 2016.
  * [15] Guoyin Li and Ting Kei Pong.  Global convergence of splitting methods for nonconvex composite optimization.  SIAM Journal on Optimization, 25(4):2434–2460, 2015.
  * [16] Fenghui Wang, Wenfei Cao, and Zongben Xu.  Convergence of multi-block bregman admm for nonconvex composite problems.  Science China Information Sciences, 61:1–12, 2018.
  * [17] Yu Yang, Qing-Shan Jia, Zhanbo Xu, Xiaohong Guan, and Costas J Spanos.  Proximal admm for nonconvex and nonsmooth optimization.  Automatica, 146:110551, 2022.
  * [18] Yu Wang, Wotao Yin, and Jinshan Zeng.  Global convergence of admm in nonconvex nonsmooth optimization.  Journal of Scientific Computing, 78:29–63, 2019.
  * [19] Yu Yang, Xiaohong Guan, Qing-Shan Jia, Liang Yu, Bolun Xu, and Costas J Spanos.  A survey of admm variants for distributed optimization: Problems, algorithms and features.  arXiv preprint arXiv:2208.03700, 2022.
  * [20] Kaizhao Sun and X Andy Sun.  A two-level distributed algorithm for nonconvex constrained optimization.  Computational Optimization and Applications, 84(2):609–649, 2023\.
  * [21] A. Engelmann, Y. Jiang, B. Houska, and T. Faulwassser.  Decomposition of non-convex optimization via bi-level distributed ALADIN.  IEEE Transactions on Control of Network Systems, 7(4):1848–1858, 2020.
  * [22] Gösta Stomberg, Alexander Engelmann, and Timm Faulwasser.  Decentralized non-convex optimization via bi-level sqp and admm.  arXiv preprint arXiv:2204.08786, 2022.
  * [23] Bin Wang, Jun Fang, Hongbin Li, Xiaojun Yuan, and Qing Ling.  Confederated learning: Federated learning with decentralized edge servers.  arXiv preprint arXiv:2205.14905, 2022.
  * [24] B. Houska and M. Diehl.  Nonlinear Robust Optimization via Sequential Convex Bilevel Programming.  Mathematical Programming, Series A, 142:539?577, 2013.
  * [25] Bernard Lemaire.  The proximal algorithm.  International series of numerical mathematics, 87:73–87, 1989.
  * [26] Alfredo N Iusem.  Augmented lagrangian methods and proximal point methods for convex optimization.  Investigación Operativa, 8(11-49):7, 1999.
  * [27] Jorge Nocedal and Stephen Wright.  Numerical optimization.  Springer Science & Business Media, New York, 2006.
  * [28] Qing Ling, Wei Shi, Gang Wu, and Alejandro Ribeiro.  Dlm: Decentralized linearized alternating direction method of multipliers.  IEEE Transactions on Signal Processing, 63(15):4051–4064, 2015\.
  * [29] Aryan Mokhtari, Wei Shi, Qing Ling, and Alejandro Ribeiro.  Dqm: Decentralized quadratically approximated alternating direction method of multipliers.  IEEE Transactions on Signal Processing, 64(19):5158–5173, 2016\.


# Search for dark matter annihilations in the center of the Earth with IceCube

The IceCube Collaboration
(a complete list of authors can be found at the end of the proceedings)

###### Abstract

The nature of Dark Matter remains one of the most important unresolved
questions of fundamental physics. Many models, including the Weakly
Interacting Massive Particles (WIMPs), assume Dark Matter to be a particle and
predict a weak coupling with Standard Model matter. If Dark Matter particles
can scatter off nuclei in the vicinity of a massive object, such as a star or
a planet, they may lose kinetic energy and become gravitationally trapped in
the center of such objects, including Earth. As Dark Matter accumulates in the
center of the Earth, self-annihilation of WIMPs into Standard Model particles
can result in an excess of neutrinos coming from the center of the Earth and
detectable at the IceCube Neutrino Observatory, situated at the geographic
South Pole. A search for excess neutrinos from these annihilations has been
performed on 10 years of IceCube data, and results have been interpreted in
the context of a number of WIMP annihilation channels
($\chi\chi\rightarrow\tau^{+}\tau^{-}/W^{+}W^{-}/b\overline{b}$) and masses
ranging from 10 GeV to 10 TeV. We present the results from this analysis and
compare the outcome with previous searches by other experiments. This analysis
yields competitive and world-leading results for masses $m_{\chi}$ > 100 GeV.

Corresponding authors: Giovanni Renzi1, J. A. Aguilar1∗
1 Université Libre de Bruxelles, Science Faculty CP230, B-1050 Brussels,
Belgium
∗ Presenter

##  1 Introduction

Over the last century, an increasing amount of compelling evidence has
emerged, indicating that approximately $27\%$ of the Universe and about $85\%$
of its matter content [1] is composed of unknown matter. This elusive form of
matter, referred to as dark matter (DM), must exhibit a low probability of
interaction with ordinary matter.

DM particles, denoted $\chi$, can be captured by the gravitational potential
of celestial bodies, such as the Sun or the Earth. This capture occurs through
repeated scattering interactions with DM particles in orbits intersecting
these celestial objects. As a result of these interactions, dark matter
particles can lose energy and fall below the gravitational escape velocity.
Over time, this will generate an over abundance of dark matter at the center
of the celestial body favoring possible process of self-annihilation and
therefore reducing the among to DM particles. The process of DM capturing in
Earth was first studied in [2]. The evolution of number density of DM
particles, $N_{\chi}$, will follow the equation that describes the interplay
between capturing and annihilation as:

$\frac{{\rm d}N_{\chi}}{{\rm d}t}=C_{\rm C}-C_{\rm A}N_{\chi}^{2}-C_{\rm
E}N_{\chi}.$ (1)

In this equation $C_{\rm C}$ refers to the capture rate and $C_{\rm
A}N_{\chi}^{2}=\Gamma_{\rm A}$ is the annihilation rate. The third process,
represented by $C_{\rm E}N_{\chi}$, is the so-called evaporation, which
describes the process in which after scattering the DM particles gains kinetic
energy above the escape velocity and therefore escapes. Evaporation can be
ignored for dark matter masses above $\sim 12\;\rm{GeV}$ [3]. The capturing
rate, depends on the local DM density, the scattering cross-section, and the
chemical composition of the target body. The local halo DM density cannot be
stated with certainty and is expected to be in the range $0.1<\rho_{\rm
DM}/{\rm GeV\>cm}^{-3}<0.5$ [4]. The value assumed conventionally and used for
the analysis presented in this work is $\rho_{\rm DM}=0.3$ GeV cm-3. Figure 1
shows the Earth’s capturing rate assuming the Earth’s composition in [5] as
function of the DM mass.

Figure 1: Capture rate as a function of the dark matter particle mass for the
spin-independent WIMP-nucleon scattering cross section value $\sigma_{\chi
N}^{\rm SI}=10^{-44}$ cm2. The peaks correspond to resonance capture with the
most abundant elements on Earth: O, Mg, Si, and Fe. Computed using DarkSusy
[6] [7] [8].

The peaks in the capturing rate are due to resonance capture with the most
abundant elements on the Earth. Most abundant elements on Earth do not carry
spin hence the most relevant scattering process is the spin-independent dark
matter-nucleon scattering, $\sigma_{\chi N}^{\rm SI}$.

As can be seen in Eq. (1), the annihilation rate depends on the dark matter
number density squared, as two dark matter particles are needed for
annihilation, as well as the self-annihilation cross-section, usually
expressed as a velocity averaged cross-section, $\langle\sigma_{\rm
A}v\rangle$. Equation (1) can be solved and has a solution that can be
expressed as:

$\Gamma_{\rm A}(t)=\frac{1}{2}C_{\rm C}\tanh^{2}\left(\frac{t}{\tau}\right),$
(2)

where $\tau=(C_{\rm C}C_{\rm A})^{-1/2}$ is the time scale for the capture and
annihilation processes to reach equilibrium. The equilibrium condition is
fulfilled when $t\geq\tau$, so that Eq. (2) becomes $C_{\rm
C}=\frac{1}{2}\Gamma_{\rm A}$. In this case the annihilation rate reaches its
maximum. For typical value of the thermal relic annihilation cross section the
Earth is far from equilibrium ($t_{\oplus}\ll\tau$). However, enhancements in
the annihilation cross-section such as a nearby clump of DM or a low-velocity
Sommerfeld effect, will increase the annihilation rate for the Earth while not
for the Sun which is, likely, already in hydrostatic equilibrium [9]. Due to
its low escape velocity, Earth capturing is very sensitive to the low tail of
the DM velocity distribution. However the low speed distribution is not
precisely known, and deviations from the standard halo model could
significantly enhance the capturing rate [10].

As DM annihilates, it will produce stable particles, such $\gamma$-rays and
neutrinos. Given the density of the Earth, only neutrinos can reach the
surface and interact in the vicinity of a neutrino telescope. This constitutes
the base of the detection principle for indirect searches of dark matter from
the center of the Earth.

##  2 The IceCube Neutrino Observatory

The IceCube Neutrino Observatory [11] is a one-cubic-kilometer neutrino
telescope situated at the geographic South Pole. The detector is composed of
5160 photo-multiplier units (DOMs) deployed on 86 strings and is buried in the
Antarctic ice-cap between the depths of 1,450 m and 2,450 m. IceCube detects
Cherenkov light induced by the passage of superluminal charged leptons and
hadrons produced in neutrino interactions in ice. The number of Cherenkov
photons and their arrival time is used to reconstruct the direction and energy
information of the incoming neutrino. Based on its geometry, IceCube can
detect neutrinos with energy ranging from a few GeV to a few PeV. A denser
sub-array, DeepCore [11], is placed between 2100 m and 2450 m of depth and
extends the detection threshold to lower energy neutrinos, between 10 GeV and
100 GeV.

##  3 Data Selection

Neutrinos produced as by-products of the DM self-annihilation at the center of
the Earth, will reach the IceCube detector with vertical up-going directions.
Hence, the data selection is designed and optimized to keep vertical tracks.
The first step of the processing focuses on reducing the muon atmospheric
background events which arrive to the detector with downward going directions.
Once data is reduced to an manageable level, we applied a random forest of
boosted decision trees (BDTs) to discriminate mis-reconstructed muon events
from neutrino events. A random forest assigns to each event a score ranging
from -1 (most background-like) to +1 (most signal-like). At this point, we can
observe differences between the highest and lowest energy signal events: low-
energy events appear more like cascade events, mostly comprised in DeepCore;
high-energy events consist of kilometer-long tracks. In a intermediate step, a
higher-performing but computationally intensive energy reconstruction is
applied. The last selection based on the random forest classification is
optimized to maximize the sensitivity of the analysis. Two BDT are applied,
one which is tuned to select high energy signal, and another one optimized for
low energy dark matter masses.

Figure 2 shows the differences in the score distributions for the LE and HE
baseline signal configurations for the two random forests.

Figure 2: Score distributions for the LE (left) and HE (right) random forests.
The baseline LE and HE signal configurations are in light purple and dark
purple, respectively. Atmospheric neutrinos and muons are in green and pink
respectively. The grey band indicates the total background estimation while
the black dots represent the data verification sample.

In this analysis we used 3,619 days of data taking over ten detector seasons
from May 2011 to May 2020. A subset of 353 days, taken sparsely over the ten
seasons, have been used as a verification dataset along the event selection
process. The verification dataset has been omitted from the analysis, leaving
a final of 3,266 available days of data.

##  4 Background and Signal Simulation

Neutrinos from the center of the Earth have an unique direction in IceCube’s
local coordinate system. There is no other direction in the sky in which
IceCube has the same detector response, therefore background estimation
techniques, like on-source/off-source methods, are not possible. This includes
the typical method used in searches for point sources in neutrinos telescopes
called right-ascension scrambling, which makes use of the uniform response of
the detector in right-ascension due to Earth’s rotation. For this reason, this
analysis relies on MonteCarlo simulation to estimate the background at a
perfect upward going direction. The atmospheric muons generated by cosmic-ray
showers are simulated using CORSIKA [12]. The diffuse atmospheric, prompt [13]
and astrophysical neutrinos [14] are simulated using two different neutrino
generators: GENIE [15] for energies up to 100 GeV and NuGen [16] beyond. The
reason for these two simulations is due to the different implementation and
energy regime of the neutrino nucleon interaction, with NuGen covering the
deep inelastic scattering and GENIE the quasi-elastic interaction. After
generating interaction vertexes, all charge particles are then propagated
through the detector and their induced Cherenkov light production and
detection also simulated.

Neutrinos from the annihilation of DM at the center of the Earth are simulated
using the WIMPsim [17]. Re-weighting of the isotropic simulation from GENIE
and NuGen becomes impractical again due to the unique direction in local
coordinates of the DM signal. Neutrinos from WIMPSim are generated as coming
from the center of the Earth ($\theta=180^{\circ}$) and their interaction
products are propagated and simulated in the same fashion as astrophysical and
atmospheric neutrinos.

We produced simulations over DM particle mass values ranging from 10 GeV to 10
TeV in three annihilation channels:
$\chi\chi\rightarrow\tau^{+}\tau^{-}$/$W^{+}W^{-}$/$b\bar{b}$. All the
simulated mass-channel configurations are resumed in Table 1. We chose two
baseline signal configurations, a low-energy one for
$\chi\chi\rightarrow\tau^{+}\tau^{-}$, $m_{\chi}=50$ GeV, and a high-energy
one for $\chi\chi\rightarrow W^{+}W^{-}$, $m_{\chi}=1$ TeV. As the event
selection selects track-like events, only muon neutrinos are considered for
these simulations.

Table 1: Summary of WIMPs simulation scenarios produced with WimpSim. Channel | Masses
---|---
$\chi\chi\rightarrow\tau^{+}\tau^{-}$ | [10, 20, 35, 50] GeV
[100, 250, 500] GeV
[1, 3, 5, 10] TeV
$\chi\chi\rightarrow W^{+}W^{-}$ | [100, 250, 500] GeV
[1, 3, 5, 10] TeV
$\chi\chi\rightarrow b\bar{b}$ | [35, 50] GeV
[100, 250, 500] GeV
[1, 3, 5, 10] TeV

##  5 Analysis Method

We used a binned Poisson likelihood method using as observables the
reconstructed energy of the event and their reconstructed zenith angle. The
expected signal from DM will arrive from almost vertical directions
($\theta_{rec}\sim 180^{\circ}$). We produce binned probability density
functions (PDFs) over these two observables based on the final distributions
obtained in the event selection process after applying a kernel density
estimation (KDE) technique to cope with statistical fluctuations due to low
statistics. For the LE analysis, we include events with energy in the range
$1\;{\rm GeV}<E_{\rm reco}<10^{4}{\rm GeV}$. We found the optimised number of
bins at $32\times 32$, for a total of $1024$ bins. For the HE analysis, the
energy range (in GeV) is $1<E_{\rm reco}<10^{5}$ and the optimised bins grid
is $100\times 100$ for a total of $10^{4}$ bins. For both analyses, the zenith
angle range is $160^{\circ}<\theta_{\rm reco}<180^{\circ}$. The PDFs for the
signal and atmospheric background are shown in Figs. 3 and 4 for LE and HE,
respectively.

Figure 3: PDFs for the LE analysis. The zenith angle, $\theta_{\rm reco}$, is
on the $x$-axis while the log-energy, $\log(E_{\rm reco}$, is on the $y$-axis.
The LE signal baseline and the atmospheric background are shown on the left
and on the right, respectively.

Figure 4: PDFs for the HE analysis. The zenith angle, $\theta_{\rm reco}$, is
on the $x$-axis while the log-energy, $\log E_{\rm reco}$, is on the $y$-axis.
The LE signal baseline and the atmospheric background are shown on the left
and on the right, respectively.

The expected fraction of events in the $i$-th bin, $\lambda_{i}$, can be
expressed as:

$\lambda_{i}=\xi S_{i}+(1-\xi)B_{i},$ (3)

where $\xi$ is the signal fraction and the parameter against with the
likelihood is minimized. The background density $B_{i}$ is composed, likewise,
of many different components: atmospheric neutrinos and muons, astrophysical
neutrinos and prompt neutrinos. Their relative normalization are kept as
nuisance parameters in the minimization and expressed as:

$\lambda_{i}=\xi
S_{i}+\sum_{j=1}^{n-1}\prod_{k=0}^{j-1}(1-n_{k})n_{j}B^{j}_{i}+\prod_{k=0}^{n-1}(1-n_{k})B^{n}_{i},$
(4)

where $n_{k}$ refers to the different background event fractions which can be
encompassed by $\overrightarrow{\eta}=(n_{1},...,n_{n-1})$ and $k_{i}$ being
the observed number of events in the $i$-th bin. We then search the
parameters, $\xi$, and $\eta$ that minimize the likelihood:

$-\log{\mathcal{L}(\xi,\overrightarrow{\eta})}=\sum_{i}^{N_{bins}}(-k_{i}\log\lambda_{i}+\lambda_{i}).$
(5)

In order to perform an hypothesis testing, we define a test-statistics as the
likelihood ratio:

$t_{\xi}=2\log\frac{\mathcal{L}(\xi,\hat{\hat{\overrightarrow{\eta}}})}{\mathcal{L}(\hat{\xi},\hat{\overrightarrow{\eta}})}=2(\log\mathcal{L}(\xi,\hat{\hat{\overrightarrow{\eta}}})-\log\mathcal{L}(\hat{\xi},\hat{\overrightarrow{\eta}})),$
(6)

where $\hat{\xi}$, $\hat{\eta}$ are the estimated parameters that minimize the
likelihood, $\mathcal{L}(\hat{\xi},\hat{\overrightarrow{\eta}})$, and
$\hat{\hat{\eta}}$ is the value of $\eta$ that minimizes the likelihood for
any given value of $\xi$. To claim a discovery we estimate the agreement with
respect to the null hypothesis $\xi=0$, and so the test-statistics is defined
as $t_{0}$. Assuming Wilk’s theorem [18, 19] and that the distribution of
$t_{0}$ under sample of null-hypothesis follows a $\chi^{2}$ distribution, we
can estimate the p-value in terms of number of sigmas $z$-score as $z-{\rm
score}=\sqrt{t_{0}}$.

##  6 Results and Conclusions

We found no significant excess over all channels and masses tested. The
highest pre-trial significance belongs to the HE analysis at $1.94\sigma$ for
the channel $\chi\chi\rightarrow b\bar{b}$ and mass $m_{\chi}=250$ GeV. The
trial factor correction gives a post-trial significance of $1.06\sigma$. Upper
limits in the signal fraction of events $\xi$, can be converted into upper
limits on the spin-independent DM-nucleon cross-section $\sigma_{\chi N}^{\rm
SI}$ after assuming a value on the DM self-annihilation cross-section.

The $90\%$ C.L. upper limits for $\sigma_{\chi N}^{\rm SI}$ are presented in
Fig. 5, compared to results from the other neutrino telescopes ANTARES [20]
and Super-Kamiokande [21].

Figure 5: Left: Spin-indenpendent DM-nucleon cross-section $90\%$ C.L. upper
limits. Limits in bold are from this analysis, ANTARES (dashed) and Super-
Kamiokande (dot-dashed) limits are included. The colour code identifies the
annihilation channels, blue for annihilation into $\tau^{+}\tau^{-}$, yellow
for $W^{+}W^{-}$ and green for $b\bar{b}$. Right: Upper limits for this
analysis (blue) for the $\chi\chi\rightarrow\tau^{+}\tau^{-}$ annihilation
channel compared to direct detection upper limits from the crystal experiments
DAMA/LIBRA [22] (grey areas) and COSINE100 [23] (green) and from XENON1T [24]
(black).

Although optimized for signals from WIMP-like particles, the results from this
analysis can also be re-casted to set limits on different DM models. An
example can be found in [25], where limits on the coupling constant of the
effective field theory of dark matter [26] were computed.

With the installation of the IceCube Upgrade [27] a significant improvement in
IceCube’s sensitivity can be expected, especially in the low-mass region.

## References

  * [1] Planck Collaboration, N. Aghanim et al. Astron. Astrophys. 641 (2020) A6.
  * [2] A. Gould Astrophys. J. 321 (Oct., 1987) 571.
  * [3] R. Garani and S. Palomares-Ruiz JCAP 2022 no. 05, (May, 2022) 042.
  * [4] L. Necib, M. Lisanti, and V. A. Belokurov Astrophys. J. 874 (2019) .
  * [5] W. McDonough Treatise on Geochemistry 2 (11, 2003) 547–568.
  * [6] T. Bringmann et al. JCAP 2018 (02, 2018) .
  * [7] P. Gondolo et al. JCAP 2004 (06, 2004) .
  * [8] J. Edsjö et al., “Darksusy, http://www.darksusy.org.”.
  * [9] C. Delaunay, P. J. Fox, and G. Perez JHEP 2009 no. 05, (May, 2009) 099.
  * [10] T. Bruch et al. Physics Letters B 674 no. 4, (2009) 250–256.
  * [11] IceCube Collaboration, M. G. Aartsen et al. JINST 12 no. 03, (2017) P03012.
  * [12] D. Heck et al. FZKA 6019 (1998) .
  * [13] R. Enberg, M. H. Reno, and I. Sarcevic Phys. Rev. D 78 (Aug, 2008) 043005.
  * [14] IceCube Collaboration, R. Abbasi et al. Astrophys. J. 928 no. 1, (2022) 50.
  * [15] C. Andreopoulos et al. arXiv e-prints (10, 2015) .
  * [16] A. Gazizov and M. Kowalski Computer Physics Communications 172 no. 3, (2005) 203–213.
  * [17] J. Edsjö, J. Elevant, and C. Niblaeus WimpSim Neutrino Monte Carlo, http://wimpsim.astroparticle.se/ .
  * [18] S. S. Wilks The Annals of Mathematical Statistics 9 (1938) 60–62.
  * [19] H. Chernoff Annals of Mathematical Statistics 25 (1954) .
  * [20] ANTARES Collaboration, A. Albert et al. Phys. Dark Univ. 16 (2017) 41–48.
  * [21] Super-Kamiokande Collaboration, P. Mijakowski J. Phys.: Conf. Ser. 1342 (Jan, 2020) 012075.
  * [22] R. Bernabei et al. Universe 4 no. 11, (2018) 116.
  * [23] G. Adhikari et al. Phys. Rev. Lett. 123 (07, 2019) .
  * [24] XENON Collaboration, E. Aprile et al. Phys. Rev. Lett. 122 no. 14, (2019) 141301.
  * [25] G. Renzi, Search for dark matter annihilation in the center of the Earth with ten years of IceCube data.  PhD thesis, Université Libre de Bruxelles, 2022.
  * [26] R. Catena JCAP 2017 no. 01, (Jan, 2017) 059–059.
  * [27] IceCube Collaboration, A. Ishihara PoS ICRC2019 (2021) 1031.

## Full Author List: IceCube Collaboration

R. Abbasi17, M. Ackermann63, J. Adams18, S. K. Agarwalla40, 64, J. A.
Aguilar12, M. Ahlers22, J.M. Alameddine23, N. M. Amin44, K. Andeen42, G.
Anton26, C. Argüelles14, Y. Ashida53, S. Athanasiadou63, S. N. Axani44, X.
Bai50, A. Balagopal V.40, M. Baricevic40, S. W. Barwick30, V. Basu40, R. Bay8,
J. J. Beatty20, 21, J. Becker Tjus11, 65, J. Beise61, C. Bellenghi27, C.
Benning1, S. BenZvi52, D. Berley19, E. Bernardini48, D. Z. Besson36, E.
Blaufuss19, S. Blot63, F. Bontempo31, J. Y. Book14, C. Boscolo Meneguolo48, S.
Böser41, O. Botner61, J. Böttcher1, E. Bourbeau22, J. Braun40, B. Brinson6, J.
Brostean-Kaiser63, R. T. Burley2, R. S. Busse43, D. Butterfield40, M. A.
Campana49, K. Carloni14, E. G. Carnie-Bronca2, S. Chattopadhyay40, 64, N.
Chau12, C. Chen6, Z. Chen55, D. Chirkin40, S. Choi56, B. A. Clark19, L.
Classen43, A. Coleman61, G. H. Collin15, A. Connolly20, 21, J. M. Conrad15, P.
Coppin13, P. Correa13, D. F. Cowen59, 60, P. Dave6, C. De Clercq13, J. J.
DeLaunay58, D. Delgado14, S. Deng1, K. Deoskar54, A. Desai40, P. Desiati40, K.
D. de Vries13, G. de Wasseige37, T. DeYoung24, A. Diaz15, J. C. Díaz-Vélez40,
M. Dittmer43, A. Domi26, H. Dujmovic40, M. A. DuVernois40, T. Ehrhardt41, P.
Eller27, E. Ellinger62, S. El Mentawi1, D. Elsässer23, R. Engel31, 32, H.
Erpenbeck40, J. Evans19, P. A. Evenson44, K. L. Fan19, K. Fang40, K. Farrag16,
A. R. Fazely7, A. Fedynitch57, N. Feigl10, S. Fiedlschuster26, C. Finley54, L.
Fischer63, D. Fox59, A. Franckowiak11, A. Fritz41, P. Fürst1, J. Gallagher39,
E. Ganster1, A. Garcia14, L. Gerhardt9, A. Ghadimi58, C. Glaser61, T.
Glauch27, T. Glüsenkamp26, 61, N. Goehlke32, J. G. Gonzalez44, S. Goswami58,
D. Grant24, S. J. Gray19, O. Gries1, S. Griffin40, S. Griswold52, K. M.
Groth22, C. Günther1, P. Gutjahr23, C. Haack26, A. Hallgren61, R. Halliday24,
L. Halve1, F. Halzen40, H. Hamdaoui55, M. Ha Minh27, K. Hanson40, J. Hardin15,
A. A. Harnisch24, P. Hatch33, A. Haungs31, K. Helbing62, J. Hellrung11, F.
Henningsen27, L. Heuermann1, N. Heyer61, S. Hickford62, A. Hidvegi54, C.
Hill16, G. C. Hill2, K. D. Hoffman19, S. Hori40, K. Hoshina40, 66, W. Hou31,
T. Huber31, K. Hultqvist54, M. Hünnefeld23, R. Hussain40, K. Hymon23, S. In56,
A. Ishihara16, M. Jacquart40, O. Janik1, M. Jansson54, G. S. Japaridze5, M.
Jeong56, M. Jin14, B. J. P. Jones4, D. Kang31, W. Kang56, X. Kang49, A.
Kappes43, D. Kappesser41, L. Kardum23, T. Karg63, M. Karl27, A. Karle40, U.
Katz26, M. Kauer40, J. L. Kelley40, A. Khatee Zathul40, A. Kheirandish34, 35,
J. Kiryluk55, S. R. Klein8, 9, A. Kochocki24, R. Koirala44, H. Kolanoski10, T.
Kontrimas27, L. Köpke41, C. Kopper26, D. J. Koskinen22, P. Koundal31, M.
Kovacevich49, M. Kowalski10, 63, T. Kozynets22, J. Krishnamoorthi40, 64, K.
Kruiswijk37, E. Krupczak24, A. Kumar63, E. Kun11, N. Kurahashi49, N. Lad63, C.
Lagunas Gualda63, M. Lamoureux37, M. J. Larson19, S. Latseva1, F. Lauber62, J.
P. Lazar14, 40, J. W. Lee56, K. Leonard DeHolton60, A. Leszczyńska44, M.
Lincetto11, Q. R. Liu40, M. Liubarska25, E. Lohfink41, C. Love49, C. J. Lozano
Mariscal43, L. Lu40, F. Lucarelli28, W. Luszczak20, 21, Y. Lyu8, 9, J.
Madsen40, K. B. M. Mahn24, Y. Makino40, E. Manao27, S. Mancina40, 48, W. Marie
Sainte40, I. C. Mariş12, S. Marka46, Z. Marka46, M. Marsee58, I. Martinez-
Soler14, R. Maruyama45, F. Mayhew24, T. McElroy25, F. McNally38, J. V. Mead22,
K. Meagher40, S. Mechbal63, A. Medina21, M. Meier16, Y. Merckx13, L. Merten11,
J. Micallef24, J. Mitchell7, T. Montaruli28, R. W. Moore25, Y. Morii16, R.
Morse40, M. Moulai40, T. Mukherjee31, R. Naab63, R. Nagai16, M. Nakos40, U.
Naumann62, J. Necker63, A. Negi4, M. Neumann43, H. Niederhausen24, M. U.
Nisa24, A. Noell1, A. Novikov44, S. C. Nowicki24, A. Obertacke Pollmann16, V.
O’Dell40, M. Oehler31, B. Oeyen29, A. Olivas19, R. Ørsøe27, J. Osborn40, E.
O’Sullivan61, H. Pandya44, N. Park33, G. K. Parker4, E. N. Paudel44, L.
Paul42, 50, C. Pérez de los Heros61, J. Peterson40, S. Philippen1, A.
Pizzuto40, M. Plum50, A. Pontén61, Y. Popovych41, M. Prado Rodriguez40, B.
Pries24, R. Procter-Murphy19, G. T. Przybylski9, C. Raab37, J. Rack-Helleis41,
K. Rawlins3, Z. Rechav40, A. Rehman44, P. Reichherzer11, G. Renzi12, E.
Resconi27, S. Reusch63, W. Rhode23, B. Riedel40, A. Rifaie1, E. J. Roberts2,
S. Robertson8, 9, S. Rodan56, G. Roellinghoff56, M. Rongen26, C. Rott53, 56,
T. Ruhe23, L. Ruohan27, D. Ryckbosch29, I. Safa14, 40, J. Saffer32, D.
Salazar-Gallegos24, P. Sampathkumar31, S. E. Sanchez Herrera24, A. Sandrock62,
M. Santander58, S. Sarkar25, S. Sarkar47, J. Savelberg1, P. Savina40, M.
Schaufel1, H. Schieler31, S. Schindler26, L. Schlickmann1, B. Schlüter43, F.
Schlüter12, N. Schmeisser62, T. Schmidt19, J. Schneider26, F. G. Schröder31,
44, L. Schumacher26, G. Schwefer1, S. Sclafani19, D. Seckel44, M. Seikh36, S.
Seunarine51, R. Shah49, A. Sharma61, S. Shefali32, N. Shimizu16, M. Silva40,
B. Skrzypek14, B. Smithers4, R. Snihur40, J. Soedingrekso23, A. Søgaard22, D.
Soldin32, P. Soldin1, G. Sommani11, C. Spannfellner27, G. M. Spiczak51, C.
Spiering63, M. Stamatikos21, T. Stanev44, T. Stezelberger9, T. Stürwald62, T.
Stuttard22, G. W. Sullivan19, I. Taboada6, S. Ter-Antonyan7, M. Thiesmeyer1,
W. G. Thompson14, J. Thwaites40, S. Tilav44, K. Tollefson24, C. Tönnis56, S.
Toscano12, D. Tosi40, A. Trettin63, C. F. Tung6, R. Turcotte31, J. P.
Twagirayezu24, B. Ty40, M. A. Unland Elorrieta43, A. K. Upadhyay40, 64, K.
Upshaw7, N. Valtonen-Mattila61, J. Vandenbroucke40, N. van Eijndhoven13, D.
Vannerom15, J. van Santen63, J. Vara43, J. Veitch-Michaelis40, M. Venugopal31,
M. Vereecken37, S. Verpoest44, D. Veske46, A. Vijai19, C. Walck54, C.
Weaver24, P. Weigel15, A. Weindl31, J. Weldert60, C. Wendt40, J. Werthebach23,
M. Weyrauch31, N. Whitehorn24, C. H. Wiebusch1, N. Willey24, D. R. Williams58,
L. Witthaus23, A. Wolf1, M. Wolf27, G. Wrede26, X. W. Xu7, J. P. Yanez25, E.
Yildizci40, S. Yoshida16, R. Young36, F. Yu14, S. Yu24, T. Yuan40, Z. Zhang55,
P. Zhelnin14, M. Zimmerman40

1 III. Physikalisches Institut, RWTH Aachen University, D-52056 Aachen,
Germany
2 Department of Physics, University of Adelaide, Adelaide, 5005, Australia
3 Dept. of Physics and Astronomy, University of Alaska Anchorage, 3211
Providence Dr., Anchorage, AK 99508, USA
4 Dept. of Physics, University of Texas at Arlington, 502 Yates St., Science
Hall Rm 108, Box 19059, Arlington, TX 76019, USA
5 CTSPS, Clark-Atlanta University, Atlanta, GA 30314, USA
6 School of Physics and Center for Relativistic Astrophysics, Georgia
Institute of Technology, Atlanta, GA 30332, USA
7 Dept. of Physics, Southern University, Baton Rouge, LA 70813, USA
8 Dept. of Physics, University of California, Berkeley, CA 94720, USA
9 Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA
10 Institut für Physik, Humboldt-Universität zu Berlin, D-12489 Berlin,
Germany
11 Fakultät für Physik & Astronomie, Ruhr-Universität Bochum, D-44780 Bochum,
Germany
12 Université Libre de Bruxelles, Science Faculty CP230, B-1050 Brussels,
Belgium
13 Vrije Universiteit Brussel (VUB), Dienst ELEM, B-1050 Brussels, Belgium
14 Department of Physics and Laboratory for Particle Physics and Cosmology,
Harvard University, Cambridge, MA 02138, USA
15 Dept. of Physics, Massachusetts Institute of Technology, Cambridge, MA
02139, USA
16 Dept. of Physics and The International Center for Hadron Astrophysics,
Chiba University, Chiba 263-8522, Japan
17 Department of Physics, Loyola University Chicago, Chicago, IL 60660, USA
18 Dept. of Physics and Astronomy, University of Canterbury, Private Bag 4800,
Christchurch, New Zealand
19 Dept. of Physics, University of Maryland, College Park, MD 20742, USA
20 Dept. of Astronomy, Ohio State University, Columbus, OH 43210, USA
21 Dept. of Physics and Center for Cosmology and Astro-Particle Physics, Ohio
State University, Columbus, OH 43210, USA
22 Niels Bohr Institute, University of Copenhagen, DK-2100 Copenhagen, Denmark
23 Dept. of Physics, TU Dortmund University, D-44221 Dortmund, Germany
24 Dept. of Physics and Astronomy, Michigan State University, East Lansing, MI
48824, USA
25 Dept. of Physics, University of Alberta, Edmonton, Alberta, Canada T6G 2E1
26 Erlangen Centre for Astroparticle Physics, Friedrich-Alexander-Universität
Erlangen-Nürnberg, D-91058 Erlangen, Germany
27 Technical University of Munich, TUM School of Natural Sciences, Department
of Physics, D-85748 Garching bei München, Germany
28 Département de physique nucléaire et corpusculaire, Université de Genève,
CH-1211 Genève, Switzerland
29 Dept. of Physics and Astronomy, University of Gent, B-9000 Gent, Belgium
30 Dept. of Physics and Astronomy, University of California, Irvine, CA 92697,
USA
31 Karlsruhe Institute of Technology, Institute for Astroparticle Physics,
D-76021 Karlsruhe, Germany
32 Karlsruhe Institute of Technology, Institute of Experimental Particle
Physics, D-76021 Karlsruhe, Germany
33 Dept. of Physics, Engineering Physics, and Astronomy, Queen’s University,
Kingston, ON K7L 3N6, Canada
34 Department of Physics & Astronomy, University of Nevada, Las Vegas, NV,
89154, USA
35 Nevada Center for Astrophysics, University of Nevada, Las Vegas, NV 89154,
USA
36 Dept. of Physics and Astronomy, University of Kansas, Lawrence, KS 66045,
USA
37 Centre for Cosmology, Particle Physics and Phenomenology - CP3, Université
catholique de Louvain, Louvain-la-Neuve, Belgium
38 Department of Physics, Mercer University, Macon, GA 31207-0001, USA
39 Dept. of Astronomy, University of Wisconsin–Madison, Madison, WI 53706, USA
40 Dept. of Physics and Wisconsin IceCube Particle Astrophysics Center,
University of Wisconsin–Madison, Madison, WI 53706, USA
41 Institute of Physics, University of Mainz, Staudinger Weg 7, D-55099 Mainz,
Germany
42 Department of Physics, Marquette University, Milwaukee, WI, 53201, USA
43 Institut für Kernphysik, Westfälische Wilhelms-Universität Münster, D-48149
Münster, Germany
44 Bartol Research Institute and Dept. of Physics and Astronomy, University of
Delaware, Newark, DE 19716, USA
45 Dept. of Physics, Yale University, New Haven, CT 06520, USA
46 Columbia Astrophysics and Nevis Laboratories, Columbia University, New
York, NY 10027, USA
47 Dept. of Physics, University of Oxford, Parks Road, Oxford OX1 3PU, United
Kingdom
48 Dipartimento di Fisica e Astronomia Galileo Galilei, Università Degli Studi
di Padova, 35122 Padova PD, Italy
49 Dept. of Physics, Drexel University, 3141 Chestnut Street, Philadelphia, PA
19104, USA
50 Physics Department, South Dakota School of Mines and Technology, Rapid
City, SD 57701, USA
51 Dept. of Physics, University of Wisconsin, River Falls, WI 54022, USA
52 Dept. of Physics and Astronomy, University of Rochester, Rochester, NY
14627, USA
53 Department of Physics and Astronomy, University of Utah, Salt Lake City, UT
84112, USA
54 Oskar Klein Centre and Dept. of Physics, Stockholm University, SE-10691
Stockholm, Sweden
55 Dept. of Physics and Astronomy, Stony Brook University, Stony Brook, NY
11794-3800, USA
56 Dept. of Physics, Sungkyunkwan University, Suwon 16419, Korea
57 Institute of Physics, Academia Sinica, Taipei, 11529, Taiwan
58 Dept. of Physics and Astronomy, University of Alabama, Tuscaloosa, AL
35487, USA
59 Dept. of Astronomy and Astrophysics, Pennsylvania State University,
University Park, PA 16802, USA
60 Dept. of Physics, Pennsylvania State University, University Park, PA 16802,
USA
61 Dept. of Physics and Astronomy, Uppsala University, Box 516, S-75120
Uppsala, Sweden
62 Dept. of Physics, University of Wuppertal, D-42119 Wuppertal, Germany
63 Deutsches Elektronen-Synchrotron DESY, Platanenallee 6, 15738 Zeuthen,
Germany
64 Institute of Physics, Sachivalaya Marg, Sainik School Post, Bhubaneswar
751005, India
65 Department of Space, Earth and Environment, Chalmers University of
Technology, 412 96 Gothenburg, Sweden
66 Earthquake Research Institute, University of Tokyo, Bunkyo, Tokyo 113-0032,
Japan

### Acknowledgements

The authors gratefully acknowledge the support from the following agencies and
institutions: USA – U.S. National Science Foundation-Office of Polar Programs,
U.S. National Science Foundation-Physics Division, U.S. National Science
Foundation-EPSCoR, Wisconsin Alumni Research Foundation, Center for High
Throughput Computing (CHTC) at the University of Wisconsin–Madison, Open
Science Grid (OSG), Advanced Cyberinfrastructure Coordination Ecosystem:
Services & Support (ACCESS), Frontera computing project at the Texas Advanced
Computing Center, U.S. Department of Energy-National Energy Research
Scientific Computing Center, Particle astrophysics research computing center
at the University of Maryland, Institute for Cyber-Enabled Research at
Michigan State University, and Astroparticle physics computational facility at
Marquette University; Belgium – Funds for Scientific Research (FRS-FNRS and
FWO), FWO Odysseus and Big Science programmes, and Belgian Federal Science
Policy Office (Belspo); Germany – Bundesministerium für Bildung und Forschung
(BMBF), Deutsche Forschungsgemeinschaft (DFG), Helmholtz Alliance for
Astroparticle Physics (HAP), Initiative and Networking Fund of the Helmholtz
Association, Deutsches Elektronen Synchrotron (DESY), and High Performance
Computing cluster of the RWTH Aachen; Sweden – Swedish Research Council,
Swedish Polar Research Secretariat, Swedish National Infrastructure for
Computing (SNIC), and Knut and Alice Wallenberg Foundation; European Union –
EGI Advanced Computing for research; Australia – Australian Research Council;
Canada – Natural Sciences and Engineering Research Council of Canada, Calcul
Québec, Compute Ontario, Canada Foundation for Innovation, WestGrid, and
Compute Canada; Denmark – Villum Fonden, Carlsberg Foundation, and European
Commission; New Zealand – Marsden Fund; Japan – Japan Society for Promotion of
Science (JSPS) and Institute for Global Prominent Research (IGPR) of Chiba
University; Korea – National Research Foundation of Korea (NRF); Switzerland –
Swiss National Science Foundation (SNSF); United Kingdom – Department of
Physics, University of Oxford.