#  Spatial Reuse in Dense Wireless Areas: A Cross-layer Optimization Approach
via ADMM

Haleh Tabrizi, Borja Peleato, Golnaz Farhadi, John M. Cioffi, and Ghadah
Aldabbagh H. Tabrizi is with the Department of Electrical Engineering,
Stanford University, Stanford, CA, 94305 USA, e-mail:
(htabrizi@stanford.edu)B. Peleato is with the Department of Electrical and
Computer Engineering, Purdue University, West Lafayette, IN, 47907 USA,
e-mail: (bpeleato@purdue.edu).J. Cioffi is with the Department of Electrical
Engineering, Stanford University, Stanford, CA, 94305 USA, and Department of
Computer Science, King Abdulaziz University, Jeddah, Saudi Arabia, e-mail:
(cioffi@stanford.edu).G. Farhadi is with Fujitsu Labs of America, Sunnyvale,
CA, 94085 USA, e-mail: (gfarhadi@us.fujitsu.com).G. Aldabbagh is with the
Department of Computer Science, KAU, Jeddah, Saudi Arabia, e-mail:
(galdabbagh@kau.edu.sa)This paper has been submitted in part to IEEE global
communications conference, Atlanta, GA, December 2013. ©2015 IEEE. Personal
use of this material is permitted. Permission from IEEE must be obtained for
all other uses, in any current or future media, including
reprinting/republishing this material for advertising or promotional purposes,
creating new collective works, for resale or redistribution to servers or
lists, or reuse of any copyrighted component of this work in other works. DOI:
10.1109/TWC.2015.2464373

###### Abstract

This paper introduces an efficient method for communication resource use in
dense wireless areas where all nodes must communicate with a common
destination node. The proposed method groups nodes based on their distance
from the destination and creates a structured multi-hop configuration in which
each group can relay its neighbor’s data. The large number of active radio
nodes and the common direction of communication toward a single destination
are exploited to reuse the limited spectrum resources in spatially separated
groups. Spectrum allocation constraints among groups are then embedded in a
joint routing and resource allocation framework to optimize the route and
amount of resources allocated to each node. The solution to this problem uses
coordination among the lower-layers of the wireless-network protocol stack to
outperform conventional approaches where these layers are decoupled.
Furthermore, the structure of this problem is exploited to obtain a semi-
distributed optimization algorithm based on the alternating direction method
of multipliers (ADMM) where each node can optimize its resources independently
based on local channel information.

###### Index Terms:

Alternating direction method of multipliers (ADMM), crosslayer optimization,
dynamic resource allocation, routing.

##  I Introduction

With the remarkable growth of wireless technology and advancement of mobile
devices, the demand for more efficient use of limited radio resources is
increasing rapidly. Optimizing performance through sophisticated physical-
layer techniques is no longer enough, so studies have considered coordinating
multiple users through scheduling and admission control, as well as joint
operation of the physical and protocol layers in wireless networks. Lin et al.
summarize such cross-layer optimization techniques in [1].

This paper studies efficient resource use in dense wireless areas where all
nodes communicate with a single destination node. Dense wireless areas are
locations populated with many radio devices that must communicate over a
shared spectrum. Some examples include concert halls and stadiums populated
with user handheld devices, which must all communicate with a common access
point or base-station (BS). Home Area Networks (HAN) are another example of
dense wireless areas, particularly when such residential local area networks
may eventually connect hundreds of digital devices within a home in the so-
called “Internet of things”.

Figure 1 suggests a configuration that can increase spectrum efficiency by
relaying data over multiple hops rather than direct communication. The node
configuration motivates a hierarchical multi-hop architecture, where spectrum
can be reused among the hops as well as among cells [2, Ch. 15]. Furthermore,
the presence of a large number of nodes in a dense wireless area offers many
possible routes that the proposed multi-hop configuration can optimize. This
approach does not require additional infrastructure, but the nodes act as
relays for their neighbor’s data.

Figure 1: Multiple hops and reuse between hops.

In this scenario all nodes communicate with a fixed destination node, hence
the reuse pattern and flow of data can be structured. The nodes that are in
close proximity of each other are grouped together to form an FDMA (frequency
division multiple access) system where each node operates on different
spectra. Each such group of nodes can relay data for neighboring groups111A
configurable MAC (e.g. platform proposed in [3]) may be used to enable a
different channel access method (random access, FDMA, etc.) for different
nodes depending on their congestion level along the path. along an optimal
routing path determined by the capacity of the links between nodes. The
capacity of each link is determined by the power and bandwidth allocated to
that link, which in turn depend on the total spectrum allocated to each group.
Finding the optimal network configuration therefore requires a joint
optimization considering all constraints.

Joint optimization and coordination of the different network layers is not a
new concept. It has been used in many different applications in wireless and
wireline networks [4], [5]. Simultaneous routing and resource allocation
methods for orthogonal multiple access systems have been proposed by [6] (FDMA
and TDMA) and [7] (CDMA). Both approaches jointly optimize network-layer
routing and physical-layer resource allocation through a convex formulation of
the problem. A joint scheduling, routing, and resource allocation for multi-
hop cellular networks with relay stations is proposed by [8]. In this case,
the general problem is non-convex, and hence an iterative algorithm is used to
coordinate the layers. Beyond resource allocation and relay selection, [9]
jointly optimizes the relay strategy that should be used between each source
and destination pair based on a cross-layer design.

This paper addresses an extended version of the simultaneous routing and
resource allocation (SRRA) framework proposed by [6], including additional
constraints on the total available spectrum resources and spectrum reuse among
the predetermined user groups. The joint optimization of routing and global
resource allocation in a dense area is a non-convex problem that can only be
solved exactly through exhaustive search methods, such as branch and bound.
Unfortunately, these methods are often slow and have exponential worst-case
performance [10]. Hence, we propose simplifying the problem into a convex form
that can be solved using ADMM (Alternating Direction Method of Multipliers),
which has $O(1/k)$ convergence rate [11].

The predetermined user groups are selected based on each node’s relative
distance to the destination and a subset of all possible links connecting
nodes from one group to a neighboring group are selected as candidate data
flow links in the optimization problem. The goal is to achieve jointly
optimized routing, _global_ spectrum assignment, and local resource allocation
for the selected links through convex optimization.

Distributed implementations of network utility maximization problems often
lead to simpler subproblems, each optimizing a subset of the decision
variables based on local information. Palomar et al. [12] propose a systematic
framework to exploit decomposition structures that lead to different
distributed algorithms. A very popular approach is using dual decomposition to
generate a highly distributed optimization problem [6], but some studies seek
faster convergence using Newton’s method [13]. Certainly, there is a tradeoff
between the amount of message-passing among computing entities, convergence
speed, and computational complexity of each distributed algorithm. This paper
exploits the structure in the proposed joint optimization framework to
generate two different distributed algorithms: The first algorithm
(centralized) separates the problem into two subproblems, allowing decoupled
optimization of network-layer and physical-layer variables. The second
algorithm (semi-distributed) employs all active nodes as computing resources
that perform parallel optimization based on local channel information.

The proposed decompositions use the alternating direction method of
multipliers (ADMM) [11] algorithm, which is more numerically stable than
conventional decomposition methods. For example, ADMM does not require strict
convexity of the objective, but the dual decomposition method proposed in [6]
does. ADMM has been widely employed in producing distributed algorithms. Some
recent examples include multicell coordinated beamforming [14] and distributed
model predictive control [15].

The main contributions of this paper are: 1) devising a hierarchical structure
of nodes that allows resource reuse among spatially separated nodes, 2)
integrating global resource allocation and reuse in a cross-layer optimization
framework of routing and resource allocation, and 3) developing a semi-
distributed optimization algorithm based on ADMM where each node optimizes the
resources based on its local channel information.

The rest of this paper is organized as follows. Section II describes the
system model, consisting of the network topology and the physical layer
constraints in a dense wireless area. Section III formulates the cross-layer
optimization problem for the previous model in a convex form. Section IV
proposes a centralized and a distributed algorithm to solve this problem, both
based on ADMM. It also provides a numerical example to illustrate the
performance and convergence of the algorithm for a small network. Finally,
Section V presents more detailed simulations evaluating the performance of the
algorithm on larger networks and Section VI concludes and summarizes the
paper.

##  II Problem Formulation

The system model and the proposed hierarchical node configuration dictate a
set of system constraints that need to be satisfied while optimizing network
performance. This section formulates and motivates these constraints.

###  II-A Network Flow Model

We consider $N$ active users/nodes randomly distributed in a given area that
need to communicate simultaneously with a common destination node, such as a
cellular base station (BS). This paper investigates uplink transmission, but
downlink communication can be formulated in a similar manner with minor
modifications. A standard directed graph is used to represent the network
topology. The graph is assumed to be connected, i.e., there is a route between
every node and the BS. The nodes are labeled as
$n\in\mathcal{N}=\left\\{1,2,...,N\right\\}$, with node $1$ representing the
BS. All nodes can transmit, receive and relay data over the existing links.
Transmission and reception occur on disjoint spectra. A link exists between
two nodes $i$ and $j$ if direct communication among the two is possible. With
a total of $L$ links in the network, denote the set of all links by
$\mathcal{L}=\left\\{1,2,...,L\right\\}$, where each link is identified with
an integer value between $1$ and $L$.

The network topology is represented by a link-node incidence matrix
$A\in\mathbb{R}^{N\times L}$, whose entries determine a link’s source and
destination nodes by $+1$ and $-1$, respectively. For example, for the
configuration shown in Figure 2, where $N=4$ and $L=3$, the link-incidence
matrix is given by: $A=\Bigg{(}\begin{smallmatrix}-1&0&0\\\ +1&-1&-1\\\
0&+1&0\\\ 0&0&+1\end{smallmatrix}\Bigg{)}$. Each column of $A$ corresponds to
a link and each row corresponds to a node. The first column indicates that
link $1$ goes from node $2$ to node $1$. The set of outgoing links from node
$n$ are denoted by $\mathcal{O}\left(n\right)$ and the set of incoming links
are denoted by $\mathcal{I}\left(n\right)$.

Figure 2: Simple node and link example.

A multi-commodity flow model [16] that satisfies flow conservation at each
node is considered. For the uplink case considered here, $x\in\mathbb{R}^{L}$
identifies the flow on each link (data rate) that is destined to node $1$
(BS). Let $A_{nl}$ denote the element on the $n$-th row and $l$-th column of
$A$, and let $r\in\mathbb{R}^{N}$ be a vector of data traffic injected in the
network at each node. Then $Ax=r$, where $r_{n}=\sum_{l}{A_{nl}x_{l}}$,
denotes the amount of data traffic injected into the network at node $n$.
Positive components will indicate flow sources, while negative ones will
indicate flow sinks. For the example in Figure 2, the flow conservation law
indicates that $r_{1}=-x_{1}$, $r_{2}=x_{1}-x_{2}-x_{3}$, $r_{3}=x_{2}$, and
$r_{4}=x_{3}$. Furthermore, the amount of flow on each link ($x_{l}$) is
constrained by that link’s physical layer capacity $c_{l}$:

$x_{l}\leq c_{l},\quad\forall l\in\mathcal{L}.$

###  II-B Physical Layer Model

An FDMA channel access method is considered, where each node’s outgoing links
are assigned disjoint frequency bands. The received signal is corrupted by
additive white Gaussian noise with power spectral density $\sigma_{l}$ and the
link capacity is given by a concave and monotone increasing function of its
allocated bandwidth $w_{l}$ and power $p_{l}$, as follows:

$c_{l}=w_{l}\log\left(1+\frac{p_{l}q_{l}}{w_{l}\sigma_{l}}\right),\quad\forall
l\in\mathcal{L}.$

Parameter $q_{l}$ is the channel gain from the source to the destination of
link $l$. The total amount of power allocated to a node $n$ is constrained to
be smaller than a pre-fixed value $P_{\max,n}$, hence

$\sum_{l\in\mathcal{O}(n)}{p_{l}\leq P_{\max,n}},\quad\forall
n\in\mathcal{N}.$

The total amount of bandwidth allocated to node $n$ is denoted by

$\sum_{l\in\mathcal{O}(n)}{w_{l}}=v_{n},\quad\forall n\in\mathcal{N}.$

If necessary, after obtaining the continuous bandwidth variables $w_{l}$, they
can be quantized based on the underlying system modulation approach. For
example, if the underlying system is an OFDMA system, the bandwidth values can
be rounded down to an integer number of subcarrier bandwidths and if the
underlying system is an LTE system, it can be rounded to the closest integer
number of resource blocks.

###  II-C Group Assignment

The objective of the proposed multi-hop configuration is to reuse as much of
spectrum as possible. Hence, direct links between far nodes and the BS are
divided into multiple hops that operate on different spectrum bands, allowing
reuse (Figure 1). To this end, nodes are grouped based on their distance to
the BS as shown in Figure 3: nodes located between radial distances $d_{g-1}$
and $d_{g}$ from the BS belong in group $\mathcal{G}_{g}$, and all nodes
within distance $d_{1}$ of the destination belong to group $\mathcal{G}_{1}$.
The coverage area is thus partitioned into $M$ disjoint sections such that
$\mathcal{G}_{g}\bigcap\mathcal{G}_{g^{\prime}}=\emptyset,\quad\mathcal{G}_{g}\neq\emptyset,\quad\forall
g,g^{\prime}$, and each node belongs to a group:
$\bigcup_{g=1}^{M}\mathcal{G}_{g}=\mathcal{N}$.

In general, the closer a link to the BS, the larger the flow it will be
carrying. The total flow on the links that terminate at the BS is the largest
and the flow decreases as the links move away from the BS as shown in Figure
3. Let $f$ denote the frequency reuse factor, such that every $(f+1)$-th group
uses the same spectrum set ($f>2$ to avoid the hidden node problem). As such,
the maximum available spectrum $W_{\max}$ is divided among groups $1$ through
$f$:

$\sum_{g=1}^{f}{W_{g}}=W_{\max},$

where $W_{g}$ denotes the amount of bandwidth allocated to group $g$. The rest
of the groups ($f+1$ through $M$) reuse the spectrum allocated to groups $1$
through $f$:

$W_{g}=W_{\left((g-1)\bmod{f}\right)+1},\;\forall g>f.$

For example, if $f=3$, the total bandwidth is divided between groups $1,2$ and
$3$, while groups $4,7,10,\ldots$ use the same spectrum as group $1$; groups
$5,8,11,\ldots$ use the same spectrum as group $2$; and groups $6,9,12,\ldots$
use the same as group $3$.

###  II-D Routing Link Assignments

In order to optimize the routing from each node to the BS, an initial set of
possible links should be determined. This subsection describes a method,
similar to that in [17], for selecting this initial pool of links (or refining
a pre-defined set of links if desired). The method assumes that the BS knows
the location of each node and uses relative distances between them to select
the initial pool of links. If this assumption does not hold, the initial pool
of links could include all those between nodes in consecutive groups.
Regardless of how the initial pool of links is selected, resource and flow
allocation are based on channel information.

The proposed method will include a link from node $i$ to node $j$ in the
initial pool if:

  1. 1.

$j\in\mathcal{G}_{g-1}$, where $i\in\mathcal{G}_{g}$.

  2. 2.

the angular distance between nodes $i$ and $j$ is less than a pre-fixed
$\theta$, as shown in Figure 3.

  3. 3.

the distance between nodes $i$ and $j$ is less than a pre-fixed $d_{th}$.

The first condition enforces a single hop between consecutive groups to allow
spectrum reuse. If $i\in\mathcal{G}_{g}$ and $j\in\mathcal{G}_{h}$ with
$h<g-1$, then the link from $i$ to $j$ traverses more than one group (spectrum
region) and defeats the purpose of reuse. On the other hand, if $h=g$ then $i$
and $j$ belong to the same group creating an extra hop within the group, and
if $h>g$, then $j$ is farther away from the BS than $i$ already is. In both
cases, $j$ has no relaying benefit for $i$.

The second and third conditions ensure that the relay node $j$ is not too far
away from the source node $i$. Parameters $d_{th}$ and $\theta$ control the
number of outgoing links from a node and allow selecting a reasonable (non-
empty) set of possible routes.

Figure 3: Link assignment conditions.

The maximum transmission power per unit bandwidth per link is limited to avoid
interference between groups operating on the same spectrum. This maximum
transmission power depends on the radial distances $d_{g}$,
$g=\left\\{1,2,...,M\right\\}$ defining the groups and on the frequency reuse
factor $f$. If every group has the same radial distance $d$, the minimum
possible distance between the transmitter of one link and a receiver suffering
its interference is equal to $d_{int}=(f-2)\times d$. The received power per
unit bandwidth has to be less than a fraction $\alpha$ of the noise power for
interference to be negligible, hence the constraint on the transmission power
is:

$p_{l}\leq w_{l}\frac{\alpha N_{0}}{K(\frac{l_{0}}{d_{int}})^{a}},$

where $a$ is the pathloss exponent, $l_{0}$ is the reference distance222This
is a conservative bound, using the shortest possible distance between
interfering nodes. Results could be improved using the actual distance for
each node, if static and known., $N_{0}$ is the noise power spectral density
at the receiver ($\sigma_{l}$), and $K$ is the attenuation factor based on the
simplified pathloss model [2]. This work denotes the bandwidth coefficient by

$\gamma=\alpha
N_{0}\left(K\left(\frac{l_{0}}{d_{int}}\right)^{a}\right)^{-1},$

such that the per link transmission power constraint becomes $p_{l}\leq
w_{l}\gamma$.

##  III Joint Routing, Global and Local Resource Allocation

Based on the network conditions, such as node density and locations, the user
groups and possible routing links are determined. Given these, the routing and
resource allocation problem can be formulated as a convex optimization problem
with the objective of maximizing a concave utility function. To maintain
fairness among all users, the utility function was chosen to be the minimum
rate among all users ($n\in\mathcal{N}\setminus\left\\{1\right\\}$). Only
uplink communication is considered here, so the set of all transmitters (all
nodes other than the BS) is
$\mathcal{\tilde{N}}=\mathcal{N}\setminus\left\\{1\right\\}$.

For each link $l\in\left\\{1,..L\right\\}$ in the initial pool, we define the
system variables $p_{l}$ and $w_{l}$, respectively denoting the amount of
power and spectrum allocated to the link, and $x_{l}$, denoting the amount of
flow on the link destined to the BS. $W_{g},\;\forall
g\in\left\\{1,...,M\right\\}$ will denote the amount of spectrum allocated to
each group. The problem then becomes:

$\displaystyle\underset{p,w,x,(r_{n},W_{g})}{\text{maximize}}$
$\displaystyle\min_{n\in\mathcal{\tilde{N}}}{r_{n}}$ (1) subject to
$\displaystyle\sum_{l}{A_{nl}x_{l}}=r_{n},\;\forall n\in\mathcal{N}$ (1a)
$\displaystyle x_{l}\leq
w_{l}\log\left(1+\frac{p_{l}q_{l}}{w_{l}N_{0}}\right),\;\forall
l\in\mathcal{L}$ (1b) $\displaystyle p_{l}\leq w_{l}\gamma,\;\forall
l\in\mathcal{L}$ (1c) $\displaystyle\sum_{l\in\mathcal{O}(n)}{p_{l}\leq
P_{\max}},\;\forall n\in\mathcal{\tilde{N}}$ (1d)
$\displaystyle\sum_{n\in\mathcal{G}_{g}}{\sum_{l\in\mathcal{O}(n)}{w_{l}}}\leq
W_{g},\;\forall g\in\left\\{1,2,...,M\right\\}$ (1e)
$\displaystyle\sum_{g=1}^{f}{W_{g}}=W_{\max}$ (1f) $\displaystyle
W_{g}=W_{\left((g-1)\bmod{f}\right)+1},\;\forall g>f$ (1g) $\displaystyle
x_{l}\geq 0,\;p_{l}\geq 0,\;w_{l}\geq 0,\;\forall l\in\mathcal{L}.$ (1h)

Constraints (1a) through (1c) are as explained in Section II. Constraint (1d)
is obtained by setting the maximum transmission power of all users equal to
$P_{\max}$. Constraints (1e) through (1g) identify the spectrum reuse
constraints introduced in Section II. The variables $W_{g}$ and $r_{n}$ are
auxiliary variables introduced to make the problem formulation easier to
understand. The value for $r_{n}$ can be derived from $x$ using equations
(1a). Similarly, it is easy to prove that constraint (1e) will be tight in the
optimal solution, so the value for $W_{g}$ can also be uniquely derived from
$w$. Consequently, both $r_{n}$ and $W_{g}$ appear between brackets in the
list of optimization variables. This notation for auxiliary variables
continues throughout the rest of the paper.

The constraints in (1) define a convex set and the objective function is
concave. Hence, the problem is convex and can be solved globally and
efficiently by interior-point methods [18]. However, there might be multiple
solutions with the same objective value. Adding a small negative factor of
power sum, $-\epsilon\sum{p_{l}}$, to the objective would guarantee that a
solution with small power values is selected among the optimal set. This
penalty term has negligible effect on the resulting minimum rate value
$r_{n}$.

The conventional direct-mode resource-allocation problem can be derived as a
specific scenario of the multi-hop optimization problem (1). In this
formulation, there is only one group of nodes, and the only possible route
from each node to the destination node is a direct link between the node and
the destination. In this case, the problem can be rewritten as:

$\displaystyle\underset{p,v,(r_{n})}{\text{maximize}}$
$\displaystyle\min_{n\in\mathcal{\tilde{N}}}{r_{n}}$ (2) subject to
$\displaystyle 0\leq r_{n}\leq
v_{n}\log\left(1+\frac{p_{n}q_{n}}{v_{n}N_{0}}\right),\;\forall
n\in\mathcal{\tilde{N}}$ $\displaystyle 0\leq p_{n}\leq P_{\max},\;\forall
n\in\mathcal{\tilde{N}}$
$\displaystyle\sum_{n\in\mathcal{\tilde{N}}}{v_{n}}=W_{\max},\;$
$\displaystyle v_{n}\geq 0,\;\forall n\in\mathcal{\tilde{N}},$

where $v_{n}$ and $p_{n}$ are equivalent to the original parameters $w_{l}$
and $p_{l}$, respectively, since there is only one outgoing link from each
node.

##  IV Distributed Optimization

Dual decomposition is an old method for solving a convex optimization problem
in a (semi-) decentralized manner333Contrary to popular belief, dual
decomposition does not provide a fully distributed solution. There still needs
to be a central entity gathering information from all nodes. [19]. Generally,
if the objective function is separable in its variables, the problem can be
split into smaller subproblems. These subproblems can then be solved
iteratively in parallel based on dual ascent (or descent). However, this
method requires some hard assumptions such as strict convexity and finiteness
of the objective function; otherwise it can result in numerical instability.

Augmented Lagrangian methods such as the method of multipliers achieve better
convergence results by adding a penalty parameter to the dual ascent
objective; however, the penalty parameter is not separable when the original
problem is. The alternating direction method of multipliers (ADMM), on the
other hand, is an algorithm that combines the decomposability property of dual
ascent with the robust convergence properties of the method of multipliers
[11].

The objective function in problem (1) is not strictly convex so dual
decomposition cannot be applied directly. A weighted quadratic regularization
term would make the objective strictly convex, but finding the appropriate
weights is a challenging problem. If the weights are too small, dual
decomposition will still suffer from numerical stability, while excessive
weights would involve an optimality cost. ADMM is in general more numerically
stable and faster in convergence than the conventional dual decomposition
method [14]. This paper thus proposes ADMM to implement distributed cross-
layer optimization algorithms for spatial reuse.

###  IV-A Disjoint Network-Layer and Physical-Layer Optimization

One common method for developing efficient algorithms in cross-layer
optimization problems is to decouple the network-layer and physical-layer
problems [6], [12]. This method divides the problem into two simpler
subproblems, the optimization of network-layer and physical-layer variables,
which can be solved separately and updated iteratively by exchanging messages
between the two layers. This section takes a similar approach in decoupling
the problems through ADMM. At each iteration $k\in\mathbb{Z}$, the network
layer and physical layer subproblems are solved in parallel and based on the
results, the dual variables are updated.

The network-layer parameters $x$ and the physical-layer parameters $w$ and $p$
are coupled by constraint (1b) in the above centralized optimization problem.
To decompose the problem, a new auxiliary variable $t\in\mathbb{R}^{L}$ and
$L$ equality constraints $t_{l}=x_{l},\;\forall l$ are introduced. ADMM
decouples constraints $t_{l}=x_{l},\forall l$ through a dual variable
$u\in\mathbb{R}^{L}$, and decomposes problem (1) into two subproblems
corresponding to the network and physical layers.

The augmented Lagrangian for the extended problem with penalty parameter
$\rho>0$ and dual variable $y$ is

$\mathfrak{L}=\min_{n\in\mathcal{\tilde{N}}}\left(r_{n}\right)-y^{T}(x-t)-\frac{\rho}{2}\|x-t\|^{2}_{2},$

subject to constraints (1a)-(1h), but it is often more convenient to scale the
dual variable $y$ and combine the linear and quadratic terms, applying ADMM to

$\tilde{\mathfrak{L}}=\min_{n\in\mathcal{\tilde{N}}}\left(r_{n}\right)-\frac{\rho}{2}\|x-t+u\|^{2}_{2},$

where $u=\rho^{-1}y$ and still subject to constraints (1a)-(1h). Both
formulations are equivalent, as shown in [11] (section 3.1.1). The network-
layer optimization problem is then:

$\displaystyle\underset{x\geq 0,(r_{n})}{\text{maximize}}$
$\displaystyle\min_{n\in\mathcal{\tilde{N}}}\left(r_{n}\right)-\frac{\rho}{2}\|x-t+u\|^{2}_{2}$
(3) subject to $\displaystyle\sum_{l}{A_{nl}x_{l}}=r_{n},\;\forall
n\in\mathcal{N},$

and the physical layer problem is:

$\displaystyle\underset{p,w,t\geq 0,(W_{g})}{\text{maximize}}$
$\displaystyle-\frac{\rho}{2}\|x-t+u\|^{2}_{2}$ (4) subject to $\displaystyle
t_{l}\leq w_{l}\log\left(1+\frac{p_{l}q_{l}}{w_{l}N_{0}}\right),\;\forall
l\in\mathcal{L}$
$\displaystyle(\ref{powerConstraint})-(\ref{individualBwConstraints}).$

Let superscripts denote iteration indices, omitted in the above subproblems
for clarity. At iteration $k=1$, the network layer parameters $x^{(1)}$ are
obtained from (3) using arbitrary initial values for the dual variables
$u^{(0)}$ and auxiliary variables $t^{(0)}$. The physical layer parameters
$p^{(1)}$, $w^{(1)}$ and $t^{(1)}$ are then obtained from (4) using the same
$u^{(0)}$ and the obtained values for $x^{(1)}$. Finally, the dual variable
$u^{(1)}$ is updated based on $u^{(0)}$, $t^{(1)}$, and $x^{(1)}$ as follows:

$u^{(k+1)}=u^{(k)}+x^{(k+1)}-t^{(k+1)}.$ (5)

ADMM has guaranteed convergence for any $\rho>0$ but in practice the speed of
convergence can change depending on its value. Unfortunately, there is no way
to know in advance which value will yield the fastest convergence. Convergence
results of this algorithm are investigated in Section IV-D through an example.

###  IV-B Semi-Distributed Device Optimization Algorithm

The authors are not aware of any fully distributed algorithm that guarantees
convergence to the global optimum using only local coordination among nodes.
In any case, such an algorithm would probably require a very large number of
iterations, since the constraints in problem (1) are not local. Information
such as data rate from each node needs to reach every potential relay down the
chain, then travel back and forth until they converge on a solution. Instead,
a semi-distributed algorithm is proposed where each node independently
optimizes its own resources based on local channel information and limited
communication with a central entity (such as the BS). The algorithm uses ADMM
to decompose problem (1) into an independent subproblem for each node and a
shared dual update. One major benefit of the semi-distributed algorithm is
that the algorithm can adapt to changing local channel conditions during the
iterations of solving the global problem. In the centralized approach,
however, the final system parameters are obtained based on the initial channel
information sent to the central unit444The addition energy consumption for
message-passing in the semi-distributed approach and the performance gain that
can be obtained by updating local channel information are very system
dependent and are out of the scope of this work..

In the semi-distributed problem there are two sets of coupling variables: the
per-link flows $x_{l}$ and the per-node bandwidths $v_{n}$. To decouple the
per-node optimization subproblems from the global problem, auxiliary sets of
variables $t_{l},\;\forall l\in\mathcal{L}$ and $b_{n},\;\forall
n\in\mathcal{\tilde{N}}$ are introduced along with the following equality
constraints:

$\displaystyle t_{l}=x_{l},\;\forall l\in\mathcal{L}$ (6a) $\displaystyle
b_{n}=v_{n},\;\forall n\in\mathcal{\tilde{N}}.$ (6b)

Even though it appears that the problem complexity increases with the
increased number of constraints and variables, it will be shown that such
decoupling generates simpler convex subproblems. The set of original
variables, $\psi\in\mathbb{R}^{L+N-1}$, and the set of new variables,
$z\in\mathbb{R}^{L+N-1}$, are defined as:
$\psi=[x_{1},...,x_{L},v_{2},...,v_{N}]$,
$z=[t_{1},...,t_{L},b_{2},...,b_{N}]$. The two variables $v_{1}$ and $b_{1}$,
which correspond to the destination node bandwidth, are eliminated.

The variable $\xi\in\mathbb{R}^{L+N-1}$ represents the set of dual variables
$\xi=(u,y)$, such that $u\in\mathbb{R}^{L}$ corresponds to the routing
variables and $y\in\mathbb{R}^{N-1}$ corresponds to the bandwidth variables.
The Lagrangian then becomes

$\tilde{\mathfrak{L}}=\min_{n\in\mathcal{\tilde{N}}}\left(r_{n}\right)-\frac{\rho}{2}\|\psi-z+\xi\|^{2}_{2},$

subject to constraints (1a)-(1h). With this notation and using ADMM, problem
(1) is decomposed into two subproblems. The first subproblem optimizes routing
and the amount of bandwidth allocated to each user:

$\displaystyle\underset{\psi=[x,v]\geq 0,(r_{n},W_{g})}{\text{maximize}}$
$\displaystyle\min_{n\in\mathcal{\tilde{N}}}\left(r_{n}\right)-\frac{\rho}{2}\|\psi-z+\xi\|^{2}_{2}$
(7) subject to $\displaystyle\sum_{l}{A_{nl}x_{l}}=r_{n},\;\forall
n\in\mathcal{N}$ $\displaystyle\sum_{n\in\mathcal{G}_{g}}{v_{n}}\leq
W_{g},\;g\in\left\\{1,2,...,M\right\\}$
$\displaystyle(\ref{bwConstraint})-(\ref{individualBwConstraints}).$

The second subproblem determines variables $z$, per link power $p_{l}$, and
bandwidth $w_{l},\;\forall l$:

$\displaystyle\underset{p,w,z=[t,b]\geq 0}{\text{maximize}}$
$\displaystyle-\frac{\rho}{2}\|\psi-z+\xi\|^{2}_{2}$ (8) subject to
$\displaystyle t_{l}\leq
w_{l}\log\left(1+\frac{p_{l}q_{l}}{w_{l}N_{0}}\right),\;\forall
l\in\mathcal{L}$
$\displaystyle\sum_{l\in\mathcal{O}(n)}{w_{l}}=b_{n},\;\forall
n\in\mathcal{\tilde{N}}$
$\displaystyle(\ref{powerConstraint})-(\ref{PmaxConstraint}).$

Subproblems (7) and (8) are then solved iteratively with dual updates:

$\xi^{(k+1)}=\xi^{(k)}+\psi^{(k+1)}-z^{(k+1)}.$ (9)

As demanded, problem (8) can now be decomposed into $N-1$ separate
subproblems, such that node $n$ solves the following problem based on its
local channel information $q_{l},\forall l\in\mathcal{O}(n)$:

$\displaystyle\underset{p_{l},w_{l},t_{l},\forall
l\in\mathcal{O}(n),b_{n}}{\text{maximize}}$
$\displaystyle-\frac{\rho}{2}\sum_{l\in\mathcal{O}(n)}{\|x_{l}-t_{l}+u_{l}\|^{2}_{2}}$
(10) $\displaystyle-\frac{\rho}{2}\|v_{n}-b_{n}+y_{n}\|^{2}_{2}$ subject to
$\displaystyle t_{l}\leq
w_{l}\log\left(1+\frac{p_{l}q_{l}}{w_{l}N_{0}}\right),\;\forall
l\in\mathcal{O}(n)$ (10a) $\displaystyle\sum_{l\in\mathcal{O}(n)}{p_{l}\leq
P_{\max}}$ (10b) $\displaystyle p_{l}\leq w_{l}\gamma,\;\forall
l\in\mathcal{O}(n)$ $\displaystyle\sum_{l\in\mathcal{O}(n)}{w_{l}}=b_{n}$
$\displaystyle p_{l},w_{l},t_{l}\geq 0,\;\forall l\in\mathcal{O}(n).$

Furthermore, subproblem (7) is separable into bandwidth $v$ and flow variables
$x$. The objective in (7) can be rewritten as
$\min_{n\in\mathcal{\tilde{N}}}\left(r_{n}\right)-\frac{\rho}{2}\|x-t+u\|^{2}_{2}-\frac{\rho}{2}\|v-b+y\|^{2}_{2}$.
Hence this subproblem can further be decomposed into routing and spectrum
assignment problems. The routing subproblem is:

$\displaystyle\underset{x\geq 0,(r_{n})}{\text{maximize}}$
$\displaystyle\min_{n\in\mathcal{\tilde{N}}}\left(r_{n}\right)-\frac{\rho}{2}\|x-t+u\|^{2}_{2}$
(11) subject to $\displaystyle\sum_{l}{A_{nl}x_{l}}=r_{n},\;\forall
n\in\mathcal{N}.$

The spectrum assignment subproblem is simply the Euclidean projection onto a
constrained set:

$\displaystyle\underset{v\geq 0,(W_{g})}{\text{maximize}}$
$\displaystyle-\frac{\rho}{2}\|v-b+y\|^{2}_{2}$ (12) subject to
$\displaystyle\sum_{n\in\mathcal{G}_{g}}{v_{n}}\leq
W_{g},\;g\in\left\\{1,2,...,M\right\\}$
$\displaystyle(\ref{bwConstraint})-(\ref{individualBwConstraints}).$

Hence problems (11) and (12) can be solved in parallel and even by two
independent entities. The dual variables $u$ and $y$ can also be updated
separately as:

$\displaystyle u^{(k+1)}$ $\displaystyle=$ $\displaystyle
u^{(k)}+t^{(k+1)}-x^{(k+1)}$ (13) $\displaystyle y^{(k+1)}$ $\displaystyle=$
$\displaystyle y^{(k)}+b^{(k+1)}-v^{(k+1)}.$ (14)

The primal residuals corresponding to link flows are
$h_{1}^{(k)}=x^{(k)}-t^{(k)}$ and the primal residuals corresponding to
bandwidth are $h_{2}^{(k)}=v^{(k)}-b^{(k)}$. The dual residuals are then
defined as $s_{1}^{(k)}=-\rho(t^{(k)}-t^{(k-1)})$ and
$s_{2}^{(k)}=-\rho(b^{(k)}-b^{(k-1)})$.

Various criteria can be considered as stopping points for the ADMM method.
According to [11], one possibility could be:

$\displaystyle\|h_{1}^{(k)}\|_{2}\leq\epsilon_{1}^{p},\;\|h_{2}^{(k)}\|_{2}\leq\epsilon_{2}^{p},$
(15)
$\displaystyle\|s_{1}^{(k)}\|_{2}\leq\epsilon_{1}^{d},\;\|s_{2}^{(k)}\|_{2}\leq\epsilon_{2}^{d}$

for small predefined $\epsilon_{1}^{p}$, $\epsilon_{2}^{p}$,
$\epsilon_{1}^{d}$, and $\epsilon_{2}^{d}$ values.

###  IV-C Implementation

This section focuses on the implementation of the semi-distributed device
optimization algorithm introduced above and investigates the required message
passing between different entities. All the information regarding channels,
network topology, and other such infrastructure can be communicated to the
nodes by the central unit (base-station) through control channels. In cellular
networks, signaling for link establishment can be done via slight
modifications to the standard radio resource control (RRC) connection
reconfiguration procedure [20]. This procedure can include signaling to the
sink nodes of each hop indicating their role as well as providing a list of
source nodes that the sink node is providing wireless broadband service to.

Based on the optimization subproblems developed in the previous section, the
entity solving problem (11) (which can be the BS) provides each user $n$ with
the flows $x_{l},\;\forall l\in\mathcal{O}(n)$ that it should carry on each of
its outgoing links and the total amount of bandwidth $v_{n}$ that it can use.
Given $v_{n}$ and $x_{l}$, each node calculates $p_{l},w_{l}$, and
$t_{l},\;\forall l\in\mathcal{O}(n)$ by solving optimization problem (10).
Each node $n$ then broadcasts $t_{l},\;\forall l\in\mathcal{O}(n)$ and
$b_{n}=\sum_{l\in\mathcal{O}(n)}{w_{l}}$.

The primal variables $x$ and $v$ can be calculated at a single central entity
in parallel or by two independent entities. Denote by CU $1$, the central unit
that collects $t_{l}$ values from all nodes and solves problem (11) to obtain
$x$. A second central unit (CU $2$) is in charge of group bandwidths: it
gathers per user calculated bandwidths $b_{n}$ and solves (12) to obtain $v$.
After each iteration, the dual variables $u$ and $y$ are updated. The dual
variable updates can be performed either at the central units or at each node,
resulting in different message passing structures. Figure 4 shows the message-
passing when the updates are performed at the CU’s, and the corresponding
distributed algorithm is summarized in Algorithm 1.

Figure 4 shows the variables that are updated by each entity and the reference
number of the problem solved to obtain the corresponding variable. In this
configuration, based on initialized values of $u^{(0)}$ and $t^{(0)}$, CU $1$
solves optimization problem (11) to obtain $x^{(1)}$. Then the sum of the dual
variable and the updated primal variable ($u^{(0)}+x^{(1)}$) is broadcast.
Similarly, based on initialized values of $y^{(0)}$ and $b^{(0)}$, CU $2$
solves optimization problem (12) to obtain $v^{(1)}$. Then ($y^{(0)}+v^{(1)}$)
is broadcast. Each node $n$ captures its required information
($y_{n}^{(1)}+v_{n}^{(1)}$) and ($u_{l}^{(0)}+x_{l}^{(1)}$), $\forall
l\in\mathcal{O}(n)$, calculates $b_{n}^{(1)}$ and $t_{l}^{(1)},\;\forall
l\in\mathcal{O}(n)$ by solving problem (10), and broadcasts the updated
variables. During the next iteration, CU $1$ calculates $u^{(1)}$, and CU $2$
calculates $y^{(1)}$ based on equations (13) and (14), respectively. In this
manner, at each iteration, all the dual and primal parameters are updated.

Figure 4: $3$-level optimization and message passing among different entities.
Algorithm 1 Proposed distributed cross-layer optimization.

1:initialize $u_{l}$, $y_{n}$, $t_{l}$, $b_{n},\small{\forall l}$ and
$\small{\forall n}$

2:repeat

3: perform in parallel:

  * •

CU $1$: given $u_{l}$ and $t_{l}$, update $x_{l}$ (11) and broadcast
$x_{l}+u_{l},\small{\forall l}$.

  * •

CU $2$: given $y_{n}$ and $b_{n}$, update $v_{n}$ (12) and broadcast
$v_{n}+y_{n},\small{\forall n}$.

4: Each node $n$: given $x_{l}+u_{l},\;\small{\forall l\in\mathcal{O}(n)},$
$v_{n}+y_{n}$ and local channel information, update

5: $t_{l},\small{\forall l\in\mathcal{O}(n)}$ and $b_{n}$ (10) and broadcast
results.

6: CU $1$ and $2$: update corresponding dual variables $u_{l},\small{\forall
l}$ and $y_{n},\small{\forall n}$ (9).

7:until stopping criterion (15) met.

Each update step in Algorithm 1 is followed by broadcasting the updated
parameters. The dual variable update is performed at the CU’s, because there
are generally more computing resources at these locations as compared to
individual nodes. According to problem (10), the only information that each
node requires from the CU’s is the sum of $x_{l}$ and $u_{l}$ for all its
outgoing links and the sum of $v_{n}$ and $y_{n}$. Hence in order to reduce
the amount of message-passing between the nodes and CU’s, the sum of the dual
and primal variables $x_{l}+u_{l},\;\forall l$ and $w_{n}+y_{n},\;\forall n$
are broadcast by the CU’s instead of the separate variable values (Figure 4).
After convergence, the optimum routing variables $x$ and physical-layer
resources $p$ and $w$ are obtained.

In summary, the semi-distributed optimization algorithm works as follows: Each
node $n$ gathers and stores the values of $x_{l}+u_{l},\;\forall
l\in\mathcal{O}(n)$, $v_{n}+y_{n}$. Then it performs an Euclidean projection
onto a set of constraints (10), which has very small complexity, obtaining
$t_{l},\;\forall l\in\mathcal{O}(n)$, and $b_{n}$. Each node then broadcasts
its updated results. Two central units, CU $1$ and CU $2$, that can be the
same or independent entities, gather the updated values of $t_{l}$ and
$b_{n}$, respectively. CU $1$ and CU $2$ then solve (11) and (12),
respectively, in parallel. These problems also have very small complexity. The
results of such optimizations are $x_{l},\;\forall l\in\mathcal{L}$, and
$v_{n},\;\forall n\in\mathcal{\tilde{N}}$. CU $1$ and CU $2$ then update the
dual variables $u$ and $y$, respectively, while storing a local copy of each
variable for the next dual variable update.

In theory, the proof of convergence for ADMM requires that all nodes update
$t_{l},\;\forall l\in\mathcal{O}(n)$ and $b_{n}$ in step-3 before the
algorithm can proceed to step-4. In most practical cases, however, the
algorithm still converges even if some nodes occasionally fail to perform the
update. A good practical approach is to allow partial updates (only a subset
of nodes update their local variables) during the first few iterations, and
impose that all nodes perform the update in subsequent ones. This approach
guarantees convergence to the same solution, while allowing faster execution
of the algorithm in cases where some nodes are significantly slower than
others. Wei and Ozdaglar proposed in [21] an asynchronous distributed version
of ADMM that converges in $O(1/k)$, but it assumes that all entities send
updates infinitely often and only guarantees convergence with probability 1,
unlike regular ADMM. The number of iterations required with this asynchronous
version can also be significantly larger than with regular ADMM.

The complexity of Algorithm 1 scales gracefully to large number of nodes.
Step-2 consists of solving a simple quadratic program (11) and performing a
Euclidean projection onto the convex set specified in (12). Both of these can
be easily scaled to large number of variables. Step-3 is performed
independently by each node, so its complexity does not grow. Finally, Step-4
computes a linear combination of optimized variables, which again is a very
simple operation whose complexity increases linearly with the number of
variables. The amount of data broadcast by the nodes also increases linearly
with the number of nodes and links. In low-power large scale networks such as
wireless sensor networks (WSN), the message-passing overhead can be comparable
to the actual data size. The proposed algorithm can then demolish the low-
power objective of WSNs; however, it can prove beneficial in dense areas
populated with devices that have higher energy resources.

###  IV-D A Numerical Example

Figure 5 illustrates the algorithm operation over a circular cellular area of
radius $280$m. A total of $11$ nodes (with destination $N=12$) are randomly
distributed in this area, creating $6$ groups separated by $d=40$m (except
$d_{1}=2d$). We assume that a total bandwidth of $W_{\max}=10$MHz is available
with a carrier frequency of $800$Mhz. The receiver noise power spectral
density is $N_{0}=10^{-11}$W/MHz. The reference distance $l_{0}$ in constraint
(1c) is set to $1$ and the pathloss exponent is $a=4$. Each user’s maximum
transmission power is set to $P_{\max}=0.5$W.

Figure 5: Top: simple node configuration (N = 12), bottom: resulting optimum
routing configuration.

Setting the angular threshold $\theta=10^{o}$ and the distance threshold
$d_{th}=1.5d$ in the link assignment step, a total of $L=15$ initial links are
generated, as shown in the top plot of Figure 5. CVX, a package for specifying
and solving convex programs [22], is used to solve the centralized
optimization problem (1), resulting in a max-min rate of $208$Kbps. For the
disjoint network layer and physical layer optimization problems, the penalty
parameter is set to $\rho=1$ (for variations refer to [11, Ch 3.4]) and
routing parameters $t^{(0)}$ are initialized to $1$. Solving the individual
subproblems (3) and (4) iteratively via ADMM results in the same max-min rate
of $208$Kbps, as shown in Figure 7. The norms of the primal and dual residuals
at each iteration are depicted in Figure 6. The results show that the
centralized algorithm converges in about $15$ iterations.

Figure 6: Norms of primal and dual residuals versus iteration for layered
optimization problem. Figure 7: Max-min rate convergence for layered
optimization problem. Figure 8: Norms of primal and dual residuals versus
iteration for the semi-distributed optimization problem.

Figure 9: Max-min rate convergence for the semi-distributed optimization
problem for a time-invariant channel (left) and $N_{0}$ changed after $50$
iterations (right).

Before implementing the subproblems (10), (11), and (12) of the semi-
distributed algorithm in Section IV-B, it is important to investigate their
nature. Optimization problem (10) is performed by all the nodes (in parallel)
and has a quadratic objective and linear constraints other than constraint
(10a). By approximating the $\log$ function with a piecewise linear function,
problem (10) can be converted to a simple Quadratic Programming (QP) problem
with a relatively small number of variables. Then CVXGEN [23], which can
generate fast custom code for small QP-representable convex optimization
problems, can be used to solve the problem.

A sequential convex programming approach [24] is used to convert subproblem
(10) into a QP: the capacity term in (10a),
$c(w_{l},p_{l})=w_{l}\log(1+\frac{p_{l}q_{l}}{w_{l}N_{0}})$, is approximated
by a piece-wise linear function:
$\tilde{c}(\beta)=\min_{(a_{1},a_{2})\in\Sigma}\left\\{a_{1}^{T}\beta+a_{2}\right\\}$,
where $\beta=\left(p_{l},w_{l}\right)$ and $\Sigma$ represents a set of planes
tangent to $c(\beta)$. For simplicity, $\Sigma$ in this case consists of only
$10$ planes. Consequently, constraint (10a) is replaced with the following
linear constraint:

$t_{l}\leq\min_{(a_{1},a_{2})\in\Sigma}\left\\{a_{1}^{T}\beta+a_{2}\right\\},\;\forall
l\in\mathcal{O}(n)$

This approximation can be repeatedly re-adjusted during the execution of the
algorithm to choose planes that are tangent at the current iterate, and hence
the approximation error is progressively reduced. Eventually, when a point is
reached where planes are tangent at the iterate points and the iterates remain
constant (reached an equilibrium point), it is ensured that the optimal
solution with zero approximation error is achieved.

The number of variables in this subproblem depends on the number of outgoing
links from each node. Each node’s unknown parameters are $t_{l}$, $p_{l}$ and
$w_{l}$, and hence if node $n$ has $L_{n}$ outgoing links, it will have
$3L_{n}$ variables. Since the nodes in Figure 5 only have $1$ or $2$ outgoing
links, each node has a $3$ or $6$ variable QP problem to solve.

Subproblem (12) is readily in QP format and subproblem (11) can be easily
converted to a QP by replacing it with its epigraph form. The problem then
becomes

$\displaystyle\underset{x\geq 0,(\nu)}{\text{maximize}}$
$\displaystyle\nu-\frac{\rho}{2}\|x-t+u\|^{2}_{2}$ (11) subject to
$\displaystyle\sum_{l}{A_{nl}x_{l}}=\nu,\;\forall n\in\mathcal{\tilde{N}},$

Both subproblems (11) and (12) with $L=15$ and $N=12$ variables, respectively,
can now be implemented using CVXGEN, as well.

Setting $\rho=0.5$ and initializing the flows $t^{(0)}$ and bandwidths
$w^{(0)}$ to $1$, a max-min rate of $208$Kbps is obtained through the above
approach. The primal and dual residual norms at each iteration are represented
in Figure 8, and Figure 9 (left) represents the max-min rate convergence
results. Simulations using CVXGEN (on an Intel core i7, 2.7 GHz processor)
took an average of $0.327$ms to solve problem (11) and $0.132$ms to solve
problem (12). Problem (10), solved in parallel by all nodes, took $0.145$ms in
average. In this scenario, the semi-distributed algorithm requires
approximately $50$ iterations to converge. Ignoring the time required for
message passing
$50\times\left(0.145+\max\left\\{0.132,0.327\right\\}\right)=23.6$ms are
required on average to obtain the optimal routing and resource allocation
solution. Ignoring the links that have negligible flows (e.g. less than
$10$Kbps in this case), the bottom plot of Figure 5 is obtained, where Link
$4$ has been eliminated due to carrying only $8$Kbps. If desired, the resource
allocation can then be re-optimized based on the updated links.

The number of iterations and time required for the algorithm to converge can
be reduced significantly through “warm start” initialization techniques [11,
ch.4.3]. The rapid convergence of ADMM when the initial point is close to the
optimal solution also allows the semi-distributed implementation to adapt to
changing network conditions. Figure 9 (right) shows the convergence of the
max-min rate when each node’s receiver noise power spectral density $N_{0}$ is
scaled by a random factor between $\frac{1}{2}$ and $\frac{5}{2}$ after the
first 50 iterations. The nodes do not need to notify the CUs of this change,
they just use the new value when solving their subproblems. ADMM took 50
iterations to converge from the initial conditions, but only 12 more
iterations to adapt to the new values for $N_{0}$.

##  V Performance Evaluation

This section evaluates the performance of the proposed hierarchical node
configuration and cross-layer optimization algorithms by using the direct
communication in (2) as the base line for comparison. Simulations consider a
single sector of a cellular network. However, the results can be easily
extended to the entire cellular area, allowing spectrum to be reused among
different sectors as well as among the groups within a sector.

###  V-A Simulation Setup

An instance of the simulation scenario is shown in Figure 10. It consists of a
single sector from a $6$-sector circular cellular network with radius $210$
meters. A total of $44$ users ($N=45$ with destination node) are randomly
distributed within this sector and the common destination node is located at
$(0,0)$. The sector is sub-divided such that the radial distance between
consecutive groups is $d=30$m. Since all nodes within group $1$ communicate
directly with the node at $(0,0)$, the radial distance to the boundary of
group $1$ ($d_{1}$) is set to double the regular distance $d$. Such value
keeps the average link lengths within group $1$ similar to other links in the
network. According to the routing-link assignment policy, with a distance
threshold of $1.5d$ and angular threshold of $15^{0}$, a total of $86$ links
were generated for the network in Figure 10.

The number of hops can be adjusted by varying $d$. For low rate requirements,
as the number of hops increases, the total required transmission power from a
source to destination decreases. However, there is a limit on the number of
hops, since there must be at least one node per group to take advantage of the
frequency reuse. Furthermore, the reduction in transmission power with
increasing number of hops comes at the expense of increasing the number of
nodes cooperating to transmit data and a larger routing table.

The total available bandwidth was assumed to be $W_{\max}=10$MHz with carrier
frequency of $800$Mhz. The receiver power spectral density is
$N_{0}=10^{-11}$W per $1$MHz bandwidth. The reference distance $l_{0}$ in
constraint (1c) is set to $1$ and the pathloss exponent $a=4$. Four different
methods are compared in terms of their average performance: multi-hop
configurations with frequency reuse of $f=3$ and $f=4$, multi-hop
configuration without reuse ($f=\infty$), and direct communication (2). The
results were averaged over $1000$ randomly generated networks, constrained to
have at least one node in each group. A negative factor of power sum,
$-\epsilon\sum{p_{l}}$, with $\epsilon=10^{-6}$ is used in the simulations, as
described in Section III.

###  V-B Results and Discussion

Figure 10: Nodes distributed randomly ($N=45$) with initial set of possible
links ($L=86$). Figure 11: Per user rate vs. per user maximum transmission
power. Figure 12: Total network transmission power vs. per user maximum
transmission power ($N=45$).

The maximum transmission power per device $P_{\max}$ is varied between $0$dBm
and $30$dBm and the minimum rate among all users is plotted in Figure 11.
Furthermore, the total transmission power employed by the network for each
$P_{\max}$ is represented in Figure 12. It is evident from these figures that
the performance of multi-hop configurations with reuse can be significantly
limited because of transmission power limits that avoid interference. To solve
convex problems (1) and (2), CVX, a package for specifying and solving convex
programs [22] is used.

The minimum rate and total power varied significantly among the $1000$
networks that were averaged, but the relative performance of each scheme
remained fairly constant. For example, when the maximum power per user is
limited to $P_{\max}=20$dBm, the standard deviation in the minimum
transmission rate for the multi-hop scheme was $8.6$ Kbps when $f=3$, $45$
Kbps when $f=4$, and $47$ Kbps with no reuse, but the rates with the latter
two were nearly identical and better than with the $f=3$, for all networks.
Similarly, the total transmission power had a standard deviation of $0.5$ dBm
for $f=3$, $2$ dBm for $f=4$, and $1.8$ dBm without frequency reuse, but the
relative positions of the curves remained constant. The standard deviations
with the direct transmission scheme were negligible in all cases.

Problem (1) contains two constraints on transmission power: constraint (1c),
which limits the transmission power on each link and constraint (1d), which
limits per node transmission power $P_{\max}$ (i.e., the total power a node
transmits over all its outgoing links). This implies that when $P_{\max}$ is
large, the transmission powers are limited by constraint (1c) and when
$P_{\max}$ is small, constraint (1d) dominates. Figure 12 shows the $P_{\max}$
threshold values for frequency reuse factors of $f=3,$ ($15$dBm) and $f=4$
($25$dBm), after which increasing $P_{\max}$ does not change the total network
transmission power and hence does not increase the achievable rate. Figure 12
also shows that total network transmission power increases linearly with
$P_{\max}$ in direct mode where constraint (1c) does not apply.

At maximum transmission power of $0$dBm, as Figure 11 depicts, multi-hopping
provides a factor of $10$ increase in data rate relative to direct
communication, while approximately $10$ times less transmission power is
required (Figure 12). However, this performance gain decreases as $P_{\max}$
increases. As theoretical studies of spatial reuse and multi-hopping [25]
suggest, using direct communication between source and destination is
preferred over spatial reuse and multi-hopping when large transmission powers
are allowed. However, for limited power sources, multi-hopping and spatial
reuse perform better than direct mode.

At moderately low transmission powers the configurations with reuse factors of
$3$ and $4$ obtain the same data rate as that without reuse, but the total
network transmission power for each case is different. In general, for the
same data rate, the required transmission power increases with the reuse
factor. When $P_{\max}=10$dBm, the total transmission power required by the
multi-hop scheme without reuse is $2.7$ times the amount of power required
when $f=3$ and $1.5$ times the amount of power required in the $f=4$ case,
despite all three provide the same minimum rate.

The above results, obtained through joint optimization of routing and resource
allocation, often require data to flow from a source node through multiple
routes to the destination node. The number of used routes is determined in
part by the initial pool of links selected during the link assignment step of
Section II-D. One method for reducing the number of used routes, albeit with a
loss in optimality, is to restrict the number of outgoing links from each node
by reducing $d_{th}$ and $\theta$ .

At the expense of further reducing network performance, a less complex
approach is to enforce a single route from each source to the destination.
This can be done by eliminating routing from the joint optimization framework
and simply allowing each node to connect to the closest neighbor in the group
above it. Such approach even though simpler to implement, performs
significantly worse than cross-layer optimization, because the nodes that are
located near the boundary of the groups will be congested with incoming
traffic. Another alternative would be to solve the problem with multiple
routes and have each node use all its resources for the most promising
outgoing link. The performance loss with each of these approaches depends on
many network parameters such as the number of nodes present and their
configuration.

##  VI Summary

The demand for better spectrum use grows dramatically as the wireless
technology and mobile devices progress rapidly. This demand is more stressed
in areas where radio devices are densely packed and need to communicate with a
common destination node. This paper proposes employing these radio devices as
both relays and computing resources for better system performance. In this
method, each device is used to relay other user’s data and perform parallel
optimization of global resource allocation. Simulation results show that the
proposed method can boost each user’s rate by a factor of $10$ while using
lower transmission power. However this gain is at the expense of higher system
complexity.

It is shown that multi-hop communication with reuse has better performance
when transmission power is limited by device power and not interference. By
further increasing spectrum reuse, the available bandwidth increases, but at
the expense of tightening transmission power limits to mitigate co-channel
interference. Pathloss can be further reduced and transmission power lowered
by increasing the number of hops between each node and the destination, but
again at the expense of higher system complexity with larger number of hops.

## Acknowledgment

The authors would like to thank Professor Stephen Boyd for his valuable
discussions. This paper was funded by the Deanship of Scientific Research
(DSR), King Abdulaziz University, under grant No. (11-15-1432 HiCi). The
authors, therefore, acknowledge with thanks DSR technical and financial
support.

## References

  * [1] X. Lin, N. Shroff, and R. Srikant, “A tutorial on cross-layer optimization in wireless networks,” _Selected Areas in Communications, IEEE Journal on_ , vol. 24, no. 8, pp. 1452–1463, 2006.
  * [2] A. Goldsmith, _Wireless Communications_. New York, NY, USA: Cambridge University Press, 2005.
  * [3] C. D. A. S. Michael Neufeld, Jeff Fifield and D. Grunwald, “Softmac - flexible wireless research platform,” _Fourth Workshop on Hot Topics in Networks (HotNets-IV)_ , 2005.
  * [4] T. S. R. Shakkottai, Sanjay and P. C. Karlsson, “Cross-layer design for wireless networks.” _Communications Magazine, IEEE 41.10_ , pp. 74–80, 2003\.
  * [5] V. Srivastava and M. Motani., “Cross-layer design: a survey and the road ahead.” _Communications Magazine, IEEE 43.12_ , pp. 112–119, 2005.
  * [6] L. Xiao, M. Johansson, and S. Boyd, “Simultaneous routing and resource allocation via dual decomposition,” _Communications, IEEE Transactions on_ , vol. 52, no. 7, pp. 1136 – 1144, july 2004.
  * [7] M. Johansson, L. Xiao, and S. Boyd, “Simultaneous routing and power allocation in cdma wireless data networks,” vol. 1, pp. 51–55 vol.1, 2003.
  * [8] B. Timus, P. Soldati, D. Kim, and J. Zander, “Cross-layer resource allocation model for cellular-relaying network performance evaluation,” _Vehicular Technology, IEEE Transactions on_ , vol. 60, no. 6, pp. 2765–2776, 2011.
  * [9] T. Ng and W. Yu, “Joint optimization of relay strategies and resource allocations in cooperative cellular networks,” _Selected Areas in Communications, IEEE Journal on_ , vol. 25, no. 2, pp. 328–339, 2007.
  * [10] S. Boyd. (2011) Branch and bound methods. [Online]. Available: http://stanford.edu/class/ee364b/lectures/bb_slides.pdf
  * [11] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, “Distributed optimization and statistical learning via the alternating direction method of multipliers,” _Foundations and Trends in Machine Learning_ , vol. 3, no. 1, pp. 1–122, 2011.
  * [12] D. Palomar and M. Chiang, “Alternative distributed algorithms for network utility maximization: Framework and applications,” _Automatic Control, IEEE Transactions on_ , vol. 52, no. 12, pp. 2254–2269, 2007.
  * [13] J. Liu and H. D. Sherali, “A distributed newton’s method for joint multi-hop routing and flow control: Theory and algorithm.” _INFOCOM, 2012 Proceedings IEEE_ , pp. 2489 – 2497, 2012.
  * [14] C. Shen, T.-H. Chang, K.-Y. Wang, Z. Qiu, and C.-Y. Chi, “Distributed robust multicell coordinated beamforming with imperfect csi: An admm approach,” _Signal Processing, IEEE Transactions on_ , vol. 60, no. 6, pp. 2988–3003, 2012.
  * [15] J. Mota, J. Xavier, P. Aguiar, and M. Puschel, “Distributed admm for model predictive control and congestion control,” in _Decision and Control (CDC), 2012 IEEE 51st Annual Conference on_ , 2012, pp. 5110–5115.
  * [16] D. P. Bertsekas, _Network Optimization: Continuous and Discrete Models_. Belmont, MA, USA: Athena Scientific, 1998.
  * [17] M. R. Sanatkar and A. Mohammadi, “Scalability analysis of wireless sensor networks using analytical techniques,” _8th Annual Communication Networks and Services Research Conference_ , vol. 2, no. 1, pp. 298–303, 2010\.
  * [18] S. Boyd and L. Vandenberghe, _Convex Optimization_. New York, NY, USA: Cambridge University Press, 2004.
  * [19] A. Nedic and A. Ozdaglar, “Distributed subgradient methods for multi-agent optimization,” _Automatic Control, IEEE Transactions on_ , vol. 54, no. 1, pp. 48–61, 2009.
  * [20] S. Sesia and et al, _LTE, The UMTS Long Term Evolution: From Theory to Practice_. Wiley and sons: 2nd ed., 2011\.
  * [21] E. Wei and A. Ozdaglar, “On the $o(1/k)$ convergence of asynchronous distributed alternating direction method of multipliers.” _arXiv preprint arXiv:1307.8254_ , 2013.
  * [22] I. CVX Research, “CVX: Matlab software for disciplined convex programming, v2.0 beta,” http://cvxr.com/cvx, Sep. 2012.
  * [23] J. Mattingley and S. Boyd, “Cvxgen: A code generator for embedded convex optimization,” _Optimization and Engineering_ , vol. 13, no. 1, pp. 1–27, 2010.
  * [24] A. Imamoto and B. Tang, “A recursive descent algorithm for finding the optimal minimax piecewise linear approximation of convex functions,” in _World Congress on Engineering and Computer Science._ IEEE, 2008, pp. 287–293.
  * [25] V. Ramaiyan and A. Kumar, “On the limits of spatial reuse and cooperative communication for dense wireless networks,” _Information Theory for Wireless Networks, 2007 IEEE Information Theory Workshop on_ , pp. 1–5, July.


∎

11institutetext: F.Sakuma, H.Asano, K.Itahashi, M.Iwasaki, K.Kanno, Y.Kato,
Y.Ma, R.Murayama, H.Outa, T.Yamaga, T.Yamazaki, Q.Zhang 22institutetext: RIKEN
Cluster for Pioneering Research, RIKEN, Wako 351-0198, Japan 22email:
<EMAIL_ADDRESS>33institutetext: S.Ajimura, T.Hiraiwa, K.Inoue,
S.Kawasaki, H.Noumi, K.Shirotori, D.Tomono 44institutetext: Research Center
for Nuclear Physics (RCNP), Osaka University, Osaka 567-0047, Japan
55institutetext: T.Akaishi, A.Sakaguchi 66institutetext: Department of
Physics, Osaka University, Osaka 560-0043, Japan 77institutetext: G.Beer
88institutetext: Department of Physics and Astronomy, University of Victoria,
Victoria, British Columbia V8W 3P6, Canada 99institutetext: H.Bhang, S.Choi
1010institutetext: Department of Physics, Seoul National University, Seoul
151-742, South Korea 1111institutetext: M.Bragadireanu, D.Sirghi, F.Sirghi
1212institutetext: National Institute of Physics and Nuclear Engineering -
IFIN HH, Bucharest, Romania 1313institutetext: P.Buehler, M.Cargnelli,
J.Marton, H.Shi, M.Simon, E.Widmann, J.Zmeskal 1414institutetext: Stefan-
Meyer-Institut für subatomare Physik, A-1020 Vienna, Austria
1515institutetext: L.Busso, O.Morra 1616institutetext: Istituto Nazionale di
Fisica Nucleare (INFN), Sezione di Torino, Torino, Italy 1717institutetext:
L.Busso 1818institutetext: Dipartimento di Fisica Generale, Universita’ di
Torino, Torino, Italy 1919institutetext: M.Bazzi, A.Clozza, C.Curceanu,
C.Guaraldo, M.Iliescu, M.Miliucci, K.Piscicchia, A.Scordo, D.Sirghi, F.Sirghi,
O.Vázquez Doce 2020institutetext: Laboratori Nazionali di Frascati dell’ INFN,
I-00044 Frascati, Italy 2121institutetext: S.Enomoto, M.Iio, S.Ishimoto,
M.Iwai, H.Noumi, K.Ozawa, M.Sato, M.Sekimoto, S.Suzuki, A.Toyoda
2222institutetext: High Energy Accelerator Research Organization (KEK),
Tsukuba 305-0801, Japan 2323institutetext: H.Fujioka, H.Kou, M.Tokuda
2424institutetext: Department of Physics, Tokyo Institute of Technology, Tokyo
152-8551, Japan 2525institutetext: Y.Fujiwara, R.S.Hayano, T.Ishikawa,
T.Suzuki, T.Yamazaki 2626institutetext: Department of Physics, The University
of Tokyo, Tokyo 113-0033, Japan 2727institutetext: T.Fukuda, Y.Mizoi
2828institutetext: Laboratory of Physics, Osaka Electro-Communication
University, Osaka 572-8530, Japan 2929institutetext: T.Hashimoto, K.Tanida
3030institutetext: ASRC, Japan Atomic Energy Agency, Ibaraki 319-1195, Japan
3131institutetext: Y.Ishiguro, K.Kato, T.Nagae 3232institutetext: Department
of Physics, Kyoto University, Kyoto 606-8502, Japan 3333institutetext: P.
Kienle†, O.Vázquez Doce 3434institutetext: Technische Universität München,
D-85748, Garching, Germany 3535institutetext: Y.Matsuda 3636institutetext:
Graduate School of Arts and Sciences, The University of Tokyo, Tokyo 153-8902,
Japan 3737institutetext: H.Ohnishi, Y.Sada, C.Yoshida 3838institutetext:
Research Center for Electron Photon Science (ELPH), Tohoku University, Sendai
982-0826, Japan 3939institutetext: S.Okada 4040institutetext: Engineering
Science Laboratory, Chubu University, Aichi 487-8501, Japan 4141institutetext:
K.Piscicchia 4242institutetext: Centro Ricerche Fermi-Museo Storico della
Fisica e Centro studi e ricerche “Enrico Fermi”, 000184 Rome, Italy
4343institutetext: H.Tatsuno 4444institutetext: Department of Physics, Tokyo
Metropolitan University, Tokyo 192-0397, Japan 4545institutetext: K.Tsukada
4646institutetext: Institute for Chemical Research, Kyoto University, Kyoto
611-0011, Japan

#  Recent results and future prospects of kaonic nuclei at J-PARC

F.Sakuma  S.Ajimura  T.Akaishi  H.Asano  M.Bazzi  G.Beer  H.Bhang
M.Bragadireanu  P.Buehler  L.Busso  M.Cargnelli  S.Choi  A.Clozza  C.Curceanu
S.Enomoto  H.Fujioka  Y.Fujiwara  T.Fukuda  C.Guaraldo  T.Hashimoto
R.S.Hayano  T.Hiraiwa  M.Iio  M.Iliescu  K.Inoue  Y.Ishiguro  T.Ishikawa
S.Ishimoto  K.Itahashi  M.Iwasaki  M.Iwai  K.Kanno  K.Kato  Y.Kato  S.Kawasaki
P.Kienle  H.Kou  Y.Ma  J.Marton  Y.Matsuda  M.Miliucci  Y.Mizoi  O.Morra
R.Murayama  T.Nagae  H.Noumi  H.Ohnishi  S.Okada  H.Outa  K.Ozawa
K.Piscicchia  Y.Sada  A.Sakaguchi  M.Sato  A.Scordo  M.Sekimoto  H.Shi
K.Shirotori  M.Simon  D.Sirghi  F.Sirghi  S.Suzuki  T.Suzuki  K.Tanida
H.Tatsuno  M.Tokuda  D.Tomono  A.Toyoda  K.Tsukada  O.Vázquez Doce  E.Widmann
T.Yamaga  T.Yamazaki  C.Yoshida  Q.Zhang  J.Zmeskal

(Received: date / Accepted: date)

###### Abstract

$\bar{K}$-nuclear bound systems, kaonic nuclei, have been widely discussed as
products of the strongly attractive $\bar{K}N$ interaction in $I=0$ channels.
Recently, we demonstrated that kaonic nuclei can be produced via in-flight
$(K^{-},N)$ reactions using the low-momentum DC kaon beam at the J-PARC E15
experiment. We observed the simplest kaonic nuclei, $K^{-}pp$, having a much
deeper binding energy than normal nuclei. For further studies, we have
proposed a series of experimental programs for the systematic investigation of
light kaonic nuclei, from $\bar{K}N$ ($\Lambda(1405)$) to $\bar{K}NNNN$. In
the new experiment approved as J-PARC E80, we will measure the $\bar{K}NNN$
($A=3$) system as a first step toward a comprehensive study.

##  1 Introduction

The study of the $\bar{K}N$ interaction is one of the most important subjects
in meson–baryon interactions in low energy quantum chromodynamics (QCD).
Extensive measurements of anti-kaonic hydrogen atoms Iwasaki:1997wf ;
Beer:2005qi ; Bazzi:2011zj  and low-energy $\bar{K}N$ scattering Martin:1980qe
have revealed the strongly attractive nature of the $\bar{K}N$ interaction in
the isospin I = 0 channel. Consequently, the possible existence of deeply
bound kaonic nuclear states (kaonic nuclei) has been widely discussed
Nogami:1963xqa ; Akaishi:2002bg ; Yamazaki:2002uh ; Shevchenko:2006xy ;
Shevchenko:2007ke ; Ikeda:2007nz ; Dote:2008in ; Ikeda:2008ub ; Wycech:2008wf
; Dote:2008hw ; Ikeda:2010tk ; Barnea:2012qa ; Bayar:2012hn ; Maeda:2013 ;
Revai:2014twa ; Dote:2014via ; Sekihara:2016vyd ; Ohnishi:2017uni ;
Dote:2017veg ; Dote:2017wkk . Kaonic nuclei are predicted to be compact due to
the strong $\bar{K}N$ attraction, suggesting that high-density nuclear matter
is realized in kaonic systems.

Among kaonic nuclei, the $\bar{K}NN$ system with $I=1/2$ and $J^{P}=0^{-}$
(symbolically denoted as $K^{-}pp$ for the $I_{z}=+1/2$ state) is of special
interest because it is the lightest $S=-1$ $\bar{K}$ nucleus. Despite
considerable experimental efforts over the past 20 years, it has been
challenging to prove the existence of $K^{-}pp$. Several groups have reported
observations of a $K^{-}pp$ candidate with a binding energy of around 100 MeV
in experiments measuring non-mesonic decay branches of $\Lambda p$ and/or
$\Sigma^{0}p$ in different reactions Agnello:2005qj ; Yamazaki:2010mu ;
Ichikawa:2014ydh . There are also contradicting reports concluding that the
reactions can be understood without a bound state Doce:2015ust ;
DelGrande:2018sbv ; Agakishiev:2014dha ; Tokiyasu:2013mwa .

Recently, we confirmed the existence of the $K^{-}pp$ bound state by using the
simplest reaction of in-flight 3He$(K^{-},N)$ at the J-PARC E15 experiment
Hashimoto:2014cri ; Sada:2016nkb ; Ajimura:2018iyx ; J-PARCE15:2020gbh . A
distinct peak structure was observed well below the mass threshold of
$K^{-}+p+p$ in the $\Lambda p$ invariant-mass (IM) spectrum, obtained from the
3He$(K^{-},\Lambda p)n$ measurement. The simplest and most natural
interpretation of this peak is $K^{-}pp$. This result is experimentally solid
compared to previously reported results.

To obtain further detailed information on kaonic nuclei, we have planned a
series of experimental programs using the $(K^{-},N/d)$ reaction on light
nuclear targets. The programs will enable a detailed study of a range of
nuclei from $\bar{K}N$ ($\Lambda(1405)$) to $\bar{K}NNNN$ using the world’s
highest intensity low-momentum kaon beam at J-PARC. The programs comprise:

  * •

Precise measurements of the $\Lambda(1405)$ state in a large momentum transfer
region via the $d(K^{-},n)$ reaction, to experimentally clarify whether it is
a baryonic state or a $\bar{K}N$ molecular state,

  * •

Investigations of the spin and parity of the $\bar{K}NN$ state via
3He$(K^{-},N)$ reactions,

  * •

A search for $\bar{K}NNN$ states via 4He$(K^{-},N)$ reactions, as a bridge to
access heavier systems, and

  * •

An advanced search for $\bar{K}NNNN$ states via the 6Li$(K^{-},d)$ reaction.

In parallel to these studies, we also intend to access the $S=-2$ kaonic
nuclei, such as the theoretically predicted $K^{-}K^{-}pp$ state. The $S=-2$
system could allow us to access even higher density systems than the $S=-1$
kaonic nuclei. As described in our Letter of Intent KKpp_LoI , one possible
approach for the measurements at J-PARC could be:

  * •

Searching for $\bar{K}\bar{K}NN$ states via $\bar{p}{{}^{3}}$He annihilation.

To ensure the measurements are systematic and precise, we are planning to
construct a totally new $4\pi$ spectrometer to measure all the particles
involved in the reactions and to reconstruct their formation and decay
exclusively. The spectrometer is designed to be highly versatile so that all
the experiments can be performed simply by changing the target materials. In
addition, for more efficient use of the high-intensity kaon beam, we have
proposed shortening the existing K1.8BR beam line for a larger kaon yield
without deteriorating the momentum resolution of the kaon beam.

In this article, we review previous $K^{-}pp$ measurements at the J-PARC E15
experiment in Sec. 2 and set out what we have learned from the experimental
investigations dedicated to kaonic nuclei. In Sec.3, we present the details of
the new experiment, approved as J-PARC E80, that focuses on searching for
$\bar{K}NNN$.

##  2 Results of the J-PARC E15 Experiment

We conducted an experimental investigation of the $K^{-}pp$ bound state using
the simplest $\bar{K}$ induced reaction of $K^{-}+{{}^{3}}$He, via the nucleon
knock-out reaction $K^{-}N\to\bar{K}n$ followed by the two-nucleon absorption
$\bar{K}+NN\to$ $K^{-}pp$. In the experiment, we used a kaon momentum of 1
GeV/$c$, around which the $K^{-}N\to\bar{K}n$ reactions have the maximum cross
section. The recoiled kaon $\bar{K}$ at a momentum $q$ behaves as an ‘off-
shell particle’ (the total energy can be lower than its intrinsic mass) within
a time range allowed by the uncertainty principle. The momentum transfer $q$
is defined as that between the incident kaon and the outgoing neutron in the
laboratory frame $q=|p^{lab}_{K^{-}}-p^{lab}_{n}|$. In the reaction, we used
the low-momentum back-scattered kaon as an off-shell kaon source and the
residual spectator nucleons $NN$ as an ‘actual target’ to form a $K^{-}pp$
state with an energy below the intrinsic mass of $M(Kpp)$ (= $m_{K}+2m_{N}$ =
2.37 GeV/$c^{2}$).

Figure 1:  (left) Efficiency and acceptance corrected data over $IM$ ($\Lambda
p$ invariant mass) and $q$ (momentum transfer). (right) $\Lambda p$ invariant
mass in the region $0.3<q<0.6$ MeV/$c$.

For the $\Lambda pn$ final states, we observed a kinematic anomaly in the
$\Lambda p$ invariant mass near the mass threshold of $M(Kpp)$ at $q\sim$ 0.4
GeV/$c$ Sada:2016nkb ; Ajimura:2018iyx ; J-PARCE15:2020gbh . As shown in Fig.
1 (left), we confirmed the existence of the bound state below $M(Kpp)$, whose
mass centroid is independent of $q$, at a binding energy as deep as $\sim$ 50
MeV. The back-scattered ‘on-shell’ kaon, whose total kaon energy is above its
intrinsic mass, can also be absorbed by the spectator nucleons without forming
a bound state. The kinematical centroid of this quasi-free absorption process
is plotted in the figure, denoted as $QF_{K}$. Along the line $QF_{K}$, there
are two event concentration points at $\theta_{n}=0$ and $\theta_{n}=\pi$, but
both are well separated from the region of interest. Figure 1 (right) shows
the $\Lambda p$ invariant mass spectrum after the acceptance and efficiency
corrections for the region $0.3<q<0.6$ GeV/$c$ where the $K^{-}pp$ bound state
is dominant. A clear peak originating from $K^{-}pp$ can be seen below
$M(Kpp)$, shown by the red line, whose binding energy reaches $\sim$ 50 MeV.

The simplest and most natural interpretation is the kaon-nuclear bound state
$K^{-}pp$. The observed large form factor of $\sim$ 400 MeV/$c$, based on a
simple plane wave impulse approximation (PWIA) Sada:2016nkb ; Ajimura:2018iyx
; J-PARCE15:2020gbh , and the large binding energy of the $K^{-}pp$ state
imply the formation of a fairly compact and dense system. The observed
structure below the mass threshold in the $\Lambda p$ spectrum has also been
interpreted as the $\bar{K}NN$ quasi-bound system: a theoretical calculation
in Ref. Sekihara:2016vyd ; Sekihara:QNP2018  demonstrated that the
experimental spectrum can be reproduced with the $\bar{K}NN$ quasi-bound
system and the quasi-free processes based on a theoretical treatment of the
3He$(K^{-},\Lambda p)n$ reaction.

Thus, the E15 experiment opened a new era of experimental research on kaonic
nuclei with the virtual kaon beam produced from in-flight $(K^{-},N)$
reactions. This has been realized by using the world’s highest intensity kaon
beam available at J-PARC. However, we encountered several issues with the
existing setup that made further investigations of the kaonic nuclei
difficult. The new series of experiments will continue the systematic research
into light kaonic nuclei using a new large acceptance detector system and the
improved K1.8BR beam line.

##  3 New Experiment (J-PARC E80 Experiment)

For the next stage, we aim to determine the mass number dependence of the
binding energy, decay width, and system size beyond $\bar{K}NN$. The mass
number dependence has been calculated with several theoretical models, as
summarized in Fig. 2. The values of the binding energy and decay width
predicted by the models vary widely due to the differences in the $\bar{K}N$
interaction models. However, almost all the predictions show that the larger
nuclei have stronger binding energies. For the width, the theoretical
calculations take into account only mesonic decay channels, such as $\pi\Sigma
N$ and $\pi\Lambda N$. The calculated width is expected to be larger if the
models adopt non-mesonic decay channels, as demonstrated in Ref. Bayar:2012hn
.

From an experimental point of view, when the mass number is large, it becomes
difficult to handle the number of particles in the final state and to deduce
the physics behind the reaction. Thus, we take a step-by-step approach. In the
new experiment, J-PARC E80, we aim to measure the $\bar{K}NNN$ ($A=3$) system
as a first step toward a comprehensive study. From the experience of E15, we
have learned that reducing the number of particles in the final state is a key
to removing ambiguity in interpreting the reaction process. Therefore, we
focus on the $K^{-}ppn\to\Lambda d$ and $\Lambda pn$ decay channels in
4He$(K^{-},\Lambda d/\Lambda pn)n$ reactions.

The new experiment will provide the mass number dependence of the kaonic
nuclei for the first time. The dependence can more clearly reveal the
$\bar{K}N$ interaction below the mass threshold, by comparing the obtained
properties of the $\bar{K}NNN$ state with those of the already reported
$K^{-}p$ ($\Lambda(1405)$) and $K^{-}pp$ states.

Figure 2: Summary of theoretical calculations of the kaonic nuclei from $A=2$
to $4$ in different models AY Akaishi:2002bg ; Yamazaki:2002uh , WG
Wycech:2008wf , BGL Barnea:2012qa , and OHHMH Ohnishi:2017uni . The result
obtained at the E15 experiment is also plotted Ajimura:2018iyx .

###  3.1 Experimental Method

To date, stopped $K^{-}$ reactions have mainly been used to search for the
$\bar{K}NNN$ state. The KEK-PS E471/E549 collaborations measured the inclusive
4He$(K^{-}_{stopped},p/n)$ reactions with a spectrometer dedicated to TOF
measurement of protons and neutrons Sato:2007sb ; Yim:2010zza . They found no
specific peak structures below the mass threshold of $M(\bar{K}NNN)$ in the
missing mass spectra from the reactions. This is due to the huge background
originating from the two-nucleon absorption processes, $\bar{K}NN\to YN$, and
quasi-free hyperon productions and its decays, $\bar{K}N\to\pi Y$, whose
production mechanisms are quite complicated. The background cannot be
discriminated kinematically from the signal with small production cross
section in the inclusive measurements.

Indeed, our inclusive analysis of the in-flight 3He$(K^{-},n)X$ measurement
showed no significant peak structures due to the huge background from quasi-
free processes and the two-nucleon absorption processes Hashimoto:2014cri . By
reducing the background by the exclusive measurement of 3He$(K^{-},\Lambda
p)n$ and identifying all the final-state particles in a wide momentum-transfer
region, the $K^{-}pp$ signal could be separated out. We found that the
$K^{-}pp$ signal has a much smaller cross section ($\sim$ 10 $\mu$b) than the
quasi-free processes ($\sim$ 10 mb) and is distributed up to $q\sim$ 600
MeV/$c$ Sada:2016nkb ; Ajimura:2018iyx ; J-PARCE15:2020gbh .

On the other hand, there are two reports of observations of the $K^{-}ppn$
candidate below the mass threshold in the $\Lambda d$ invariant mass spectrum
with stopped $K^{-}$ reactions (FINUDA) and heavy-ion collisions (FOPI). The
FINUDA collaboration reported that the candidate has a binding energy of
$\sim$ 60 MeV and a width of $\sim$ 40 MeV Agnello:2007aa . The candidate
reported by the FOPI collaboration has a much deeper binding energy of $\sim$
150 MeV with a broader width of $\sim$ 100 MeV Herrmann:2005 . However,
because these measurements were performed only inclusively, possible
contributions from multi-nucleon absorption processes and $N^{*}/Y^{*}$ decays
to the peak structure cannot be excluded. Furthermore, the statistics reported
were very limited and thus the results remain speculative.

Therefore, the key to the experimental search is to adopt a simple reaction
and to measure it exclusively. Adopting a simple reaction, such as in-flight
$\bar{K}$ induced reactions with light target nuclei, enables us to specify
the reaction channel using the momentum-transfer dependence. Exclusive
measurements are crucial for distinguishing small and broad signals from the
large and widely distributed quasi-free and multi-nucleon absorption
background.

In the new experiment, we will perform exclusive measurements of the
production and decay of the $K^{-}ppn$ state using the in-flight reaction

$\displaystyle K^{-}+^{4}{\rm He}$ $\displaystyle\to$ $\displaystyle
K^{-}ppn+n$

followed by the expected no-mesonic decays

$\displaystyle K^{-}ppn$ $\displaystyle\to$ $\displaystyle\Lambda+d{\rm,}$
$\displaystyle K^{-}ppn$ $\displaystyle\to$ $\displaystyle\Lambda+p+n.$

We aim to determine the binding energy and width from the invariant mass
reconstruction of the decays. The invariant mass will be obtained as a
function of the momentum transfer to distinguish the bound-state production
from the quasi-free processes and multi-nucleon absorption processes by the
event kinematics as demonstrated in the E15 analysis. In the
$K^{-}+{{}^{4}}{\rm He}$ reaction, it is possible to measure the isospin
partner of $K^{-}ppn$, i.e., the $K^{-}pnn$ state via

$\displaystyle K^{-}+^{4}{\rm He}$ $\displaystyle\to$ $\displaystyle
K^{-}pnn+p,$ $\displaystyle K^{-}pnn$ $\displaystyle\to$
$\displaystyle\Lambda+n+n.$

Comparing the properties of the isospin partners is of special importance for
investigating the internal composition of the kaonic nuclei. On the other
hand, the measurement is challenging, because two-neutron detection is
required to identify the $K^{-}pnn$ decay.

###  3.2 Apparatus

We aim to produce the $\bar{K}NNN$ state using the $(K^{-},N)$ reaction at the
K1.8BR beam line, as successfully performed in the previous experiment E15; a
recoiled virtual kaon ($\bar{K}$) generated by $K^{-}N\to\bar{K}N$ processes
can be directly induced into residual nucleons within the strong interaction
range. We will utilize 1.0 GeV/$c$ incident kaons to maximize the $\bar{K}N$
reaction rate at around zero degrees. The incoming $K^{-}$ beam will be
identified and its momentum analyzed by the beam-line spectrometer. The beam
kaon will irradiate a liquid 4He target located at the final focus point, and
all the particles generated from the reactions will be identified with a
cylindrical detector system (CDS) surrounding the target system. A conceptual
design of the CDS is shown in Fig. 3. It is mainly composed of a large
superconducting solenoid magnet, a cylindrical wire drift chamber, and a
cylindrical neutron detector. The kaonic nuclei will then be identified via
invariant-mass reconstruction of the decay particles. By detecting the nucleon
coming from the initial $(K^{-},N)$ reaction, or by identifying it with the
missing mass technique, we will realize exclusive measurement of the kaonic
nuclei. The details of the apparatus used for the experiment can be found in
Ref. E80_proposal .

Figure 3:  Conceptual design of the CDS.

###  3.3 Expected Spectrum

Figure 4:  Expected spectra of (top) $\Lambda d$ and (bottom) $\Lambda pn$ in
the $K^{-}{{}^{4}}$He $\to\Lambda d+n$ and $\Lambda pn+n$ final state with
100-G kaons on target, respectively. We assume the cross section of the
$K^{-}ppn$, quasi-free, and background to be 10 $\mu$b, 10 $\mu$b, and 20
$\mu$b, respectively, based on the E15 results.

Figure 4 shows the expected spectra of the $K^{-}$ 4He $\to\Lambda dn$ and
$\Lambda pnn$ final states with 100-G kaons on target, corresponding to three
weeks data taking under 90-kW beam power of the J-PARC Main Ring accelerator
E80_proposal . For the $K^{-}ppn$ state in the figures, we assume a similar
distribution to that of $K^{-}pp$ observed in E15, with a binding energy of
$\sim$ 50 MeV, a width of $\sim$ 100 MeV, and a Gaussian form factor of $\sim$
400 MeV/$c$ with 10 $\mu$b ($\sigma\cdot BR$). For the quasi-free process and
broad background, we also used similar parameters to those obtained for the
$K^{-}{{}^{3}}$He $\to\Lambda pn$ final state in E15. The spectra for each
process are generated using the Geant4 simulation, in which we assume the same
resolution of the existing CDS to reconstruct the simulated tracks:
$5.3\%\times p_{t}\oplus 0.5\%$ for charged particles Sada:2016nkb  and $\sim
10\%\times p\oplus 7\%$ for a neutron 111based on the preliminary result of
$K^{-}$${}^{3}{\rm He}\to\pi\Sigma pn$ analysis Sakuma:2020dnc .

In the assumed conditions, the bound states are clearly identified not only in
the invariant mass spectra but also on the 2-dimensional plane of the momentum
transfer versus the invariant mass ($q$–$IM$) distribution, which is the most
important feature of the experiment. By investigating the specific structures
in the $q$–$IM$ distribution, we can kinematically identify the signal of the
bound state. Figure 5 shows a demonstration of the signal enhancement by
selecting a momentum transfer of $0.3<q<0.6$ GeV/$c$. In this region, we
expect the bound state to be clearly separated from the quasi-free process, as
in the case of the E15 experiment shown in Fig. 1.

Figure 5:  Expected invariant-mass spectra of (left) $\Lambda d$ and (light)
$\Lambda pn$ in the region $0.3<q<0.6$ GeV/$c$, obtained from Fig. 4.

##  4 Summary

We demonstrated that kaonic nuclei can be produced via in-flight $(K^{-},N)$
reactions using the low-momentum DC kaon beam at the J-PARC E15 experiment. We
observed that the simplest kaonic nuclei, $K^{-}pp$, has a much deeper binding
energy than normal nuclei. We also found that the large form factor obtained
in a PWIA analysis implies the possible formation of a compact and dense
system. For the next stage, we have proposed a series of experimental programs
for the systematic investigation of light kaonic nuclei, from $\bar{K}N$
($\Lambda(1405)$) to $\bar{K}NNNN$. Through the experiments, we will determine
the features of kaonic nuclei depending on the mass number $A$, i.e., nuclear
density, which is related to spontaneous and explicit chiral symmetry breaking
in QCD. In the new experiment approved as J-PARC E80 we will measure the
$\bar{K}NNN$ ($A=3$) system as a first step toward a comprehensive study.

###### Acknowledgements.

We are grateful to all the staff members of J-PARC/KEK/JAEA for their
extensive efforts in the successful operation of the facility. This work is
supported in part by MEXT Grants-in-Aid 26800158, 17K05481, 26287057,
24105003, 14102005, 17070007, 18H05402, 18H01237, and 20K04006. Part of this
work is also supported by the Ministero degli Affari Esteri e della
Cooperazione Internazionale, Direzione Generale per la Promozione del Sistema
Paese (MAECI), StrangeMatter project.

## References

  * (1) M. Iwasaki, et al., Phys. Rev. Lett. 78, 3067 (1997)
  * (2) G. Beer, et al., Phys. Rev. Lett. 94, 212302 (2005)
  * (3) M. Bazzi, et al., Phys. Lett. B704, 113 (2011)
  * (4) A.D. Martin, Nucl. Phys. B179, 33 (1981)
  * (5) Y. Nogami, Phys. Lett. 7, 288 (1963)
  * (6) Y. Akaishi, T. Yamazaki, Phys. Rev. C65, 044005 (2002)
  * (7) T. Yamazaki, Y. Akaishi, Phys. Lett. B535, 70 (2002)
  * (8) N.V. Shevchenko, A. Gal, J. Mares, Phys. Rev. Lett. 98, 082301 (2007)
  * (9) N.V. Shevchenko, A. Gal, J. Mares, J. Revai, Phys. Rev. C76, 044004 (2007)
  * (10) Y. Ikeda, T. Sato, Phys. Rev. C76, 035203 (2007)
  * (11) A. Doté, T. Hyodo, W. Weise, Nucl. Phys. A804, 197 (2008)
  * (12) Y. Ikeda, T. Sato, Phys. Rev. C79, 035201 (2009)
  * (13) S. Wycech, A.M. Green, Phys. Rev. C79, 014001 (2009)
  * (14) A. Doté, T. Hyodo, W. Weise, Phys. Rev. C79, 014003 (2009)
  * (15) Y. Ikeda, H. Kamano, T. Sato, Prog. Theor. Phys. 124, 533 (2010)
  * (16) N. Barnea, A. Gal, E.Z. Liverts, Phys. Lett. B712, 132 (2012)
  * (17) M. Bayar, E. Oset, Phys. Rev. C88, 044003 (2013)
  * (18) S. Maeda, Y. Akaishi, T. Yamazaki, Proc. Jpn. Acad. B89, 418 (2013)
  * (19) J. Révai, N.V. Shevchenko, Phys. Rev. C90, 034004 (2014)
  * (20) A. Doté, T. Inoue, T. Myo, Prog. Theor. Exp. Phys. 2015, 043D02 (2015)
  * (21) T. Sekihara, E. Oset, A. Ramos, Prog. Theor. Exp. Phys. 2016, 123D03 (2016)
  * (22) S. Ohnishi, et al., Phys. Rev. C95, 065202 (2017)
  * (23) A. Doté, T. Inoue, T. Myo, Phys. Rev. C95, 062201 (2017)
  * (24) A. Doté, T. Inoue, T. Myo, Phys. Lett. B784, 405 (2018)
  * (25) M. Agnello, et al., Phys. Rev. Lett. 94, 212303 (2005)
  * (26) T. Yamazaki, et al., Phys. Rev. Lett. 104, 132502 (2010)
  * (27) Y. Ichikawa, et al., Prog. Theor. Exp. Phys. 2015, 021D01 (2015)
  * (28) O. Vzquez Doce, et al., Phys. Lett. B758, 134 (2016)
  * (29) R. Del Grande, et al., Eur. Phys. J. C79, 190 (2019)
  * (30) G. Agakishiev, et al., Phys. Lett. B742, 242 (2015)
  * (31) A.O. Tokiyasu, et al., Phys. Lett. B728, 616 (2014)
  * (32) T. Hashimoto, et al., Prog. Theor. Exp. Phys. 2015, 061D01 (2015)
  * (33) Y. Sada, et al., Prog. Theor. Exp. Phys. 2016, 051D01 (2016)
  * (34) S. Ajimura, et al., Phys. Lett. B789, 620 (2019)
  * (35) T. Yamaga, et al., Phys. Rev. C 102, 044002 (2020)
  * (36) Letter of Intent for J-PARC, 2009, ”Double Anti-kaon Production in Nuclei by Stopped Anti-proton Annihilation”, http://j-parc.jp/researcher/Hadron/en/Proposal_e.html
  * (37) T. Sekihara, E. Oset, A. Ramos, JPS Conf. Proc. 26, 023009 (2019)
  * (38) M. Sato, et al., Phys. Lett. B 659, 107 (2008)
  * (39) H. Yim, et al., Phys. Lett. B 688, 43 (2010)
  * (40) M. Agnello, et al., Phys. Lett. B 654, 80 (2007)
  * (41) N. Herrmann, et al., Proc. EXA05 Conference ISBN 3-7001-3616-1, 73 (2005)
  * (42) Proposal for J-PARC, 2020, ”Systematic investigation of the light kaonic nuclei”, http://j-parc.jp/researcher/Hadron/en/Proposal_e.html
  * (43) F. Sakuma, et al., AIP Conf. Proc. 2249, 020005 (2020)