# Worst-Case Per-User Error Bound for Asynchronous Unsourced Multiple Access

Jyun-Sian Wu, Pin-Hsun Lin, Marcel A. Mross and Eduard A. Jorswieck  Institute
for Communications Technology, Technische Universität Braunschweig, Germany
{jyun-sian.wu, p.lin, m.mross<EMAIL_ADDRESS>
###### Abstract

This work considers an asynchronous $\textsf{K}_{\text{a}}$-active-user
unsourced multiple access channel (AUMAC) with the worst-case asynchronicity.
The transmitted messages must be decoded within $n$ channel uses, while some
codewords are not completely received due to asynchronicities. We consider a
constraint of the largest allowed delay of the transmission. The AUMAC lacks
the permutation-invariant property of the synchronous UMAC since different
permutations of the same codewords with a fixed asynchronicity are
distinguishable. Hence, the analyses require calculating all
$2^{\textsf{K}_{\text{a}}}-1$ combinations of erroneously decoded messages.
Moreover, transmitters cannot adapt the corresponding codebooks according to
asynchronicity due to a lack of information on asynchronicities. To overcome
this challenge, a uniform bound of the per-user probability of error (PUPE) is
derived by investigating the worst-case of the asynchronous patterns with the
delay constraint. Numerical results show the trade-off between the energy-per-
bit and the number of active users for different delay constraints. In
addition, although the asynchronous transmission reduces interference, the
required energy-per-bit increases as the receiver decodes with incompletely
received codewords, compared to the synchronous case.

##  I Introduction

Internet-of-things (IoT), sensor networks, and ultra-reliable low latency
massive machine-type communications have attracted attention for 6G
communications and beyond. The main challenges of the codebook designs for
these systems are: 1) short-blocklength codewords and 2) a large number of
devices that an access point has to serve. Classical information theory uses
the multiple-access channel (MAC) to analyze these systems. The classical MAC
considers individual codebooks for all devices. However, the dramatically
increasing number of devices prohibits using individual codebooks practically.
In [1], the author proposes a new system model, called unsourced multiple-
access channel (UMAC). For UMAC systems, all transmitters share an identical
codebook, and the amount of data transmitted at each transmitter is the same.

There are several aspects to investigating UMAC. Authors in [2] investigate
the first-order capacity when the numbers of users are some functions of the
blocklength, and users apply individual codebooks for identification and an
identical codebook for transmitting information. The second-order asymptotic
achievable rates of the grant-free random access system, where users access
the channel without any prior request, are analyzed in [3, 4]. However, the
achievable rates vanish if the number of transmitters increases
asymptotically. Therefore, authors in [1] investigate the energy efficiency of
synchronous UMAC with per-user error probability (PUPE) constraint. Authors in
[5] propose the T-fold ALOHA and a low-complexity coding scheme for the grant-
free Gaussian random access channel, where the coding scheme is based on the
compute-and-forward [6] scheme and coding for a binary adder channel. Authors
in [5] also analyze the energy efficiency of the T-fold ALOHA and the low-
complexity coding scheme.

Asynchronous systems are worth investigating due to the difficulty of
synchronizing a large number of devices. For asynchronous classical MAC, the
capacity is the same as the synchronous MAC [7], assuming the ratio of delay
to blocklength asymptotically vanishes. For asynchronous UMAC (AUMAC), authors
in [8, 9] utilize the T-fold ALOHA [5] and the orthogonal frequency-division
multiplexing (OFDM), transforming the time-asynchronous problem to a
frequency-shift problem. The maximum delay in [9] must be smaller than the
length of the cyclic prefix. Authors in [10] apply a sparse orthogonal
frequency-division multiple access (OFDMA) scheme and compressed sensing-based
algorithms to reliably identify arbitrarily asynchronous devices and decode
messages.

We consider the AUMAC system with a bounded delay, i.e., maximum delay
$\textsf{D}_{\text{m}}\in\mathbb{Z}^{+}\cup 0$, and
$\frac{\textsf{D}_{\text{m}}}{n}$ is a constant w.r.t. the blocklength $n$.
Transmitters transmit a fixed payload size with an identical finite-length $n$
codebook. The delays of active users are smaller than $\textsf{D}_{\text{m}}$.
In our considered model, the messages have to be decoded within $n$ channel
uses. Receivers decoding without completely receiving codewords are
investigated in broadcast channels [11, 12]. We analyze the PUPE of AUMAC with
decoding from incompletely received codewords while assuming the blocklength
is finite. To provide a more precise analysis than the typically used Berry-
Esseen theorem (BET) in finite blocklength analyses [13], we apply the
saddlepoint approximation [14]. In the synchronous UMAC, for any $1\leq
k\leq\textsf{K}_{\text{a}}$, all combinations that $k$ out of
$\textsf{K}_{\text{a}}$ messages are decoded erroneously have identical tail
probabilities due to the permutation-invariant property. However, the
permutation-invariant property is invalid due to the asynchronicity. In
particular, each $k$ out of $\textsf{K}_{\text{a}}$ combination of the
erroneously decoded messages has a different tail probability, while
$k\in[\textsf{K}_{\text{a}}]$. Therefore, the analysis requires the sum of
$2^{\textsf{K}_{\text{a}}}-1$ different tail probabilities. To overcome this
computational challenge, we derive a uniform upper bound of PUPE for our
considered AUMAC. This bound allows us to: 1) analyze the PUPE without
calculating every combination of the erroneously decoded messages and 2)
evaluate the required energy to satisfy the PUPE constraint and transmit the
payload. Analyses show that even though the AUMAC has less interference than
synchronous UMAC, the reduction of PUPE due to the increasing number of
received symbols is more significant than the increment of the PUPE due to the
interference. Numerical results compare achievable energy efficiencies for the
proposed AUMAC to synchronous UMAC, which can be considered a special case of
AUMAC with $\textsf{D}_{\text{m}}=0$. Numerical results show that compared to
synchronous UMAC, transmitters in AUMAC require more energy to reliably
transmit messages with a constant $\frac{\textsf{D}_{\text{m}}}{n}$.

Notation: We will denote $f^{(i)}(t)$ as the $i$-th derivative of $f(x)$ at
the point $x=t$ and $f_{1,y}^{(i)}(x,t)$ as the $i$-th partial derivative of
$f_{1}(x,y)$ w.r.t. $y$ at the point $y=t$. We use the indicator function
$\mathbbm{1}(\cdot)$, the natural logarithm $\log(\cdot)$, and the Landau
symbol $O(\cdot)$. The binomial coefficient of $n$ out of $k$ is represented
by $\binom{n}{k}$. The number of permutations of $k$ is denoted as $k!$. We
define $j=\sqrt{-1}$. We denote $[k]=\\{1,2,...,k\\}$ and
$\mathcal{F}\setminus\mathcal{T}=\\{x:x\in\mathcal{F},x\not\in\mathcal{T}\\}$,
where $\mathcal{F}$ and $\mathcal{T}$ are two sets. We also denote
$\mathbb{Z}^{+}_{0}=\mathbb{Z}^{+}\cup 0$. For any set
$\mathcal{F}=\\{F_{1},F_{2},...,F_{|\mathcal{F}|}\\}$, we denote
$\\{X_{m}\\}_{m\in\mathcal{F}}=\\{X_{F_{1}},X_{F_{2}},...,X_{F_{|\mathcal{F}|}}\\}$.

##  II System model and preliminaries

We consider an AUMAC, which has additive white Gaussian noise (AWGN), one
receiver, and multiple transmitters, where the number of active transmitters
is denoted by a positive integer $\textsf{K}_{\text{a}}$. All transmitters
utilize the same codebook with the same maximal power constraint,
$\mathsf{P}^{\prime}$, to transmit the same (and fixed) size of payloads,
i.e., ${\log\textsf{M}}$ nats, to the receiver. The codewords are independent
and identically distributed (i.i.d.) generated from a Gaussian distribution
with mean zero and variance $\mathsf{P}$, where
$\mathsf{P}<\mathsf{P}^{\prime}$ due to the power backoff. The power backoff
reduces the probability that the maximal power constraint violations occur.

###### Definition 1

We define the asynchronicity in terms of the vector of time shifts (delay) as

$\displaystyle D^{\textsf{K}_{\text{a}}}\\!:=[d_{1},\ d_{2},\ ...,\
d_{\textsf{K}_{\text{a}}}]\in\\{\mathbb{Z}^{+}_{0}\\}^{\textsf{K}_{\text{a}}},\vspace{-2.5mm}$

where $0=d_{1},\;d_{i}\leq\textsf{D}_{\text{m}}$ and $d_{i}\leq
d_{\ell},\;\forall\ell>i$ for all $i\in[\textsf{K}_{\text{a}}]$. The $i$-th
entry, $d_{i}$, represents the delay of the $i$-th received codeword relative
to the first received codeword, and $\textsf{D}_{\text{m}}$ denotes the delay
constraint. We define $\alpha:=\frac{\textsf{D}_{\text{m}}}{n}\in[0,1)$, which
is constant w.r.t. the blocklength $n$, and $\bar{\alpha}=1-\alpha$.

We assume that the receiver has perfect knowledge of the asynchronicity [15]
and jointly detects the transmitted messages. Asynchronous communication
systems may result from asynchronous clocks between transmitters and
receivers, different idle times among transmitters, or channel delays.

###### Remark 1

We consider that every transmitter transmits with the same codebook, and the
receiver is not interested in identifying the senders of the received
codewords. Therefore, $d_{i}$ indicates the delay of the $i$-th received
codeword but does not indicate the identification of the transmitter.

In the asynchronous model, the number of transmitted codewords symbols of each
channel use can be different. For a given delay $D^{\textsf{K}_{\text{a}}}\\!$
and the set of erroneously decoded messages
$\mathcal{S}\subseteq[\textsf{K}_{\text{a}}]$, we define a vector

$\displaystyle
a^{n}(\mathcal{S},\,D^{\textsf{K}_{\text{a}}}\\!)\\!:=\\![a_{1}(\mathcal{S},\,D^{\textsf{K}_{\text{a}}}\\!),a_{2}(\mathcal{S},\,D^{\textsf{K}_{\text{a}}}\\!),\\!...,a_{n}(\mathcal{S},\,D^{\textsf{K}_{\text{a}}}\\!)],$
(1)

where $a_{i}(\mathcal{S},D^{\textsf{K}_{\text{a}}}\\!)\leq
a_{\ell}(\mathcal{S},D^{\textsf{K}_{\text{a}}}\\!)$, $\forall\ell>i$,
$i\in[n]$ and
$a_{i}(\mathcal{S},D^{\textsf{K}_{\text{a}}}\\!)\in\mathbb{Z}^{+}_{0}$,
$\forall i\in[n]$. For a given $D^{\textsf{K}_{\text{a}}}\\!$ and a given
$i\in[n]$, the $i$-th entry of
$a^{n}(\mathcal{S},D^{\textsf{K}_{\text{a}}}\\!)$, i.e.,
$a_{i}(\mathcal{S},D^{\textsf{K}_{\text{a}}}\\!)$, indicates the number of
simultaneously received symbols, which belong to $\mathcal{S}$, at the $i$-th
channel use. To simplify notations, we use $a^{n}:=[a_{1},a_{2},...,a_{n}]$ to
represent $a^{n}(\mathcal{S},D^{\textsf{K}_{\text{a}}}\\!)$. For example,
considering a $\textsf{K}_{\text{a}}$-active-user AUMAC with
$D^{\textsf{K}_{\text{a}}}\\!=[0,1,3,5,...,5]$ in Fig. 1, for the set
$\mathcal{S}\\!=\\!\\{1,2\\}$,
$a^{n}\\!=\\![a_{1}\\!=\\!1,a_{[n]\setminus[1]}\\!=\\!2]$; for the set
$\mathcal{S}\\!=\\!\\{2,3,4\\}$,
$a^{n}\\!=\\![a_{1}\\!=\\!0,a_{2}\\!=\\!1,a_{3}\\!=\\!1,a_{4}\\!=\\!2,a_{5}\\!=\\!2,a_{[n]\setminus[5]}\\!=\\!3]$.
Note that for a given $\mathcal{S}$ and $D^{\textsf{K}_{\text{a}}}\\!$,
$a_{[n]\setminus[\alpha n]}=|\mathcal{S}|$.

Figure 1: A $\textsf{K}_{\text{a}}$-active-user AUMAC with
$D^{\textsf{K}_{\text{a}}}\\!=[0,1,3,5,...,5]$.

For any $\ell\in[n]$, we define a shift function
$\tau_{d_{i}}(X^{n}_{i},\ell):=X_{i,\ell-d_{i}}$, where $X_{i,\ell-d_{i}}$ is
the $(\ell-d_{i})$-th element of $X^{n}_{i}$, and if $\ell-d_{i}\not\in[n]$,
$X_{i,\ell-d_{i}}=0,\ \forall i\in[\textsf{K}_{\text{a}}]$. The received
symbol at the receiver at time $\ell\in[n]$ is

$\displaystyle{Y}_{\ell}=\sum_{i=1}^{\textsf{K}_{\text{a}}}\tau_{d_{i}}\left(X^{n}_{i},\ell\right)+Z_{\ell},\vspace{-2mm}$
(2)

where the channel input
$X^{n}_{i}\\!\in\\!\chi^{n}\\!\subset\\!\mathbb{R}^{n}$, where
$\chi^{n}\\!:=\\!\\{x^{n}:\\!x^{n}\\!\in\mathbb{R}^{n},\;\|x^{n}\|^{2}\\!\leq\\!n\mathsf{P}^{\prime}\\}$
is the channel input satisfying the maximal power constraint and
$Z_{\ell}\sim\mathcal{N}(0,1)$ is an i.i.d. AWGN, $\forall\ell\in[n]$.

###### Definition 2

An
$(n,{\textsf{M}},\epsilon,\textsf{K}_{\text{a}},\alpha,D^{\textsf{K}_{\text{a}}}\\!)-$code,
$\mathcal{C}_{1}$, for an AUMAC described by
$P_{Y|X_{[\textsf{K}_{\text{a}}]}}$, consists of

  * •

one message set $\mathcal{M}=\\{1,2,...,{\textsf{M}}\\}$,

  * •

one encoder $f:\mathcal{M}\to\chi^{n}$,

  * •

one decoder $g:\mathbb{R}^{n}\to\binom{[\mathcal{M}]}{\textsf{K}_{\text{a}}},$
where $\binom{[\mathcal{M}]}{\textsf{K}_{\text{a}}}$ is a set containing
$\textsf{K}_{\text{a}}$ distinct elements from the set $\mathcal{M}$,

and the delay $D^{\textsf{K}_{\text{a}}}\\!$ fulfills the delay constraint
$\alpha n$ in Def. 1, the PUPE satisfies

$\displaystyle{\mbox{P}_{\text{PUPE}|D^{\textsf{K}_{\text{a}}}\\!}}:=\frac{1}{\textsf{K}_{\text{a}}}\sum_{i=1}^{\textsf{K}_{\text{a}}}\mbox{Pr}(\tilde{\mathcal{E}}_{i}|D^{\textsf{K}_{\text{a}}}\\!)\leq\epsilon,\vspace{-2.5mm}$
(3)

where $\tilde{\mathcal{E}}_{i}:=\\{\;\cup_{\ell\neq
i}\\{M_{i}=M_{\ell}\\}\cup\\{M_{i}\not\in
g({Y}^{n})\\}\cup\\{\|f(M_{i})\|^{2}>n\mathsf{P}^{\prime}\\}\\}$,
$i\in[\textsf{K}_{\text{a}}]$, and $M_{i}\sim\mbox{Unif}(\mathcal{M})$ is the
$i$-th transmitted message.

##  III Main Results

There are several achievable schemes in the synchronous multiple-access
models. For synchronous UMAC models, shell codes achieve better second-order
asymptotic rates than the i.i.d. Gaussian codes [4, 3]. However, in our
considered asynchronous model, the receiver decodes the messages solely based
on the first $n$ received symbols, i.e., some codewords are incompletely
received when the receiver starts to decode. The decoding based on
incompletely received codewords does not satisfy the definition of shell
codes. This fact prevents us from analyzing the models by the uniform
distribution on a power shell. Therefore, we consider an achievable scheme
that all transmitters share the same codebook i.i.d. generated by a Gaussian
distribution $P_{X}$. The receiver performs maximum information density
decoding with knowledge of $D^{\textsf{K}_{\text{a}}}\\!$, which is defined by

$\displaystyle
g(Y^{n})=\underset{X^{n}_{[\textsf{K}_{\text{a}}]}\in\mathcal{C}_{1}}{\arg\max}\
\sum_{\ell=1}^{n}i\left(\big{\\{}\tau_{d_{m}}(X^{n}_{m},\ell)\big{\\}}_{m\in[\textsf{K}_{\text{a}}]};Y_{\ell}\right),\vspace{-1mm}$
(4)

where

$\displaystyle i\Big{(}\big{\\{}\tau_{d_{m}}(X^{n}_{m},\ell)$
$\displaystyle\big{\\}}_{m\in[\textsf{K}_{\text{a}}]};Y_{\ell}\Big{)}:=$
$\displaystyle\log\left(\frac{dP_{Y|X_{[\textsf{K}_{\text{a}}]}}\left(Y_{\ell}|\big{\\{}\tau_{d_{m}}(X^{n}_{m},\ell)\big{\\}}_{m\in[\textsf{K}_{\text{a}}]}\right)}{dP_{Y}(Y_{\ell})}\right).$

The corresponding finite-blocklength (FBL) analysis results are summarized in
the following.

###### Theorem 1

Fix $0<\mathsf{P}<\mathsf{P}^{\prime}$. There exists an
$(n,{\textsf{M}},\epsilon,\textsf{K}_{\text{a}},\alpha,D^{\textsf{K}_{\text{a}}}\\!)-$code
for an AUMAC such that the PUPE can be upper bounded by the following:

$\displaystyle\sum\limits_{\mathcal{S}\subseteq[\textsf{K}_{\text{a}}]}\\!\\!\frac{|\mathcal{S}|g_{1}(\\!a^{n}\\!,t_{0}(a^{n}))}{\textsf{K}_{\text{a}}\sqrt{2\pi}}\\!\Big{[}g_{2}(a^{n}\\!,\\!t_{0}(a^{n}))\\!\\!+\\!\xi(a^{n},\\!t_{0}(a^{n}))\\!\Big{]}\\!\\!+\\!p_{0}\\!\\!\leq\\!\epsilon,\vspace{-1mm}$
(5)

if there exists a $t_{0}(a^{n})\in(0,1)$ such that
$E^{(1)}_{t}(a^{n},t_{0}(a^{n}))=|\mathcal{S}|{\log\textsf{M}}$, where

$\displaystyle g_{1}(a^{n},t):=$
$\displaystyle\exp(t|\mathcal{S}|{\log\textsf{M}}-E(a^{n},t)),$ (6)
$\displaystyle g_{2}(a^{n},t):=$
$\displaystyle\left(t(1-t)\sqrt{-E_{t}^{(2)}(a^{n},t)}\right)^{-1},$ (7)
$\displaystyle E(a^{n},t):=$
$\displaystyle\frac{1}{2}\sum\limits_{i=1}^{n}\left(\\!t\log(1\\!+\\!a_{i}\mathsf{P})\\!+\\!\log\left(1\\!-\\!\frac{a_{i}\mathsf{P}t^{2}}{1\\!+\\!a_{i}\mathsf{P}}\right)\right),$
(8) $\displaystyle\xi(a^{n},t):=$ $\displaystyle\frac{1}{2\pi
j}\int^{t+j\infty}_{t-j\infty}\exp\left(-\frac{E^{(2)}_{t}(a^{n},t)}{2}(\rho-t)^{2}\right)$
$\displaystyle\hskip
42.67912pt\cdot\frac{1}{\rho(1-\rho)}\sum_{m=1}^{\infty}\frac{\bar{\xi}(a^{n},t,\rho)^{m}}{m!}d\rho,\\!\\!$
(9) $\displaystyle\bar{\xi}(a^{n},t,\rho):=$
$\displaystyle-\sum_{i=3}^{\infty}E^{(i)}_{t}(a^{n},t)\frac{(\rho-t)^{i}}{i!},\vspace{-1mm}$
(10)

and
$p_{0}:=\frac{\textsf{K}_{\text{a}}(\textsf{K}_{\text{a}}-1)}{2{\textsf{M}}}+\sum\limits_{i=1}^{\textsf{K}_{\text{a}}}\mbox{Pr}(\|X^{n}_{i}\|^{2}>n\mathsf{P}^{\prime}).$

In our PUPE analysis, two main tools are the Taylor expansion and the inverse
Laplace transform. The former tool is used to expand the exponent of
$g_{1}(a^{n},t)$ at $t=t_{0}(a^{n})$ and the latter one is used to derive the
probability density function (PDF) from the cumulant-generating function
(CGF), where CGF is the logarithm of the moment-generating function (MGF). The
sum of higher order terms of the Taylor expansion at $t=t_{0}(a^{n})$ is
represented by $\bar{\xi}(a^{n},t_{0}(a^{n}),t)$. The proof is relegated to
Appendix A.

Theorem 1 can evaluate the PUPE for any delay $D^{\textsf{K}_{\text{a}}}\\!$
satisfying the delay constraint $\alpha n$. However, evaluating (5) requires
calculating all $\mathcal{S}\subseteq[\textsf{K}_{\text{a}}]$, which is
infeasible if $\textsf{K}_{\text{a}}$ is sufficiently large. Additionally,
even though we can design different codebooks for different
$D^{\textsf{K}_{\text{a}}}\\!$ satisfying the PUPE constraints, the
transmitters cannot select the codebook corresponding to a particular
$D^{\textsf{K}_{\text{a}}}\\!$ since they have no information on delays.
Therefore, a codebook that satisfies the PUPE constraint regardless of
$D^{\textsf{K}_{\text{a}}}\\!$ is required. In the following, we derive a
uniform upper bound of the PUPE for all $D^{\textsf{K}_{\text{a}}}\\!$’s
satisfying delay constraint $\alpha n$.

###### Definition 3

An $(n,{\textsf{M}},\epsilon,\textsf{K}_{\text{a}},\alpha)-$code,
$\mathcal{C}_{2}$, for an AUMAC described by
$P_{Y|X_{[\textsf{K}_{\text{a}}]}}$ consists of one message set $\mathcal{M}$,
one encoder $f$, and one decoder $g$ defined by

$\displaystyle
g(Y^{n})=\underset{X^{n}_{[\textsf{K}_{\text{a}}]}\in\mathcal{C}_{2}}{\arg\max}\
\sum_{\ell=1}^{n}i\left(\big{\\{}\tau_{d_{m}}(X^{n}_{m},\ell)\big{\\}}_{m\in[\textsf{K}_{\text{a}}]};Y_{\ell}\right),\vspace{-1mm}$
(11)

such that for the power constraint $\mathsf{P}^{\prime}$ and any
$D^{\textsf{K}_{\text{a}}}\\!$ satisfying the maximum delay constraint, the
PUPE satisfies

$\displaystyle{\mbox{P}_{\text{PUPE}}}:=\max\limits_{D^{\textsf{K}_{\text{a}}}\\!:\
d_{\textsf{K}_{\text{a}}}\leq\alpha
n}\sum_{i=1}^{\textsf{K}_{\text{a}}}\frac{1}{\textsf{K}_{\text{a}}}\mbox{Pr}(\tilde{\mathcal{E}}_{i}|D^{\textsf{K}_{\text{a}}}\\!)\leq\epsilon,\vspace{-1mm}$
(12)

where $\tilde{\mathcal{E}}_{i}$ is defined in Definition 2.

Based on the PUPE defined in Def. 3, we find the $a^{*n}_{\iota}$ that leads
to the uniform upper bound of the PUPE, where $\iota\in\\{0,1\\}$.

###### Theorem 2

Fix $0<\mathsf{P}<\mathsf{P}^{\prime}$. There exists an
$(n,{\textsf{M}},\epsilon,\textsf{K}_{\text{a}},\alpha)-$code for AUMAC, such
that the PUPE can be upper bounded by the following:

$\displaystyle\frac{1}{\textsf{K}_{\text{a}}\sqrt{2\pi}}\sum_{|\mathcal{S}|=1}^{\textsf{K}_{\text{a}}}\left(\binom{\textsf{K}_{\text{a}}-1}{|\mathcal{S}|}\frac{|\mathcal{S}|g_{1}(a_{0}^{n*},t_{0}(a_{0}^{n*}))}{T_{0}^{*}\sqrt{-E^{(2)}_{t}(a_{0}^{n*},\underline{t}_{0})}}\right.$
$\displaystyle\\!\\!+\\!\\!\left.\binom{\textsf{K}_{\text{a}}-1}{|\mathcal{S}|-1}\frac{|\mathcal{S}|g_{1}(a_{1}^{n*},t_{0}(a_{1}^{n*}))}{T_{1}^{*}\sqrt{-E^{(2)}_{t}(a_{1}^{n*},\underline{t}_{1})}}\\!\right)\\!\\!+\\!p_{0}\\!+\\!O\left(\frac{\exp(-n)}{\sqrt{n}}\right)\\!\leq\epsilon,\vspace{-1mm}$
(13)

if $t_{0}(a^{n*}_{\iota})\in\mathcal{A}\cap\mathcal{B},\
\bar{t}_{\iota}\in\mathcal{A}\cap\bar{\mathcal{B}}$, $\
\underline{t}_{\iota}\in\mathcal{A}\cap\underline{\mathcal{B}}$, and
$\underline{t}_{\iota}\leq t_{0}(a^{n})\leq\bar{t}_{\iota}$, where
$a_{\iota}^{n*}=[\iota^{\alpha n},|\mathcal{S}|^{n-\alpha n}]$,
$T^{*}_{\iota}:=\min\\{\underline{t}_{\iota}-\underline{t}_{\iota}^{2},\
\bar{t}_{\iota}-\bar{t}_{\iota}^{2}\\}$,
$\iota:=\mathbbm{1}(1\in\mathcal{S})$,
$\mathcal{A}\\!:=\\!\left\\{t\\!:\\!t\in(0,1)\right\\}$,

$\displaystyle\mathcal{B}\\!:=\\!\left\\{t\\!:\\!E^{(1)}_{t}(a^{n*}_{\iota},t)=|\mathcal{S}|{\log\textsf{M}}\right\\},$
(14)
$\displaystyle\bar{\mathcal{B}}\\!:=\\!\left\\{\\!t\\!:\\!\\!\sum_{i=1}^{n}\\!\frac{a^{*}_{\iota,i}\mathsf{P}t}{\\!1\\!\\!+\\!a^{*}_{\iota,i}\mathsf{P}(1\\!-\\!t^{2})\\!}\\!\\!=\\!\\!\frac{n}{2}\\!\log(1\\!\\!+\\!|\mathcal{S}|\mathsf{P})\\!\\!-\\!|\mathcal{S}|{\log\textsf{M}}\\!\right\\},\\!\\!\\!\\!\\!$
(15)
$\displaystyle\underline{\mathcal{B}}\\!:=\\!\left\\{t\\!:\\!\frac{|\mathcal{S}|n\mathsf{P}t}{\\!1\\!\\!+\\!|\mathcal{S}|\mathsf{P}(1\\!-\\!t^{2})\\!}\\!=\\!\\!\sum_{i=1}^{n}\\!\frac{\log(1\\!\\!+\\!a^{*}_{\iota,i}\mathsf{P})}{2}\\!-\\!|\mathcal{S}|{\log\textsf{M}}\\!\right\\},\\!\\!\\!\\!$
(16) $\displaystyle\text{and }a^{*}_{\iota,i}\text{ is the }i\text{-th element
of }a^{n*}_{\iota}.$

The proof is relegated to Appendix B.

The benefit of having a uniform PUPE upper bound is that it allows us to
analyze the performance of an
$(n,{\textsf{M}},\epsilon,\textsf{K}_{\text{a}},\alpha)-$code without
calculating all tail probabilities of the corresponding possible $\mathcal{S}$
but scaling the uniform PUPE upper bound by a binomial coefficient.

The term $g_{1}(a^{n},t_{0}(a^{n}))\xi(a^{n},\\!t_{0}(a^{n}))$ in Theorem 1 is
expressed by $O\left(\frac{\exp(-n)}{\sqrt{n}}\right)$ in Theorem 2 since
$\bar{\xi}(a^{n},t,\rho)$ behaves as $O(n^{-\frac{1}{2}})$ [14]. We refer to
both the approximations obtained by ignoring
$g_{1}(a^{n},t_{0}(a^{n}))\xi(a^{n},t_{0}(a^{n}))$ and
$O\left(\frac{\exp(-n)}{\sqrt{n}}\right)$ in (5) and (2), respectively, as
saddlepoint approximations. The saddlepoint approximation has an exponentially
decreasing approximation error w.r.t. $n$, allowing us to obtain sufficiently
precise approximations of the FBL PUPE compared to the BET.

###### Remark 2

The upper bound in Theorem 1 decreases as $a_{i}$ increases for any
$D^{\textsf{K}_{\text{a}}}\\!$ and $\mathcal{S}$ because
$\frac{\partial}{\partial a_{i}}g_{1}(a^{n},t)g_{2}(a^{n},t)\leq 0$ for all
$t\in(0,1)$ and $i\in[\alpha n]$. In fact, having more overlap in the
transmission leads to more interference. A larger number of overlapping
symbols has one positive and one negative effect on the receiver: it leads to
more received energy but, meanwhile, more interference. By our analysis, we
found that the positive effect is dominant.

##  IV Numerical Results

Based on the results of Theorem 1 and Theorem 2 without considering
$g_{1}(a^{n},t_{0}(a^{n}))\xi(a^{n},t_{0}(a^{n}))$ and
$O\left(\frac{\exp(-n)}{\sqrt{n}}\right)$, we numerically evaluate the energy-
per-bit versus the number of active users. We define the energy-per-bit as
$\frac{\text{E}_{\text{b}}}{\text{N}_{\text{0}}}:=\frac{n\mathsf{P}^{\prime}}{{\log\textsf{M}}}$.
The PUPE upper bounds from Theorem 2 are compared to two UMAC schemes under
different scenarios with the following parameters: ${\log\textsf{M}}=100$,
$n=4000$, $\epsilon=10^{-3}$ and $\textsf{K}_{\text{a}}\in[50,160]$. In Fig.
2, the purple dash curve of [1] is evaluated by numerically optimizing
$\mathsf{P}<\mathsf{P}^{\prime}$. We also evaluate the
$\frac{\text{E}_{\text{b}}}{\text{N}_{\text{0}}}$ for a $16$-fold ALOHA [5]
with Theorem 1. The $16$-fold ALOHA splits the transmission into $V$ subblocks
such that the collision probability is less than $0.9\epsilon$ [5]. Each
subblock has blocklength $\tilde{n}=n/V$ and the delay constraint for each
subblock is $\alpha\tilde{n}$. We assume that the messages have to be decoded
within $\tilde{n}$ channel uses. The black dot-dash curve shows that the
required $\frac{\text{E}_{\text{b}}}{\text{N}_{\text{0}}}$ for $16$-fold ALOHA
from Theorem 1. The synchronous UMAC can be considered to be a special case of
AUMAC with $\alpha=0$. Additionally, analyzing synchronous UMAC does not
require calculating the number of permutations of erroneously decoded
messages, i.e., $|\mathcal{S}|!$, so we modify (6) in Theorem 1 as follows:

$\displaystyle g_{1}(a^{n},t):=$
$\displaystyle\exp\left(t\log\left(\binom{{\textsf{M}}-\textsf{K}_{\text{a}}}{|\mathcal{S}|}\right)-E(a^{n},t)\right).$
(17)

The yellow solid curve shows the required
$\frac{\text{E}_{\text{b}}}{\text{N}_{\text{0}}}$ for the synchronous UMAC. It
is computed by numerically optimizing $\mathsf{P}<\mathsf{P}^{\prime}$ in
Theorem 1 with
$D^{\textsf{K}_{\text{a}}}\\!=\\{d_{1}=d_{2}=...=d_{\textsf{K}_{\text{a}}}=0\\}$
and adapting (17).

We numerically optimize $\mathsf{P}$ in Theorem 2 with $\alpha=0.2$ and
$\alpha=0.4$ and compare the $\frac{\text{E}_{\text{b}}}{\text{N}_{\text{0}}}$
of AUMAC and that of synchronous UMAC. Numerical results show that the AUMAC
that has larger $\alpha$ causes the transmitters to consume more energy to
transmit in the worst case of delay. Observing the curves of $\alpha=0.2$,
$\alpha=0.4$ and $\alpha=0$ (synchronous), we can conclude that for the AUMAC
with larger $\alpha$, which means fewer interference for the first $\alpha n$
channel uses, the PUPE increases. It is because the receiver decodes the
messages based on fewer transmitted codewords symbols, which is equivalently
based on less received energy. This effect is illuminated in Remark 2. Thus,
codebooks of our considered model require more energy to achieve the same PUPE
constraint.

Figure 2: $\frac{\text{E}_{\text{b}}}{\text{N}_{\text{0}}}$ of AUMAC compared
to synchronous UMAC for different numbers of active users.

##  V Conclusions

In this work, we analyze the FBL performance of the asynchronous UMAC system
with bounded and non-vanishing delay constraints $\alpha n$. The derivations
based on the saddlepoint approximation provide FBL performance bounds. We also
investigate a uniform upper bound of the PUPE, which highly simplifies the
analysis to multiply the uniform upper bound with the corresponding binomial
coefficient instead of calculating tail probabilities of all error events. The
numerical results show the trade-off between
$\frac{\text{E}_{\text{b}}}{\text{N}_{\text{0}}}$ and delay constraint $\alpha
n$. Although asynchronous transmissions have less interference, reducing the
error probability of the first few codewords, it increases PUPE as the
receiver decodes shorter codewords, which is analytically shown in Theorem 2
and is numerically shown in our numerical results. Compared to the synchronous
case, the achievable energy-per-bit
$\frac{\text{E}_{\text{b}}}{\text{N}_{\text{0}}}$ for the asynchronous case
shows that the required $\frac{\text{E}_{\text{b}}}{\text{N}_{\text{0}}}$
increases as the receiver decodes shorter codewords, even though interference
reduces.

## References

  * [1] Y. Polyanskiy, “A perspective on massive random-access,” in _2017 IEEE International Symposium on Information Theory (ISIT)_ , 2017, pp. 2523–2527.
  * [2] X. Chen, T.-Y. Chen, and D. Guo, “Capacity of Gaussian many-access channels,” _IEEE Transactions on Information Theory_ , vol. 63, no. 6, pp. 3516–3539, 2017.
  * [3] R. C. Yavas, V. Kostina, and M. Effros, “Random access channel coding in the finite blocklength regime,” _IEEE Transactions on Information Theory_ , vol. 67, no. 4, pp. 2115–2140, 2021.
  * [4] ——, “Gaussian multiple and random access channels: Finite-blocklength analysis,” _IEEE Transactions on Information Theory_ , vol. 67, no. 11, pp. 6983–7009, 2021.
  * [5] O. Ordentlich and Y. Polyanskiy, “Low complexity schemes for the random access Gaussian channel,” in _2017 IEEE International Symposium on Information Theory (ISIT)_ , 2017, pp. 2528–2532.
  * [6] B. Nazer and M. Gastpar, “Compute-and-forward: Harnessing interference with structured codes,” in _2008 IEEE International Symposium on Information Theory_ , 2008, pp. 772–776.
  * [7] T. Cover, R. McEliece, and E. Posner, “Asynchronous multiple-access channel capacity,” _IEEE Transactions on Information Theory_ , vol. 27, no. 4, pp. 409–413, 1981.
  * [8] K. Andreev, S. S. Kowshik, A. Frolov, and Y. Polyanskiy, “Low complexity energy efficient random access scheme for the asynchronous fading MAC,” in _2019 IEEE 90th Vehicular Technology Conference (VTC2019-Fall)_ , 2019, pp. 1–5.
  * [9] S. S. Kowshik, K. Andreev, A. Frolov, and Y. Polyanskiy, “Short-packet low-power coded access for massive MAC,” in _2019 53rd Asilomar Conference on Signals, Systems, and Computers_ , 2019, pp. 827–832.
  * [10] X. Chen, L. Liu, D. Guo, and G. W. Wornell, “Asynchronous massive access and neighbor discovery using OFDMA,” _IEEE Transactions on Information Theory_ , vol. 69, no. 4, pp. 2364–2384, 2023.
  * [11] P.-H. Lin, S.-C. Lin, P.-W. Chen, M. A. Mross, and E. A. Jorswieck, “Second order rate regions of Gaussian broadcast channels under heterogeneous blocklength constraints,” _IEEE Transactions on Communications_ , pp. 1–1, 2023.
  * [12] M. Mross, P. Lin, and E. A. Jorswieck, “Second-order performance of early decoding with shell codes in Gaussian broadcast channels,” in _IEEE International Symposium on Information Theory, ISIT 2023, Taipei, Taiwan, June 25-30, 2023_. IEEE, 2023, pp. 2123–2128. [Online]. Available: https://doi.org/10.1109/ISIT54713.2023.10206483
  * [13] A. Martinez and A. G. i Fàbregas, “Saddlepoint approximation of random-coding bounds,” in _2011 Information Theory and Applications Workshop_ , 2011, pp. 1–6.
  * [14] J. Jensen, _Saddlepoint approximations_. Clarendon Press Oxford, 1995, oxford Science Publications.
  * [15] I. Csiszár, L. Farkas, and T. Kói, “Error exponents for asynchronous multiple access channels, controlled asynchronism may outperform synchronism,” _IEEE Transactions on Information Theory_ , vol. 67, no. 12, pp. 7684–7707, 2021.
  * [16] Y. Polyanskiy, H. V. Poor, and S. Verdu, “Channel coding rate in the finite blocklength regime,” _IEEE Trans. Inf. Theory_ , vol. 56, no. 5, pp. 2307–2359, 2010.
  * [17] Y. Polyanskiy and Y. Wu, _Information Theory: From Coding to Learning_. Cambridge University Press, 2023.
  * [18] G. D. Finn and D. Mugglestone, “Tables of the line broadening function $h(a,v)$,” _Monthly Notices Roy. Astronom. Soc._ , vol. 129, p. 221–235, 1965.

##  Appendix A Proof of Theorem 1

Theorem 1 is derived by the maximal information density decoder with the
random coding union (RCU) bound [16] to express the per-user probability of
error (PUPE) as a sum of tail probabilities. In the following, we first show
the expressions of PUPE regarding tail probabilities. Then, we apply the
Taylor expansion and the inverse Laplace transform to derive the tail
probabilities.

We define
$\tilde{\mathcal{E}}_{\ell}:=\\{\;\cup_{i\neq\ell}\\{M_{\ell}=M_{i}\\}\cup\\{M_{\ell}\not\in
g({Y}^{n})\\}\cup\\{\|f(M_{\ell})\|^{2}>n\mathsf{P}^{\prime}\\}\\}$, which
represents the $\ell$-th user’s error event of the PUPE,
${\mathcal{E}}_{\ell}:=\\{\\{M_{\ell}\neq M_{i},\;\forall
i\neq\ell\\}\cap\\{\|f(M_{i})\|^{2}\leq n\mathsf{P}^{\prime},\;\forall
i\in[\textsf{K}_{\text{a}}]\\}\\}$, $\ell\in[\textsf{K}_{\text{a}}]$, which
represents the event that the other transmitted messages are distinct to the
$i$-th transmitted message and transmitted codewords fulfill the power
constraint, and
$p_{0}:=\frac{\textsf{K}_{\text{a}}(\textsf{K}_{\text{a}}-1)}{2{\textsf{M}}}+\sum\limits_{i=1}^{\textsf{K}_{\text{a}}}\mbox{Pr}(\|X^{n}_{i}\|^{2}>n\mathsf{P}^{\prime})$
is the upper bound of the probability that collisions or power constraint
violations occur.

The PUPE of an
$(n,{\textsf{M}},\epsilon,\textsf{K}_{\text{a}},\alpha,D^{\textsf{K}_{\text{a}}}\\!)$-code
can be upper bounded by the union bound as follows:

$\displaystyle{\mbox{P}_{\text{PUPE}|D^{\textsf{K}_{\text{a}}}\\!}}:$
$\displaystyle=\sum_{\ell=1}^{\textsf{K}_{\text{a}}}\frac{1}{\textsf{K}_{\text{a}}}\mbox{Pr}\Big{(}\tilde{\mathcal{E}}_{\ell}|D^{\textsf{K}_{\text{a}}}\\!\Big{)}$
(18) $\displaystyle\leq
p_{0}+\sum_{\ell=1}^{\textsf{K}_{\text{a}}}\frac{1}{\textsf{K}_{\text{a}}}\mbox{Pr}\Big{(}M_{\ell}\not\in
g(Y^{n})|D^{\textsf{K}_{\text{a}}}\\!,{\mathcal{E}}_{\ell}\Big{)}.$ (19)

To simplify the notation, we omit the condition $D^{\textsf{K}_{\text{a}}}\\!$
in the following derivation. For any subset
$\mathcal{S}\subseteq[\textsf{K}_{\text{a}}]$, we define

$\displaystyle\tilde{\gamma}(\bar{X}^{n}_{\mathcal{S}},X^{n}_{[\textsf{K}_{\text{a}}]\setminus\mathcal{S}})$
$\displaystyle\hskip
22.76228pt\\!\\!:=\\!\\!\sum_{\ell=1}^{n}i(\\{\tau_{d_{m}}(\bar{X}^{n}_{m},\ell)\\}_{m\in\mathcal{S}},\\!\\{\tau_{d_{m}}(X^{n}_{m},\ell)\\}_{m\in[\textsf{K}_{\text{a}}]\setminus\mathcal{S}};Y_{\ell})$
and
$\displaystyle\gamma(\bar{X}^{n}_{\mathcal{S}})\\!\\!:=\\!\\!\sum_{\ell=1}^{n}\\!i(\\{\tau_{d_{m}}(\bar{X}^{n}_{m},\ell)\\}_{m\in\mathcal{S}};\\!Y_{\ell}|\\{\tau_{d_{m}}(X^{n}_{m},\ell)\\}_{m\in[\textsf{K}_{\text{a}}]\setminus\mathcal{S}}).$

We define a set

$\displaystyle\Sigma(\ell):=\\{\mathcal{S}:\mathcal{S}\subseteq[\textsf{K}_{\text{a}}],\ell\in\mathcal{S}\\},$
(20)

which contains all possible subsets $\mathcal{S}$ of the error event
$\\{M_{\ell}\not\in g(Y^{n})\\}$. Substitute the definition of the maximal
information density decoder into $\mbox{Pr}(M_{\ell}\not\in
g(Y^{n})|{\mathcal{E}}_{\ell})$, we have

$\displaystyle\mbox{Pr}(M_{\ell}\not\in g(Y^{n})|{\mathcal{E}}_{\ell})$
$\displaystyle=\mbox{Pr}\left(\bigcup\limits_{\begin{subarray}{c}{\mathcal{S}\in\Sigma(\ell),}\\\
{\bar{X}^{n}_{\mathcal{S}}\neq
X^{n}_{\mathcal{S}}}\end{subarray}}\left\\{\tilde{\gamma}(\bar{X}^{n}_{\mathcal{S}},X^{n}_{[\textsf{K}_{\text{a}}]\setminus\mathcal{S}})>\tilde{\gamma}(X^{n}_{[\textsf{K}_{\text{a}}]})\right\\}\bigg{|}{\mathcal{E}}_{\ell}\right)$
(21)
$\displaystyle=\mbox{Pr}\left(\bigcup_{\begin{subarray}{c}{\mathcal{S}\in\Sigma(\ell),}\\\
{\bar{X}^{n}_{\mathcal{S}}\neq
X^{n}_{\mathcal{S}}}\end{subarray}}\left\\{\gamma(\bar{X}^{n}_{\mathcal{S}})>\gamma(X^{n}_{\mathcal{S}})\right\\}\bigg{|}{\mathcal{E}}_{\ell}\right)$
(22)
$\displaystyle=\mathds{E}\left[\mbox{Pr}\left(\bigcup_{\begin{subarray}{c}{\mathcal{S}\in\Sigma(\ell),}\\\
{\bar{X}^{n}_{\mathcal{S}}\neq
X^{n}_{\mathcal{S}}}\end{subarray}}\left\\{\gamma(\bar{X}^{n}_{\mathcal{S}})>\gamma(X^{n}_{\mathcal{S}})\right\\}\bigg{|}X^{n}_{[\textsf{K}_{\text{a}}]},Y^{n},{\mathcal{E}}_{\ell}\right)\right]$
(23)
$\displaystyle\leq\mathds{E}\bigg{[}\min\bigg{\\{}1,\sum_{\mathcal{S}\in\Sigma(\ell)}\binom{{\textsf{M}}-\textsf{K}_{\text{a}}}{|\mathcal{S}|}|\mathcal{S}|!$
$\displaystyle\hskip
68.28644pt\left.\cdot\mbox{Pr}\left(\gamma(\bar{X}^{n}_{\mathcal{S}})>\gamma(X^{n}_{\mathcal{S}})\bigg{|}X^{n}_{[\textsf{K}_{\text{a}}]},Y^{n},{\mathcal{E}}_{\ell}\right)\right]$
(24)
$\displaystyle\leq\mathds{E}\bigg{[}\min\bigg{\\{}1,\sum_{\mathcal{S}\in\Sigma(\ell)}{\textsf{M}}^{|\mathcal{S}|}\exp(-\gamma(X^{n}_{\mathcal{S}}))\bigg{\\}}\bigg{]}$
(25)
$\displaystyle\leq\sum_{\mathcal{S}\in\Sigma(\ell)}\mbox{E}\bigg{[}\min\bigg{\\{}1,{\textsf{M}}^{|\mathcal{S}|}\exp(-\gamma(X^{n}_{\mathcal{S}}))\bigg{\\}}\bigg{]}$
(26)
$\displaystyle\leq\sum_{\mathcal{S}\in\Sigma(\ell)}\mbox{Pr}\bigg{(}{\textsf{M}}^{|\mathcal{S}|}\exp(-\gamma(X^{n}_{\mathcal{S}}))\geq
U\bigg{)}$ (27)
$\displaystyle=\sum_{\mathcal{S}\in\Sigma(\ell)}\mbox{Pr}\bigg{(}\log\Big{(}{\textsf{M}}^{|\mathcal{S}|}\exp(-\gamma(X^{n}_{\mathcal{S}}))\Big{)}\\!-\\!\log(U)\\!\geq\\!0\bigg{)}\\!\\!$
(28)
$\displaystyle=\sum_{\mathcal{S}\in\Sigma(\ell)}\mbox{Pr}(W_{\mathcal{S}}\geq
0),$ (29)

where (21) is due to the definition of the maximum information density
decoder, (22) is due to the chain rule of information density. The random
coding scheme and union bound are used in (23) and (24), respectively. Note
that the asynchronous model does not have the permutation-invariant property.
Therefore, the number of permutations of the erroneously decoded messages,
$|\mathcal{S}|!$, is summed up. The inequality (25) follows from the fact that
$\binom{{\textsf{M}}-\textsf{K}_{\text{a}}}{|\mathcal{S}|}\cdot|\mathcal{S}|!\leq{\textsf{M}}^{|\mathcal{S}|}$
and

$\mbox{Pr}(\gamma(\bar{X}^{n}_{\mathcal{S}})>\gamma(X^{n}_{\mathcal{S}}))\leq\exp(-\gamma(X^{n}_{\mathcal{S}})),$

where $\bar{X}^{n}_{\mathcal{S}}$ is an independent copy of
$X^{n}_{\mathcal{S}}$[17, Corollary 18.4]. The inequality (26) follows from
$\min\\{1,\beta_{1}+\beta_{2}\\}\leq\min\\{1,\beta_{1}\\}+\min\\{1,\beta_{2}\\}$
for $\beta_{1},\beta_{2}\in\mathbb{R}$. The inequality (27) follows from
$\mathds{E}[\min\\{1,V\\}]=\mbox{Pr}(V\geq U)$ [16, eq.(77)] for a non-
negative random variable $V$, where $U\sim\mbox{Unif}(0,1)$ is independent of
$V$. The equality (29) follows from defining

$W_{\mathcal{S}}:=\log\left({\textsf{M}}^{|\mathcal{S}|}\exp(-\gamma(X_{\mathcal{S}}^{n}))\right)-\log(U).$

We apply the CGF, the Taylor expansion, and the inverse Laplace transform to
derive $\mbox{Pr}(W_{\mathcal{S}}\geq 0)$. We denote by
$\psi_{W_{\mathcal{S}}}(t)=\log(\mathds{E}[\exp(tW_{\mathcal{S}})])$ the CGF
of the random variable $W_{\mathcal{S}}$ with parameter $t$.

$\displaystyle\psi_{W_{\mathcal{S}}}(t)$
$\displaystyle=\log\bigg{(}\mathds{E}\left[\exp\left(t\log\left({\textsf{M}}^{|\mathcal{S}|}\exp(-\gamma(X_{\mathcal{S}}^{n}))\right)\right.\right.-t\log(U)\Big{)}\Big{]}\bigg{)}$
(30)
$\displaystyle=t|\mathcal{S}|{\log\textsf{M}}-\log(1-t)+\log(\mathds{E}[\exp(-t\cdot\gamma(X^{n}_{\mathcal{S}}))])$
(31) $\displaystyle=t|\mathcal{S}|{\log\textsf{M}}-\log(1-t)-E(a^{n},t)$ (32)
$\displaystyle=\tilde{\psi}_{W_{\mathcal{S}}}(t)-\log(1-t),$ (33)

where $t\in(0,1)$,
$\tilde{\psi}_{W_{\mathcal{S}}}(t):=t|\mathcal{S}|{\log\textsf{M}}-E(a^{n},t)$,
and (32) is due to the following definition in Theorem 1,

$\displaystyle E(a^{n},t)$
$\displaystyle:=-\log(\mathds{E}[\exp(-t\gamma(X^{n}_{\mathcal{S}}))])$
$\displaystyle=\frac{1}{2}\sum\limits_{i=1}^{n}\left(t\log(1+a_{i}\mathsf{P})+\log\left(1-\frac{a_{i}\mathsf{P}t^{2}}{1+a_{i}\mathsf{P}}\right)\right),$

where

$\displaystyle\exp(\\!-t\\!\cdot\\!\gamma(X^{n}_{\mathcal{S}}))\\!\\!=\\!\prod_{\ell=1}^{n}\\!\left(\\!\frac{dP_{Y|X_{[\textsf{K}_{\text{a}}]}}\\!\left(Y_{\ell}|\big{\\{}\tau_{d_{m}}(X^{n}_{m},\ell)\big{\\}}_{m\in[\textsf{K}_{\text{a}}]}\right)}{dP_{Y}(Y_{\ell})}\\!\right)^{\\!t}\\!.$

For $t\in(0,1)$, the CGF converges, which is proved as follows. Since the CGF
is the summation of the logarithm of the following $n$ terms,

$\displaystyle\mathds{E}[\exp(t\cdot
i(\\{\tau_{d_{m}}(X^{n}_{m},\ell)\\}_{m\in\mathcal{S}};Y_{\ell}|\\{\tau_{d_{m}}(X^{n}_{m},\ell)\\}_{m\in[\textsf{K}_{\text{a}}]\setminus\mathcal{S}}))],$
(34)

$\ell=1,2,...,n$, for a CGF to converge, a sufficient condition is that (34)
converges in term of $t$ for all $\ell\in[n]$. We apply the Gaussian integral
to derive (34). The corresponding range of convergence for any $\ell\in[n]$ is
$t\in\left(-\frac{1+a_{\ell}\mathsf{P}}{a_{\ell}\mathsf{P}},\sqrt{\frac{1+a_{\ell}\mathsf{P}}{a_{\ell}\mathsf{P}}}\right)$.
When $t<0$, it is possible that
$|\mathcal{S}|{\log\textsf{M}}>\sum_{\ell=1}^{n}\frac{1}{2}\log(1+a_{\ell}\mathsf{P})$,
which means that the corresponding error probability approaches $1$. For all
$\ell\in[n]$ and $a_{\ell}\in\mathbb{Z}^{+}$,
$\sqrt{\frac{1+a_{\ell}\mathsf{P}}{a_{\ell}\mathsf{P}}}\geq 1$. If
$a_{\ell}=0$, $\ell\in[n]$, it represents that no codeword symbol is
transmitted at the $\ell$-th channel use. Thus, the information density is $0$
and the corresponding (34) must converge. Therefore, in Theorem 1 and Theorem
2, we choose $t=t_{0}(a^{n})\in(0,1)$, which fulfills

$\displaystyle|\mathcal{S}|{\log\textsf{M}}=E^{(1)}_{t}(a^{n},t_{0}(a^{n})),$
(35)

to guarantee the convergence.

The PDF of $W_{\mathcal{S}}$ is obtained by the inverse Laplace transform:

$\displaystyle f_{W_{\mathcal{S}}}(w)=\frac{1}{2\pi
j}\int^{c+j\infty}_{c-j\infty}\exp(\psi_{W_{\mathcal{S}}}(t)-tw)dt,$ (36)

where $c\in(0,1).$ The probability, $\mbox{Pr}(W_{\mathcal{S}}\geq 0)$, is
obtained by changing the order of integration, i.e.,

$\displaystyle\mbox{Pr}(W_{\mathcal{S}}\geq 0)$ $\displaystyle=\frac{1}{2\pi
j}\int^{\infty}_{0}\left\\{\int^{c+j\infty}_{c-j\infty}\exp(\psi_{W_{\mathcal{S}}}(t)-tw)dt\right\\}dw$
(37) $\displaystyle=\frac{1}{2\pi
j}\int^{c+j\infty}_{c-j\infty}\exp(\psi_{W_{\mathcal{S}}}(t))\frac{dt}{t}$
(38) $\displaystyle=\frac{1}{2\pi
j}\int^{c+j\infty}_{c-j\infty}\exp(\tilde{\psi}_{W_{\mathcal{S}}}(t))\frac{dt}{t(1-t)}.$
(39)

The last equality follows from (33). By applying the Taylor expansion to
$\tilde{\psi}_{W_{\mathcal{S}}}(t)$ at the point $t=t_{0}(a^{n})$, which
fulfills (35), we have

$\displaystyle\tilde{\psi}_{W_{\mathcal{S}}}(t)=$ $\displaystyle
t_{0}(a^{n})|\mathcal{S}|{\log\textsf{M}}-E(a^{n},t_{0}(a^{n}))$
$\displaystyle+\big{[}|\mathcal{S}|{\log\textsf{M}}-E^{(1)}_{t}(a^{n},t_{0}(a^{n}))\big{]}(t-t_{0}(a^{n}))$
$\displaystyle-E^{(2)}_{t}(a^{n},t_{0}(a^{n}))\frac{(t\\!-\\!t_{0}(a^{n}))^{2}}{2}\\!+\\!\bar{\xi}(a^{n},t_{0}(a^{n}),t),$
(40)

where
$\big{[}|\mathcal{S}|{\log\textsf{M}}-E^{(1)}_{t}(a^{n},t_{0}(a^{n}))\big{]}(t-t_{0}(a^{n}))=0$
due to (35),

$\bar{\xi}(a^{n},t_{0}(a^{n}),t):=\sum_{i=3}^{\infty}-E^{(i)}_{t}(a^{n},t_{0}(a^{n}))\frac{(t-t_{0}(a^{n}))^{i}}{i!}$

is the sum of higher order terms of Taylor expansion, and $t_{0}(a^{n})$
satisfies (35). Substitute (40) and $c=t_{0}(a^{n})$ into (39), we have

$\displaystyle\frac{1}{2\pi
j}\int^{t_{0}(a^{n})+j\infty}_{t_{0}(a^{n})-j\infty}\exp(\tilde{\psi}_{W_{\mathcal{S}}}(t))\frac{dt}{t(1-t)}$
$\displaystyle=\frac{\eta}{j}\\!\int^{t_{0}(a^{n})\\!+j\infty}_{t_{0}(a^{n})-j\infty}\\!\\!\exp\\!\\!\bigg{(}\\!\beta\frac{(t\\!-\\!t_{0}(a^{n}))^{2}\\!}{2}\\!+\\!\bar{\xi}(a^{n},t_{0}(a^{n}),t)\\!\bigg{)}\frac{dt}{t(1\\!-\\!t)}$
$\displaystyle=\\!\\!\eta\bigg{\\{}\frac{1}{j}\int^{t_{0}(a^{n})\\!+j\infty}_{t_{0}(a^{n})-j\infty}\\!\exp\bigg{(}\beta\frac{(t\\!-\\!t_{0}(a^{n}))^{2}}{2}\bigg{)}\frac{dt}{t(1\\!-\\!t)}$
$\displaystyle\hskip 142.26378pt+2\pi\xi(a^{n},t_{0}(a^{n}))\bigg{\\}}\\!\\!$
(42)
$\displaystyle\\!=\\!\eta\bigg{\\{}\\!\int^{\infty}_{-\infty}\\!\\!\exp\\!\bigg{(}-\beta\frac{\rho^{2}}{2}\\!\bigg{)}\frac{d\rho}{t_{0}(a^{n})\\!\\!+\\!j\rho}$
$\displaystyle\quad\\!+\\!\int^{\infty}_{-\infty}\\!\exp\bigg{(}\\!\\!-\beta\frac{\rho^{2}}{2}\bigg{)}\frac{d\rho}{1\\!-\\!t_{0}(a^{n})\\!-\\!j\rho}\\!+\\!2\pi\xi(a^{n},t_{0}(a^{n}))\bigg{\\}},$
(43)

where $\eta\\!:=\\!\frac{g_{1}(a^{n}\\!,t_{0}(a^{n}))}{2\pi}$,
$\beta\\!:=\\!-E^{(2)}_{t}(a^{n},t_{0}(a^{n}))$,
$\rho\\!:=\\!\frac{t-t_{0}(a^{n})}{j}$, and

$g_{1}(a^{n},t):=\exp(t|\mathcal{S}|{\log\textsf{M}}-E(a^{n},t)).$

The equality (42) follows from $e^{x}=1+\sum_{i=1}^{\infty}\frac{x^{i}}{i!}$
and by letting $x=\bar{\xi}(a^{n},t_{0}(a^{n}),t)$,

$\displaystyle\xi(a^{n},t_{0}(a^{n})):=\\!\frac{1}{2\pi
j}\int^{t_{0}(a^{n})+j\infty}_{t_{0}(a^{n})-j\infty}\\!\\!\exp\left(\frac{\beta}{2}(t\\!-\\!t_{0}(a^{n}))^{2}\\!\right)$
$\displaystyle\hskip
85.35826pt\cdot\frac{1}{t(1-t)}\sum_{m=1}^{\infty}\frac{\bar{\xi}(a^{n},t_{0}(a^{n}),t)^{m}}{m!}dt.$
(44)

By multiplying $\frac{t_{0}(a^{n})-j\rho}{t_{0}(a^{n})-j\rho}$ to the first
integral in (43), we have

$\displaystyle\int^{\infty}_{-\infty}\exp\bigg{(}-\beta\frac{\rho^{2}}{2}\bigg{)}\frac{d\rho}{t_{0}(a^{n})+j\rho}$
$\displaystyle=\int^{\infty}_{-\infty}\exp\bigg{(}-\beta\frac{\rho^{2}}{2}\bigg{)}\frac{t_{0}(a^{n})d\rho}{t_{0}(a^{n})^{2}+\rho^{2}}$
$\displaystyle\hskip
56.9055pt-\int^{\infty}_{-\infty}\exp\bigg{(}-\beta\frac{\rho^{2}}{2}\bigg{)}\frac{j\rho
d\rho}{t_{0}(a^{n})^{2}+\rho^{2}}$ (45)
$\displaystyle=\int^{\infty}_{-\infty}\exp\bigg{(}-\beta\frac{\rho^{2}}{2}\bigg{)}\frac{t_{0}(a^{n})d\rho}{t_{0}(a^{n})^{2}+\rho^{2}}$
(46)
$\displaystyle=2\pi\exp\left(\frac{t_{0}(a^{n})^{2}\beta}{2}\right)Q\left(t_{0}(a^{n})\sqrt{\beta}\right)$
(47) $\displaystyle\leq\frac{\sqrt{2\pi}}{t_{0}(a^{n})}\frac{1}{\sqrt{\beta}}$
(48)
$\displaystyle=\frac{\sqrt{2\pi}}{t_{0}(a^{n})}\frac{1}{\sqrt{-E^{(2)}_{t}(a^{n},t_{0}(a^{n}))}},$
(49)

where the second integral in (45) is the integral of an odd function, which
equals $0$. By applying the Voigt function [18] to the integral in (46), we
have (47). The inequality (49) follows from the upper bound of the Gaussian
Q-function, $Q(x)\leq\frac{1}{x\sqrt{2\pi}}\exp(-\frac{x^{2}}{2})$. The last
equality follows from the definition:
$\beta\\!:=\\!-E^{(2)}_{t}(a^{n},t_{0}(a^{n}))$. By multiplying
$\frac{1-t_{0}(a^{n})+j\rho}{1-t_{0}(a^{n})+j\rho}$ with the same steps used
in deriving (49), the second integral in (43) is bounded by

$\displaystyle\\!\\!\\!\\!\int^{\infty}_{-\infty}\\!\\!\exp\bigg{(}\\!\\!-\beta\frac{\rho^{2}}{2}\bigg{)}\frac{d\rho}{1\\!-\\!t_{0}(a^{n})\\!-\\!j\rho}$
$\displaystyle\hskip
71.13188pt\leq\frac{\sqrt{2\pi}}{1-t_{0}(a^{n})}\frac{1}{\sqrt{-E^{(2)}_{t}(a^{n},t_{0}(a^{n}))}}.\\!\\!\\!$
(50)

Consequently, we can upper bound the sum of the two integrations in (43) as
follows:

$\displaystyle\eta\bigg{\\{}\\!\int^{\infty}_{-\infty}\\!\\!\exp\\!\bigg{(}-\beta\frac{\rho^{2}}{2}\\!\bigg{)}\frac{d\rho}{t_{0}(a^{n})\\!\\!+\\!j\rho}$
$\displaystyle\quad\\!+\\!\int^{\infty}_{-\infty}\\!\exp\bigg{(}\\!\\!-\beta\frac{\rho^{2}}{2}\bigg{)}\frac{d\rho}{1\\!-\\!t_{0}(a^{n})\\!-\\!j\rho}\\!+\\!2\pi\xi(a^{n},t_{0}(a^{n}))\bigg{\\}}$
$\displaystyle\leq\frac{g_{1}(a^{n},t_{0}(a^{n}))}{\sqrt{2\pi}(1-t_{0}(a^{n}))t_{0}(a^{n})}\frac{1}{\sqrt{-E^{(2)}_{t}(a^{n},t_{0}(a^{n}))}}$
$\displaystyle\hskip
88.2037pt+g_{1}(a^{n},t_{0}(a^{n}))\xi(a^{n},t_{0}(a^{n})).$ (51)

By combining (19), (29), (39), (43), and (51), the PUPE of the AUMAC system
for a given $D^{\textsf{K}_{\text{a}}}\\!$ is

$\displaystyle\sum_{\ell=1}^{\textsf{K}_{\text{a}}}\frac{1}{\textsf{K}_{\text{a}}}\mbox{Pr}(M_{\ell}\not\in
g(Y^{n})|D^{\textsf{K}_{\text{a}}}\\!,{\mathcal{E}}_{\ell})+p_{0}$
$\displaystyle\leq\sum_{\ell=1}^{\textsf{K}_{\text{a}}}\\!\frac{1}{\textsf{K}_{\text{a}}}\\!\sum_{\begin{subarray}{c}\mathcal{S}\in\Sigma(\ell)\end{subarray}}\\!\\!\bigg{\\{}\\!\frac{g_{1}(a^{n},t_{0}(a^{n}))}{(1\\!-\\!t_{0}(a^{n}))t_{0}(a^{n})}\frac{1}{\sqrt{\\!-2\pi
E^{(2)}_{t}\\!(a^{n},t_{0}(a^{n}))}}$ $\displaystyle\hskip
65.44142pt+g_{1}(a^{n},t_{0}(a^{n}))\xi(a^{n},t_{0}(a^{n}))\bigg{\\}}+p_{0}$
(52)
$\displaystyle=\sum_{\mathcal{S}\subseteq[\textsf{K}_{\text{a}}]}\frac{|\mathcal{S}|}{\textsf{K}_{\text{a}}}\bigg{\\{}\frac{g_{1}(a^{n},t_{0}(a^{n}))}{(1-t_{0}(a^{n}))t_{0}(a^{n})}\frac{1}{\sqrt{-2\pi
E^{(2)}_{t}(a^{n},t_{0}(a^{n}))}}$ $\displaystyle\hskip
62.59596pt+g_{1}(a^{n},t_{0}(a^{n}))\xi(a^{n},t_{0}(a^{n}))\bigg{\\}}+p_{0},$
(53)

where $\Sigma(\ell)$ is defined in (20).

##  Appendix B Proof of Theorem 2

In the following, in addition to Theorem 1, we derive a uniform upper bound of
the PUPE of an $(n,{\textsf{M}},\epsilon,\textsf{K}_{\text{a}},\alpha)$-code
as indicated in Theorem 2. In particular, we will find the worst-case
asynchronicity, which implies finding the worst-case of $a^{n}$ and
$t_{0}(a^{n})$ in Theorem 1. To simplify the derivation, we denote
$\iota:=\mathbbm{1}(1\in\mathcal{S})$ and all possible $a^{n}$’s w.r.t.
$\iota$ by the set $\mathcal{F}_{k,\iota}:=\\{a^{n}:\iota,|\mathcal{S}|=k\\}$,
where $a^{n}$ is defined in (1) as a function of $\mathcal{S}$ and
$D^{\textsf{K}_{\text{a}}}\\!$.

We will show that for all $t\in(0,1)$, there exists an $a_{\iota}^{n*}$
resulting in a uniform upper bound of PUPE for all
$a^{n}\in\mathcal{F}_{|\mathcal{S}|,\iota}$, such that the upper bound of the
PUPE in (5) has the following property

$\displaystyle
g_{1}(a^{n}\\!,t)g_{2}(a^{n}\\!,t)\leq\\!g_{1}(a_{\iota}^{n*}\\!,t)g_{2}(a_{\iota}^{n*}\\!,t).$

However, for any $a^{n}\\!\in\\!\mathcal{F}_{|\mathcal{S}|,\iota}$, the order
between

$g_{1}(a^{n}\\!,t_{0}(a^{n}))g_{2}(a^{n}\\!,t_{0}(a^{n}))$

and

$g_{1}(a_{\iota}^{n*}\\!,t_{0}(a_{\iota}^{n*}))g_{2}(a_{\iota}^{n*}\\!,t_{0}(a_{\iota}^{n*}))$

is not fixed, since the sign of $\frac{\partial}{\partial
t}g_{1}(a^{n},t)g_{2}(a^{n},t)$ is not the same for all $t\in(0,1)$.
Therefore, for fixed $a_{\iota}^{n*}$, we will show that the choices of
$T^{*}_{0}$, $T^{*}_{1}$, $\underline{t}_{0}$, and $\underline{t}_{1}$
uniformly upper bound the PUPE regardless $D^{\textsf{K}_{\text{a}}}\\!$. We
start from (5) restated as follows

$\displaystyle\sum_{\mathcal{S}\subseteq[\textsf{K}_{\text{a}}]}\frac{|\mathcal{S}|}{\textsf{K}_{\text{a}}\sqrt{2\pi}}g_{1}(a^{n},t_{0}(a^{n}))g_{2}(a^{n},t_{0}(a^{n})),$
(54)

while omitting the term $p_{0}$ and also the approximation error term
$\xi(a^{n},t_{0}(a^{n}))$ since we do not bound these terms.

To proceed, we use the following lemma.

###### Lemma 1

Let $g_{1}(a^{n},t)=\exp(f_{1}(a^{n},t))$ and
$g_{2}(a^{n},t)=(f_{2}(a^{n},t))^{-\frac{1}{2}}$, where
$a^{n}\in\\{\mathbb{Z}^{+}_{0}\\}^{n}$, $t\in(0,1)$,
$f_{1}(a^{n},t)\in\mathbb{R}$ and $f_{2}(a^{n},t)>0$. Then
$g_{1}(a^{n},t)g_{2}(a^{n},t)$ is a non-increasing function w.r.t.
$a_{i},\,\forall i\\!\in\\![n]$ if $f_{1,a_{i}}^{(1)}(a_{i},t)\\!\leq\\!0$ and
$f_{2,a_{i}}^{(1)}(a_{i},t)\\!\geq\\!0$.

The proof of Lemma 1 is relegated to Appendix C.

Then, we apply Lemma 1 by defining
$f_{1}(a^{n},t):=t|\mathcal{S}|{\log\textsf{M}}-E(a^{n},t)$ and
$f_{2}(a^{n},t):=-(t-t^{2})^{2}E_{t}^{(2)}(a^{n},t)$. The first derivatives of
$f_{1}(a^{n},t)$ and $f_{2}(a^{n},t)$ w.r.t. $a_{i}$ are expressed as follows,
respectively

$\displaystyle f_{1,a_{i}}^{(1)}(a_{i},t)=$
$\displaystyle\frac{\mathsf{P}(t^{2}-t)+a_{i}\mathsf{P}^{2}(t^{3}-t)}{2(1+a_{i}\mathsf{P})(1+a_{i}\mathsf{P}-a_{i}\mathsf{P}t^{2})}$
(55)

and

$\displaystyle f_{2,a_{i}}^{(1)}(a_{i},t)=$
$\displaystyle\frac{\mathsf{P}(1-t)^{2}t^{2}(1+a_{i}\mathsf{P}+3a_{i}\mathsf{P}t^{2})}{(1+a_{i}\mathsf{P}-a_{i}\mathsf{P}t^{2})^{3}}.$
(56)

For $t\in(0,1)$, it is clear that $f^{(1)}_{1,a_{i}}(a_{i},t)\leq 0$ and
$f^{(1)}_{2,a_{i}}(a_{i},t)\geq 0$. We then conclude that
$g_{1}(a^{n},t)g_{2}(a^{n},t)$ is a non-increasing function w.r.t.
$a_{i},\;i\in[\alpha n]$ according to Lemma 1. It implies that the PUPE of any
given $\mathcal{S}$ decreases with increasing $a_{i},\;i\in[\alpha n]$.
Namely, reducing $a_{i},\;i\in[\alpha n]$ will upper bound the error
probability. Therefore, to upper bound the PUPE, we can consider the following
case, where the number of transmitted symbols that belong to $\mathcal{S}$ at
the first $\alpha n$ channel use, $a_{[\alpha n]}$, are reduced to the
minimum, which is $a^{*}_{\iota,[\alpha n]}=\iota$. Namely
$a^{n*}_{\iota}=[\iota^{\alpha n},|\mathcal{S}|^{n-\alpha n}]$. Consequently,
for all $a^{n}\in{\mathcal{F}}_{|\mathcal{S}|,\iota}$ and a given
$t=t_{0}(a^{n})$, we have

$\displaystyle g_{1}(a^{n},t_{0}(a^{n}))g_{2}(a^{n},t_{0}(a^{n}))$
$\displaystyle\hskip 71.13188pt\leq
g_{1}(a^{n*}_{\iota},t_{0}(a^{n}))g_{2}(a^{n*}_{\iota},t_{0}(a^{n})).$ (57)

We have shown that the error probability is non-increasing w.r.t. $a_{i}$.
However, the sign of $\frac{\partial}{\partial t}g_{1}(a^{n},t)g_{2}(a^{n},t)$
w.r.t. $t$ changes for $t\in(0,1)$. To solve it, we can show that given
$a^{n*}_{\iota}$, if $t_{0}(a^{n*}_{\iota})\in\mathcal{A}\cap\mathcal{B}$,
$\bar{t}_{\iota}\in\mathcal{A}\cap\bar{\mathcal{B}}$,
$\underline{t}_{\iota}\in\mathcal{A}\cap\underline{\mathcal{B}}$, and
$\underline{t}_{\iota}\leq t_{0}(a^{n})\leq\bar{t}_{\iota}$, then there exist
a uniform upper bound of the error probability for all
$D^{\textsf{K}_{\text{a}}}\\!$ satisfying delay constraint $\alpha n$, where

$\displaystyle\mathcal{A}\\!:=\left\\{t\\!:\\!t\in(0,1)\right\\},$ (58)
$\displaystyle\mathcal{B}\\!:=\left\\{t\\!:\\!E^{(1)}_{t}(a^{n*}_{\iota},t)=|\mathcal{S}|{\log\textsf{M}}\right\\},$
(59)
$\displaystyle\bar{\mathcal{B}}\\!:=\left\\{t\\!:\\!\sum_{i=1}^{n}\frac{a^{*}_{\iota,i}\mathsf{P}t}{1\\!+\\!a^{*}_{\iota,i}\mathsf{P}(1\\!-\\!t^{2})}\\!=\\!\frac{n}{2}\log(1\\!+\\!|\mathcal{S}|\mathsf{P})\\!-\\!|\mathcal{S}|{\log\textsf{M}}\right\\},$
(60)
$\displaystyle\underline{\mathcal{B}}\\!:=\left\\{t\\!:\\!\frac{|\mathcal{S}|n\mathsf{P}t}{1\\!+\\!|\mathcal{S}|\mathsf{P}\\!-\\!|\mathcal{S}|\mathsf{P}t^{2}}\\!=\\!\sum_{i=1}^{n}\frac{1}{2}\log(1\\!\\!+\\!a^{*}_{\iota,i}\mathsf{P})\\!-\\!|\mathcal{S}|{\log\textsf{M}}\right\\},$
(61)

and $a^{*}_{\iota,i}$ is the $i$-th element of $a^{n*}_{\iota}$.

To proceed, we find upper bounds of $g_{1}(a^{n*}_{\iota},t_{0}(a^{n}))$ and
$g_{2}(a^{n*}_{\iota},t_{0}(a^{n}))$ as $u_{1}$ and $u_{2}$. Then we upper
bound $g_{1}(a^{n*}_{\iota},t_{0}(a^{n}))g_{2}(a^{n*}_{\iota},t_{0}(a^{n}))$
by $u_{1}u_{2}$. Since the second partial derivative w.r.t. $t$,

$\displaystyle
f^{(2)}_{1,t}(a^{n},t)=\sum_{i=1}^{n}\frac{(a_{i}\mathsf{P})+(a_{i}\mathsf{P})^{2}+(a_{i}\mathsf{P})^{2}t^{2}}{(1+a_{i}\mathsf{P}-a_{i}\mathsf{P}t^{2})^{2}},$
(62)

is positive for $t\in(0,1)$, $f_{1}(a^{n},t)$ is a convex function regarding
$t$. Moreover, $f^{(1)}_{1,t}(a^{n},t_{0}(a^{n}))=0$ by (59). Namely,
$f_{1}(a^{n},t)$ achieves minimum at $t=t_{0}(a^{n})$. Therefore, for any
$a^{n}\in\mathcal{F}_{|\mathcal{S}|,\iota}$, we have

$\displaystyle g_{1}(a^{n},t_{0}(a^{n}))\leq
g_{1}(a^{n},t_{0}(a^{n*}_{\iota}))\leq
g_{1}(a^{n*}_{\iota},t_{0}(a^{n*}_{\iota})),$ (63)

where the first inequality is because $g_{1}(a^{n},t)$ achieves the minimum at
$t=t_{0}(a^{n})$. If $a^{n}=a^{n*}_{\iota}$, the equalities hold. The second
inequality follows from the fact that $g_{1}(a^{n},t)$ is a non-decreasing
function for a given $t\in(0,1)$ w.r.t. $a_{i},\;\forall i\in[\alpha n]$,
since

$g_{1,a_{i}}^{(1)}(a_{i},t)=\exp(f_{1}(a^{n},t))\cdot
f_{1,a_{i}}^{(1)}(a_{i},t)\leq 0,$

for $t\in(0,1)$, where $f_{1,a_{i}}^{(1)}(a_{i},t)$ is given in (55).

We define $f_{2}(a^{n},t):=(f_{3}(t))^{2}f_{4}(a^{n},t)$, where
$f_{3}(t):=t-t^{2}$ and $f_{4}(a^{n},t):=-E^{(2)}_{t}(a^{n},t)$. Since the
first partial derivative w.r.t. $t$ of $f_{4,t}(a^{n},t)$ is as follows

$\displaystyle
f^{(1)}_{4,t}(a^{n},t)=\sum_{i=1}^{n}2(a_{i}\mathsf{P})^{2}t\frac{3+3a_{i}\mathsf{P}+a_{i}\mathsf{P}t^{2}}{(1+a_{i}\mathsf{P}-a_{i}\mathsf{P}t^{2})^{3}},$
(64)

which is positive for $t\in(0,1)$, $f_{4}(a^{n},t)$ is a non-decreasing
function of $t$. Then, we have

$\displaystyle f_{4}(a^{n*}_{\iota},t_{0}(a^{n*}_{\iota}))\geq
f_{4}(a^{n*}_{\iota},\underline{t}_{\iota}).$ (65)

By the condition, $\underline{t}_{\iota}\leq t_{0}(a^{n})\leq\bar{t}_{\iota}$,
there exists a $\lambda\in[0,1]$ such that
$t_{0}(a^{n})=\lambda\underline{t}_{\iota}+\bar{\lambda}\bar{t}_{\iota}$,
where $\bar{\lambda}=1-\lambda$. Since $f_{3}(t)$ is concave, it satisfies

$\displaystyle f_{3}(t_{0}(a^{n}))\geq$ $\displaystyle\lambda
f_{3}(\underline{t}_{\iota})+\bar{\lambda}f_{3}(\bar{t}_{\iota})$ (66)
$\displaystyle\geq$
$\displaystyle\lambda\min\\{f_{3}(\bar{t}_{\iota}),f_{3}(\underline{t}_{\iota})\\}+\bar{\lambda}\min\\{f_{3}(\bar{t}_{\iota}),f_{3}(\underline{t}_{\iota})\\}$
(67) $\displaystyle=$
$\displaystyle\min\\{f_{3}(\bar{t}_{\iota}),f_{3}(\underline{t}_{\iota})\\}=:T^{*}_{\iota}.$
(68)

Consequently, for all $a^{n}\in\mathcal{F}_{|\mathcal{S}|,\iota}$, we have

$\displaystyle g_{1}(a^{n*}_{\iota},t_{0}(a^{n}))\leq
g_{1}(a^{n*}_{\iota},t_{0}(a^{n*}_{\iota})),$ (69)

stated in (63), and

$\displaystyle g_{2}(a^{n*}_{\iota},t_{0}(a^{n}))$
$\displaystyle=\frac{1}{f_{3}(t_{0}(a^{n}))\sqrt{f_{4}(a^{n*}_{\iota},t_{0}(a^{n}))}}$
(70)
$\displaystyle\leq\frac{1}{T^{*}_{\iota}\sqrt{f_{4}(a^{n*}_{\iota},t_{0}(a^{n}))}}$
(71)
$\displaystyle\leq\frac{1}{T^{*}_{\iota}\sqrt{f_{4}(a^{n*}_{\iota},\underline{t}_{\iota})}},$
(72)

where (70) is by definition, (71) follows from (68) and (72) follows from
(65).

Consequently, we have

$\displaystyle g_{1}(a^{n},t_{0}(a^{n}))g_{2}(a^{n},t_{0}(a^{n}))$
$\displaystyle\leq
g_{1}(a^{n*}_{\iota},t_{0}(a^{n}))g_{2}(a^{n*}_{\iota},t_{0}(a^{n}))$ (73)
$\displaystyle\leq
g_{1}(a^{n*}_{\iota},t_{0}(a^{n*}_{\iota}))g_{2}(a^{n*}_{\iota},t_{0}(a^{n}))$
(74)
$\displaystyle\leq\frac{g_{1}(a^{n*}_{\iota},t_{0}(a^{n*}_{\iota}))}{T^{*}_{\iota}\sqrt{-E^{(2)}_{t}(a^{n*}_{\iota},\underline{t}_{\iota})}},$
(75)

which completes the proof of Theorem 2.

##  Appendix C Proof of Lemma 1

Let $g_{1}(a^{n},t)=\exp(f_{1}(a^{n},t))$ and
$g_{2}(a^{n},t)=(f_{2}(a^{n},t))^{-\frac{1}{2}}$, where
$a^{n}\in\\{\mathbb{Z}^{+}_{0}\\}^{n}$, $t\in(0,1)$ and
$f_{1}(a^{n},t)\in\mathbb{R}$, and $f_{2}(a^{n},t)>0$. Then the first partial
derivative of $g_{1}(a^{n},t)g_{2}(a^{n},t)$ w.r.t. $a_{i}$ is

$\displaystyle\frac{\partial}{\partial a_{i}}g_{1}(a^{n},t)g_{2}(a^{n},t)$
$\displaystyle=g_{2}(a^{n},t)\frac{\partial}{\partial
a_{i}}g_{1}(a^{n},t)+g_{1}(a^{n},t)\frac{\partial}{\partial
a_{i}}g_{2}(a^{n},t)$ (76)
$\displaystyle=g_{2}(a^{n},t)\exp(f_{1}(a^{n},t))\frac{\partial}{\partial
a_{i}}f_{1}(a^{n},t)$ $\displaystyle\hskip
42.67912pt-\frac{1}{2}g_{1}(a^{n},t)(f_{2}(a^{n},t))^{-\frac{3}{2}}\frac{\partial}{\partial
a_{i}}f_{2}(a^{n},t).$ (77)

Therefore, $g_{1}(a^{n},t)g_{2}(a^{n},t)$ is a non-increasing function w.r.t.
$a_{i}$, if $\frac{\partial}{\partial a_{i}}f_{1}(a^{n},t)\leq 0$ and
$\frac{\partial}{\partial a_{i}}f_{2}(a^{n},t)\geq 0$.


# Fast Metasurface Hybrid Lens Design using a Semi-Analytical Model

A. Cléroux Cuillerier∗, J. Borne and S. Thibault
Centre d’optique, photonique et laser, Université Laval
Université Laval
Quebec City, Canada
∗<EMAIL_ADDRESS>

###### Abstract

We propose a new method for integrating metasurfaces in optical design using
semi-analytical modelling of dielectric nanostructures. The latter computes
the output phase of an electric field incident on the metasurface, allowing
their use with ray-tracing software. This tool provides a method to use
metasurfaces in optical systems while using built-in optimization processes to
avoid time-consuming computation. To demonstrate the applicability and
versatility of our method, we present variations of a triplet composed of
refractive elements and a metasurface. For each of the systems, similar
optical performances are achieved. Our unique and innovative approach to
joining metasurfaces and ray tracing has the potential to promote the design
of innovative systems by exploiting the richness of metasurfaces and the
functionality of conventional lens design software.

_K_ eywords Optics $\cdot$ Metasurface $\cdot$ Lens Design

##  1 Introduction

Research on metasurfaces has been motivated by their broad potential for
manipulating light. Exotic applications exploiting their increased control on
the wavefront, their polarization sensitivity or their unique chromatic
dispersion were proposed in the past years [4, 1, 3, 17, 12]. Many of these
technologies were demonstrated by wave optics computation, sometimes followed
by experimental work. Indeed, high nonlinear responses [16], multifunctional
metalens [23, 11] or non-Hermitian metasurfaces [9] can be found in the
literature, to cite a few.

The utility of metasurfaces has been praised many times; some even say they
have the potential to replace traditional optics [20, 26, 14]. However, their
integration in conventional lens design remains a little-explored avenue.
Considering metasurfaces in a complete optical system is inconceivable if
common time-consuming numerical methods are used. The macroscopic scale of
elements and the large number of iterations required to design a suitable
system [24] would require too much computing power. Therefore, a method to
efficiently represent metasurfaces in a ray-tracing process is needed. Such a
method could originate from modelling the optical properties of meta-atoms.

Thus, we propose an approach for using metasurfaces in ray-tracing
applications thanks to the efficient characterization of their optical
properties provided by a model. Section 2 describes the ray-tracing equations
adapted for propagation through metasurfaces. Such equations will later be
implemented in ray-tracing software like Zemax${}^{\textrm{\textregistered}}$.
Then, in section 3, a semi-analytical model computing the optical properties
of a metasurface composed of dielectric nanocylinders is given. This modelling
is a continuation of our previous work, where rectangular structures were
studied [2, 6]. Results obtained by this model are then compared to FDTD
simulations in section 4. Finally, we present optical designs using different
refractive lenses and metasurfaces targeting the same application to display
the method’s versatility in section 5.

##  2 Ray-tracing through a metasurface

As the ray-tracing method is based upon calculating the angle of refraction at
each optical surface, it is natural to start our implementation with the
generalized law of refraction for metasurfaces proposed by Yu et al. [28].
Here is the case for the general refraction

$\displaystyle
n_{t}\sin\theta_{t}-n_{i}\sin\theta_{i}=\frac{\lambda}{2\pi}\partialderivative{\varphi}{r},$
(1)

where the metalens implements a phase $\varphi$ on the wavefront, and $n_{i}$,
$\theta_{i}$ and $n_{t}$, $\theta_{t}$ are the refractive index and the
incidence angle before and after the metasurface.

In optical design software such as Zemax${}^{\textrm{\textregistered}}$ or
CodeV${}^{\textrm{\textregistered}}$, rays are represented by a three-
dimensional vector
$\mathbf{R}=x\hat{\mathbf{x}}+z\hat{\mathbf{z}}+z\hat{\mathbf{z}}$. The
orientation of such a ray is defined by its direction cosines $l$, $m$ and
$n$.

$\displaystyle l=\cos{u_{x}},\quad m=\cos{u_{y}},\quad n=\cos{u_{z}},$ (2)

where $u_{x}$, $u_{y}$ et $u_{z}$ are the angles formed by the vector
$\mathbf{R}$ and the axes $\hat{\mathbf{x}}$, $\hat{\mathbf{y}}$ and
$\hat{\mathbf{z}}$. Considering metasurfaces are flat components placed
perpendicular to the optical axis in most cases, we can establish that
$\theta_{i}=90^{\circ}-u\implies\sin{\theta_{i}}=\cos{u}$. Thus, by
rearranging Eq. 1 with $n_{t}=n_{i}=1$ using the direction cosines, the
following set of equations for the coordinates
$(l^{\prime},m^{\prime},n^{\prime})$ of the refracted beam can be obtained

$\displaystyle l^{\prime}$
$\displaystyle=l+\frac{\lambda}{2\pi}\partialderivative{\varphi(x,y)}{x},$ (3)
$\displaystyle m^{\prime}$
$\displaystyle=m+\frac{\lambda}{2\pi}\partialderivative{\varphi(x,y)}{y},$ (4)
$\displaystyle n^{\prime}$
$\displaystyle=\sqrt{1-(l^{\prime})^{2}-(m^{\prime})^{2}}.$ (5)

The normalization condition $l^{2}+m^{2}+n^{2}=1$ was used for the
$n^{\prime}$-coordinate. If one wishes to perform ray tracing through a
metasurface, one must know the phase shift induced for a given incident beam

##  3 Semi-analytical model

As stated previously, obtaining the phase shift implemented by the metasurface
via conventional full-wave simulations is inconceivable for lens design
purposes considering the standard optimization process. Therefore, a model
that can efficiently output the optical properties used for ray-tracing is
advisable. Therefore, this section concerns the development of an expression
for the phase shift of cylindrical nanostructures, an essential parameter for
performing the ray tracing as previously determined.

Related to our earlier work about nanofins [2], the guided propagation through
a nanocylindrical waveguide can be described as a phase accumulation of the
electrical field given by the optical path of an equivalent wave in a medium
of length $H$ as a first-order approximation.

$\displaystyle\varphi=\arg\quantity{\sum_{j}\tau_{j}\exp(-i\beta_{j}H)},$ (6)

where $\tau_{j}$ and $\beta_{j}$ are respectively the complex amplitude and
the propagation constant of the waveguide mode $j$. For the needs of our
application, we consider the fundamental mode, leaving only to obtain the
complex amplitude and the propagation constant of the latter.

###  3.1 Propagation constant

The propagation constant $\beta$ of a single nanocylindrical structure can be
obtained by solving the following dispersion relation [27, 5]:

$\displaystyle\left(\frac{\beta
n}{a}\right)^{2}\left(\frac{1}{p_{1}^{2}}+\frac{1}{q_{2}^{2}}\right)^{2}$
$\displaystyle=\left[\frac{J^{\prime}_{n}(p_{1}a)}{p_{1}J_{n}(p_{1}a)}+\frac{K^{\prime}_{n}(q_{2}a)}{q_{2}K_{n}(q_{2}a)}\right]$
$\displaystyle\qquad\cdot\left[\frac{k_{1}^{2}J^{\prime}_{n}(p_{1}a)}{p_{1}J_{n}(p_{1}a)}+\frac{k_{2}^{2}K_{n}^{\prime}(q_{2}a)}{q_{2}K_{n}(q_{2}a)}\right]$
(7)

with $p_{1}^{2}=\omega^{2}\mu_{0}\epsilon_{1}-\beta^{2}$,
$q_{2}^{2}=\beta^{2}-\omega^{2}\mu_{0}\epsilon_{2}$. $\omega$, $\mu_{0}$,
$\epsilon_{1}$ and $\epsilon_{2}$ are the common physical constants. $J_{n}$
and $K_{n}$ are the Bessel function of the first kind and the modified Bessel
function of the second kind order $n$ respectfully with prime as the
derivative of those function.

From the propagation constant of an isolated structure, it is necessary to add
a correction to consider the coupling effects, knowing that a metasurface is
composed of an array of equidistant nanostructures. To do so, Coupled Mode
Theory (CMT) is often used because of its simplicity, but multiple
approximations should be met for its computation to be accurate [10]. To
consider the periodicity, we will apply a correction to the propagation
constant following the Minot method [19]. Using the slowly-varying envelope
and assuming a coupling $C$ which is much smaller than the propagation
constant, a long-established approximation of the effective propagation
constant of an array of waveguides is expressed as

$\displaystyle\beta_{\text{array}}=\beta_{\text{isolated}}+\frac{4\pi}{\lambda}C\cos(k_{x}S),$
(8)

where $S$ the distance between waveguides, $k_{x}$ is the transverse Floquet-
Bloch waves and $C$ the coupling. Here, we will assume that
$k_{x}=\frac{p+1}{N+1}\pi$ for the $p$ structure within an $N$-size array
[19]. The coupling $C$ can be computed following

$\displaystyle
C=\frac{k_{0}^{2}}{2\beta_{\text{isolated}}}\left(n^{2}-1\right)\frac{\langle
001\rangle}{\langle 00\rangle},$ (9)

with $\langle 00\rangle$ and $\langle 001\rangle$ expressed as

$\displaystyle\langle m\ m^{\prime}\rangle$
$\displaystyle=\iint\limits_{\text{whole array}}M_{m}^{*}M_{m^{\prime}}dA,$
(10) $\displaystyle\langle m\ m^{\prime}\ n\rangle$
$\displaystyle=\iint\limits_{\text{n section
cylinder}}M_{m}^{*}M_{m^{\prime}}dA,$ (11)

where $M_{p}$ correspond to the transverse field distribution of the $p$
structure of the array. For our implementation, $N=2$ and $p=0$ are used,
corresponding to the correction of the closest neighbour. For more details on
the calculations applied for cylinder structures, see appendix A.

###  3.2 Complex amplitude

Using the Snell-Descartes law, one obtains the transmission angles at the
first and second sides of the nanocylinder. Because of the high contrast of
the refractive indices often used in dielectric metasurfaces [3, 13], multiple
scattering events might occur. To take this phenomenon into account, we
consider the Fabry-Perot effect of the nanostructure by introducing a simple
geometric sum formulated as

$\displaystyle\tau=\frac{t_{1}t_{2}}{1-r_{1}r_{2}e^{-i2\beta H}},$ (12)

where $t_{1,2}$, $r_{1,2}$ correspond to the transmission and reflection
coefficients at the first or second interface of the nanofin, and $\beta$ is
the propagation constant inside of the nanocylinder of height $H$.

Finally, from Eq. 6, the output phase-shift $\varphi$ produced by a nanorod of
radius $a$ considering only the fundamental mode is given by

$\displaystyle\varphi(a)=\arg\quantity{\frac{t_{1}t_{2}e^{-i\beta_{\text{array}}(a)H}}{1-r_{1}r_{2}e^{-i2\beta_{\text{array}}(a)H}}}.$
(13)

##  4 Meta-atoms properties

The semi-analytical model allows for computing the phase shift introduced by
those nanostructures under various illuminations. Fig. 1(a) presents the phase
shift $\varphi$ for 600 nm high TiO2 nanocylinders on a SiO2 substrate for
wavelengths $\lambda$ from 568 nm to 628 nm. The structures at a varying
radius between 45 and 110 nm are located 350 nm apart. These parameters were
used since they correspond to metasurfaces often represented in the literature
[18, 7, 22]. As one can observe, the phase can be accurately modelled for
various wavelengths when compared to FDTD simulations using Lumerical software
[25]. We can also extend our analysis to illuminations at various incidence
angles, as shown in Fig. 1(b) where the incident beam of light comes at a
$20^{\circ}$ field of view.

(a) On-axis illumination

(b) Illumination at $20^{\circ}$ field of view

Figure 1: Output phase shift at various wavelengths for our model and
Lumerical simulations. TiO2 meta-atoms are 600 nm high on a SiO2 substrate
with a spacing of 350 nm between structures (n${}_{\text{TiO}_{2}}$=2.6139 and
nsub=1.4584 at $\lambda=588$ nm) with TM plane wave illumination beam.

To assess the magnitude of the error in our predictions relative to numerical
simulations and the impact this would have in a lens design context, we refer
to the work of Patoux et al. [21]. Through their research regarding
nanofabrication errors and their effect on metasurfaces performances, they
concluded that the fabrication error must be limited to 10 nm for a metalens
to maintain acceptable quality. We have identified that the deviations
observed in Fig. 1 are equivalent to a variation of 1-4 nm in radius on
average which is well below the specified limit of 10 nm. Furthermore, a
numerical study on the variation of the Strehl ratio showed that errors
similar to those shown in Fig. 1 still lead to a metalens with a Strehl ratio
above 0.9. Such results confirm that the accuracy of our model compared to
FDTD simulations is sufficient to represent metasurfaces with performance
limited by nanofabrication errors. Therefore, our approach is relevant to
computing the properties of the nanostructures and shows the interest of a
semi-analytical model to perform a quick initial and optimized design.

##  5 Optical design with metasurfaces

To effectively integrate metasurfaces into the design process, it is of utmost
importance that the designer can modify and optimize its characteristics just
like any other optical component of the system. Therefore, the model’s
implementation in the design process was carried out with this objective as a
priority. In this section, a design methodology is proposed and then applied
to present "meta-refractive" hybrid optical systems as a demonstration.

###  5.1 Design process

Figure 2: Flowchart of the ray-tracing implementation in Zemax using a custom
DLL to compute the nanocylinder’s distribution, the output wavefront and the
incident rays’ refraction angle.

The optical design software Zemax has been selected for our first proof of
concept. Its ability to process user-defined surfaces via a dynamic link
library (DLL) is an adequate solution to handle calculations related to the
semi-analytical model. Its operating principle is shown schematically in Fig.
2. First, Zemax initiates rays and propagates them until the metasurface is
reached. Once intercepted, the DLL starts by defining the metasurface from the
characteristics of the meta-atoms and the target phase shift profile
(subsection 5.1.1). The model then computes the wavefronts so that the
refraction angles of the rays are given as output by the DLL (subsection
5.1.2). Exiting the DLL, traditional ray-tracing is resumed across the rest of
the optical system. Ultimately, the lens designer can use performance analysis
tools to evaluate the system (subsection 5.1.3).

####  5.1.1 Defining the metasurface

Directly from the Zemax GUI, a lens designer can choose suitable materials and
dimensions for the nanostructures. The first step of the design process starts
by defining the physical characteristics of the structures to be used:
material, dimensions and spacing between structures. For a given height, it is
necessary to define an interval of radius $(a_{min},a_{max})$ allowing an
adequate phase covering, i.e.
$\absolutevalue{\varphi(a_{max})-\varphi(a_{min})}\geq 2\pi$. A lens designer
can find such information with simple preliminary analyzes using the model
corroborated by numerical simulations.

The design phase shift can be expressed in two ways. First, if one desires to
produce metalenses specifically designed for focusing, a hyperbolic phase
profile is best suited (Eq. 14). Such phase shift leads to a perfect focal
spot on the optical axis at a distance $f$ for a plane incident wavefront.
However, for a metasurface located between multiple elements, a polynomial
phase profile offers much-needed additional degrees of freedom. The phase
relation should be defined as Eq. 15.

$\displaystyle\varphi_{1}(r)$
$\displaystyle=\frac{2\pi}{\lambda}(f-\sqrt{r^{2}+f^{2}}),$ (14)
$\displaystyle\varphi_{2}(r)$
$\displaystyle=\sum_{i}a_{i}\quantity(\frac{r}{R})^{2i}.$ (15)

####  5.1.2 Computing the output wavefront

Once the metasurface target phase shift is defined, the model is then put to
work. Starting from the design phase profile
$\varphi(\mathbf{r},\lambda_{d})$, the nanorods’ distribution
$a(\mathbf{r},\lambda_{d})$ that reproduces exactly $\varphi(\mathbf{r})$ at
the design wavelength $\lambda_{d}$ is found. To do so, the phase relation
developed in section 3 should be inverted in order to obtain a distribution of
propagation constant as $\beta=\varphi(r)/H$. Since the only unknown of this
relation is now the radius $a$ of the nanocylinder, the radius distribution
needed to obtain a target phase distribution is computed.

For the sake of this numerical application, a continuous distribution of
structures of varying dimensions is assumed (see inset of Fig. 2). From the
spatial distribution of the structures, the model is again used to determine
the resulting wavefront for each wavelength and field of view. Finally, the
ray-tracing equations developed previously (Eq. 3 to 5) define the angles of
refraction of the incident beam from the phase gradient of the metasurface.
Finally, the ray’s updated information is given back to Zemax so that ray
tracing can continue throughout the system.

####  5.1.3 Coming up with a design

One of the advantages of our approach is to be able to exploit directly the
optimization algorithms included in the design software. Thus, it is possible
to optimize for the most efficient optical system at several wavelengths and
fields of view, with a proper merit function, while considering the unique
behaviour of metasurfaces. Moreover, a lens designer can study how a slight
variation in the dimensions of the structures can impact the performance and
may enhance the system.

The optimization completed and the performance deemed sufficient, one might
want to perform a confirmation with full-wave simulations. The wavefront
obtained by simulations, possibly slightly different from the one predicted by
the model, can be substituted in Zemax${}^{\textrm{\textregistered}}$. If
necessary, fine adjustments can be made to downstream elements to compensate
for these differences and regain the target optical performances. This way,
optimizing performances using the model in the first place provides an
adequately accurate starting point that only requires realistic adjustments
when moving to numerical FDTD results.

Ultimately, the resulting metasurface is fabricated using conventional
nanofabrication methods. Specialized software such as MetaOptics [8] makes it
possible to obtain a technical drawing of the distribution of nanostructures
to be designed.

###  5.2 Hybrid refractive and metalens designs

To demonstrate the efficiency and versatility of our method for obtaining
metasurface-based optical systems, we now present several variations of a
triplet composed of two refractive elements and a metasurface (see Fig. 3).
The decision to design a triplet was taken considering that this type of
system is the simplest that can correct all of the seven Seidel aberrations
over a wide field of view [15].

According to our process, the first step is defining the structures used in
the design. For continuity purposes, we have chosen the nanocylinders
presented in section 4 for the metasurface, knowing that these respect the
phase coverage condition and that the model agrees with numerical FDTD
simulations. We then arranged the metasurface as either the last (Fig. 3a),
second (Fig. 3b) or first (Fig. 3c) element among a combination of two lenses.
The arbitrary goal being to design an imaging system with an $f/\\#=6$
suitable for a $\pm 10^{\circ}$ field of view (FoV), we optimized the
metasurface’s phase shift profile in addition to the characteristics of the
other lenses to obtain a system with almost diffraction-limited performance.
For reasons mentioned earlier, we chose a polynomial phase profile over its
hyperbolic counterpart to achieve this task.

Figure 3: Variations of a $f/\\#=6$ triplet using two refractive lenses and a
metasurface ($\lambda=588$ nm, FoV = $\pm 10^{\circ}$). The metasurface is
either the last (a), the second (b) or the first element (c). For each
configuration, the spot diagrams and the modulation transfer function (MTF)
were chosen as a criterion for performance evaluation.

All three systems display good optical performances, a fact suggested by the
focal spots of a smaller size than 10 microns. The Airy spot being 8.6 microns
in diameter for these systems means that almost diffraction-limited
performance is achieved for the whole field of view. An analysis of the
modulation transfer function (MTF) supports these conclusions. Indeed, high
contrast is achieved at resolutions suitable for most imaging requirements. As
mentioned previously, we can expect a slight variation of the predicted phase
compared to the FDTD simulations generating a degradation in the optical
quality of the metasurface. However, we have determined that for errors listed
in Fig. 2, the greatest degradation on the RMS spot size was 0.1 $\mu$m for
system 3a), 0.2 $\mu$m for 3b) and 0.4 $\mu$m for 3c). Such variations do not
represent a considerable loss of optical quality, highlighting the adequate
accuracy of our model.

The optimization process made it possible to adapt the metasurface, just like
conventional optics, to the system and thus lead to similar overall
performance regardless of its positioning. It demonstrates one of the main
advantages of our method: considering metasurfaces as a standard optical
element in the optimization process while obtaining an accurate performance
evaluation.

##  6 Conclusion

To conclude, we present a new method to implement metasurfaces in the lens
design process via a semi-analytical model used in conjunction with ray-
tracing software. We had to adapt the ray tracing equations to consider a
phase shift induced by metasurfaces in the first place. Instead of going
through lengthy numerical simulations to predict the phase shift, the latter
is provided by a semi-analytical model that we have developed. Based on an
effective medium theory, the Fabry-Perot effect and a closed-form
approximation of the coupled-wave theory, our model predicts the phase shift
introduced by cylindrical nanostructures. We then proceeded to its
implementation in the optical design software Zemax through a DLL. Starting
from a target phase shift for the metasurface, the model integrated into the
DLL renders the distribution of nanostructures and calculates the output
wavefront so that a ray of a given wavelength and orientation can continue
through the system. It led to the development of a 3-step feedback process
exploiting the ability to include the metasurface in the optimization routine
of an optical system. For demonstration purposes, we have designed hybrid
refractive and metasurface systems composed of two refractive elements and a
metasurface. The latter, being located at distinct positions for each system,
was able to adapt just like conventional elements to render a system with
near-diffraction-limited performance across the entire field of view.

Further developments in the generalization of the method to other types of
meta-atoms are to be expected. For example, using intrinsically achromatic
structures or efficient structures with a large field of view would make it
possible to exploit the full potential of metasurfaces in optical design.
Ultimately, experimental work could corroborate the accuracy of the semi-
analytical model and the relevance of our approach.

## APPENDIX A : CMT for nanocylinders

The sum of mode for the radial component inside (1) and outside (2) of the
waveguide can be expressed as [27]

$\displaystyle E_{r}^{(1)}(r,\theta)$
$\displaystyle=\frac{1}{p_{1}^{2}}\sum\left[-i\beta
p_{1}A_{n}J_{n}^{\prime}(p_{1}r)+\frac{\omega\mu_{0}n}{r}B_{n}J_{n}(p_{1}r)\right]e^{in\theta}$
(16) $\displaystyle E_{r}^{(2)}(r,\theta)$
$\displaystyle=-\frac{1}{q_{2}^{2}}\sum\left[-i\beta
q_{2}C_{n}K^{\prime}_{n}(q_{2}r)+\frac{\omega\mu_{0}n}{r}D_{n}K_{n}(q_{2}r)\right]e^{in\theta},$
(17)

where $A_{n},B_{n},C_{n}$ and $D_{n}$ are arbitrary constants. However, to be
solutions of the dispersion relation of the waveguide, the constants should
follow these relations

$\displaystyle B_{n}$
$\displaystyle=-\frac{-\frac{i\omega\epsilon_{1}J_{n}^{\prime}(p_{1}a)}{p_{1}}+\frac{i\omega\epsilon_{2}K_{n}^{\prime}(q_{2}a)}{q_{2}}\cdot\frac{J_{n}(p_{1}a)}{-K_{n}(q_{2}a)}}{\frac{\beta
J_{n}(p_{1}a)}{ap_{1}^{2}}-\frac{\beta
K_{n}(q_{2}a)}{aq_{2}^{2}}\cdot\frac{J_{n}(p_{1}a)}{-K_{n}(q_{2}a)}}A_{n}$
(18)

or

$\displaystyle B_{n}$ $\displaystyle=-\frac{\frac{\beta
J_{n}(p_{1}a)}{ap_{1}^{2}}-\frac{\beta
K_{n}(q_{2}a)}{aq_{2}^{2}}\cdot\frac{J_{n}(p_{1}a)}{-K_{n}(q_{2}a)}}{\frac{i\omega\mu_{0}J_{n}^{\prime}(p_{1}a)}{p1}-\frac{i\omega\mu_{0}K_{n}^{\prime}(q_{2}a)}{q_{2}}\cdot\frac{J_{n}(p_{1}a)}{-K_{n}(q_{2}a)}}A_{n}$
(19) $\displaystyle C_{n}$
$\displaystyle=\frac{J_{n}(p_{1}a)}{K_{n}(q_{2}a)}A_{n}$ (20) $\displaystyle
D_{n}$ $\displaystyle=\frac{J_{n}(p_{1}a)}{K_{n}(q_{2}a)}B_{n}$ (21)

where $a$ is the radius of the nanocylinder under consideration. If we only
considering the fundamental mode $HE_{11}$, the $\theta$ dependency is not
present. Thus, the surface integral in one dimension can be computed with the
following equations, where $[f(x)]^{*}$ is the complex conjugate of the
function $f(x)$.

$\displaystyle\langle 00\rangle$
$\displaystyle=\int_{-\infty}^{-a}\left[E_{r}(|x|)^{(2)}\right]^{*}E_{r}(|x|)^{(2)}dx$
$\displaystyle\qquad+\int_{-a}^{a}\left[E_{r}(|x|)^{(1)}\right]^{*}E_{r}(|x|)^{(1)}dx$
$\displaystyle\qquad\qquad\cdot\int_{a}^{\infty}\left[E_{r}(|x|)^{(2)}\right]^{*}E_{r}(|x|)^{(2)}dx$
(22) $\displaystyle\langle 001\rangle$
$\displaystyle=\int_{-a}^{a}\left[E_{r}(|x|)^{(1)}\right]^{*}E_{r}(|x\pm
2\pi|)^{(2)}dx$ (23)

##  7 Funding Information

Natural Sciences and Engineering Research Council of Canada
(NSERC)(RGPIN-2016-05962).

Disclosures. The authors declare no conflicts of interest.

## References

  * [1] K Achouri and C Caloz.  Design, concepts, and applications of electromagnetic metasurfaces.  Nanophotonics, 7(6):1095–1116, 2018.
  * [2] Jeck Borne, Denis Panneton, Michel Piché, and Simon Thibault.  Semi-analytical model of the optical properties of a metasurface composed of nanofins.  Journal of the Optical Society of America B, 38(11):3509, 2021.
  * [3] Shengyuan Chang, Xuexue Guo, and Xingjie Ni.  Optical metasurfaces: Progress and applications.  Annual Review of Materials Research, 48(1):279–302, 2018.
  * [4] Hou-Tong Chen, Antoinette J Taylor, and Nanfang Yu.  A review of metasurfaces: physics and applications.  Reports on Progress in Physics, 79(7):076401, 2016.
  * [5] D K Cheng.  Field and Wave Electromagnetics.  Addison-Wesley series in electrical engineering. Addison-Wesley, 1989\.
  * [6] Alexandre Cléroux Cuillerier, Jeck Borne, Xavier Dallaire, and Simon Thibault.  Toward hybrid refractive and metalens design.  In OSA Optical Design and Fabrication 2021 (Flat Optics, Freeform, IODC, OFT), page 120781W. Optica Publishing Group, 2021.
  * [7] Manuel Decker, Wei Ting Chen, Thomas Nobis, Alexander Y Zhu, Mohammadreza Khorasaninejad, Zameer Bharwani, Federico Capasso, and Jörg Petschulat.  Imaging Performance of Polarization-Insensitive Metalenses.  ACS Photonics, 6(6), 2019.
  * [8] Raghu Dharmavarapu, Soon Hock Ng, Fatima Eftekhari, Saulius Juodkazis, and Shanti Bhattacharya.  Metaoptics: opensource software for designing metasurface optical element gdsii layouts.  Opt. Express, 28(3):3505–3516, Feb 2020.
  * [9] Shaohua Dong, Guangwei Hu, Qiang Wang, Yuxiang Jia, Qing Zhang, Guangtao Cao, Jiafu Wang, Shuqing Chen, Dianyuan Fan, Weixiang Jiang, Ying Li, Andrea Alù, and Cheng Wei Qiu.  Loss-Assisted Metasurface at an Exceptional Point.  ACS Photonics, 7(12):3321–3327, 2020.
  * [10] Leonard Eyges and Peter Wintersteiner.  Modes of an Array of Dielectric Waveguides.  Journal of the Optical Society of America, 71(11):1351–1360, 1981\.
  * [11] Aaron L. Holsteen, Ahmet Fatih Cihan, and Mark L. Brongersma.  Temporal color mixing and dynamic beam shaping with silicon metasurfaces.  Science, 365(6450):257–260, 2019.
  * [12] Jie Hu, Sankhyabrata Bandyopadhyay, Yu-hui Liu, and Li-yang Shao.  A review on metasurface : From principle to smart metadevices.  Frontier in Physics, 8(January):1–20, 2021.
  * [13] SM Kamali, E Arbabi, A Arbabi, and A Faraon.  A review of dielectric optical metasurfaces for wavefront control.  Nanophotonics, 7(6):1041–1068, 2018.
  * [14] Mohammadreza Khorasaninejad and Federico Capasso.  Metalenses: Versatile multifunctional photonic components.  Science, 358(6367), 2017.
  * [15] M. J. Kidger.  Fundamental Optical Design.  SPIE Press, Bellingham, Wash, 2001.
  * [16] Alexander Krasnok, Mykhailo Tymchenko, and Andrea Alù.  Nonlinear metasurfaces: a paradigm shift in nonlinear optics.  Materials Today, 21(1):8–21, 2018.
  * [17] Haowen Liang, Augusto Martins, Ben-Hur V. Borges, Jianying Zhou, Emiliano R. Martins, Juntao Li, and Thomas F. Krauss.  High performance metalenses: numerical aperture, aberrations, chromaticity, and trade-offs.  Optica, 6(12):1461, 2019.
  * [18] Yaoyao Liang, Hongzhan Liu, Faqiang Wang, Hongyun Meng, Jianping Guo, Jinfeng Li, and Zhongchao Wei.  High-efficiency, near-diffraction limited, dielectric metasurface lenses based on crystalline titanium dioxide at visible wavelengths.  Nanomaterials, 8(5), 2018.
  * [19] Christophe Minot, Nadia Belabas, Juan Ariel Levenson, and Jean-Marie Moison.  Analytical first-order extension of coupled-mode theory for waveguide arrays.  Optics Express, 18(7):7157, 2010.
  * [20] Seong-Won Moon, Yeseul Kim, Gwanho Yoon, and Junsuk Rho.  Recent progress on ultrathin metalenses for flat optics.  iScience, 23(12):101877, 2020.
  * [21] Adelin Patoux, G. Agez, Christian Girard, Vincent Paillard, Peter R. Wiecha, Aurélie Lecestre, Franck Carcenac, Guilhem Larrieu, and A. Arbouet.  Challenges in nanofabrication for efficient optical metasurfaces.  Scientific reports, 11 1:5620, 2021.
  * [22] Z. Shen, R. Li, Y. Z. Xue, Z. Y. Qiu, Z. Y. Xiang, J. Y. Zhou, and B. F. Zhang.  Generation of Optical Vortices with Polarization-Insensitive Metasurfaces.  IEEE Photonics Journal, 12(4), 2020.
  * [23] Zhujun Shi, Mohammadreza Khorasaninejad, Yao Wei Huang, Charles Roques-Carmes, Alexander Y. Zhu, Wei Ting Chen, Vyshakh Sanjeev, Zhao Wei Ding, Michele Tamagnone, Kundan Chaudhary, Robert C. Devlin, Cheng Wei Qiu, and Federico Capasso.  Single-Layer Metasurface with Controllable Multiwavelength Functions.  Nano Letters, 18(4):2420–2427, 2018.
  * [24] Warren J Smith.  Modern optical engineering; the design of optical systems.  McGraw-Hill, New York, 1966.
  * [25] Lumerical Inc.  http://www.lumerical.com/tcad-products/fdtd/.
  * [26] Ming Lun Tseng, Hui-hsin Hsiao, Cheng Hung Chu, Mu Ku Chen, Greg Sun, and Ai-qun Liu.  Metalenses : Advances and Applications.  Advanced Optical Materials, 6(18):1–16, 2018.
  * [27] Cavour W. H Yeh and F. I Shimabukuro.  The essence of dielectric waveguides.  Springer, New York, 2008.  Bibliogr.
  * [28] Nanfang Yu, Patrice Genevet, Mikhail A. Kats, Francesco Aieta, Jean-Philippe Tetienne, Federico Capasso, and Zeno Gaburro.  Light propagation with phase discontinuities: Generalized laws of reflection and refraction.  Science, 334(6054):333–337, 2011.