# Capacity Analysis and Sum Rate Maximization for the SCMA Cellular Network
Coexisting with D2D Communications

Yukai Liu  Wen<EMAIL_ADDRESS>Shanghai Institute of Advanced
Communications and Data Sciences, Department of Electronic Engineering,
Shanghai Jiao Tong University, Shanghai 200240, China

###### Abstract

Sparse code multiple access (SCMA) is the most concerning scheme among non-
orthogonal multiple access (NOMA) technologies for 5G wireless communication
new interface. Another efficient technique in 5G aimed to improve spectral
efficiency for local communications is device-to-device (D2D) communications.
Therefore, we utilize the SCMA cellular network coexisting with D2D
communications for the connection demand of the Internet of things (IOT), and
improve the system sum rate performance of the hybrid network. We first derive
the information-theoretic expression of the capacity for all users and find
the capacity bound of cellular users based on the mutual interference between
cellular users and D2D users. Then we consider the power optimization problem
for the cellular users and D2D users jointly to maximize the system sum rate.
To tackle the non-convex optimization problem, we propose a geometric
programming (GP) based iterative power allocation algorithm. Simulation
results demonstrate that the proposed algorithm converges fast and well
improves the sum rate performance.

###### keywords:

SCMA; D2D; cellular capacity; geometric programming

##  1 Introduction

The large connection and high data rate constitute the core goals for the 5G
wireless networks. Some key technologies of 5G systems include device-to-
device (D2D) communications, non-orthogonal multiple access techniques (NOMA)
along with massive multiple-input multiple-output (MIMO), ultra-dense radio
networking, all-spectrum access, and so on [1]. The conventional orthogonal
multiple access (OMA) schemes allocate orthogonal resource blocks (RBs) either
in time, frequency, or code domains to different users [2, 3], which is
underloaded as such schemes occupy more RBs than the users. However, NOMA can
support massive user access via non-orthogonal RB allocation [4], which is an
overloaded system. Two main approaches of the power-domain NOMA and the code-
domain NOMA are widely investigated. When multiple users transmit at the same
RB by power-domain NOMA, users are allocated with different power levels,
while the code-domain NOMA transmits at the same RB in distinguishing
spreading codes [5, 6].

Among the currently proposed code domain NOMA, sparse code multiple access
(SCMA) has been commonly considered as a promising method. SCMA is proposed by
Nikopour and Baligh [7], in which, coded bits are directly mapped into the
multi-dimensional complex lattice point (called codeword) and the codewords
are designed to be sparse. The sparsity of the SCMA codewords enables massive
connectivity and the use of suboptimal message passing algorithm (MPA) to
detect multiple users. For general SCMA system, the users’ SCMA codeword in
each dimension will be modulated to an Orthogonal Frequency Division Multiple
(OFDM) subcarrier before the air interface. The number of users is larger than
the number of subcarriers, which lead to the non-orthogonal feature. On the
other side, D2D communication is a promising technology for local
communications, which allows nearby devices to communicate without base
station (BS) or with limited BS involvement, and improves the link
reliability, spectral efficiency and system capacity [8]. It is expected that
the social network structure can be further expanded through D2D
communications, which supports the development of the Internet of things
(IOT). Many recent works have investigated D2D communications and SCMA
cellular communications hybrid network [9, 10, 11, 12, 13, 14, 15, 16, 17, 18,
19], and often focus on optimization problem [9, 11, 12, 13, 14, 15, 16, 17,
18, 19, 20].

When D2D communications are considered in uplink cellular network, two main
modes have been utilized. In simple mode, a proportion of available resources
are allocated to D2D users while both cellular users and D2D users utilize
SCMA. In another mode, SCMA is only employed by cellular users and some
subcarriers can be reused by D2D users’ transmission. Liu $et$ $al.$ [9] have
considered the analytical model for the SCMA enhanced cellular network
coexisting with D2D. Based on the approximate system area spectral efficiency
(ASE), the performance between two different modes has been compared and the
system design guidance along with SCMA codebook allocation scheme have been
proposed. In [10], the distance between D2D transmitter and the base station
has been restricted by specific methods, which improves the block error rate
(BER) performance as well as the convergence behavior.

Both simple mode and complex mode have been considered into the research of
optimization problem [11, 12, 15, 16]. A graph-based joint mode selection and
resource allocation algorithm with predefined threshold has been proposed to
maximize the system sum rate in [11], and [12] has further considered more
issues like partner assignment and admission control. Kim $et$ $al.$ [13] have
considered joint power and resource allocation by using heuristic and inner
approximation method, and they have proposed a two-phase algorithm with low
complexity. The Opportunistic SCMA codebook allocation (OSCA) [14] and the
interference-aware hyper-graph based codebook allocation (IAHCA) [15, 16] have
been proposed to obtain the suboptimal solution and reduce the cross-tier
interference between cellular users and D2D users. In OSCA, codebooks are
opportunistically assigned based on the channel conditions of the cellular
links and D2D links. While in IAHCA, hypergraph model is used to characterize
the interference, and each RB is allowed to be shared by one cellular user and
more than one D2D users. Moreover, a channel gain based mode selection
criterion as well as a greedy style partner assignment scheme have been
proposed in [16], while the total transmit power minimization has been
considered in [17]. Besides, the game based resource allocation scheme has
been proposed in [18], and the block-chain-based transaction flow has been
analyzed for resource allocation in [19]. Authors in [20] have utilized
Lagrangian relaxation (LR) method and the Dinkelbach method successively for
energy-efficient maximization problem.

Recent works for SCMA and NOMA capacity analysis are summarized as in [5, 21,
22, 23, 24, 25, 26]. Le $et$ $al.$ [21] have considered the low-density
spreading (LDS) channels with fading and analyzed the spectral efficiency. In
[22], the authors have analyzed and compared several NOMA schemes in code
domain and proposed the theoretical derivations. Besides, Shental $et$ $al.$
[5] have analyzed the closed form expression of low density code domain NOMA.

Moreover, Han $et$ $al.$ [23] have considered the downlink SCMA systems with
user grouping character. The power allocation scheme with maximum capacity has
been proposed for both intragroup and intergroup case. Then, the Lagrange dual
decomposition based algorithm has been proposed to find the near-optimal
solution. The downlink SCMA model is also utilized in [24], where the authors
have briefly analyzed the system capacity and separated the codebook
assignment and power allocation. Chen $et$ $al.$ [25] have derived the
capacity region of uplink SCMA system based on the theory of multiple access
channels and analyzed the influence between data rate, power allocation
function, and resource assignment matrix. Further, the common and individual
outage probability regions have been proposed and optimized with Lagrangian
dual method and adaptive algorithm. In [26], the authors have compared and
analyzed the channel capacities of SCMA and low-density signature multiple
access (LDSMA) schemes.

In recent related work, when sum rate maximization problems are considered,
concrete analysis for the expression of cellular users’ capacity is lacked
[12, 20, 23]. Besides, the complex mode which brings more mutual interference
between cellular users and D2D users is worth optimization for system sum
rate. This paper considers cellular uplink transmission simultaneously with
D2D communications, where SCMA is employed only for cellular users and D2D
users work in complex mode.

In this paper, our main contribution contains the derivation of SCMA capacity
bound for general SCMA system and the sum rate maximization problem. We first
try to derive the information-theoretic expression of the capacity for all
users in the hybrid network model. Since it is hard to derive a close form
expression of cellular capacity, we find the capacity bound for cellular users
based on general SCMA codeword structure. After the theoretic derivation and
analysis, we focus on the sum rate maximization problem by power optimization
allocation for different users. Since this is a non-convex optimization
problem, we propose an algorithm based on geometric programming (GP), which
transforms parts of the objective function and rewrite the problem in standard
GP form to get the optimal solution with convex optimization theory. By
iteration, the solution of GP problem will converge to the optimal solution of
the original problem. The convergence and gain of the proposed algorithm are
conducted by simulation part. The performance shows that this algorithm is
appropriate and significant.

The rest of the paper is organized as follows. The network model and SCMA are
described in Section 2. The theoretic analysis of system capacity for cellular
users and D2D users are presented in Section 3. In Section 4, the optimization
problem with power allocation is formulated and the GP based iterative
algorithm is proposed and introduced in detail. Section 5 devotes to the
numerical results, and the paper finally concludes in Section 6.

Throughout this paper, the following notations will be used [27]. The bold
lowercase letter $\bm{x}$ denotes a column vector, and a matrix is represented
by a bold uppercase letter $\bm{X}$. $\bm{I}$ denotes the identity matrix. The
superscript $(\cdot)^{T}$ denotes matrix transpose and $(\cdot)^{*}$ denotes
conjugate transpose matrix. $diag(\bm{x})$ denotes a diagonal matrix with the
diagonal entries being vector $\bm{x}$. $E[x]$ denotes the expectation of a
random variable $x$. For a set $A$, $|A|$ denotes the number of elements, and
$A\setminus n$ denotes the set $A$ in which the element $n$ is excluded.

##  2 System Model

In this section, we introduce the hybrid network with cellular users and D2D
users. The factor graph and mapping matrix are also introduced for SCMA, which
is essential for capacity analysis. Then, the expression of SCMA received
signal is proposed.

###  2.1 Hybrid Network

Figure 1 shows a hybrid network example of a cellular coexisting with two D2D
pairs [27]. In our proposed network, we consider $J$ cellular users and
$J_{D}$ D2D pairs. The $i$th cellular user is denoted as CUi, and the D2D
transmitter as well as the receiver of the $i$th D2D pair are denoted as DTi
and DRi, respectively. Cellular users are allocated with SCMA codewords spread
over OFDM tones and transmitted to the air interface. While for each D2D
transmitter, one of the same OFDM subcarriers is occupied.

Figure 1: The cellular and D2D hybrid network.

###  2.2 SCMA Structure

The SCMA encoder can be defined as a mapping from $\log_{2}(M)$ coded bits to
a $K$-dimensional complex codebook of size $M$. We select $N$ dimensions from
$K$ for the sparse codebook design, in other words, each $K$-dimensional
complex codeword has only $N$-dimensional non-zero elements. In SCMA, $J$
cellular users are allocated with $K$ OFDM subcarriers. In this model, we
consider that each D2D pair will employ one of the $K$ resource blocks (OFDM
subcarriers) so that $J_{D}\leq K$. For each CUj, we introduce an indicator
vector $\bm{f}_{j}=(f_{1j},f_{2j},\dots,f_{Kj})^{T}$ of size $K$, which $k$th
elements ${f}_{kj}$ is defined as

$f_{kj}=\left\\{\begin{aligned} 0,&\quad x_{jk}=0,\\\ 1,&\quad x_{jk}\neq
0,\end{aligned}\right.$ (1)

for $j=1,2,...,J$, where ${x}_{jk}$ is the data of user $j$ transmitting on
the subcarrier $k$. $\bm{F}=(\bm{f}_{1},\bm{f}_{2},...,\bm{f}_{J})$ denotes
the corresponding indicator matrix, where the set of non-zero entries in a
column corresponds to the subcarriers that a cellular user will transmit over,
and the set of nonzero entries in a row denotes the users who collide on the
same subcarrier, The Forney factor graph can depict the structure of a
indicator matrix [28]. Figure 2 shows a factor graph with $J=6$, $K=4$, and
$N=2$, and Eq. (2) is the corresponding indicator matrix. The size of the
matrix determines the number of layer nodes and resource nodes, and when
${f}_{kj}=1$, the corresponding nodes are connected in the factor graph.
Define two sets $\xi_{k}=\\{j|{f}_{kj}\neq 0\\}$ for $k=1,...,K$ and
$\zeta_{j}=\\{k|{f}_{kj}\neq 0\\}$ for $j=1,...,J$. Besides,
$|\xi_{k}|\triangleq d_{c}$ and $|\zeta_{j}|\triangleq d_{f}$ denote that each
subcarrier supports $d_{c}$ users and each user occupies $d_{f}$ subcarriers,
respectively. In the proposed network model, $J_{D}$ of $K$ resource nodes are
also used by D2D communications.

Figure 2: Factor graph with $J=6$, $K=4$, and $N=2$.

$\bm{F}=\left[\begin{array}[]{cccccc}1&1&1&0&0&0\\\ 1&0&0&1&1&0\\\
0&1&0&1&0&1\\\ 0&0&1&0&1&1\end{array}\right].$ (2)

The design of SCMA codewords has been regarded as the optimization problem for
the mapping matrix and the SCMA encoder [7, 29]. Let a non-zero SCMA codeword
$\bm{x}_{N,j}$ be $N$-dimensional for CUj, and $\bm{u}_{2N,j}$ be the
equivalent $2N$-dimensional QAM constellation [27]. The codeword can be
written as

$\bm{x}_{N,j}=\Delta_{j}\cdot(\bm{E}_{r}+i\cdot\bm{E}_{i})\cdot\bm{M}\cdot\bm{u}_{2N,j},$
(3)

where $\Delta_{j}$ denotes the constellation operators for CUj, i.e. phase
rotation, $\bm{M}$ denotes the $2N\times 2N$ unitary rotation matrix,
$\bm{E}_{r}$ and $\bm{E}_{i}$ are the $N\times 2N$ matrices which can select
components from vector $\bm{u}_{j}$ that corresponds to the real part and
imaginary part of QAM symbols, respectively, and $i=\sqrt{-1}$ is the
imaginary unit.

###  2.3 SCMA Transmission

In the SCMA enhanced cellular uplink transmission, the received signal
$\bm{y}$ at BS contains cellular users’ data as well as D2D users’ data,
expressed as

$\bm{y}=\sum_{j=1}^{J}diag(\bm{h}_{j}^{cb})\bm{x}_{j}+\sum_{j=1}^{J_{D}}h^{db}_{j}\bm{x}^{{}^{\prime}}_{j}+\bm{n},$
(4)

where $\bm{h}_{j}^{cb}=(h_{j1}^{cb},h_{j2}^{cb},...,h_{jK}^{cb})^{T}$ is the
channel vector between CUj and BS, and
$\bm{x}_{j}=(x_{j1},x_{j2},...,x_{jK})^{T}$ is the SCMA codeword vector of
CUj. Besides, $h^{db}_{j}$ is the channel state information (CSI) from DTj to
BS, and
$\bm{x}^{{}^{\prime}}_{j}=(x_{j1}^{{}^{\prime}},x_{j2}^{{}^{\prime}},...,x_{jK}^{{}^{\prime}})^{T}$
is the D2D codeword vector which is only non-zero in a certain row. $\bm{n}$
is the additive white Gaussian noise vector with mean $\bm{0}$ and covariance
matrix $N_{0}\bm{I}$, denoted as
$\bm{n}\sim\bm{\textit{CN}}(\bm{0},N_{0}\bm{I})$.

##  3 Capacities of Cellular Users and D2D Users

In this section, we focus on the achievable rates for cellular users and D2D
users and derive the information-theoretic expression, respectively. Besides,
we provide specific analysis for the capacity bound of cellular users and for
D2D users the closed form expression is given.

###  3.1 Capacity of Cellular Users

Since $\bm{n}$ is a complex Gaussian random vector, we can get the
distribution for the noise in subcarrier $k$, that is,
$n_{k}\sim\bm{\textit{CN}}(0,N_{0})$. Assume that subcarrier $k$ is occupied
by DTk, and define $E[x_{k}^{{}^{\prime}}(x_{k}^{{}^{\prime}})^{*}]\triangleq
P_{k}^{{}^{\prime}}$. Then according to Eq. (4), we have

$x_{k}^{{}^{\prime}}\sim\bm{\textit{CN}}(0,P_{k}^{{}^{\prime}}).$ (5)

For cellular users’ communications, the equivalent noise term includes D2D
interference and noise. Thus, the equivalent noise $\tilde{n}_{k}\triangleq
h_{k}^{db}x_{k}^{{}^{\prime}}+n_{k}$ and

$\tilde{n}_{k}\sim\bm{\textit{CN}}(0,N_{0}+|h_{k}^{db}|^{2}P_{k}^{{}^{\prime}}).$
(6)

Let $\bm{K}_{\tilde{\bm{n}}}$ be the covariance matrix of
$\tilde{\bm{n}}=(\tilde{n}_{1},\dots,\tilde{n}_{K})^{T}$, and define
$\tilde{N}_{k}\triangleq N_{0}+|h_{k}^{db}|^{2}P_{k}^{{}^{\prime}}$ for
simplicity. Then we have
$\bm{K}_{\tilde{\bm{n}}}=diag\\{\tilde{N}_{1},\tilde{N}_{2},...,\tilde{N}_{K}\\}$.
Rewrite Eq. (4) as

$\bm{y}=\bm{H}\bm{x}+\tilde{\bm{n}},$ (7)

where $\bm{H}=(\bm{H}_{1},...,\bm{H}_{J})$,
$\bm{H}_{j}=diag\\{h_{j1}^{cb},...,h_{jK}^{cb}\\}$, for $j=1,...,J$, and
$\bm{x}=(\bm{x}_{1}^{T},\bm{x}_{2}^{T},...,\bm{x}_{J}^{T})^{T}$. Let
$\bm{K}_{\bm{x}}$ be the covariance matrix of $\bm{x}$. Then the covariance
$\bm{K}_{\bm{y}}$ of $\bm{y}$ is

$\bm{K}_{\bm{y}}=\bm{K}_{\tilde{\bm{n}}}+\bm{H}\bm{K}_{\bm{x}}\bm{H}^{*}.$ (8)

By Shannon’s theorem, the capacity is calculated based on the mutual
information. The channel in Eq. (7) is considered as random and time-
invariant. By [30], for a specific realization, we have the following analysis
process as shown in Eq. (9),

$\displaystyle I(\bm{x};\bm{y}|\bm{H}=\textup{H})$ (9)
$\displaystyle=h(\bm{y}|\bm{H}=\textup{H})-h(\tilde{\bm{n}})$
$\displaystyle\leq\log\det(\pi e\bm{K}_{\bm{y}})-\log\\{(\pi
e)^{K}\prod_{k=1}^{K}\tilde{N}_{k}\\}$
$\displaystyle=\log\det(\bm{K}_{\bm{y}})-\sum_{k=1}^{K}\log(N_{0}+|h_{k}^{db}|^{2}P_{k}^{{}^{\prime}}),$

where the equality in Eq. (9) holds when the input $\bm{x}$ follows Gaussian
input. The maximum mutual information is not a closed form expression as the
determinant is not explicit. Let
$\bm{K}_{\bm{x}_{j}}=E[\bm{x}_{j}\bm{x}_{j}^{*}]$. Then it is easy to show
that $\bm{K}_{\bm{x}_{j}}$ is a Hermitian matrix and
$\bm{K}_{\bm{x}}=diag\\{\bm{K}_{\bm{x}_{1}},\cdots,\bm{K}_{\bm{x}_{J}}\\}$ is
a block diagonal matrix.

Notice that $\bm{K}_{\bm{x}_{j}}$ is not a diagonal matrix in general, it is
valuable to derive the capacity bounds with Eq. (9). We introduce two lemmas
for further analysis.

######  Lemma 3.1.

Let the $N$-dimensional square matrices $\bm{Q_{1}}$, $\bm{Q_{2}}$ be
Hermitian. The eigenvalues of $\bm{Q_{1}}$, $\bm{Q_{2}}$, and
$\bm{Q_{1}}+\bm{Q_{2}}$ are defined as ${\lambda_{n}(\bm{Q_{1}})}_{n=1}^{N}$,
${\lambda_{n}(\bm{Q_{2}})}_{n=1}^{N}$, and
${\lambda_{n}(\bm{Q_{1}}+\bm{Q_{2}})}_{n=1}^{N}$, which are ordered as
$\lambda_{1}\leq\lambda_{2}\leq...\leq\lambda_{N}$. Then for $n=1,...,N$

$\displaystyle\lambda_{n}(\bm{Q_{1}})+\lambda_{1}(\bm{Q_{2}})$
$\displaystyle\leq\lambda_{n}(\bm{Q_{1}}+\bm{Q_{2}})$ (10)
$\displaystyle\leq\lambda_{n}(\bm{Q_{1}})+\lambda_{N}(\bm{Q_{2}}).$

######  Lemma 3.2.

Let the $N$-dimensional square matrices $\bm{Q}$, $\bm{S}$ be Hermitian and
nonsingular, respectively. The eigenvalues of $\bm{Q}$ and
$\bm{S}\bm{Q}\bm{S}^{*}$ are defined as ${\lambda_{n}(\bm{Q})}_{n=1}^{N}$ and
${\lambda_{n}(\bm{S}\bm{Q}\bm{S}^{*})}_{n=1}^{N}$. The singular values of
$\bm{S}$ are defined as $0<\sigma_{1}\leq\sigma_{2}\leq...\leq\sigma_{N}$. For
each $n=1,...,N$, there is a real number $\theta_{n}$ $\in$
$[\sigma_{1}^{2},\sigma_{N}^{2}]$ such that

$\lambda_{n}(\bm{S}\bm{Q}\bm{S}^{*})=\theta_{n}\lambda_{n}(\bm{Q}).$ (11)

Notice that Lemma 3.1 is the classic Weyl’s Theorem and Lemma 3.2 can be
easily derived by Lemma 3.1 in [31], the proof is omitted here.

In order to find the upper bound and the lower bound of
$\log\det(\bm{K}_{\tilde{\bm{n}}}+\bm{H}\bm{K}_{\bm{x}}\bm{H}^{*})$ by Eq.
(8), (9), we need to analyze the eigenvalues. Define
${\lambda_{k}(\bm{Q})}_{k=1}^{K}$ as the eigenvalues of an $K$-dimensional
square matrix $\bm{Q}$ in nondecreasing order, then we produce the analysis in
three aspects.

####  3.1.1 Upper Bound

By using Lemma 3.1, we have

$\displaystyle\lambda_{k}(\bm{K}_{\tilde{\bm{n}}}$
$\displaystyle+\bm{H}\bm{K}_{\bm{x}}\bm{H}^{*})$ (12)
$\displaystyle\leq\lambda_{k}(\bm{K}_{\tilde{\bm{n}}})+\lambda_{K}(\bm{H}\bm{K}_{\bm{x}}\bm{H}^{*})$
$\displaystyle=\tilde{N}_{k}+\lambda_{K}\left(\sum_{j=1}^{J}\bm{H}_{j}\bm{K}_{\bm{x}_{j}}\bm{H}_{j}^{*}\right)$
$\displaystyle\leq\tilde{N}_{k}+\sum_{j=1}^{J}\lambda_{K}(\bm{H}_{j}\bm{K}_{\bm{x}_{j}}\bm{H}_{j}^{*}).$

Notice that the first equality in Eq. (12) occurs if and only if there is
nonzero vector $\bm{\alpha}$ such that
$\bm{K}_{\tilde{\bm{n}}}\bm{\alpha}=\lambda_{k}(\bm{K}_{\tilde{\bm{n}}})\bm{\alpha}$,
$\bm{H}\bm{K}_{\bm{x}}\bm{H}^{*}\bm{\alpha}=\lambda_{K}(\bm{H}\bm{K}_{\bm{x}}\bm{H}^{*})\bm{\alpha}$,
and
$(\bm{K}_{\tilde{\bm{n}}}+\bm{H}\bm{K}_{\bm{x}}\bm{H}^{*})\bm{\alpha}=\lambda_{k}(\bm{K}_{\tilde{\bm{n}}}+\bm{H}\bm{K}_{\bm{x}}\bm{H}^{*})\bm{\alpha}$.
In general condition, $\bm{K}_{\tilde{\bm{n}}}$ and
$\bm{H}\bm{K}_{\bm{x}}\bm{H}^{*}$ should have no common eigenvector, then the
inequalities in Eq. (12) are strict inequalities.

Since $\bm{K}_{\bm{x}_{j}}$ is Hermitian, $\bm{H}_{j}$ is nonsingular, and
$\bm{H}_{j}\bm{H}_{j}^{*}=diag\\{|h_{j1}^{cb}|^{2},|h_{j2}^{cb}|^{2}...,|h_{jK}^{cb}|^{2}\\}$,
by using Lemma 3.2, for CUj we have

$\lambda_{K}(\bm{H}_{j}\bm{K}_{\bm{x}_{j}}\bm{H}_{j}^{*})\leq\lambda_{K}(\bm{K}_{\bm{x}_{j}})\cdot\max\limits_{k}(|h_{jk}^{cb}|^{2}).$
(13)

According to the SCMA codebook design in Eq. (3), we introduce a $K\times N$
matrix $\bm{V}_{j}$ to transform non-zero $N$-dimensional codeword into
$K$-dimensional codeword with $K-N$ zero elements. Then we have

$\displaystyle\bm{x}_{j}$ $\displaystyle=\bm{V}_{j}\cdot
e^{i\theta_{j}}\cdot(\bm{E}_{r}+i\cdot\bm{E}_{i})\cdot\bm{M}\cdot\bm{u}_{j}$
(14)
$\displaystyle=e^{i\theta_{j}}\cdot\bm{V}_{j}\cdot\bm{M^{{}^{\prime}}}\cdot\bm{u}_{j},$

where $e^{i\theta_{j}}$ denotes the phase rotation here and
$\bm{M}^{{}^{\prime}}\triangleq(\bm{E}_{r}+i\cdot\bm{E}_{i})\cdot\bm{M}$
denotes a $N\times 2N$ complex matrix. Then

$\displaystyle\bm{K}_{\bm{x}_{j}}$
$\displaystyle=E[e^{i\theta_{j}}(\bm{V}_{j}\bm{M^{{}^{\prime}}}\bm{u}_{j})(\bm{V}_{j}\bm{M^{{}^{\prime}}}\bm{u}_{j})^{*}e^{-i\theta_{j}}]$
(15) $\displaystyle=(\bm{V}_{j}\bm{M^{{}^{\prime}}})\cdot
E[\bm{u}_{j}(\bm{u}_{j})^{T}]\cdot(\bm{V}_{j}\bm{M^{{}^{\prime}}})^{*}.$

Let $\bm{V}_{j}\bm{M^{{}^{\prime}}}\triangleq(\bm{M}_{j1},\bm{M}_{j2})$ where
$\bm{M}_{j1}$ and $\bm{M}_{j2}$ are $K\times N$ matrices. Besides, let
$E[\bm{u}_{j}(\bm{u}_{j})^{T}]\triangleq diag[\bm{A}_{j1},\bm{A}_{j2}]$ where
$\bm{A}_{j1}$ and $\bm{A}_{j2}$ are $N\times N$ diagonal matrices in Gaussian
input. Then we have

$\bm{K}_{\bm{x}_{j}}=\bm{M}_{j1}\bm{A}_{j1}\bm{M}_{j1}^{*}+\bm{M}_{j2}\bm{A}_{j2}\bm{M}_{j2}^{*}.$
(16)

By using Lemma 3.1 repeatedly, we have

$\displaystyle\lambda_{K}(\bm{K}_{\bm{x}_{j}})$ (17)
$\displaystyle=\lambda_{K}(\bm{M}_{j1}\bm{A}_{j1}\bm{M}_{j1}^{*}+\bm{M}_{j2}\bm{A}_{j2}\bm{M}_{j2}^{*})$
$\displaystyle\leq\lambda_{K}(\bm{M}_{j1}\bm{A}_{j1}\bm{M}_{j1}^{*})+\lambda_{K}(\bm{M}_{j2}\bm{A}_{j2}\bm{M}_{j2}^{*})$
$\displaystyle\triangleq\tau_{K,j}.$

Combing Eq. (12), (13), and (17), we conclude that

$\displaystyle\lambda_{k}(\bm{K}_{\tilde{\bm{n}}}$
$\displaystyle+\bm{H}\bm{K}_{\bm{x}}\bm{H}^{*})$ (18)
$\displaystyle<\tilde{N}_{k}+\sum_{j=1}^{J}\tau_{K,j}\cdot\max\limits_{k}(|h_{jk}^{cb}|^{2}).$

This formula reveals the theoretic upper bound of each eigenvalue.

####  3.1.2 Lower Bound

Similarly, we can analyze the lower bound for
$\lambda_{k}(\bm{K}_{\tilde{\bm{n}}}+\bm{H}\bm{K}_{\bm{x}}\bm{H}^{*})$. By
using Lemma 3.1 and Lemma 3.2, we have

$\displaystyle\lambda_{k}(\bm{K}_{\tilde{\bm{n}}}$
$\displaystyle+\bm{H}\bm{K}_{\bm{x}}\bm{H}^{*})$ (19)
$\displaystyle\geq\lambda_{k}(\bm{K}_{\tilde{\bm{n}}})+\lambda_{1}(\bm{H}\bm{K}_{\bm{x}}\bm{H}^{*})$
$\displaystyle\geq\tilde{N}_{k}+\sum_{j=1}^{J}\lambda_{1}(\bm{H}_{j}\bm{K}_{\bm{x}_{j}}\bm{H}_{j}^{*})$
$\displaystyle\geq\tilde{N}_{k}+\sum_{j=1}^{J}[\lambda_{1}(\bm{K}_{\bm{x}_{j}})\cdot\min\limits_{k}(|h_{jk}^{cb}|^{2})].$

By Lemma 3.1 and Eq. (16), we have

$\displaystyle\lambda_{1}(\bm{K}_{\bm{x}_{j}})$ (20)
$\displaystyle=\lambda_{1}(\bm{M}_{j1}\bm{A}_{j1}\bm{M}_{j1}^{*}+\bm{M}_{j2}\bm{A}_{j2}\bm{M}_{j2}^{*})$
$\displaystyle\geq\lambda_{1}(\bm{M}_{j1}\bm{A}_{j1}\bm{M}_{j1}^{*})+\lambda_{1}(\bm{M}_{j2}\bm{A}_{j2}\bm{M}_{j2}^{*})$
$\displaystyle\triangleq\tau_{1,j}.$

Thus, we conclude in general condition that

$\displaystyle\lambda_{k}(\bm{K}_{\tilde{\bm{n}}}$
$\displaystyle+\bm{H}\bm{K}_{\bm{x}}\bm{H}^{*})$ (21)
$\displaystyle>\tilde{N}_{k}+\sum_{j=1}^{J}\tau_{1,j}\cdot\min\limits_{k}(|h_{jk}^{cb}|^{2}).$

####  3.1.3 Capacity Bound

We have

$\displaystyle C_{c}$
$\displaystyle=\max\limits_{\bm{K}_{\bm{x}}}I(\bm{x};\bm{y}|\bm{H}=\textup{H})$
(22)
$\displaystyle=\log\det(\bm{K}_{\bm{y}})-\sum_{k=1}^{K}\log(\tilde{N}_{k})$
$\displaystyle=\sum_{k=1}^{K}\log\lambda_{k}(\bm{K}_{\bm{y}})-\sum_{k=1}^{K}\log(\tilde{N}_{k}),$

where the derivation is based on $\bm{x}\sim\textit{CN}(0,\bm{K}_{\bm{x}})$.

According to Eq. (14), the arbitrariness of $\bm{V}_{j}$ and
$\bm{M}^{{}^{\prime}}$ makes the $\bm{K}_{\bm{x}_{j}}$ non-diagonal, however,
with Gaussian input we can define
$\bm{u}_{j}\sim\textit{CN}(0,\bm{K}_{\bm{u}_{j}})$ where
$\bm{K}_{\bm{u}_{j}}=diag\\{P_{j,1}^{u},P_{j,2}^{u},...,P_{j,2N}^{u}\\}$. Then
according to Eq. (17) we further analyze that

$\displaystyle\tau_{K,j}$ (23)
$\displaystyle=\lambda_{K}(\bm{M}_{j1}\bm{A}_{j1}\bm{M}_{j1}^{*})+\lambda_{K}(\bm{M}_{j2}\bm{A}_{j2}\bm{M}_{j2}^{*})$
$\displaystyle=\lambda_{K}(\bm{M}_{j1}diag\\{P_{j,1}^{u},...,P_{j,N}^{u}\\}\bm{M}_{j1}^{*})$
$\displaystyle+\lambda_{K}(\bm{M}_{j2}diag\\{P_{j,N+1}^{u},...,P_{j,2N}^{u}\\}\bm{M}_{j2}^{*}).$

To maximize the rate of cellular network, the upper bound of SCMA capacity
should be figured out based on Eq. (18), (22), and (23), which is summarized
in follows

$\displaystyle
C_{c}<\sum_{k=1}^{K}\log\left\\{1+\frac{\sum_{j=1}^{J}\tau_{K,j}\cdot\max\limits_{k}(|h_{jk}^{cb}|^{2})}{N_{0}+|h_{k}^{db}|^{2}P_{k}^{{}^{\prime}}}\right\\},$
(24)

where we should notice that, the bound for proposed network with general SCMA
system is also determined by SCMA codebook design and factor graph structure.
For more practical situation, SCMA codebook is designed based on optimal
rotation matrix [32]. Thus, the analysis of optimization problem formulation
in section IV needs some simplification.

###  3.2 Capacity of D2D users

For D2D users, each D2D pair is allocated with one of the $K$ OFDM tones. For
example, we assume that the $\ell$th D2D pair selects the $\ell$th OFDM tone
for $\ell=1,2,...,J_{D}$. The received signal of DRℓ is disturbed by some
cellular users’ data, which is shown as

$y_{\ell}=h_{\ell}^{dd}x_{\ell}^{{}^{\prime}}+\sum_{j\in\xi_{\ell}}h_{j\ell}^{cd}x_{j\ell}+n_{\ell},$
(25)

where $h_{\ell}^{dd}$ is the CSI between $\ell$th D2D pair, $h_{j\ell}^{cd}$
is the CSI between CUj and DRℓ, and the noise is still defined as
$n_{\ell}\sim\bm{\textit{CN}}(0,N_{0})$. Assume that

$x_{j\ell}\sim\bm{\textit{CN}}(0,P_{j\ell}),$ (26)

where $P_{j\ell}\triangleq E[x_{j\ell}x_{j\ell}^{*}]$. Let the equivalent
noise term
$n_{\ell}^{{}^{\prime\prime}}\triangleq\sum_{j\in\xi_{\ell}}h_{j\ell}^{cd}x_{j\ell}+n_{\ell}$
with the interference term. Then we have

$n_{\ell}^{{}^{\prime\prime}}\sim\bm{\textit{CN}}(0,N_{0}+\sum_{j\in\xi_{\ell}}|h_{j\ell}^{cd}|^{2}P_{j\ell}).$
(27)

We can rewrite the received signal at DRℓ as

$y_{\ell}=h_{\ell}^{dd}x_{\ell}^{{}^{\prime}}+n_{\ell}^{{}^{\prime\prime}}.$
(28)

Similarly, we can calculate the mutual information for D2D pair $\ell$ as

$\displaystyle
I(x_{\ell}^{{}^{\prime}};y_{\ell}|h_{\ell}^{dd})=h(y_{\ell}|h_{\ell}^{dd})-h(n_{\ell}^{{}^{\prime\prime}})$
(29) $\displaystyle\leq\log\left\\{\pi
e\left(|h_{\ell}^{dd}|^{2}P_{\ell}^{{}^{\prime}}+N_{0}+\sum_{j\in\xi_{\ell}}|h_{j\ell}^{cd}|^{2}P_{j\ell}\right)\right\\}$
$\displaystyle-\log\left\\{\pi
e\left(N_{0}+\sum_{j\in\xi_{\ell}}|h_{j\ell}^{cd}|^{2}P_{j\ell}\right)\right\\}$
$\displaystyle=\log\left\\{1+\frac{|h_{\ell}^{dd}|^{2}P_{\ell}^{{}^{\prime}}}{N_{0}+\sum_{j\in\xi_{\ell}}|h_{j\ell}^{cd}|^{2}P_{j\ell}}\right\\}.$

Thus, when the $\ell$th OFDM tone is chosen for the $\ell$th D2D pair,
according to Eq. (29), we have

$\displaystyle C_{d}$
$\displaystyle=\sum_{\ell=1}^{J_{D}}\log\left\\{1+\frac{|h_{\ell}^{dd}|^{2}P_{\ell}^{{}^{\prime}}}{N_{0}+\sum_{j\in\xi_{\ell}}|h_{j\ell}^{cd}|^{2}P_{j\ell}}\right\\}$
(30)
$\displaystyle=\sum_{\ell=1}^{J_{D}}\log\left(1+\gamma_{\ell}^{d}\right),$

where $\gamma_{\ell}^{d}$ denotes the signal-to-interference-plus-noise ratio
(SINR) for the $\ell$th D2D receiver.

##  4 Sum Rate Maximization

In this section, we focus on the sum rate maximization issue. We formulate the
optimization problem as the power allocation of cellular users and D2D users.
The GP based iterative algorithm is introduced to find the numerical solution.

###  4.1 Problem Formulation

In section 3, the analysis is based on general SCMA codebook. However,
according to the sparsity of SCMA codeword, most non-diagonal entries of
$\bm{K}_{\bm{x}_{j}}$ are zero. Actually, most SCMA system models utilize
optimal codebook design principles to obtain the diagonal covariance matrix
[28, 29, 33]. In this section, we consider the general situation as

$\displaystyle\bm{K}_{\bm{x}_{j}}$ $\displaystyle=E[\bm{x}_{j}\bm{x}_{j}^{*}]$
(31) $\displaystyle=diag\\{E[x_{j1}x_{j1}^{*}],\dots,E[x_{jK}x_{jK}^{*}]\\}$
$\displaystyle\triangleq diag\\{P_{j1},P_{j2},...,P_{jK}\\}.$

Notice that the maximum mutual information in Eq. (9) can get closed form
expression, we have the simplified derivation as

$\lambda_{k}(\bm{K}_{y})=\tilde{N}_{k}+\sum_{j\in\xi_{k}}|h_{jk}^{cb}|^{2}P_{jk},k=1,2,...,K,$
(32)

where we should notice that $P_{jk}=0$ for $j\notin\xi_{k}$ with the
corresponding factor graph. Then we have

$\displaystyle\log\det(\bm{K}_{y})$ (33)
$\displaystyle=\sum_{k=1}^{K}\log\left(N_{0}+|h_{k}^{db}|^{2}P_{k}^{{}^{\prime}}+\sum_{j\in\xi_{k}}|h_{jk}^{cb}|^{2}P_{jk}\right).$

Finally, the capacity of cellular users is

$C_{c}=\sum_{k=1}^{K}\log\left\\{1+\frac{\sum_{j\in\xi_{k}}|h_{jk}^{cb}|^{2}P_{jk}}{N_{0}+|h_{k}^{db}|^{2}P_{k}^{{}^{\prime}}}\right\\}.$
(34)

Besides, for decoding the codewords of the CUj at the BS, we define the SINR
in subcarrier $k\in\zeta_{j}$ as

$\gamma_{jk}^{c}=\frac{|h_{jk}^{cb}|^{2}P_{jk}}{N_{0}+|h_{k}^{db}|^{2}P_{k}^{{}^{\prime}}}.$
(35)

Accordingly, the sum rate maximization problem is formulated as

$\displaystyle\bm{P1:}$
$\displaystyle\max_{\bm{P},\bm{P^{{}^{\prime}}}}(C_{c}+C_{d})$ (36)
$\displaystyle s.t.$ $\displaystyle C1:\gamma_{jk}^{c}\geq\gamma_{0}^{c},\quad
j=1,\dots,J,\quad k\in\zeta_{j},$ $\displaystyle
C2:\gamma_{\ell}^{d}\geq\gamma_{0}^{d},\quad\ell=1,\dots,J_{D},$
$\displaystyle C3:0\leq P_{j}=\sum_{k\in\zeta_{j}}P_{jk}\leq P_{0},\quad
j=1,\dots,J,$ $\displaystyle C4:0\leq P_{\ell}^{{}^{\prime}}\leq
P_{0}^{{}^{\prime}},\quad\ell=1,\dots,J_{D},$

where the constraints $C1$ and $C2$ denote the minimum QoS requirements of
cellular users and D2D users, determined by the SINR threshold
$\gamma_{0}^{c}$ and $\gamma_{0}^{d}$, respectively, $C3$ and $C4$ mean the
transmitted power limitation. Since the objective function is nonlinear and
nonconvex, P1 is a nonconvex optimization problem. We decide to rely on the
tractable Geometric Programming problem to find the appropriate algorithm.

###  4.2 Geometric Programming

Geometric Programing (GP) is a special optimization problem with useful
properties [34]. The standard form of GP is nonlinear and nonconvex, while the
convex form of GP can make use of the standard convex optimization theory.
Define a monomial $g$ as

$g(\bm{x})=cx_{1}^{a_{1}}x_{2}^{a_{2}}...x_{n}^{a_{n}},$ (37)

where $c\geq 0$ is the multiplicative constant, $a_{\ell}\in\bm{R}$ for
$\ell=1,2,...,n$ is the exponents, and the argument $x_{\ell}>0$ for
$\ell=1,2,...,n$. Then a posynomial is defined as a summation of monomials
shown as

$G(\bm{x})=\sum_{m=1}^{M}c_{m}x_{1}^{a_{1}^{m}}x_{2}^{a_{2}^{m}}...x_{n}^{a_{n}^{m}},$
(38)

where $a_{\ell}^{m}\in\bm{R}$ and $c_{m}\geq 0$ for $\ell=1,2,\dots,n$ and
$m=1,2,\dots,M$. The standard form for a GP is as follows.

$\displaystyle\min_{\bm{x}}\quad G_{0}(\bm{x})$ (39) $\displaystyle s.t.$
$\displaystyle G_{s}(\bm{x})\leq 1,\quad s=1,\dots,S,$ $\displaystyle
h_{t}(\bm{x})=1,\quad t=1,\dots,T,$

where the posynomial $G_{s}(\cdot)$ and the monomial $h_{t}(\cdot)$ for
$s=0,1,\dots,S$ and $t=1,\dots,T$ are defined as

$G_{s}(\bm{x})=\sum_{m=1}^{M_{s}}c_{sm}x_{1}^{a_{1}^{sm}}x_{2}^{a_{2}^{sm}}...x_{n}^{a_{n}^{sm}},$
(40)

and

$h_{t}(\bm{x})=c_{t}x_{1}^{a_{1}^{t}}x_{2}^{a_{2}^{t}}...x_{n}^{a_{n}^{t}}.$
(41)

In Eq. (40) and Eq. (41), we have $c_{sm}>0$ and $c_{t}>0$ for
$m=1,2,...,M_{s}$, $s=0,1,...,S$, and $t=1,2,...,T$.

For the standard form of GP, the posynomials are not convex functions. In
order to simplify the problem, the logarithmic conversion as $y_{\ell}=\log
x_{\ell}$, $b_{sm}=\log c_{sm}$, and $b_{t}=\log c_{t}$ are utilized. Then the
equivalent problem in term of $\bm{y}$ is

$\displaystyle\min_{\bm{y}}\quad
G_{0}^{{}^{\prime}}(\bm{y})=\log\sum_{m=1}^{M_{0}}\exp(\bm{a}_{0m}^{T}\bm{y}+b_{0m})$
(42) $\displaystyle s.t.\quad
G_{i}^{{}^{\prime}}(\bm{y})=\log\sum_{m=1}^{M_{s}}\exp(\bm{a}_{sm}^{T}\bm{y}+b_{sm})\leq
0,$ $\displaystyle\quad\quad s=1,2,\dots,S,$ $\displaystyle\quad
h_{t}^{{}^{\prime}}(\bm{y})=\bm{a}_{t}^{T}\bm{y}+b_{t}=0,\quad t=1,2,\dots,T,$

where $\bm{a}_{sm}=(a_{1}^{sm},\cdots,a_{n}^{sm})^{T}$ and
$\bm{a}_{t}=(a_{1}^{t},\cdots,a_{n}^{t})^{T}$. The problem in (42) is the
convex form of GP.

###  4.3 The Proposed GP Based Iterative Algorithm

The objective function in $\bm{P1}$ is obviously neither a posynomial nor a
convex function. We consider the deformation of the objective function to
construct the GP problem. For the transmitted power of cellular users, we
further improve the constraint $C3$ as $P_{jk}\leq\frac{P_{0}}{d_{f}},\forall
k\in\zeta_{j}$, $\forall j$. According to Eq. (30), (34), we have

$\displaystyle C_{c}+C_{d}$ (43)
$\displaystyle=\sum_{k=1}^{K}\log(1+\sum_{j\in\xi_{k}}{\gamma_{jk}^{c}})+\sum_{\ell=1}^{J_{D}}\log\left\\{1+\gamma_{\ell}^{d}\right\\}$
$\displaystyle=\log\left\\{\prod_{k=1}^{K}\frac{g_{k}(\bm{p},\bm{p^{{}^{\prime}}})}{f_{k}(\bm{p^{{}^{\prime}}})}\cdot\prod_{\ell=1}^{J_{D}}\frac{g_{\ell}^{{}^{\prime}}(\bm{p},\bm{p^{{}^{\prime}}})}{f_{\ell}^{{}^{\prime}}(\bm{p})}\right\\},$

where

$\displaystyle f_{k}(\bm{p^{{}^{\prime}}})\triangleq
N_{0}+|h_{k}^{db}|^{2}P_{k}^{{}^{\prime}},$ (44) $\displaystyle
f_{\ell}^{{}^{\prime}}(\bm{p})\triangleq
N_{0}+\sum_{j\in\xi_{\ell}}|h_{j\ell}^{cd}|^{2}P_{j\ell},$ $\displaystyle
g_{k}(\bm{p},\bm{p^{{}^{\prime}}})\triangleq
N_{0}+|h_{k}^{db}|^{2}P_{k}^{{}^{\prime}}+\sum_{j\in\xi_{k}}|h_{jk}^{cb}|^{2}P_{jk},$
$\displaystyle g_{\ell}^{{}^{\prime}}(\bm{p},\bm{p^{{}^{\prime}}})\triangleq
N_{0}+\sum_{j\in\xi_{\ell}}|h_{j\ell}^{cd}|^{2}P_{j\ell}+|h_{\ell}^{dd}|^{2}P_{\ell}^{{}^{\prime}}.$

Note that the logarithmic function is a monotonic function, we can change the
maximization of the objective function to the minimization of its reciprocal.
Then problem $\bm{P1}$ can be transformed into

$\displaystyle\bm{P2:}$
$\displaystyle\min_{\bm{p},\bm{p^{{}^{\prime}}}}\quad\prod_{k=1}^{K}\prod_{\ell=1}^{J_{D}}\frac{f_{k}(\bm{p^{{}^{\prime}}})f_{\ell}^{{}^{\prime}}(\bm{p})}{g_{k}(\bm{p},\bm{p^{{}^{\prime}}})g_{\ell}^{{}^{\prime}}(\bm{p},\bm{p^{{}^{\prime}}})},$
(45) $\displaystyle s.t.$ $\displaystyle
C1:(\gamma_{jk}^{c})^{-1}\gamma_{0}^{c}\leq 1,j=1,\dots,J,k\in\zeta_{j},$
$\displaystyle C2:(\gamma_{\ell}^{d})^{-1}\gamma_{0}^{d}\leq
1,\ell=1,\dots,J_{D},$ $\displaystyle C3:P_{jk}(\frac{P_{0}}{d_{f}})^{-1}\leq
1,j=1,\dots,J,k\in\zeta_{j},$ $\displaystyle
C4:P_{\ell}^{{}^{\prime}}(P_{0}^{{}^{\prime}})^{-1}\leq 1,\ell=1,\dots,J_{D},$

where all the constraints are rewritten as posynomials with upper bounds.

The transformed objective function is a ratio between two posynomials, which
belongs to Complementary GP and is an intractable NP-hard problem [34].
However, the denominator posynomial can be approximated by a monomial such
that the ratio is a posynomial as shown in Lemma 4.1.

######  Lemma 4.1.

Let the $g(\bm{x})=\sum_{\ell}u_{\ell}(\bm{x})$ is a posynomial where
$u_{\ell}(\bm{x})$ is a monomial. The argument $\bm{x}$ is with all positive
elements. There must have a monomial $\tilde{g}(\bm{x})$ as

$g(\bm{x})\geq\tilde{g}(\bm{x})\triangleq\prod_{\ell}\left(\frac{u_{\ell}(\bm{x})}{\beta_{\ell}}\right)^{\beta_{\ell}},$
(46)

where $\beta_{\ell}>0$ and $\sum_{\ell}\beta_{\ell}=1$, and that becomes an
equality when $\beta_{\ell}=\frac{u_{\ell}(\bm{x})}{g(\bm{x})}$ for
$\forall\ell$.

######  Proof 4.2.

See Appendix.

By using Lemma 4.1, we can rely on GP to find the globally optimal solution.
$\bar{\bm{p}}\triangleq(\bm{p}^{T},(\bm{p}^{{}^{\prime}})^{T})^{T}$ denotes
the complete optimization variables in our algorithm. Specifically, we focus
on the expansion for the denominator of the objective function in $\bm{P2}$
and set an initial power vector $\bar{\bm{p}}_{0}$ to calculate a group of
coefficients. Then, the denominator is transformed into a monomial
$\tilde{g}(\bar{\bm{p}})$ with calculated coefficients. This problem is GP in
standard form, and we can modify it to a GP in convex form and use the
interior point method to solve. With the optimal solution $\bar{\bm{p}^{*}}$,
we solve different GP constantly to improve the accuracy for the approximation
of original problem. Notice that the Karush-Kuhn-Tucker (KKT) conditions is
satisfied obviously as we take derivatives of $g(\bar{\bm{p}})$ and
$\tilde{g}(\bar{\bm{p}})$, this is verified as the best local monomial
approximation [35]. The algorithm is convergent to the optimal solution of
$\bm{P2}$, and we will show numerical results in simulation section. The
proposed algorithm is summarized in Algorithm 1, and we define

$\displaystyle\prod_{k=1}^{K}\prod_{\ell=1}^{J_{D}}g_{k}(\bar{\bm{p}})g_{\ell}^{{}^{\prime}}(\bar{\bm{p}})\triangleq\sum_{s}u_{s}(\bar{\bm{p}})$
(47) $\displaystyle\triangleq\sum_{s}c_{s}(P_{1}^{a_{s}^{1}}\cdots
P_{J}^{a_{s}^{J}})[(P^{{}^{\prime}}_{1})^{b_{s}^{1}}\cdots(P^{{}^{\prime}}_{J_{D}})^{b_{s}^{J_{D}}}],$

where $c_{s}>0$, and $a_{s}^{j},b_{s}^{\ell}\in\bm{R}$.

Algorithm 1 Iterative GP power allocation algorithm.

1: Initialize $\bar{\bm{p}}=\bar{\bm{p}}_{0}$ and set the maximum iterative
number $T^{\max}$.

2: Set $I=1$.

3: while $I\leq T^{\max}$ do

4: Compute coefficients $\beta_{s}$
by$\beta_{s}=u_{s}(\bar{\bm{p}})/[\prod_{k=1}^{K}\prod_{\ell=1}^{J_{D}}g_{k}(\bar{\bm{p}})g_{\ell}^{{}^{\prime}}(\bar{\bm{p}})]$.

5: Substitute the denominator of the objective function in $\bm{P2}$ by
$\tilde{g}(\bar{\bm{p}})$,
where$\tilde{g}(\bm{p},\bm{p^{{}^{\prime}}})=\prod_{s}\left[\frac{u_{s}(\bm{p},\bm{p^{{}^{\prime}}})}{\beta_{s}}\right]^{\beta_{s}}$.

6: Solve the new GP optimization problem in convex optimization methods with
the objective function
$\frac{\prod_{k=1}^{K}\prod_{\ell=1}^{J_{D}}f_{k}(\bm{p^{{}^{\prime}}})f_{\ell}^{{}^{\prime}}(\bm{p})}{\tilde{g}(\bm{p},\bm{p^{{}^{\prime}}})}$
and the same constraints as in $\bm{P2}$.

7: Update $\bar{\bm{p}}={\bar{\bm{p}}^{*}}$, where ${\bar{\bm{p}}^{*}}$ is the
optimal solution of the new GP problem.

8: $I=I+1$.

9: end while

10: return The optimal solution ${\bar{\bm{p}}^{*}}$ of $\bm{P2}$.

The main complexity in Algorithm 1 for each iteration reflects on the
computation of coefficients $\bm{\beta}$ and the process of solving GP.
Actually, the computation of new function
$\tilde{g}(\bm{p},\bm{p^{{}^{\prime}}})$ requires
$(K+J_{D})(d_{c}+1)+\sum_{s}(J+J_{D}+3)$ multiplications and
$\sum_{s}(J+J_{D}+1)$ exponentiations. Besides, GP is usually solved by
interior point methods, which is proved to have polynomial time complexity
[34, 36].

##  5 Simulation Results

In this section, we perform simulations to demonstrate the performance of our
power allocation algorithm. Notice that uplink SCMA system with large users
can be resolved in small-scale SCMA structure with same overloading feature,
the SCMA scheme with $J=6$, $K=4$, and $N=2$ is claimed to be the basic
regular SCMA system. Our simulation scenario utilizes the basic SCMA scheme
with the same factor graph in figure 2. Considering the decoding performance
of D2D communications in proposed network [27], we set no more than two D2D
pairs in this scenario. Besides, notice that the small modulation and simple
channel coding are preferable as the devices in massive machine communications
are low-complexity [37], we adopt Quadrature Phase Shift Keying (QPSK) for
cellular users here. Thus, each cellular user’s total power should be equally
allocated to two subcarriers according to the SCMA codebook structure. The
channels for both cellular communications and D2D communications are generated
with a normalized Rayleigh fading component and a distance-dependent path loss
model. Table 1 summarizes the simulation parameters [16, 33]. Different scenes
and methods are compared to prove that this algorithm is appropriate for this
specific optimization problem.

Table 1: Simulation parameters. No. CU, $J$ | 6
---|---
No. subcarriers, $K$ | 4
No. non-zero dimension, $N$ | 2
Noise power, $N_{0}$ | -174dBm/Hz
Power limitation of CU, $P_{0}$ | 30dBm
Power limitation of DU, $P_{0}^{{}^{\prime}}$ | 30dBm
SINR of CU, $\gamma_{0}^{c}$ | 0dB
SINR of DU, $\gamma_{0}^{d}$ | 10dB
Bandwidth, $B$ | 180kHz
Initial Power of CU, $p_{0}$ | $P_{0}/2$
Initial Power of DU, $p_{0}^{{}^{\prime}}$ | $P_{0}^{{}^{\prime}}/2$
Cell radius | 500m
Distance of each D2D pair | 1$\sim$20m
Path loss Model of CU | 37.6log10($d_{C}$[km])+128.1
Path loss Model of D2D | 40log10($d_{D}$[km])+148

###  5.1 Convergence

The convergence of this algorithm is shown from figure 3 to figure 5. Figure 3
illustrates the power allocation versus number of iterations by Algorithm 1,
where CU and DU denote cellular user and D2D user, respectively. The number of
D2D pairs $J_{D}$ is $1$ in figure 3 with a specific numerical result from the
tests. The figure shows that the power allocation algorithm converges very
fast and achieves the optimal power only after three iterations. We can
analyze that after few series of approximation, the solutions converge to a
point satisfying the KKT conditions of the origin problem in $\bm{P2}$. Thus,
transforming complementary GP into standard GP by this algorithm often
performs a fast convergence rate. Besides, two users are allocated with
maximal power from all seven users as they may suffer from relatively small
interference.

Moreover, we consider the scenarios with two pairs of D2D users, which of
course will improve the computational complexity of the optimization problem
as the variables increase. We can observe from figure 4 that only one CU is
allocated with maximal power as the interference among cellular users and D2D
users increases. Besides, the power allocation converges very close to the
optimal power allocation by the $2$th iterations, which still proves that this
algorithm performs a fast convergence rate.

Figure 3: Convergence of power allocation, $J_{D}=1$. Figure 4: Convergence of
power allocation, $J_{D}=2$. Figure 5: Convergence of system sum rate.

The sum rate versus number of iterations is shown in figure 5, which compares
the scenarios with different D2D pairs. The algorithm always converges within
five iterations while the rate achieves maximum after three iterations from
figure 5 in both scenarios. From this aspect, the algorithm is also proved to
have a better convergence rate. Besides, the increase of D2D pairs makes the
performance of sum rate rises about $1.2\%$ while other parameters remain
same.

Figure 6: Algorithm performance comparison versus cellular power limitation
with $J_{D}=1$. Figure 7: Algorithm performance comparison versus D2D’s power
limitation with $J_{D}=1$. Figure 8: Algorithm performance comparison versus
cellular power limitation with $J_{D}=2$. Figure 9: Algorithm performance
comparison versus D2D’s power limitation with $J_{D}=2$.

###  5.2 Algorithm Performance

We compare the performance of the proposed algorithm with random power
allocation from figure 6 to figure 9. We set different power limitation values
for both cellular users and D2D users to observe the performance of final
system sum rate. The sum rate versus the power limitation of cellular users
$P_{0}$ is shown in figure 6, while figure 7 illustrates system sum rate
versus power limitation of D2D user, $P_{0}^{{}^{\prime}}$. Both the results
in figure 6 and figure 7 utilize the scenario of $J_{D}=1$. The two figures
prove that the proposed iterative GP based algorithm achieves better sum rate
performance than the random power allocation, which demonstrates the advantage
of the proposed algorithm. That is because the iterative algorithm always
converges to the optimal solution to the original optimization problem.
Besides, it can be founded that the increase of threshold for cellular users
will improve the performance of sum rate. However, when the threshold for D2D
users increases, the sum rate will decrease. The reason is that the rate of
cellular users always has larger effect on the system sum rate, and the less
D2D’s power allocation could give improvement to the whole system while the
communication of D2D users must be guaranteed.

Figure 8 and figure 9 show the similar comparison for the system sum rate
performance with the scenario of $J_{D}=2$. We can observe that the gap of
performances between proposed algorithm and random scheme becomes larger
compared to the scenario of $J_{D}=1$. It means that the effectiveness of our
proposed algorithm is more prominent with much mutual interference case.
Besides, the variety of either cellular or D2D’s threshold doesn’t provide
obvious change to the performance of system sum rate. The probable reason is
that cellular communications and D2D communications have made equal effect on
the performance of system sum rate with much mutual interference.

##  6 Conclusion

We investigate the hybrid cellular network with D2D communications, where SCMA
schemes is employed for cellular users. To improve the performance of system
sum rate region, we first derive the capacity for cellular users as well as
D2D users. Then we consider the joint power optimization problem to maximize
the system sum rate and propose a GP based iterative algorithm. Simulation
results prove the fast convergence and show the much improvement of proposed
algorithm compared with random power allocation. With the increasing number of
D2D users, the proposed algorithm performs prominent for the hybrid network
model.

## ACKNOWLEDGEMENT

This paper is supported by National key project 2018YFB1801102 and
2020YFB1807700, by NSFC 62071296, and STCSM 20JC1416502, 22JC1404000.

## Appendix

Here is the proof of Lemma 4.1. The arithmetic-geometric mean (AM-GM)
inequality is shown as

$\sum_{\ell}\beta_{\ell}x_{\ell}\geq\prod_{\ell}x_{\ell}^{\beta_{\ell}},$ (48)

where $x_{\ell}>0$, $\beta_{\ell}\geq 0$, and $\sum_{\ell}\beta_{\ell}=1$, for
$\forall\ell$. We utilize $y_{\ell}=\beta_{\ell}x_{\ell}$ and rewrite Eq. (48)
as

$\sum_{\ell}y_{\ell}\geq\prod_{\ell}(\frac{y_{\ell}}{\beta_{\ell}})^{\beta_{\ell}},$
(49)

with equality for $\beta_{\ell}=y_{\ell}/\sum_{\ell}y_{\ell}$, for
$\forall\ell$. This equality can be extended to the case of general monomial
and polynomial [35].

## References

  * Ltd. [2013] LTD. H T C.  5g: A technology vision[R].  [S.l.]: Huawei Technol. Co.. Ltd., Shenzhen, China, 2013.
  * Wu et al. [2015] WU Q, CHEN W, TAO M, et al.  Resource allocation for joint transmitter and receiver energy efficiency maximization in downlink ofdma systems[J].  IEEE Transactions on Communications, 2015, 63(2):416-430.
  * Dai et al. [2015] DAI L, WANG B, YUAN Y, et al.  Non-orthogonal multiple access for 5g: solutions, challenges, opportunities, and future research trends[J].  IEEE Journal on Selected Areas in Communications, 2015, 16(9):74-81.
  * Wu et al. [2018] WU Q, CHEN W, NG D, et al.  Spectral and energy-efficient wireless powered iot networks: Noma or tdma?[J].  IEEE Transactions on Vehicular Technology, 2018, 67(7):6663-6667.
  * Shental et al. [2017] SHENTAL O, ZAIDEL B, SHITZ S.  Low-density code-domain noma: better be regular[C]//IEEE Int. Symp. Inf. Theory (ISIT).  [S.l.]: IEEE, 2017: 2628-2632.
  * Chen et al. [2018a] CHEN Y, BAYESTEH A, WU Y, et al.  Toward the standardization of non-orthogonal multiple access for next generation wireless networks[J].  IEEE Commun. Magazine, 2018, 56(3):19-27.
  * Nikopour et al. [2013] NIKOPOUR H, BALIGH H.  Sparse code multiple access[C]//IEEE 24th Int. Symp. Pers. Indoor Mobile Radio Commun. (PIMRC).  [S.l.]: IEEE, 2013: 332-336.
  * Tehrani et al. [2014] TEHRANI M, UYSAL M, YANIKOMEROGLU H.  Device-to-device communication in 5g cellular networks: Challenges, solutions, and future directions[J].  IEEE Commun. Magazine, 2014, 52(5):86-92.
  * Liu et al. [2017] LIU J, SHENG M, LIU L, et al.  Modeling and analysis of scma enhanced d2d and cellular hybrid network[J].  IEEE Transactions on Communications, 2017, 65(1):173-185.
  * Li et al. [2016a] LI Y, SHA X, JIANG F, et al.  Sparse codes allocation in d2d communications underlaying cellular networks[C]//IEEE Int. Wireless Commun. Mobile Comput. Conf. (IWCMC).  [S.l.]: IEEE, 2016a: 948-952.
  * Zhao et al. [2016] ZHAO K, SHENG M, LIU J, et al.  Graph-based joint mode selection and resource allocation scheme for d2d and cellular hybrid network using scma[C]//IEEE Wireless Commun. Signal Process. (WCSP).  [S.l.]: IEEE, 2016: 1-5.
  * Li et al. [2016b] LI Y, SHENG M, ZHU Y, et al.  Sum rate maximization in underlay scma device-to-device networks[C]//IEEE Global Commun. (GLOBECOM).  [S.l.]: IEEE, 2016b: 1-6.
  * Kim et al. [2019] KIM H, NGUYEN H, KANG G, et al.  Device-to-device communication underlaying an uplink scma system[J].  IEEE Access, 2019, 7:21756-21768.
  * Zhao et al. [2015] ZHAO K, SHI Y, DAI Y, et al.  Resource allocation in device-to-device communication underlaid cellular network using scma: An opportunistic approach[C]//IEEE Int. Conf. Commun. China (ICCC).  [S.l.]: IEEE, 2015: 1-6.
  * Dai et al. [2016] DAI Y, SHENG M, ZHAO K, et al.  Interference-aware resource allocation for d2d underlaid cellular network using scma: A hypergraph approach[C]//IEEE Wireless Commun. Network. Conf. (WCNC).  [S.l.]: IEEE, 2016: 1-6.
  * Zhai et al. [2017] ZHAI D, SHENG M, WANG X, et al.  Energy-saving resource management for d2d and cellular coexisting networks enhanced by hybrid multiple access technologies[J].  IEEE Transactions on Wireless Communications, 2017, 16(4):2678-2692.
  * Sun et al. [2016] SUN Z, SHENG M, ZHAI D, et al.  Resource management for d2d underlaying cellular network with hybrid multiple access technologies[C]//IEEE Wireless Commun. Signal Process. (WCSP).  [S.l.]: IEEE, 2016: 1-5.
  * Liu et al. [2016] LIU Y, FANG X, YANG H, et al.  A mapping scheme of users to scma layers for d2d communication[C]//IEEE Veh. Technol. Conf. (VTC Spring).  [S.l.]: IEEE, 2016: 1-5.
  * Kong et al. [2019] KONG L, LING L, ZHANG X.  Scma based resource management of d2d communications for maximum sum-revenue[C]//Wireless Optical Commun. Conf. (WOCC).  [S.l.]: IEEE, 2019: 1-5.
  * Chen et al. [2019a] CHEN R, SHI F, ZHAO C.  A novel energy-efficient power allocation for d2d communications underlaying cellular network[C]//IEEE Wireless Commun. Signal Process. (WCSP).  [S.l.]: IEEE, 2019a: 1-6.
  * Le et al. [2018a] LE M, FERRANTE G, QUEK T, et al.  Fundamental limits of low-density spreading noma with fading[J].  IEEE Transactions on Wireless Communications, 2018, 17(7):4648-4659.
  * Le et al. [2018b] LE M, FERRANTE G, CASO G, et al.  On information-theoretic limits of code-domain noma for 5g[J].  IET. Commun., 2018, 12(15):1864-1871.
  * Han et al. [2019] HAN S, HUANG Y, MENG W, et al.  Optimal power allocation for scma downlink systems based on maximum capacity[J].  IEEE Transactions on Communications, 2019, 67(2):1480-1489.
  * Zhang et al. [2018] ZHANG Y, WANG X, WANG D, et al.  Capacity analysis and resource allocation of layered multicast in scma networks[C]//IEEE Wireless Commun. Network. Conf. (WCNC).  [S.l.]: IEEE, 2018: 1-6.
  * Chen et al. [2018b] CHEN J, WANG Z, XIANG W, et al.  Outage probability region and optimal power allocation for uplink scma systems[J].  IEEE Transactions on Communications, 2018, 66(10):4965-4980.
  * Chen et al. [2018c] CHEN S, PENG K, ZHANG Y, et al.  A comparative study on the capacity of scma and ldsma[C]//IEEE Int. Wireless Commun. Mobile Comput. Conf. (IWCMC).  [S.l.]: IEEE, 2018c: 1-5.
  * Liu et al. [2019] LIU Y, CHEN W, WEI F, et al.  Scma receiver design for cellular uplink transmission coexisted with d2d[C]//IEEE Wireless Commun. Signal Process. (WCSP).  [S.l.]: IEEE, 2019: 1-6.
  * Wei et al. [2017] WEI F, CHEN W.  Low complexity iterative receiver design for sparse code multiple access[J].  IEEE Transactions on Communications, 2017, 65(2):621-634.
  * Taherzadeh et al. [2014] TAHERZADEH M, NIKOPOUR H, BAYESTEH A, et al.  Scma codebook design[C]//IEEE Veh. Technol. Conf. (VTC Fall).  [S.l.]: IEEE, 2014: 1-5.
  * Tse et al. [2005] TSE D, VISWANATH P.  Fundamentals of wireless communication[M].  [S.l.]: Cambridge university press, 2005.
  * Horn et al. [2013] HORN R, JOHNSON C.  Matrix analysis[M].  [S.l.]: Cambridge university press, 2013.
  * Boutros et al. [1998] BOUTROS J, VITERBO E.  Signal space diversity: A power- and bandwidth-efficient diversity technique for the rayleigh fading channel[J].  IEEE Transactions on Information Theory, 1998, 44(4):1453-1467.
  * Li et al. [2016c] LI Z, CHEN W, WEI F, et al.  Joint codebook assignment and power allocation for scma based on capacity with gaussian input[C]//IEEE/CIC Int. Conf. Commun. China (ICCC).  [S.l.]: IEEE, 2016c: 1-6.
  * Chiang et al. [2007] CHIANG M, TAN C, PALOMAR D, et al.  Power control by geometric programming[J].  IEEE Transactions on Wireless Communications, 2007, 6(7):2640-2651.
  * Boyd et al. [2005] BOYD S, KIM S, VANDENBERGHE L, et al.  A tutorial on geometric programming[J].  Optimization and Engineering, 2005.
  * Nesterov et al. [1994] NESTEROV Y, NEMIROVSKY A.  Interior point polynomial methods in convex programming[M].  [S.l.]: SIAM Press, 1994.
  * Chen et al. [2019b] CHEN G, CUI Y, CHENG H, et al.  Analysis and optimization of successful symbol transmission rate for grant-free massive access with massive mimo[J].  IEEE Commun. Letters, 2019, 23(12):2381-2385.

lyk.epsYukai Liu received the B.S. degree from Shanghai Jiao Tong University
(SJTU) in 2016. He is currently pursuing the Ph.D. degree on Information and
Communication Engineering from SJTU. His current research interests include
wireless communications, with particular focus on multiple access and device-
to-device communications.

cw.epsWen Chen is a tenured professor at Shanghai Jiao Tong University. He is
a fellow of Chinese Institute of Electronics, and a distinguished lecturer of
IEEE Communications Society and IEEE Vehicular Technology Society. He is an
Editor of IEEE Transactions on Wireless Communications, IEEE Transactions on
Communications, IEEE Access and IEEE Open Journal of Vehicular Technology, and
the Shanghai chapter chair of IEEE Vehicular Technology Society. His research
interests include reconfigurable meta-surface, multiple access, wireless AI,
and green networks. He has published 100+ articles in IEEE journals and more
than 120 papers in IEEE Conferences, with citations more than 8000 in Google
scholar.


# Flatfield Ultrafast Imaging with Single-Shot Non-Synchronous Array
Photography

Matthew Sheinman  1,* Shyamsunder Erramilli  1,2 Lawrence Ziegler  2,3 Mi K.
Hong  1 and Jerome Mertz3,4  1Department of Physics, Boston University, 590
Commonwealth Avenue, Boston, Massachusetts 02215, USA  2Photonics Center,
Boston University, 8 St. Mary’s St., Boston, Massachusetts 02215, USA
3Department of Chemistry, Boston University, 590 Commonwealth Avenue, Boston,
Massachusetts 02215, USA  4Department of Biomedical Engineering, Boston
University, 44 Cummington Mall, Boston, Massachusetts 02215, USA
<EMAIL_ADDRESS>
###### Abstract

We present a method for acquiring a sequence of time-resolved images in a
single shot, called Single-Shot Non-Synchronous Array Photography (SNAP). In
SNAP, a pulsed laser beam is split by a diffractive optical element into an
array of angled beamlets whose illumination fronts remain perpendicular to the
optical axis. Different time delays are imparted to each beamlet by an
echelon, enabling them to probe a scene in rapid succession. The beamlets are
then imaged onto different regions of a camera by a lenslet array. Because the
illumination fronts remain flat (head-on) independently of beamlet angle, the
temporal resolution of SNAP is fundamentally limited only by the laser pulse
duration, akin to a "global shutter" in conventional imaging. We demonstrate
SNAP by capturing the evolution of a laser induced plasma filament for 20
frames at an average rate of 4.2 Tfps.

††journal: osac††articletype: Research Article

##  1 Introduction

Ultrafast imaging has the potential to reveal fundamental mechanisms in
chemistry, physics, and biology[1, 2, 3]. Often, ultrafast dynamics are
measured by repeated acquisitions. [4, 5, 6, 7]. However, this approach is not
feasible for events that are difficult or impossible to reproduce, such as
optical rogue waves[8], irreversible phase transitions[9], or laser-plasma
interactions[10], thus motivating the development of techniques that are
single shot.

Compressed imaging techniques have been successful at capturing ultrafast
scenes in a single shot. Compressed Ultrafast Photography (CUP) has the
advantage of being completely passive, that is, it can capture light emitted
from an ultrafast scene without the need of a structured probe beam [11].
However, the compressed nature of CUP imposes a sparsity requirement for image
reconstruction. In contrast, a variety of approaches allow the direct probing
of ultrafast events in a single shot without the need for reconstruction[12,
13, 14].

One approach involves spectral encoding, which associates portions of an
ultrafast pulse’s spectrum to different times. This approach enables a number
of temporally resolved images to be captured in a single acquisition with high
spatial resolution. A representative such technique is sequentially timed all-
optical mapping photography (STAMP), or its successor, SF-STAMP, which can
capture 25 frames in a single shot [15, 16]. However, because each image (time
point) is probed only by a portion of the original spectrum, temporal
resolution degrades as the number of images increases.

Another approach associates a unique spatial frequency to each time. This
disassociates temporal resolution from the number of frames, since each image
is acquired with the full pulse spectrum. In FRAME, a scene is probed by a
train of pulses with discrete spatial frequencies [17]. Other approaches
utilize digital holography [18, 19]. In particular, Time and Spatial-Frequency
Multiplexing (TSFM) has been used to capture up to 14 images in a single shot,
and could potentially capture more images relatively easily. TSFM uses a
diffractive optical element (DOE) to produce an array of reference beams for
off-axis holography, where each reference beam is assigned a unique time delay
by an echelon. However, because spatial frequency techniques involve
multiplexing, SNR becomes degraded as the number of images increases.

Yet another approach involves associating incidence angle to time, which has
the advantage of avoiding loss of temporal resolution while also not requiring
multiplexing. Angular encoding via pairs of two-by-two mirror arrays has been
demonstrated[20], though this technique is difficult to generalize to larger
frame numbers.

Here, we demonstrate that angular encoding by SNAP can easily capture a large
number of frames [21]. As in TSFM, we use a DOE and an echelon to produce an
array of time delayed beamlets. However, instead of using these beamlets as
reference beams for holography, we probe the scene directly. Because the DOE
preserves flat illumination fronts, each beamlet is incident on the entire
scene simultaneously[19]. We then project the images from each beamlet onto
different regions of a camera, where they are recorded in a single shot [22].
As a demonstration, we captured 20 shadowgraph images of the formation of a
laser induced plasma filament at an average framerate of 4.2 Tfps.

We note that while finalizing our manuscript for submission, we became aware
of very recent work by Zhu et al. based on a similar concept [23], though
without the benefit of flatfield (head-on) illumination provided by a DOE.

##  2 Principle of SNAP

Figure 1: Working principle of SNAP. (a) Schematic of setup. (b) The echelon
as viewed head on. The black indices (i,j) identify individual sections; the
yellow text denotes their thickness. (c) Fourier-plane array of beamlets (red
circles). Dashed circle represents the objective frequency support. To
maximize spatial resolution, the FPA should nearly fill the objective. (d)
Probe pulses generated by the DOE maintain flat illumination fronts
perpendicular to the optical axis. This prevents any "rolling shutter" effect,
even at the large incidence angles required to maximize spatial resolution.

The layout of our system is illustrated in Fig. 1a. SNAP resolves time delays
by probing a scene with an array of angles and then recording each of these
angles individually. This can be understood through well known principals in
optics[24]. A femtosecond pulse transmitted through a DOE produces an array of
pulsed beamlets $a_{i,j}(t)e^{i\vec{k}_{i,j}\cdot\vec{\rho}}$. Here, $(i,j)$
are indices denoting diffraction order, and $a_{i,j}$ are amplitudes for
spatially homogeneous ultrashort fields, but with intensities varying due to
the efficiency of the DOE. The transverse vector $\vec{k}_{i,j}$ corresponds
to the DOE diffraction angle $\theta_{i,j}$.

As shown in Fig. 1a, the DOE is imaged onto a scene by a 4f relay formed by
lenses $L_{1}$ and $L_{2}$. An echelon is placed at the Fourier plane of this
relay (fig. 1b)[9]. Each beamlet is focused through a unique section of this
echelon, imparting the unique time delay

$\tau_{i,j}=\frac{d_{i,j}}{c}(n_{g}-1)$ (1)

where $d_{i,j}$ is the section thickness, $n_{g}$ is the group velocity of the
echelon material, and $c$ is the speed of light in vacuum.

Immediately after the scene of the ultrafast event being imaged, the train of
beamlets is given by

$E_{sample}(\vec{\rho},t)=T(\vec{\rho},t)\sum_{i,j}a_{i,j}(t-\tau_{i,j})e^{i\vec{k}_{i,j}\cdot\vec{\rho}}$
(2)

where $T(\vec{\rho},t)$ is the complex amplitude transmittance of the
ultrafast scene.

From the Fourier shift theorem the field at the objective back focal plane is
then given by

$E_{FPA}(\rho,t)=\sum_{i,j}a_{i,j}(t-\tau_{i,j})\tilde{T}(\frac{\rho}{\lambda
f_{obj}}-\frac{k_{i,j}}{2\pi},t)$ (3)

That is, our array of angles has been transformed into an array of
displacements in the objective’s focal plane. The pitch $p_{FP}$ of the
beamlet array in this place is

$p_{FP}=\frac{\lambda
f_{obj}}{2\pi}\Delta{k}=f_{obj}\frac{f_{1}}{f_{2}}\theta_{0}$ (4)

where $\Delta{k}$ is the magnitude of the difference between adjacent wave
vectors and we have made use of the small-angle approximation.

To image each beamlet separately, we place a microlens array (MLA) in a 4f
configuration with our objective. If the pitch of the MLA is chosen to closely
match $p_{FP}$, the field at the camera plane becomes

$E_{camera}(\rho,t)=\sum_{i,j}a_{i,j}(t-\tau_{i,j})T(\frac{f_{obj}}{f_{MLA}}\rho-\frac{\lambda
f_{obj}}{2\pi}k_{i,j},t)$ (5)

Finally, because our camera exposure time is orders of magnitude longer than
our probe pulse duration, our signal is given by the time integral

$I_{camera}(\rho)=\int_{-\infty}^{\infty}\sum_{i,j}|a_{i,j}(t-\tau_{i,j})|^{2}|T(\frac{f_{obj}}{f_{MLA}}\rho-\frac{\lambda
f_{obj}}{2\pi}k_{i,j},t)|^{2}dt$ (6)

Here we have assumed that our beamlets are well separated on the camera
sensor, allowing us to omit cross-terms. Equation 6 describes an array of
images of the ultrafast scene intensity transmittance at different times. The
temporal resolution associated with each sub-image is simply the pulse
duration $\Delta t_{i,j}$ of each corresponding probe beamlet, which differ
slightly owing to dispersion from the echelon and the lenses $L_{1}$ and
$L_{2}$.

Equations 5 and 6 are approximations that do not account for the finite
aperture of each microlens. A more rigorous version of Eq. 5 includes a
frequency cutoff defined by the MLA pitch $p$

$P(\xi)=\begin{cases}1,&|\xi|<p/2\\\ 0,&|\xi|>p/2\end{cases}$ (7)

Applying equation 7 to equation 3, we can infer the system’s spatial
resolution from the maximum spatial frequency $f_{max}$ associated with each
sub-image

$f_{max}=\frac{p}{2\lambda f_{obj}}$ (8)

In effect, the objective’s frequency support is divided between each sub-image
(Fig. 1c). This introduces a trade-off between spatial resolution and number
of sub-images, as observed in earlier angle encoding techniques [20]. To
maximize spatial resolution then, the maximum probe beamlet tilt should
approach the objective numerical aperture.

The advantage of using a DOE to impart angular tilts to the beamlets can now
be understood (Fig 1c). The illumination fronts of the beamlets remain
perpendicular to the optical axis, regardless of tilt angle, meaning that each
beamlet is incident on the ultrafast scene head on. Without this advantage,
the system would suffer from a "rolling shutter" effect which would become
more pronounced for higher tilt angles. Indeed, with our existing setup, the
most tilted probe beamlets generated without a DOE would require a full on the
order of 1 ps to sweep through the FOV.

We note that Eq. 8 fails to take into account the vignetting of the FPA by our
MLA. Vignetting can cause a loss of spatial resolution and temporal cross-talk
between sub-images. We expect these effects to become more pronounced farther
from the center of our FOV. In principle, such effects can be mitigated by
departing from the 4f geometry and placing the MLA directly in the objective
back focal plane, at the expense of increased aberrations. In our current
setup, we elected to minimize aberrations and adopt the 4f geometry described
above.

##  3 Methods

###  3.1 Plasma Sample and Shadowgraphy

Figure 2: Imaging with SNAP. (a) A more detailed schematic. The camera can be
focused on the sample or displaced to produce a shadowgram. (b) Simultaneously
acquired raw images of the plasma shadowgram at different time points. (c)
Comparison of the resolution in sub-images (-1,-1) and (2,-1) when $\Delta$ is
near 0. (d) Same sub-images when $\Delta$ is about 15 mm.

To demonstrate SNAP, we image the formation of a laser induced plasma in air.
A diagram of the system’s pump arm for inducing the plasma formation is shown
in Fig. 2a. To better highlight SNAP’s ability to image features in both space
and time, the pump arm was designed to produce a plasma that exhibited spatial
structure. Specifically, high energy pulses focused with a NA greater than 0.3
are known to produce plasma filaments with two distinct foci [25].
Accordingly, a diverging lens is used to expand the pump beam and fill the
focusing lens, thus increasing NA. In addition, the pump arm includes a delay
stage to allow temporal calibration, as detailed below. Our pulse has a peak
wavelength of 800 nm with a measured pulse duration of 50 fs.

To enhance the signal contrast produced by the plasma filament, we perform
shadowgraph imaging by intentionally translating our camera approximately 3.5
mm beyond the nominal image plane (Fig. 2a). In this way, we image the
diffraction pattern (shadowgram) caused by the filament rather than the
filament itself.

To maximize SNR in as many sub-images as possible, we overexpose and discard
the (0,0) image. Further discarding the four sub-images in the corners, we
thus capture a total of 20 images, corresponding to 20 time points.

###  3.2 Fourier Plane Tuning

A more complete schematic of the system’s imaging arm is shown in Fig. 2a.
Additional relay lenses are placed between the objective and MLA. This serves
the practical purpose of imaging the objective’s focal plane, which would
otherwise be inaccessible, onto the MLA’s focal plane. In addition, this relay
can be used to tune the value of $p_{FP}$ in front of the MLA.

We use a DOE with a diffraction angle of $\theta_{0}=9.3^{\circ}$ at 808 nm
(HOLOEYE DOE 804). Our lenses have focal lengths $f_{1}$ = 50 mm, $f_{2}=$ 100
mm, $f_{3}=$ 100 mm, and $f_{4}=$ 70 mm. We use a $10\times$ objective and a
MLA with pitch $p=1$ mm, focal length $f_{MLA}=4.8$ mm (Advanced Microoptic
Systems APO-Q-P1000-R2.2). The displacement of the lens $L_{1}$ from its
nominal 4f configuration by an amount $\Delta$ is explained below.

Equation 5 imposes a requirement on the relation between DOE angle $\theta$,
MLA pitch $p$, and the optics chosen for the system. The above relay
introduces a new degree of freedom to this relation:

$p_{FP}=f_{obj}\frac{f_{1}f_{4}}{f_{2}f_{3}}\theta$ (9)

In our case $p_{FP}=$ 1.02 mm, which is slightly larger than our MLA pitch,
and may be responsible for irregularities in the resolution between sub-images
(Fig. 2b). To further tune our system, we displaced the lens $L_{2}$ and DOE
by an distance $\Delta$. With this displacement, a high degree of resolution
uniformity can be achieved (Fig 2c). We speculate that this improvement is due
to more refined tuning of the probe beam angles, and thus the pitch $p_{FP}$.

We note that, in addition to helping tune the pitch $p_{FP}$, displacing
$L_{2}$ causes the probe beams to become divergent. For this reason, we expect
the plasma shadowgram to expand more rapidly than the plasma itself.

###  3.3 Temporal Calibration

Figure 3: Temporal calibration of SNAP. (a) A representative sub-image after
background subtraction and normalization where the plasma has only partially
formed. (b) A line profile (averaged over 8 rows) is used to identify the
plasma leading edge. (c) Three of the 20 fits made for each sub-image, with
the dashed line representing a reference point near the filament center. (d)
Time stamps obtained from the temporal displacement of each of the 20 fits.

Our echelon is designed to provide an increase in glass thickness of 100 $\mu
m$ per step and a group index $n_{g}=1.53$. With these parameters, Eq. 1
predicts a time step of about 177 fs and frame rate of 5.6 Tfps. However, in
the actual experiment, deviations from these values are observed, which we
attribute to spherical abberations in the imaging relay formed by $L_{1}$ and
$L_{2}$.

To correct for these deviations, we calibrate our system by measuring the
relative arrival times of our probe beamlets. We do this by creating a
conventional, many-acquisition ultrafast video for each sub-image. First, we
acquire a background image by positioning our delay stage such that all
beamlets arrive at the sample plane before the plasma begins to form. We then
translate our delay stage in 25 $\mu m$ steps, each causing the plasma to form
167 fs later, and acquire the corresponding images. An example of one such
image (background subtracted) is shown in Fig. 3a. The pump energy is about 2
mJ per pulse.

To extract time delays from this data, we identify the leading edges of the
shadowgrams based on a threshold (Fig. 3b). The position of the leading edges
is then plotted against the probe beamlet arrival times. Example linear fits
are shown in Fig. 3c. The intersections of these fits with a location in the
center of the FOV then provide a measure of the probe beam relative arrival
times. These are plotted in Fig. 3d, enabling us to correct for the deviations
from a linear trend.

For this method of temporal calibration to be effective, a common spatial axis
must be found for each subimage. To achieve this, we place an obstruction in
the plane of the plasma and image its shadowgram. We then perform normalized
cross correlations between each of these shadowgrams, using subimage (0,1) as
a reference image. With this choice of reference, the mean peak cross
correlation value is .95.

##  4 Results

Figure 4: Evolution of a plasma filament. (a) Four frames taken from our 20
sub-images, with associated line-scans (averaged over 8 rows). The sequential
formation of two foci is evident. (b) A plot of the shadowgram leading edge
progression based on calibrated time stamps.

A pump energy of about 1 mJ per pulse was found to produce laser plasmas with
two clearly distinguishable foci, one leading and the other trailing. Example
line scans from four of our 20 frames are shown in Fig. 4a. At early times
when the plasma just begins to form and is highly localized, a sharp drop
behind the plasma’s leading edge is observed. As the plasma continues to form,
a second focus becomes apparent. Finally, at long times, the sharp drop
becomes a gradual roll-off at the edge of the second focus.

The single-shot data is fit to a line and compared with the data from the
many-acquisition temporal calibration (Fig. 4b). The single-shot expansion
rate is observed to be 310$\pm$20 um/ps. The largest sources of uncertainty in
our estimate are the shadow decay at the edge of our linescans and the
uncertainty in spatial calibration by cross correlation.

To compare with conventional many-acquisition imaging, we take the average of
the expansion rates measured during temporal calibration for the 8 inner
subimages, where abberations should be at a minimum. We find a value of
318$\pm$5 um/ps, consistent with the single-shot value. Discrepancies with the
known value of the speed of light are not examined, but may be due to the
previously mentioned divergence in the probe beam.

##  5 Conclusion

In this work, we demonstrated that SNAP is a viable single-shot alternative to
more conventional multi-shot ultrafast imaging techniques. There is a
fundamental tradeoff in SNAP between spatial resolution and number of temporal
samples. However, the temporal resolution of SNAP is limited only by the
duration of the probe pulse. Importantly, the use of a DOE to generate the
beamlet array means that the entire FOV of the ultrafast scene is probed
simultaneously, analogous to a "global shutter".

We imaged the evolution of a laser induced plasma filament both with SNAP and
with conventional many-acquisition imaging, and verified that the expansion
rates of the resultant plasma shadowgrams were consistent. Moreover, SNAP was
able to resolve a double focus structure in the filament, illustrating its
effectiveness at capturing spatiotemporal data.

SNAP uses readily available optical components that are easy to assemble, and
its generalization to a larger number of temporal samples is straightforward.
As such, SNAP can be useful for a wide variety of ultrafast imaging
experiments where a high temporal resolution is required.

## Funding

National Institute of Health R21GM128020.

## References

  * [1] T. Gorkhover, S. Schorb, R. Coffee, M. Adolph, L. Foucar, D. Rupp, A. Aquila, J. D. Bozek, S. W. Epp, B. Erk _et al._ , “Femtosecond and nanometre visualization of structural dynamics in superheated nanoparticles,” Nature photonics 10, 93–97 (2016).
  * [2] M. Imada, A. Fujimori, and Y. Tokura, “Metal-insulator transitions,” Reviews of modern physics 70, 1039 (1998).
  * [3] V. Tuchin, “Methods and algorithms for the measurement of the optical parameters of tissues,” Tissue Optics pp. 143–192 (2007).
  * [4] D. Papazoglou and S. Tzortzakis, “In-line holography for the characterization of ultrafast laser filamentation in transparent media,” Applied Physics Letters 93, 041120 (2008).
  * [5] Z. Wang, B. Zeng, G. Li, H. Xie, W. Chu, F. He, Y. Liao, W. Liu, H. Gao, and Y. Cheng, “Time-resolved shadowgraphs of transient plasma induced by spatiotemporally focused femtosecond laser pulses in fused silica glass,” Optics letters 40, 5726–5729 (2015).
  * [6] H. Hu, T. Liu, and H. Zhai, “Comparison of femtosecond laser ablation of aluminum in water and in air by time-resolved optical diagnosis,” Optics express 23, 628–635 (2015).
  * [7] N. Šiaulys, A. Melninkaitis, and A. Dubietis, “In situ study of two interacting femtosecond filaments in sapphire,” Optics letters 40, 2285–2288 (2015).
  * [8] D. R. Solli, C. Ropers, P. Koonath, and B. Jalali, “Optical rogue waves,” Nature 450, 1054–1057 (2007).
  * [9] P. R. Poulin and K. A. Nelson, “Irreversible organic crystalline chemistry monitored in real time,” Science 313, 1756–1760 (2006).
  * [10] Z. Li, H.-E. Tsai, X. Zhang, C.-H. Pai, Y.-Y. Chang, R. Zgadzaj, X. Wang, V. Khudik, G. Shvets, and M. C. Downer, “Single-shot visualization of evolving plasma wakefields,” in _AIP Conference Proceedings,_ vol. 1777 (AIP Publishing LLC, 2016), p. 040010.
  * [11] L. Gao, J. Liang, C. Li, and L. V. Wang, “Single-shot compressed ultrafast photography at one hundred billion frames per second,” Nature 516, 74–77 (2014).
  * [12] J. Liang and L. V. Wang, “Single-shot ultrafast optical imaging,” Optica 5, 1113–1127 (2018).
  * [13] D. Faccio and A. Velten, “A trillion frames per second: the techniques and applications of light-in-flight photography,” Reports on Progress in Physics 81, 105901 (2018).
  * [14] H. Mikami, L. Gao, and K. Goda, “Ultrafast optical imaging technology: principles and applications of emerging methods,” Nanophotonics 5, 497–509 (2016).
  * [15] K. Nakagawa, A. Iwasaki, Y. Oishi, R. Horisaki, A. Tsukamoto, A. Nakamura, K. Hirosawa, H. Liao, T. Ushida, K. Goda _et al._ , “Sequentially timed all-optical mapping photography (stamp),” Nature Photonics 8, 695–700 (2014).
  * [16] T. Suzuki, R. Hida, Y. Yamaguchi, K. Nakagawa, T. Saiki, and F. Kannari, “Single-shot 25-frame burst imaging of ultrafast phase transition of ge2sb2te5 with a sub-picosecond resolution,” Applied Physics Express 10, 092502 (2017).
  * [17] A. Ehn, J. Bood, Z. Li, E. Berrocal, M. Aldén, and E. Kristensson, “Frame: femtosecond videography for atomic and molecular dynamics,” Light: Science & Applications 6, e17045–e17045 (2017).
  * [18] J. Moon, S. Yoon, Y.-S. Lim, and W. Choi, “Single-shot imaging of microscopic dynamic scenes at 5 thz frame rates by time and spatial frequency multiplexing,” Optics express 28, 4463–4474 (2020).
  * [19] H.-Y. Huang, Z.-J. Cheng, Y. Yang, Q.-Y. Yue, and C.-S. Guo, “Single-shot ultrafast sequential holographic imaging with high temporal resolution and a large field of view,” Optics letters 44, 4885–4888 (2019).
  * [20] S. Yeola, D. Kuk, and K.-Y. Kim, “Single-shot ultrafast imaging via spatiotemporal division of femtosecond laser pulses,” J. Opt. Soc. Am. B 35, 2822–2827 (2018).
  * [21] M. Sheinman, S. Erramilli, and J. Mertz, “Trillion frames per second single shot imaging with angular encoding,” Bulletin of the American Physical Society (2021).
  * [22] C. Guo, W. Liu, X. Hua, H. Li, and S. Jia, “Fourier light-field microscopy,” Optics express 27, 25573–25594 (2019).
  * [23] Q. Zhu, Y. Cai, X. Zeng, H. Long, L. Zeng, Y. Zhu, X. Lu, and J. Li, “Single-shot framing integration photography with high spatial resolution at 5.3* 1012 frames per second by an inversed 4f system,” arXiv preprint arXiv:2110.01941 (2021).
  * [24] J. Mertz, _Introduction to optical microscopy_ (Cambridge University Press, 2019).
  * [25] Y. Yu, L. Jiang, Q. Cao, X. Shi, Q. Wang, G. Wang, and Y. Lu, “Ultrafast imaging the light-speed propagation of a focused femtosecond laser pulse in air and its ionized electron dynamics and plasma-induced pulse reshaping,” Applied Physics A 122, 205 (2016).