# Neveu’s exchange formula for analysis of
wireless networks with hotspot clusters

Naoto Miyoshi The author is with Department of Mathematical and Computing
Science, Tokyo Institute of Technology, Tokyo 152-8552, Japan. E-mail:
<EMAIL_ADDRESS>author’s work was supported by the Japan Society
for the Promotion of Science (JSPS) Grant-in-Aid for Scientific Research (C)
19K11838.

###### Abstract

Theory of point processes, in particular Palm calculus within the stationary
framework, plays a fundamental role in the analysis of spatial stochastic
models of wireless communication networks. Neveu’s exchange formula, which
connects the respective Palm distributions for two jointly stationary point
processes, is known as one of the most important results in the Palm calculus.
However, its use in the analysis of wireless networks seems to be limited so
far and one reason for this may be that the formula in a well-known form is
based upon the Voronoi tessellation. In this paper, we present an alternative
form of Neveu’s exchange formula, which does not rely on the Voronoi
tessellation but includes the one as a special case. We then demonstrate that
our new form of the exchange formula is useful for the analysis of wireless
networks with hotspot clusters modeled using cluster point processes.
Keywords: Stationary point processes, Palm calculus, Neveu’s exchange formula,
cluster point processes, device-to-device networks, hotspot clusters, coverage
probability, device discovery.

##  1 Introduction

Spatial stochastic models have been widely accepted in the literature as
mathematical models for the analysis of wireless communication networks, where
irregular locations of wireless nodes, such as base stations (BSs) and user
devices, are modeled using spatial point processes on the Euclidean plane
(see, e.g., [7, 8, 24, 23, 36, 12] for monographs and [6, 20, 26, 32] for
recent survey and tutorial articles). In such analysis of wireless networks,
the theory of point processes, in particular Palm calculus within the
stationary framework, plays a fundamental role. Neveu’s exchange formula,
which connects the respective Palm distributions for two jointly stationary
point processes, is known as one of the most important results in the Palm
calculus. However, its use in the analysis of wireless networks seems to be
limited so far and one reason for this may be that the formula in a well-known
form is based upon the Voronoi tessellation (see, e.g, [9, Section 6.3]). In
this paper, we present an alternative form of Neveu’s exchange formula, which
does not rely on the Voronoi tessellation but includes the one as a special
case, and then demonstrate that it is useful for the analysis of spatial
stochastic models based on cluster point processes.

A cluster point process represents a state such that there exist a large
number of clusters consisting of multiple points and is used to model the
locations of wireless nodes in an (urban) area with a number of hotspots.
Indeed, many researchers have adopted the cluster point processes in their
models of various wireless networks such as ad hoc networks [21],
heterogeneous networks [41, 16, 38, 39, 1, 40, 44], device-to-device (D2D)
networks [2], wireless powered networks [14], unmanned aerial vehicle assisted
networks [43], and so on. In this paper, we focus on so-called stationary
Poisson-Poisson cluster processes (PPCPs) (see, e.g., [13, 33]) and apply the
new form of the exchange formula to the analysis of stochastic models based on
them.

We first use the exchange formula for the Palm characterization, where we
derive the intensity measure, the generating functional and the nearest-
neighbor distance distribution for a stationary PPCP under its Palm
distribution. Although these results are known in the literature (see, e.g,
[11, 21]), we here give them simple and unified proofs using the new form of
the exchange formula. We next consider some applications to wireless networks
modeled using stationary PPCPs, where we examine the problems of coverage and
device discovery in a D2D network. The coverage analysis of a D2D network
model based on a cluster point process was considered in [2], where a device
communicates with another device in the same cluster. In contrast to this, we
assume here that a device receives messages from the nearest transmitting
device, which is possibly in a different cluster because clusters may overlap
in space. For this model, we derive the coverage probability using the
exchange formula. On the other hand, in the problem of device discovery,
transmitting devices transmit broadcast messages and a receiving device can
detect the transmitters if it can successfully decode the broadcast messages.
Such a problem was studied in [25, 10, 27] when the devices are located
according to a homogeneous Poisson point process (PPP) and in [28] when the
devices are located according to a Ginibre point process (see, e.g., [34, 35]
for the Ginibre point process and its applications to wireless networks). We
consider the case where the devices are located according to a stationary
PPCPs and derive the expected number of transmitting devices discovered by a
receiving device. We should note that Neveu’s exchange formula is also
introduced in a more general form in [31, 29, 22], so that the form presented
in the paper may be within its scope. Nevertheless, we see in the rest of the
paper that our new form would be valuable and could spread the application
fields of the exchange formula.

The rest of the paper is organized as follows. The new form of Neveu’s
exchange formula is derived in the next section, where the relations with the
existing forms are also discussed. In Section 3, the exchange formula is
applied to the Palm characterization of a stationary PPCP, where alternative
proofs of the intensity measure, the generating functional and the nearest-
neighbor distance distribution under the Palm distribution are given. In
section 4, some applications to wireless network models are examined, where
for a D2D network model based on a stationary PPCP, the coverage probability
and the expected number of discovered devices are derived using the exchange
formula. The results of numerical experiments are also presented there.
Concluding remarks are provided in Section 5.

##  2 Neveu’s exchange formula

In this section, we discuss point processes on the $d$-dimensional Euclidean
space $\mathbb{R}^{d}$ within the stationary framework (see, e.g., [9, Chapter
6] for details on the stationary framework). In what follows,
$\mathcal{B}(\mathbb{R}^{d})$ denotes the Borel $\sigma$-field on
$\mathbb{R}^{d}$ and $\delta_{x}$ denotes the Dirac measure with mass at
$x\in\mathbb{R}^{d}$. Let $(\Omega,\mathcal{F},\mathbb{P})$ denote a
probability space. On $(\Omega,\mathcal{F})$, a flow
$\\{\theta_{t}\\}_{t\in\mathbb{R}^{d}}$ is defined such that $\theta_{t}$:
$\Omega\to\Omega$ is $\mathcal{F}$-measurable and bijective satisfying
$\theta_{t}\circ\theta_{u}=\theta_{t+u}$ for $t,u\in\mathbb{R}^{d}$, where
$\theta_{0}$ is the identity for $0=(0,0,\ldots,0)\in\mathbb{R}^{d}$; so that
$\theta_{t}^{-1}=\theta_{-t}$ for $t\in\mathbb{R}^{d}$. We assume that the
probability measure $\mathbb{P}$ is invariant to the flow
$\\{\theta_{t}\\}_{t\in\mathbb{R}^{d}}$ (in other words,
$\\{\theta_{t}\\}_{t\in\mathbb{R}^{d}}$ preserves $\mathbb{P}$) in the sense
that $\mathbb{P}\circ\theta_{t}^{-1}=\mathbb{P}$ for any $t\in\mathbb{R}^{d}$,
where $\theta_{t}^{-1}(A)=\\{\omega\in\Omega:\theta_{t}(\omega)\in A\\}$ for
$A\in\mathcal{F}$. A point process $\Phi=\sum_{n=1}^{\infty}\delta_{X_{n}}$ on
$\mathbb{R}^{d}$ is said to be compatible with the flow
$\\{\theta_{t}\\}_{t\in\mathbb{R}^{d}}$ if it holds that
$\Phi(B)\circ\theta_{t}=\Phi(\theta_{t}(\omega),B)=\Phi(\omega,B+t)=\Phi(B+t)$
for $\omega\in\Omega$, $B\in\mathcal{B}(\mathbb{R}^{d})$ and
$t\in\mathbb{R}^{d}$, where $B+t=\\{x+t\in\mathbb{R}^{d}:x\in B\\}$; that is,
for $t\in\mathbb{R}^{d}$ and $n\in\mathbb{N}=\\{1,2,\ldots\\}$, there exists
an $n^{\prime}\in\mathbb{N}$ such that
$X_{n}\circ\theta_{t}=X_{n^{\prime}}-t$. Under the assumption of the
$\\{\theta_{t}\\}_{t\in\mathbb{R}^{d}}$-invariance of $\mathbb{P}$, a point
process $\Phi$ compatible with $\\{\theta_{t}\\}_{t\in\mathbb{R}^{d}}$ is
stationary in $\mathbb{P}$ and furthermore, two point processes $\Phi$ and
$\Psi$, both of which are compatible with
$\\{\theta_{t}\\}_{t\in\mathbb{R}^{d}}$, are jointly stationary in
$\mathbb{P}$.

Let $\Phi=\sum_{n=1}^{\infty}\delta_{X_{n}}$ and
$\Psi=\sum_{m=1}^{\infty}\delta_{Y_{m}}$ denote point processes on
$\mathbb{R}^{d}$, which are both simple, compatible with
$\\{\theta_{t}\\}_{t\in\mathbb{R}^{d}}$ and have positive and finite
intensities $\lambda_{\Phi}$ and $\lambda_{\Psi}$, respectively. Thus, $\Phi$
and $\Psi$ are jointly stationary in probability $\mathbb{P}$ and the
respective Palm probabilities $\mathbb{P}_{\Phi}^{0}$ and
$\mathbb{P}_{\Psi}^{0}$ are well-defined. Note that
$\mathbb{P}_{\Phi}^{0}(\Phi(\\{0\\})=1)=\mathbb{P}_{\Psi}^{0}(\Psi(\\{0\\})=1)=1$.
In this paper, when we consider the event
$\\{\Phi(\\{0\\})=1\\}\in\mathcal{F}$, we assign index $0$ to the point at the
origin; that is, $X_{0}=0$ on $\\{\Phi(\\{0\\})=1\\}$, and this is also the
case for $\Psi$; that is, $Y_{0}=0$ on $\\{\Psi(\\{0\\})=1\\}$. To present an
alternative form of Neveu’s exchange formula, we introduce a family of shift
operators $S_{t}$, $t\in\mathbb{R}^{d}$, on the set of measures $\eta$ on
$(\mathbb{R}^{d},\mathcal{B}(\mathbb{R}^{d}))$ by $S_{t}\eta(B)=\eta(B+t)$ for
$B\in\mathcal{B}(\mathbb{R}^{d})$. For example, operating $S_{t}$ on the point
process $\Psi=\sum_{m=1}^{\infty}\delta_{Y_{m}}$, we have
$S_{t}\Psi=\sum_{m=1}^{\infty}\delta_{Y_{m}-t}=\Psi\circ\theta_{t}$. The shift
operators $S_{t}$, $t\in\mathbb{R}^{d}$, also work on a function $h$ on
$\mathbb{R}^{d}$ such as $S_{t}h(x)=h(x+t)$ for $x\in\mathbb{R}^{d}$.

###### Theorem 1

For the two jointly stationary point processes
$\Phi=\sum_{n=1}^{\infty}\delta_{X_{n}}$ and
$\Psi=\sum_{m=1}^{\infty}\delta_{Y_{m}}$ described above, we assume that a
family of point processes $\Psi_{n}=\sum_{k=1}^{\kappa_{n}}\delta_{Y_{n,k}}$,
$n\in\mathbb{N}$, can be constructed such that

  1. 1.

$S_{-X_{n}}\Psi_{n}=\sum_{k=1}^{\kappa_{n}}\delta_{X_{n}+Y_{n,k}}$,
$n\in\mathbb{N}$, form a partition of $\Psi$; that is,
$\Psi=\sum_{n=1}^{\infty}S_{-X_{n}}\Psi_{n}$.

  2. 2.

$\widetilde{\Phi}=\sum_{n=1}^{\infty}\delta_{(X_{n},\Psi_{n})}$ is a
stationary marked point process with the set of counting measures on
$\mathbb{R}^{d}$ as its mark space.

Then, for any nonnegative random variable $W$ defined on
$(\Omega,\mathcal{F})$,

$\lambda_{\Psi}\,\mathbb{E}_{\Psi}^{0}[W]=\lambda_{\Phi}\,\mathbb{E}_{\Phi}^{0}\biggl{[}\int_{\mathbb{R}^{d}}W\circ\theta_{y}\,\Psi_{0}(\mathrm{d}y)\biggr{]}=\lambda_{\Phi}\,\mathbb{E}_{\Phi}^{0}\Biggl{[}\sum_{k=1}^{\kappa_{0}}W\circ\theta_{Y_{0,k}}\Biggr{]},$
(1)

where $\mathbb{E}_{\Phi}^{0}$ and $\mathbb{E}_{\Psi}^{0}$ denote the
expectations with respect to the Palm probabilities $\mathbb{P}_{\Phi}^{0}$
and $\mathbb{P}_{\Psi}^{0}$, respectively, and
$\Psi_{0}=\sum_{k=1}^{\kappa_{0}}\delta_{Y_{0,k}}$ denotes the mark associated
with the point $X_{0}=0$ on $\\{\Phi(\\{0\\})=1\\}$.

  * Proof:

As with the proof of the exchange formula in [9, Theorem 6.3.7], we start our
proof with the mass transport formula (see, e.g., [9, Theorem 6.1.34]); that
is, for any measurable function $\xi$:
$\Omega\times\mathbb{R}^{d}\to\overline{\mathbb{R}}_{+}$,

$\lambda_{\Phi}\,\mathbb{E}_{\Phi}^{0}\biggl{[}\int_{\mathbb{R}^{d}}\xi(y)\,\Psi(\mathrm{d}y)\biggr{]}=\lambda_{\Psi}\,\mathbb{E}_{\Psi}^{0}\biggl{[}\int_{\mathbb{R}^{d}}\xi(-x)\circ\theta_{x}\,\Phi(\mathrm{d}x)\biggr{]}.$
(2)

Let $\xi(y)=W\circ\theta_{y}\,\Psi_{0}(\\{y\\})$ on $\\{\Phi(\\{0\\})=1\\}$.
Then, the left-hand side of (2) becomes

$\lambda_{\Phi}\,\mathbb{E}_{\Phi}^{0}\biggl{[}\int_{\mathbb{R}^{d}}W\circ\theta_{y}\,\Psi_{0}(\\{y\\})\,\Psi(\mathrm{d}y)\biggr{]}=\lambda_{\Phi}\,\mathbb{E}_{\Phi}^{0}\biggl{[}\int_{\mathbb{R}^{d}}W\circ\theta_{y}\,\Psi_{0}(\mathrm{d}y)\biggr{]}.$

On the other hand, the right-hand side of (2) is reduced to

$\lambda_{\Psi}\,\mathbb{E}_{\Psi}^{0}\biggl{[}\int_{\mathbb{R}^{d}}\bigl{(}W\circ\theta_{-x}\,\Psi_{0}(\\{-x\\})\bigr{)}\circ\theta_{x}\,\Phi(\mathrm{d}x)\biggr{]}=\lambda_{\Psi}\,\mathbb{E}_{\Psi}^{0}\biggl{[}W\sum_{n=1}^{\infty}S_{-X_{n}}\Psi_{n}(\\{0\\})\biggr{]}=\lambda_{\Psi}\,\mathbb{E}_{\Psi}^{0}[W],$

where the first equality follows from $\Psi_{0}\circ\theta_{X_{n}}=\Psi_{n}$
for $n\in\mathbb{N}$ and the fact that the point of $\Psi$ at the origin on a
sample $\omega\in\\{\Psi(\\{0\\})=1\\}$ is shifted to location $-x$ on the
shifted sample $\theta_{x}(\omega)$ for $x\in\mathbb{R}^{d}$, and the last
equality follows since there exists exactly one point, say $X_{n}$, of $\Phi$
such that its mark $\Psi_{n}=\sum_{k=1}^{\kappa_{n}}\delta_{Y_{n,k}}$ has a
point, say $Y_{n,k}$, satisfying $X_{n}+Y_{n,k}=0$ on $\\{\Psi(\\{0\\})=1\\}$.
The proof is completed.

###### Remark 1

Let $W\equiv 1$ in (1). Then, we have
$\mathbb{E}_{\Phi}^{0}[\kappa_{0}]=\lambda_{\Psi}/\lambda_{\Phi}$ and
therefore, each $\Psi_{n}$ in Theorem 1 has finite points. In the form of the
exchange formula in [9, Theorems 6.3.7 & 6.3.19], the point process $\Psi$ is
partitioned by the Voronoi tessellation for $\Phi$, which corresponds to a
special case of (1) such that $S_{-X_{n}}\Psi_{n}(\cdot)=\Psi(\cdot\cap
V_{\Phi}(X_{n}))$ for $n\in\mathbb{N}$, where $V_{\Phi}(X_{n})$ denotes the
Voronoi cell of point $X_{n}$ of $\Phi$. The condition in [9] such that there
are no points of $\Psi$ on the boundaries of Voronoi cells $V_{\Phi}(X_{n})$,
$n\in\mathbb{N}$, is covered by our Condition 1 in Theorem 1, where
$S_{-X_{n}}\Psi_{n}$, $n\in\mathbb{N}$, form a partition of $\Psi$ and have no
common points. On the other hand, our Theorem 1 considers only the case where
the point process $\Phi$ is simple unlike [9, Theorem 6.3.7]. However, this
would be enough for applications to wireless networks and it could be extended
to the non-simple case. It should also be noted that, in [31, 29, 22], a more
general formula is introduced under the name of Neveu’s exchange formula, from
which the mass transport formula (2) is derived. In that sense, our form (1)
may be within its scope. Nevertheless, we can see in the following sections
that Theorem 1 is valuable in the sense that it is tractable and can spread
the application fields of the exchange formula.

##  3 Applications to cluster point processes

In this section, we demonstrate that Neveu’s exchange formula (1) in Theorem 1
is useful to characterize the Palm distribution of stationary cluster point
processes. A cluster point process is, roughly speaking, constructed by
placing point processes (usually with finite points), called offspring
processes, around respective points of another point process, called a parent
process, and is used to represent a state such that there exist a large number
of clusters consisting of multiple points (see Figure 1). In particular, we
focus here on a stationary PPCP described next.

Figure 1: A sample of a $2$-dimensional cluster point process ($\bullet$
represents the points of the cluster point process and $\mathop{\times\hbox
to0.0pt{\kern-7.47925pt$+$\hss\hss}}$ represents the points of the parent
process)

###  3.1 Poisson-Poisson cluster processes

Let $\Phi=\sum_{n=1}^{\infty}\delta_{X_{n}}$ denote a homogeneous PPP on
$\mathbb{R}^{d}$, which works as the parent process, and let
$\Psi_{n}=\sum_{k=1}^{\kappa_{n}}\delta_{Y_{n,k}}$, $n\in\mathbb{N}$, denote a
family of finite (therefore inhomogeneous) and mutually independent PPPs on
$\mathbb{R}^{d}$, which are also independent of $\Phi$ and work as the
offspring processes. Then, PPCP $\Psi=\sum_{m=1}^{\infty}\delta_{Y_{m}}$ is
given as

$\Psi=\sum_{n=1}^{\infty}S_{-X_{n}}\Psi_{n}=\sum_{n=1}^{\infty}\sum_{k=1}^{\kappa_{n}}\delta_{X_{n}+Y_{n,k}}.$
(3)

The PPCP $\Psi$ constructed as above is stationary since the parent process
$\Phi$ is stationary and the offspring processes $\Psi_{n}$, $n\in\mathbb{N}$,
are independent and identically distributed (see, e.g., [9, Example 2.3.18]).
We assume that $\Phi$ has a positive and finite intensity $\lambda_{\Phi}$,
and $\Psi_{n}$, $n\in\mathbb{N}$, have an identical intensity measure
$\Lambda_{\mathrm{o}}=\mu\,Q$, where $\mu$ is a positive constant and $Q$ is a
probability distribution on $(\mathbb{R}^{d},\mathcal{B}(\mathbb{R}^{d}))$.
Thus, the number of points in each offspring process follows a Poisson
distribution with mean $\mu$, so that the intensity of $\Psi$ is equal to
$\lambda_{\Psi}=\lambda_{\Phi}\mu$, and offspring points are scattered on
$\mathbb{R}^{d}$ according to $Q$ independently of each other. We further
assume that $Q$ is diffuse; that is, $Q(\\{x\\})=0$ for any
$x\in\mathbb{R}^{d}$, to make $\Psi$ simple. We refer to $S_{-X_{n}}\Psi_{n}$
in (3) as the cluster associated with $X_{n}$ for $n\in\mathbb{N}$. Two main
examples of the PPCPs are the (modified) Thomas point process and the Matérn
cluster process (see, e.g., [15, Example 5.5]). When $Q$ is an isotropic
normal distribution, then the obtained PPCP is called the Thomas point
process. On the other hand, when $Q$ is the uniform distribution on a fixed
ball centered at the origin, then the result is called the Matérn cluster
process. Note that the PPCP $\Psi$ and its parent process $\Phi$ satisfy the
conditions of Theorem 1.

###  3.2 Characterization of Palm distribution

For a stationary point process $\Psi$, let $\Psi^{!}:=\Psi-\delta_{0}$ on the
event $\\{\Psi(\\{0\\})=1\\}$, which is referred to as the reduced Palm
version of $\Psi$.

###### Lemma 1

For the stationary PPCP $\Psi$ described in Section 3.1, the intensity measure
of the reduced Palm version $\Psi^{!}$ (with respect to the Palm distribution)
is given by

$\Lambda_{\Psi}^{0}(B):=\mathbb{E}_{\Psi}^{0}[\Psi^{!}(B)]=\lambda_{\Phi}\,\mu\,|B|+\mu\int_{\mathbb{R}^{d}}Q(B-y)\,Q^{-}(\mathrm{d}y),\quad
B\in\mathcal{B}(\mathbb{R}^{d}),$ (4)

where $|\cdot|$ denotes the Lebesgue measure on
$(\mathbb{R}^{d},\mathcal{B}(\mathbb{R}^{d}))$ and $Q^{-}(B)=Q(-B)$ with
$-B=\\{-x:x\in B\\}$ for $B\in\mathcal{B}(\mathbb{R}^{d})$.

  * Proof:

Since the offspring processes $\Psi_{n}$, $n\in\mathbb{N}$, are PPPs, the PPCP
$\Psi$ is a Cox point process; that is, once the parent process
$\Phi=\sum_{n=1}^{\infty}\delta_{X_{n}}$ is given, $\Psi$ is conditionally an
inhomogeneous PPP with a conditional intensity measure
$\mu\sum_{n=1}^{\infty}S_{-X_{n}}Q$ (see, e.g., [9, Example 2.3.13]). Since
the reduced Palm version of a PPP is identical in distribution to its original
version (not conditioned on $\\{\Psi(\\{0\\})=1\\}$) by Slivnyak’s theorem
(see, e.g., [18, Proposition 13.1.VII] or [9, Theorem 3.2.4]), we have

$\mathbb{E}_{\Psi}^{0}[\Psi^{!}(B)\mid\Phi]=\mathbb{E}[\Psi(B)\mid\Phi]=\mu\sum_{n=1}^{\infty}Q(B-X_{n}),\quad
B\in\mathcal{B}(\mathbb{R}^{d}).$

Taking the expectation with respect to $\mathbb{P}_{\Psi}^{0}$ and then
applying Theorem 1, we obtain

$\displaystyle\mathbb{E}_{\Psi}^{0}[\Psi^{!}(B)]$
$\displaystyle=\mu\,\mathbb{E}_{\Psi}^{0}\Biggl{[}\sum_{n=1}^{\infty}Q(B-X_{n})\Biggr{]}=\mathbb{E}_{\Phi}^{0}\Biggl{[}\int_{\mathbb{R}^{d}}\Biggl{(}\sum_{n=1}^{\infty}Q(B-X_{n})\Biggr{)}\circ\theta_{y}\,\Psi_{0}(\mathrm{d}y)\Biggr{]}$
$\displaystyle=\mathbb{E}_{\Phi}^{0}\Biggl{[}\int_{\mathbb{R}^{d}}\sum_{n=0}^{\infty}Q(B-X_{n}+y)\,\Psi_{0}(\mathrm{d}y)\Biggr{]}$
$\displaystyle=\mu\int_{\mathbb{R}^{d}}\Biggl{(}Q(B+y)+\mathbb{E}_{\Phi}^{0}\Biggl{[}\sum_{n=1}^{\infty}Q(B-X_{n}+y)\Biggr{]}\Biggr{)}\,Q(\mathrm{d}y),$

where $\lambda_{\Psi}=\lambda_{\Phi}\mu$ is used in the second equality, the
third equality follows because, for any $n\in\mathbb{N}$ and
$y\in\mathbb{R}^{d}$, there exists an $n^{\prime}\in\mathbb{N}\cup\\{0\\}$
such that $X_{n}\circ\theta_{y}=X_{n^{\prime}}-y$ on $\\{\Phi(\\{0\\})=1\\}$,
and in the last equality, we apply Campbell’s formula (see, e.g., [30,
Proposition 2.7] or [9, Theorem 1.2.5]) for $\Psi_{0}$. For the expectation in
the last expression above, Slivnyak’s theorem, Campbell’s formula for $\Phi$
and then Fubini’s theorem yield

$\displaystyle\mathbb{E}_{\Phi}^{0}\Biggl{[}\sum_{n=1}^{\infty}Q(B-X_{n}+y)\Biggr{]}$
$\displaystyle=\lambda_{\Phi}\int_{\mathbb{R}^{d}}Q(B-x)\,\mathrm{d}x=\lambda_{\Phi}\int_{\mathbb{R}^{d}}\int_{B-x}Q(\mathrm{d}z)\,\mathrm{d}x$
$\displaystyle=\lambda_{\Phi}\int_{\mathbb{R}^{d}}\int_{B-z}\mathrm{d}x\,Q(\mathrm{d}z)=\lambda_{\Phi}\,|B|,$

which completes the proof.

###### Remark 2

The second term on the right-hand side of (4) is of course equal to $\mu\int
Q(B+y)\,Q(\mathrm{d}y)$. We adopt the form in Lemma 1 due to its
interpretability. Since $Q$ is the distribution for the position of an
offspring point viewed from its parent, $Q^{-}$ represents the distribution
for the location of the parent of the offspring point at the origin on the
event $\\{\Psi(\\{0\\})=1\\}$. On the other hand, $\mu\,Q(B-y)$ gives the
expected number of offspring points falling in
$B\in\mathcal{B}(\mathbb{R}^{d})$ among a cluster whose parent is shifted to
$y\in\mathbb{R}^{d}$. In other words, the second term on the right-hand side
of (4) represents the expected number of offspring points falling in $B$ among
the cluster which is given to have one point at the origin. Since the first
term on the right-hand side of (4) is equal to
$\Lambda_{\Psi}(B)=\mathbb{E}[\Psi(B)]$, Lemma 1 states that the intensity
measure for the reduced Palm version of a stationary PPCP is given as the sum
of the intensity measure of the stationary version and that of a cluster which
has one point at the origin. Lemma 1 is also a slight generalization of the
result in [42, Section 2.2].

The observation in Remark 2 is further enhanced by the following proposition.

###### Proposition 1

For the stationary PPCP $\Psi=\sum_{m=1}^{\infty}\delta_{Y_{m}}$ described in
Section 3.1, the generating functional of the reduced Palm version $\Psi^{!}$
(with respect to the Palm distribution) is given by

$\mathcal{G}_{\Psi}^{0}(h):=\mathbb{E}_{\Psi}^{0}\Biggl{[}\prod_{m=1}^{\infty}h(Y_{m})\Biggr{]}=\mathcal{G}_{\Psi}(h)\int_{\mathbb{R}^{d}}\widetilde{h}(z)\,Q^{-}(\mathrm{d}z),$
(5)

for any measurable function $h$: $\mathbb{R}^{d}\to[0,1]$, where
$\mathcal{G}_{\Psi}$ is the generating functional of the stationary version of
$\Psi$ given as

$\mathcal{G}_{\Psi}(h):=\mathbb{E}\Biggl{[}\sum_{m=1}^{\infty}h(Y_{m})\Biggr{]}=\mathcal{G}_{\Phi}(\widetilde{h})=\exp\biggl{(}-\lambda_{\Phi}\int_{\mathbb{R}^{d}}[1-\widetilde{h}(x)]\,\mathrm{d}x\biggr{)},$
(6)

and $\widetilde{h}(x)$ denotes the generating functional of an offspring
process $\Psi_{1}$ whose parent is shifted to $x\in\mathbb{R}^{d}$;

$\widetilde{h}(x)=\mathcal{G}_{\Psi_{1}}(S_{x}h)=\exp\biggl{(}-\mu\int_{\mathbb{R}^{d}}\bigl{[}1-h(x+y)\bigr{]}\,Q(\mathrm{d}y)\biggr{)}.$
(7)

Note that in Proposition 1 above, $\mathcal{G}_{\Phi}$ is the generating
functional of the parent process $\Phi$. The relation
$\mathcal{G}_{\Psi}(h)=\mathcal{G}_{\Phi}(\widetilde{h})$ with
$\widetilde{h}(x)=\mathcal{G}_{\Psi_{1}}(S_{x}h)$ in (6) and (7) is known to
hold for more general cluster point processes (see, e.g., [17, Example 6.3(a)]
or [9, Proposition 2.3.12 & Lemma 2.3.20]), whereas the last equalities in (6)
and (7) follow because $\Phi$ and $\Psi_{1}$ are PPPs, respectively (see,
e.g., [30, Exercise 3.6] or [9, Corollary 2.1.5]). The relation (5) is derived
in [21, Lemma 1], to which we give another proof using the exchange formula in
Theorem 1.

  * Proof:

As stated in the proof of Lemma 1, once the parent process
$\Phi=\sum_{n=1}^{\infty}\delta_{X_{n}}$ is given, the PPCP $\Psi$ is
conditionally an inhomogeneous PPP with the conditional intensity measure
$\mu\sum_{n=1}^{\infty}S_{-X_{n}}Q$. Since the reduced Palm version of a PPP
is identical in distribution to its original (not conditioned) version, we
have

$\displaystyle\mathbb{E}_{\Psi}^{0}\Biggl{[}\prod_{m=1}^{\infty}h(Y_{m})\Biggm{|}\Phi\Biggr{]}$
$\displaystyle=\mathbb{E}\Biggl{[}\prod_{m=1}^{\infty}h(Y_{m})\biggm{|}\Phi\Biggr{]}$
$\displaystyle=\exp\biggl{(}-\mu\sum_{n=1}^{\infty}\int_{\mathbb{R}^{d}}(1-h(y))\,Q(\mathrm{d}y-X_{n})\biggr{)}=\prod_{n=1}^{\infty}\widetilde{h}(X_{n}),$

where the generating functional of a PPP is applied in the second equality.
Taking the expectation with respect to $\mathbb{P}_{\Psi}^{0}$ and then
applying Theorem 1, we obtain

$\displaystyle\mathcal{G}_{\Psi}^{0}(h)$
$\displaystyle=\mathbb{E}_{\Psi}^{0}\Biggl{[}\prod_{n=1}^{\infty}\widetilde{h}(X_{n})\Biggr{]}=\frac{1}{\mu}\,\mathbb{E}_{\Phi}^{0}\Biggl{[}\int_{\mathbb{R}^{d}}\prod_{n=0}^{\infty}\widetilde{h}(X_{n}-z)\,\Psi_{0}(\mathrm{d}z)\Biggr{]}$
$\displaystyle=\int_{\mathbb{R}^{d}}\widetilde{h}(-z)\,\mathbb{E}_{\Phi}^{0}\Biggl{[}\prod_{n=1}^{\infty}\widetilde{h}(X_{n}-z)\Biggr{]}\,Q(\mathrm{d}z),$

where Campbell’s formula for $\Psi_{0}$ is applied in the last equality. By
Slivnyak’s theorem and the stationarity for $\Phi$, we have
$\mathbb{E}_{\Phi}^{0}\bigl{[}\prod_{n=1}^{\infty}\widetilde{h}(X_{n}-z)\bigr{]}=\mathbb{E}\bigl{[}\prod_{n=1}^{\infty}\widetilde{h}(X_{n})\bigr{]}=\mathcal{G}_{\Phi}(\widetilde{h})$,
which completes the proof.

###### Remark 3

The right-hand side of (5) is given as the generating functional
$\mathcal{G}_{\Psi}(h)$ of the stationary version of $\Psi$ multiplied by the
integral term $\int\widetilde{h}(z)\,Q^{-}(\mathrm{d}z)$. Since
$\widetilde{h}(z)$ represents the generating functional of an offspring
process whose parent is shifted to $z\in\mathbb{R}^{d}$ and $Q^{-}$ is the
distribution of the location of the parent point of the offspring at the
origin on the event $\\{\Psi(\\{0\\})=1\\}$, this integral term represents the
generating functional of the cluster which is given to have a point at the
origin. In other words, Proposition 1 implies that, for a stationary PPCP, its
Palm version is obtained by the superposition of the original stationary
version and an additional independent offspring process whose parent is placed
such that it has an offspring point at the origin. This observation is already
found in, e.g., [40] and is also interpreted such that a point
$z\in\mathbb{R}^{d}$ is first sampled from the distribution $Q^{-}$ and the
Palm version of $\Phi$ at $z$ is then obtained as $\Phi+\delta_{z}$ by
Slivnyak’s theorem, which works as a parent process of the Palm version of
$\Psi$. Proposition 1 indeed supports this interpretation.

###  3.3 Nearest-neighbor distance distributions

For a stationary point process $\Psi$ on $\mathbb{R}^{d}$, let
$\Psi^{!}=\sum_{m=1}^{\infty}\delta_{Y_{m}}$ be its reduced Palm version on
$\\{\Psi(\\{0\\})=1\\}$ and let $Y_{*}$ denote the nearest point of $\Psi^{!}$
from the origin. Then, the nearest-neighbor distance distribution for $\Psi$
is defined as the probability distribution for $\|Y_{*}\|$ with respect to
$\mathbb{P}_{\Psi}^{0}$, where $\|\cdot\|$ denotes the Euclidean distance. We
show below that the nearest-neighbor distance distribution for a stationary
PPCP is obtained in a similar way to Proposition 1.

###### Proposition 2

For the stationary PPCP $\Psi$ described in Section 3.1, the complementary
nearest-neighbor distance distribution is given by

$\mathbb{P}_{\Psi}^{0}(\|Y_{*}\|>r)=\mathcal{G}_{\Phi}(h^{*}_{r})\int
h^{*}_{r}(t)\,Q^{-}(\mathrm{d}t),\quad r\geq 0,$ (8)

where $h^{*}_{r}(x)=e^{-\mu Q(b_{0}(r)-x)}$ and $b_{0}(r)$ denotes a
$d$-dimensional ball centered at the origin with radius $r$.

  * Proof:

As with the proof of Proposition 1, we consider the conditional probability
given the parent process $\Phi=\sum_{n=1}^{\infty}\delta_{X_{n}}$ and obtain

$\displaystyle\mathbb{P}_{\Psi}^{0}(\|Y_{*}\|>r\mid\Phi)$
$\displaystyle=\mathbb{P}_{\Psi}^{0}\bigl{(}\Psi^{!}(b_{0}(r))=0\mid\Phi\bigr{)}=\mathbb{P}\bigl{(}\Psi(b_{0}(r))=0\mid\Phi\bigr{)}$
$\displaystyle=\prod_{n=1}^{\infty}e^{-\mu
Q(b_{0}(r)-X_{n})}=\prod_{n=1}^{\infty}h^{*}_{r}(X_{n}),$ (9)

where the second equality follows from Slivnyak’s theorem and the third does
because $\Psi^{!}$ is conditionally an inhomogeneous PPP with the intensity
measure $\mu\sum_{n=1}^{\infty}S_{-X_{n}}Q$ when $\Phi$ is given. The rest of
the proof is similar to that of Proposition 1.

###### Remark 4

In (8), the term $\mathcal{G}_{\Phi}(h^{*}_{r})$ is the complementary contact
distance distribution and that is obtained by taking the expectation of
(Proof:) with respect to $\mathbb{P}$, instead of $\mathbb{P}_{\Psi}^{0}$
(see, e.g., [33]). The result of Proposition 2 is consistent with the existing
ones in, e.g., [11, 3, 4, 37] and gives a unified approach to derive the
nearest-neighbor distance distributions for stationary PPCPs.

##  4 Applications to wireless networks with hotspot clusters

In this section, we apply Theorem 1 to the analysis of a D2D network with
hotspot clusters modeled using a stationary PPCP. We here suppose $d=2$, but
unless otherwise specified, the discussion holds for $d\geq 2$ theoretically.

###  4.1 Model of a Device-to-Device Network

Wireless devices are distributed on $\mathbb{R}^{d}$ according to a stationary
point process $\Psi=\sum_{m=1}^{\infty}\delta_{Y_{m}}$. At each time slot,
each device is in transmission mode with probability $p\in(0,1)$ or in
receiving mode with probability $1-p$ independently of the others (half duplex
with random access). Devices in the transmission mode transmit signals but can
not receive ones, whereas the devices in the receiving mode can receive
signals but can not transmit ones. We assume that all transmitting devices
transmit signals with identical transmission power (normalized to one) and
share a common frequency spectrum. The path-loss function representing
attenuation of signals with distance is given by $\ell$ satisfying
$\ell(r)\geq 0$, $r>0$, and
$\int_{\epsilon}^{\infty}\ell(r)\,r^{d-1}\,\mathrm{d}r<\infty$ for
$\epsilon>0$. We further assume that all wireless links receive Rayleigh
fading effects while we ignore shadowing effects. We focus on the device at
the origin, referred to as the typical device, under the condition of
$\\{\Psi(\\{0\\})=1\\}$ and examine whether the typical device can decode
messages from other transmitting devices. Let
$\Psi_{\mathsf{Tx}}=\sum_{m=1}^{\infty}\delta_{Y^{\prime}_{m}}$ denote the
sub-process of $\Psi$ representing the locations of devices in the
transmission mode and for each $m\in\mathbb{N}$, let $H_{m}$ denote a random
variable representing the fading effect on signals transmitted from the device
at $Y^{\prime}_{m}$, where $H_{m}$, $m\in\mathbb{N}$, are mutually
independent, independent of $\Psi_{\mathsf{Tx}}$ and exponentially distributed
with unit mean due to the Rayleigh fading. With this setup, the received
signal power by the typical device amounts to $H_{m}\,\ell(Y^{\prime}_{m})$
when it receives signals from the device at $Y^{\prime}_{m}$. Hence, if the
typical device is in the receiving mode and communicates with the transmitting
device at $Y^{\prime}_{m}$, the signal-to-interference-plus-noise ratio (SINR)
is given as

$\mathsf{SINR}_{m}=\frac{H_{m}\,\ell(\|Y^{\prime}_{m}\|)}{\sum_{\begin{subarray}{c}j=1\\\
j\neq m\end{subarray}}^{\infty}H_{j}\,\ell(\|Y^{\prime}_{j}\|)+N},$ (10)

where $N$ denotes a constant representing noise at the origin. We suppose that
the typical device can successfully decode a message from the device at
$Y^{\prime}_{m}$ if the typical device is in the receiving mode and
$\mathsf{SINR}_{m}$ in (10) exceeds a predefined threshold $\theta>0$.

###  4.2 Coverage analysis

We here suppose that a device in the receiving mode communicates with the
nearest device in transmission mode. The probability that the typical device
can successfully decode a message from its partner is called the coverage
probability and is given by

$\mathsf{CP}(\theta)=(1-p)\sum_{m=1}^{\infty}\mathbb{P}_{\Psi}^{0}\bigl{(}\mathsf{SINR}_{m}>\theta,\|Y^{\prime}_{m}\|\leq\|Y^{\prime}_{j}\|,j\in\mathbb{N}\bigr{)},$
(11)

where $1-p$ on the right-hand side indicates that the typical device must be
in the receiving mode and the sum over $m\in\mathbb{N}$ represents the
probability that the SINR from the nearest transmitting device exceeds the
threshold $\theta$. We now suppose that the point process $\Psi$ representing
the locations of devices is given as a stationary PPCP studied in Section 3.
Then, $\Psi_{\mathsf{Tx}}=\sum_{m=1}^{\infty}\delta_{Y^{\prime}_{m}}$
representing the locations of devices in the transmission mode is also a
stationary PPCP, where the parent process remains the same as the homogeneous
PPP $\Phi$ with intensity $\lambda_{\Phi}$, whereas the offspring processes
$\Psi^{\prime}_{n}=\sum_{k=1}^{\kappa^{\prime}_{n}}\delta_{Y^{\prime}_{n,k}}$,
$n\in\mathbb{N}$, are finite PPPs with the intensity measure $p\mu Q$.

###### Theorem 2

For the model of a D2D network described in Section 4.1 with the devices
deployed according to a stationary PPCP in Section 3.1, the coverage
probability is given by

$\mathsf{CP}(\theta)=(1-p)p\mu\int_{\mathbb{R}^{d}}\bigl{(}I_{1,\theta}(t)+I_{2,\theta}(t)\bigr{)}\,Q^{-}(\mathrm{d}t),$
(12)

where $Q^{-}$ is given in Lemma 1 and

$\displaystyle I_{1,\theta}(t)$ $\displaystyle=\int_{\mathbb{R}^{d}}e^{-\theta
N/\ell(\|y\|)}\,C_{\theta}(y,t)E_{\theta}(y)\,Q(\mathrm{d}y-t),$
$\displaystyle I_{2,\theta}(t)$
$\displaystyle=\lambda_{\Phi}\int_{\mathbb{R}^{d}}\\!\int_{\mathbb{R}^{d}}e^{-\theta
N/\ell(\|y\|)}\,C_{\theta}(y,t)\,C_{\theta}(y,x)\,E_{\theta}(y)\,Q(\mathrm{d}y-x)\,\mathrm{d}x,$
$\displaystyle E_{\theta}(y)$
$\displaystyle=\exp\biggl{(}-\lambda_{\Phi}\int_{\mathbb{R}^{d}}\bigl{[}1-C_{\theta}(y,w)\bigr{]}\,\mathrm{d}w\biggr{)},$
$\displaystyle C_{\theta}(y,x)$
$\displaystyle=\exp\biggl{(}-p\mu\,\biggl{[}1-\int_{\|z\|>\|y\|}\Bigl{(}1+\theta\,\frac{\ell(\|z\|)}{\ell(\|y\|)}\Bigr{)}^{-1}\,Q(\mathrm{d}z-x)\biggr{]}\biggr{)}.$
(13)

Before proceeding on the proof of Theorem 2, we give an intuitive
interpretation to the result of it. First, as stated in the preceding section,
$Q^{-}$ denotes the distribution for the location of the parent point of the
typical device at the origin. Thus, $p\mu I_{1,\theta}(t)$ and $p\mu
I_{2,\theta}(t)$ in (12) represent the cases where the typical device, whose
parent is located at $t\in\mathbb{R}^{d}$, communicates with the transmitting
device in the same cluster and in a different cluster, respectively; that is,
the location of the communication partner is sampled from a finite PPP with
the intensity measure $p\mu Q(\mathrm{d}y-t)$ in $I_{1,\theta}(t)$ and is from
one with $p\mu Q(\mathrm{d}y-x)$ in $I_{2,\theta}(t)$, where $x$ is also
sampled from a homogeneous PPP with intensity $\lambda_{\Phi}$. Moreover,
$E_{\theta}(y)$ represents the effect from other clusters which are neither
the one having the typical device nor the one having its communication partner
at $y$. Finally, $C_{\theta}(y,x)$ represents the effect of the cluster with
the parent point at $x\in\mathbb{R}^{d}$ when the typical device communicates
with the transmitting device at $y$.

  * Proof:

Similar to the proof of Proposition 1, once the parent process
$\Phi=\sum_{n=1}^{\infty}\delta_{X_{n}}$ is given, the point process
$\Psi_{\mathsf{Tx}}$ representing the locations of devices in the transmission
mode is conditionally an inhomogeneous PPP with the conditional intensity
measure $p\mu\sum_{n=1}^{\infty}S_{-X_{n}}Q$. Thus, we can use the
corresponding approach to that obtaining the coverage probability for a
cellular network with BSs deployed according to a PPP (see, e.g,, [5] or [12,
Section 5.2]). Since $H_{m}$, $m\in\mathbb{N}$, are mutually independent,
exponentially distributed, and also independent of $\Phi$, we have from (11),

$\displaystyle\mathbb{P}_{\Psi}^{0}\bigl{(}\mathsf{SINR}_{m}>\theta,\|Y^{\prime}_{m}\|\leq\|Y^{\prime}_{j}\|,j\in\mathbb{N}\bigm{|}\Phi\bigr{)}$
$\displaystyle=\mathbb{P}_{\Psi}^{0}\Biggl{(}H_{m}>\frac{\theta}{\ell(\|Y^{\prime}_{m}\|)}\,\Biggl{(}\sum_{\begin{subarray}{c}j=1\\\
j\neq
m\end{subarray}}^{\infty}H_{j}\,\ell(\|Y^{\prime}_{j}\|)+N\Biggr{)},\>\|Y^{\prime}_{m}\|\leq\|Y^{\prime}_{j}\|,j\in\mathbb{N}\Biggm{|}\Phi\Biggr{)}$
$\displaystyle=\mathbb{E}_{\Psi}^{0}\Biggl{[}e^{-\theta
N/\ell(\|Y^{\prime}_{m}\|)}\prod_{\begin{subarray}{c}j=1\\\ j\neq
m\end{subarray}}^{\infty}\biggl{(}1+\theta\,\frac{\ell(\|Y^{\prime}_{j}\|)}{\ell(\|Y^{\prime}_{m}\|)}\,\biggr{)}^{-1}\,\bm{1}_{\\{\|Y^{\prime}_{j}\|>\|Y^{\prime}_{m}\|\\}}\Biggm{|}\Phi\Biggr{]},$

where $\bm{1}_{A}$ denotes the indicator function for set $A$ and we use
$\mathbb{P}(H_{m}>a)=e^{-a}$ for $a\geq 0$ and
$\mathbb{E}[e^{-sH_{j}}]=(1+s)^{-1}$ in the last equality. Summing the above
expression over $m\in\mathbb{N}$, we have from Slivnyak’s Theorem for $\Psi$
conditioned on $\Phi$ and the refined Campbell formula (see, e.g., [18,
Theorem 13.2.III], [30, Theorem 9.1] or [9, Theorem 3.1.9]),

$\displaystyle\sum_{m=1}^{\infty}\mathbb{P}_{\Psi}^{0}\bigl{(}\mathsf{SINR}_{m}>\theta,\|Y^{\prime}_{m}\|\leq\|Y^{\prime}_{j}\|,j\in\mathbb{N}\bigm{|}\Phi\bigr{)}$
$\displaystyle=p\mu\sum_{n=1}^{\infty}\int_{\mathbb{R}^{d}}e^{-\theta
N/\ell(\|y\|)}\,\mathbb{E}\Biggl{[}\prod_{j=1}^{\infty}\biggl{(}1+\theta\,\frac{\ell(\|Y^{\prime}_{j}\|)}{\ell(\|y\|)}\,\biggr{)}^{-1}\,\bm{1}_{\\{\|Y^{\prime}_{j}\|>\|y\|\\}}\Biggm{|}\Phi\Biggr{]}\,Q(\mathrm{d}y-X_{n}).$
(14)

Furthermore, the generating functional of a PPP applying to the above
expectation yields

$\displaystyle\mathbb{E}\Biggl{[}\prod_{j=1}^{\infty}\biggl{(}1+\theta\,\frac{\ell(\|Y^{\prime}_{j}\|)}{\ell(\|y\|)}\,\biggr{)}^{-1}\,\bm{1}_{\\{\|Y^{\prime}_{j}\|>\|y\|\\}}\Biggm{|}\Phi\Biggr{]}$
$\displaystyle=\exp\biggl{(}-p\mu\sum_{i=1}^{\infty}\int_{\mathbb{R}^{d}}\biggl{[}1-\biggl{(}1+\theta\,\frac{\ell(\|z\|)}{\ell(\|y\|)}\,\biggr{)}^{-1}\,\bm{1}_{\\{\|z\|>\|y\|\\}}\biggr{]}\,Q(\mathrm{d}z-X_{i})\biggr{)}$
$\displaystyle=\prod_{i=1}^{\infty}\exp\biggl{(}-p\mu\,\biggl{[}1-\int_{\|z\|>\|y\|}\biggl{(}1+\theta\,\frac{\ell(\|z\|)}{\ell(\|y\|)}\,\biggr{)}^{-1}\,Q(\mathrm{d}z-X_{i})\biggr{]}\biggr{)}=\prod_{i=1}^{\infty}C_{\theta}(y,X_{i}).$

Plugging this into (Proof:), taking the expectation with respect to
$\mathbb{P}_{\Psi}^{0}$ and then applying Neveu’s exchange formula in Theorem
1, we have

$\displaystyle\mathsf{CP}(\theta)$
$\displaystyle=(1-p)p\mu\,\mathbb{E}_{\Psi}^{0}\Biggl{[}\sum_{n=1}^{\infty}\int_{\mathbb{R}^{d}}e^{-\theta
N/\ell(\|y\|)}\prod_{i=1}^{\infty}C_{\theta}(y,X_{i})\,Q(\mathrm{d}y-X_{n})\Biggr{]}$
$\displaystyle=(1-p)p\,\mathbb{E}_{\Phi}^{0}\Biggl{[}\int_{\mathbb{R}^{d}}\sum_{n=0}^{\infty}\int_{\mathbb{R}^{d}}e^{-\theta
N/\ell(\|y\|)}\prod_{i=0}^{\infty}C_{\theta}(y,X_{i}-t)\,Q(\mathrm{d}y-X_{n}+t)\,\Psi_{0}(\mathrm{d}t)\Biggr{]}$
$\displaystyle=(1-p)p\mu\,\int_{\mathbb{R}^{d}}\mathbb{E}_{\Phi}^{0}\Biggl{[}\sum_{n=0}^{\infty}\int_{\mathbb{R}^{d}}e^{-\theta
N/\ell(\|y\|)}\prod_{i=0}^{\infty}C_{\theta}(y,X_{i}-t)\,Q(\mathrm{d}y-X_{n}+t)\,\Biggr{]}\,Q(\mathrm{d}t),$
(15)

where we note the existence of $X_{0}=0$ on $\\{\Phi(\\{0\\})=1\\}$ in the
second equality and apply Campbell’s formula in the third equality. Noting
that $X_{0}=0$ on $\\{\Phi(\\{0\\})=1\\}$, we separate the expectation in
(Proof:) into

$\displaystyle\mathbb{E}_{\Phi}^{0}\Biggl{[}\sum_{n=0}^{\infty}\int_{\mathbb{R}^{d}}e^{-\theta
N/\ell(\|y\|)}\prod_{i=0}^{\infty}C_{\theta}(y,X_{i}-t)\,Q(\mathrm{d}y-X_{n}+t)\,\Biggr{]}$
$\displaystyle=\int_{\mathbb{R}^{d}}e^{-\theta
N/\ell(\|y\|)}\,C_{\theta}(y,-t)\,\mathbb{E}_{\Phi}^{0}\Biggl{[}\prod_{i=1}^{\infty}C_{\theta}(y,X_{i}-t)\Biggr{]}\,Q(\mathrm{d}y+t)$
$\displaystyle\quad\mbox{}+\mathbb{E}_{\Phi}^{0}\Biggl{[}\sum_{n=1}^{\infty}\int_{\mathbb{R}^{d}}e^{-\theta
N/\ell(\|y\|)}\,C_{\theta}(y,-t)\,C_{\theta}(y,X_{n}-t)\,\prod_{\begin{subarray}{c}i=1\\\
i\neq
n\end{subarray}}^{\infty}C_{\theta}(y,X_{i}-t)\,Q(\mathrm{d}y-X_{n}+t)\,\Biggr{]},$
(16)

and consider the two terms on the right-hand side of (Proof:) one by one. For
the first term, the generating functional of a PPP yields

$\displaystyle(\text{1st term of~{}\eqref{eq:cover-prf3}})$
$\displaystyle=\int_{\mathbb{R}^{d}}e^{-\theta
N/\ell(\|y\|)}\,C_{\theta}(y,-t)\,\exp\biggl{(}-\lambda_{\Phi}\int_{\mathbb{R}^{d}}[1-C_{\theta}(y,w)]\,\mathrm{d}w\biggr{)}\,Q(\mathrm{d}y+t)$
$\displaystyle=I_{1,\theta}(-t).$ (17)

On the other hand, applying Campbell’s formula and the generating functional
for $\Phi$ to the second term on the right-hand side of (Proof:), we have

$\displaystyle(\text{2nd term of~{}\eqref{eq:cover-prf3}})$
$\displaystyle=\lambda_{\Phi}\int_{\mathbb{R}^{d}}\int_{\mathbb{R}^{d}}e^{-\theta
N/\ell(\|y\|)}\,C_{\theta}(y,-t)\,C_{\theta}(y,x)\,\mathbb{E}_{\Phi}^{0}\Biggl{[}\prod_{i=1}^{\infty}C_{\theta}(y,X_{i}-t)\Biggr{]}Q(\mathrm{d}y-x)\,\mathrm{d}x$
$\displaystyle=\lambda_{\Phi}\int_{\mathbb{R}^{d}}\int_{\mathbb{R}^{d}}e^{-\theta
N/\ell(\|y\|)}\,C_{\theta}(y,-t)\,C_{\theta}(y,x)\,\exp\biggl{(}-\lambda_{\Phi}\int_{\mathbb{R}^{d}}[1-C_{\theta}(y,w)]\,\mathrm{d}w\biggr{)}\,Q(\mathrm{d}y-x)\,\mathrm{d}x$
$\displaystyle=I_{2,\theta}(-t).$ (18)

Finally, plugging (Proof:) and (Proof:) into (Proof:), and then into (Proof:),
we have (12) and the proof is completed.

When $d=2$ and the distribution $Q$ for the locations of offspring points
depends only on the distance; that is,
$Q(\mathrm{d}y)=f_{\mathrm{o}}(\|y\|)\,\mathrm{d}y$ for $y\in\mathbb{R}^{2}$,
we obtain a numerically computable form of the coverage probability.

###### Corollary 1

When $d=2$ and $Q(\mathrm{d}y)=f_{\mathrm{o}}(\|y\|)\,\mathrm{d}y$,
$y\in\mathbb{R}^{2}$, the coverage probability in Theorem 2 is reduced to

$\mathsf{CP}(\theta)=2\pi(1-p)p\mu\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}e^{-\theta
N/\ell(s)}\,\widehat{E}_{\theta}(s)\,\widehat{C}_{\theta}(s,u)\,\widehat{I}_{\theta}(s,u)\,\mathrm{d}s\,f_{\mathrm{o}}(u)\,u\,\mathrm{d}u,$
(19)

where

$\displaystyle\widehat{I}_{\theta}(s,u)$ $\displaystyle=g(s\mid
u)+2\pi\lambda_{\Phi}\int_{0}^{\infty}\widehat{C}_{\theta}(s,r)\,g(s\mid
r)\,r\,\mathrm{d}r,$ $\displaystyle\widehat{E}_{\theta}(s)$
$\displaystyle=\exp\biggl{(}-2\pi\lambda_{\Phi}\int_{0}^{\infty}\bigl{[}1-\widehat{C}_{\theta}(s,v)\bigr{]}\,v\,\mathrm{d}v\biggr{)},$
$\displaystyle\widehat{C}_{\theta}(s,r)$
$\displaystyle=\exp\biggl{(}-p\mu\biggl{[}1-\int_{s}^{\infty}\Bigl{(}1+\theta\,\frac{\ell(q)}{\ell(s)}\Bigr{)}^{-1}g(q\mid
r)\,\mathrm{d}q\biggr{]}\biggr{)},$ (20) $\displaystyle g(s\mid r)$
$\displaystyle=2s\int_{0}^{\pi}f_{\mathrm{o}}(\sqrt{s^{2}+r^{2}-2sr\cos\varphi})\,\mathrm{d}\varphi.$

  * Proof:

Since the distribution $Q$ depends only on the distance, it holds that
$Q^{-}(\mathrm{d}t)=Q(\mathrm{d}t)=f_{\mathrm{o}}(\|t\|)\,\mathrm{d}t$,
$t\in\mathbb{R}^{2}$, and (12) is reduced to

$\displaystyle\mathsf{CP}(\theta)$
$\displaystyle=(1-p)p\mu\int_{\mathbb{R}^{2}}\bigl{(}I_{1,\theta}(t)+I_{2,\theta}(t)\bigr{)}\,f_{\mathrm{o}}(\|t\|)\,\mathrm{d}t$
$\displaystyle=2\pi(1-p)p\mu\int_{0}^{\infty}\bigl{(}\widehat{I}_{1,\theta}(u)+\widehat{I}_{2,\theta}(u)\bigr{)}\,f_{\mathrm{o}}(u)\,u\,\mathrm{d}u,$
(21)

where the polar coordinate conversion is applied in the second equality and

$\displaystyle\widehat{I}_{1,\theta}(u)$
$\displaystyle=\int_{0}^{\infty}e^{-\theta
N/\ell(s)}\,\widehat{C}_{\theta}(s,u)\,\widehat{E}_{\theta}(s)\,g(s\mid
u)\,\mathrm{d}s,$ $\displaystyle\widehat{I}_{2,\theta}(u)$
$\displaystyle=2\pi\lambda_{\Phi}\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}e^{-\theta
N/\ell(s)}\,\widehat{C}_{\theta}(s,u)\,\widehat{C}_{\theta}(s,r)\,\widehat{E}_{\theta}(s)\,g(s\mid
r)\,\mathrm{d}s\,r\,\mathrm{d}r.$

Therefore, we have

$\widehat{I}_{1,\theta}(u)+\widehat{I}_{2,\theta}(u)=\int_{0}^{\infty}e^{-\theta
N/\ell(s)}\,\widehat{C}_{\theta}(s,u)\,\widehat{E}_{\theta}(s)\,\widehat{I}_{\theta}(s,u)\,\mathrm{d}s.$

Plugging this into (21), we have (19) and the proof is completed.

###  4.3 Device discovery

We next consider the problem of device discovery. Devices in the transmission
mode transmit broadcast messages, whereas a device in the receiving mode can
discover the transmitters if it can successfully decode the broadcast
messages. When a device in the receiving mode receives the signal from one
transmitting device, the signals from all other transmitting devices work as
interference. Then, the expected number of transmitting devices discovered by
the typical device is represented by

$\mathcal{N}(\theta)=(1-p)\,\mathbb{E}_{\Psi}^{0}\Biggl{[}\sum_{m=1}^{\infty}\bm{1}_{\\{\mathsf{SINR}_{m}>\theta\\}}\Biggr{]}.$
(22)

###### Proposition 3

Consider the D2D network model described in Section 4.1 with the devices
deployed according to a stationary PPCP given in Section 3.1. Then, the
expected number $\mathcal{N}(\theta)$ of transmitting devices discovered by
the typical device is obtained by (12) in Theorem 2 replacing the integral
range $\|z\|>\|y\|$ in (2) by $\mathbb{R}^{d}$. Moreover, when $d=2$ and
$Q(\mathrm{d}y)=f_{\mathrm{o}}(\|y\|)\,\mathrm{d}y$ for $y\in\mathbb{R}^{2}$,
$\mathcal{N}(\theta)$ is reduced to (19) in Corollary 1 replacing the integral
range $(s,\infty)$ in (20) by $(0,\infty)$.

  * Proof:

Since
$\mathbb{E}_{\Psi}^{0}\bigl{[}\sum_{m=1}^{\infty}\bm{1}_{\\{\mathsf{SINR}_{m}>\theta\\}}\bigr{]}=\sum_{m=1}^{\infty}\mathbb{P}_{\Psi}^{0}(\mathsf{SINR}_{m}>\theta)$,
the difference between (11) and (22) is only the event
$\\{\|Y^{\prime}_{m}\|\leq\|Y^{\prime}_{j}\|,j\in\mathbb{N}\\}$. This leads to
the difference of the integral ranges in $C_{\theta}(y,x)$ in (2) and in
$\widehat{C}_{\theta}(s,r)$ in (20).

###### Remark 5

Since
$\mathbb{P}_{\Psi}^{0}\bigl{(}\bigcup_{m=1}^{\infty}\\{\mathsf{SINR}_{m}>\theta\\}\bigr{)}\leq\sum_{m=1}^{\infty}\mathbb{P}_{\Psi}^{0}(\mathsf{SINR}_{m}>\theta)=\mathcal{N}(\theta)$,
Proposition 3 also gives an upper bound for the coverage probability with the
max-SINR association policy, where a device in the receiving mode receives a
message with the strongest SINR. This upper bound is known to be exact for
$\theta>1$ since
$\sum_{m=1}^{\infty}\bm{1}_{\\{\mathsf{SINR}_{m}>\theta\\}}\leq 1+\theta^{-1}$
(see [19] or [12, Lemma 5.1.2]).

###  4.4 Numerical experiments

We present the results of numerical experiments for the analytical results
obtained in Sections 4.2 and 4.3. We set $d=2$ and the distribution $Q$ for
the location of the offspring points as
$Q(\mathrm{d}y)=f_{\mathrm{o}}(\|y\|)\,\mathrm{d}y$ and
$f_{\mathrm{o}}(s)=e^{-s^{2}/(2\sigma^{2})}/(2\pi\sigma^{2})$, $s\geq 0$; that
is, $Q$ is the isotropic normal distribution with variance $\sigma^{2}$, so
that the resulting PPCP $\Psi$ is the Thomas point process. Furthermore, the
path-loss function is set as $\ell(r)=r^{-\beta}$, $r>0$, with $\beta>2$.

The numerical results for the coverage probability are given in Figure 2,
where the values of $\mathsf{CP}(\theta)$ with different values of $\theta$
and $\sigma^{2}$ are plotted. The other parameters are fixed at
$\lambda_{\Phi}=\pi^{-1}$, $\mu=10$, $p=0.5$, $\beta=4$ and $N=0$. For
comparison, the values when the devices are located according to a homogeneous
PPP are also displayed in the figure with the label “$\sigma^{2}\to\infty$.”
From Figure 2, we can see that, as the value of $\sigma^{2}$ increases, the
coverage probability decreases and is closer to that with the homogeneous PPP.
This is contrary to the case of cellular networks, where the coverage
probability increases and is closer to that with the homogeneous PPP from
below as the variance of the locations of offspring points increases (see
[33]). This difference is thought to be due to the fact that the locations of
a receiving device and its communication partner are near to each other in the
PPCP-deployed D2D network since they are both points of the same PPCP, whereas
the location of a receiver is likely far from that of the associated BS in the
PPCP-deployed cellular network since their locations are independent of each
other.

Figure 2: Coverage probability as a function of SINR threshold
($\lambda_{\Phi}=\pi^{-1}$, $\mu=10$, $p=0.5$, $\beta=4$ and $N=0$)

The results of the device discovery is given in Figure 3, where we know that
the closed form expression of the expected number of discovered devices is
obtained as
$\mathcal{N}^{(\mathsf{PPP})}(\theta)=(1-p)(\beta/2\pi)\sin(2\pi/\beta)\,\theta^{-2/\beta}$
for the case of the homogeneous PPP with $N\equiv 0$ (see, e.g., [25]). The
figure shows similar features to the coverage probability.

Figure 3: Expected number of discovered devices as a function of SINR
threshold ($\lambda_{\Phi}=\pi^{-1}$, $\mu=10$, $p=0.5$, $\beta=4$ and $N=0$)

##  5 Conclusion

In this paper, we have presented an alternative form of Neveu’s exchange
formula for jointly stationary point processes on $\mathbb{R}^{d}$ and then
demonstrated that it is useful for the analysis of spatial stochastic models
given based on stationary PPCPs. We have first applied it to the Palm
characterization for a stationary PPCP and then to the analysis of a D2D
network modeled using a stationary PPCP. Although we have only considered some
fundamental problems, we expect that the new form of the exchange formula will
be utilized for the analysis of more sophisticated models leading up to the
development of 5G and beyond networks.

## References

  * [1] M. Afshang and H. S. Dhillon.  Poisson clustered process based analysis of HetNets with correlated user and base station locations.  IEEE Transactions on Wireless Communications, 17:2417–2431, 2018\.
  * [2] M. Afshang, H. S. Dhillon, and P. H. J. Chong.  Modeling and performance analysis of clustered device-to-device networks.  IEEE Transactions on Wireless Communications, 15:4957–4972, 2016\.
  * [3] M. Afshang, C. Saha, and H. S. Dhillon.  Nearest-neighbor and contact distance distributions for Matérn cluster process.  IEEE Communications Letters, 21:2686–2689, 2017.
  * [4] M. Afshang, C. Saha, and H. S. Dhillon.  Nearest-neighbor and contact distance distributions for Thomas cluster process.  IEEE Wireless Communications Letters, 6:130–133, 2017.
  * [5] J. G. Andrews, F. Baccelli, and R. K. Ganti.  A tractable approach to coverage and rate in cellular networks.  IEEE Transactions on Communications, 59:3122–3134, 2011.
  * [6] J. G. Andrews, A. K. Gupta, and H. S. Dhillon.  A primer on cellular network analysis using stochastic geometry.  arXiv:1604.03183 [cs.IT], 2016.
  * [7] F. Baccelli and B. Błaszczyszyn.  Stochastic geometry and wireless networks, Volume I: Theory.  Foundations and Trends in Networking, 3:249–449, 2009.
  * [8] F. Baccelli and B. Błaszczyszyn.  Stochastic geometry and wireless networks, Volume II: Applications.  Foundations and Trends in Networking, 4:1–312, 2009.
  * [9] F. Baccelli, B. Błaszczyszyn, and M. Karray.  Random measures, point processes, and stochastic geometry.  https://hal.inria.fr/hal-02460214, 2020.
  * [10] F. Baccelli, N. Khud, R. Laroia, J. Li, T. Richardson, S. Shakkottai, S. Tavildar, and X. Wu.  On the design of device-to-device autonomous discovery.  In 2012 Fourth International Conference on Communication Systems and Networks (COMSNETS 2012), pages 1–9, 2012.
  * [11] M. Baudin.  Likelihood and nearest-neighbor distance properties of multidimensional Poisson cluster processes.  Journal of Applied Probability, 18:879–888, 1981.
  * [12] B. Błaszczyszyn, M. Haenggi, P. Keeler, and S. Mukherjee.  Stochastic Geometry Analysis of Cellular Networks.  Cambridge University Press, 2018.
  * [13] B. Błaszczyszyn and D. Yogeshwaran.  Directionally convex ordering of random measures, shot noise fields, and some applications to wireless communications.  Advances in Applied Probability, 41:623–646, 2009.
  * [14] L. Chen, W. Wang, and C. Zhang.  Stochastic wireless powered communication networks with truncated cluster point process.  IEEE Transactions on Vehicular Technology, 66:11286–11294, 2017\.
  * [15] S. N. Chiu, D. Stoyan, W. S. Kendall, and J. Mecke.  Stochastic Geometry and its Applications.  Wiley, 3rd edition, 2013.
  * [16] Y. J. Chun, M. O. Hasna, and A. Ghrayeb.  Modeling heterogeneous cellular networks interference using Poisson cluster processes.  IEEE Journal on Selected Areas in Communications, 33:2182–2195, 2015.
  * [17] D. J. Daley and D. Vere-Jones.  An Introduction to the Theory of Point Processes: Volume I: Elementary Theory and Methods.  Springer, 2nd edition, 2003.
  * [18] D. J. Daley and D. Vere-Jones.  An Introduction to the Theory of Point Processes: Volume II: General Theory and Structure.  Springer, 2nd edition, 2008.
  * [19] H. S. Dhillon, R. K. Ganti, F. Baccelli, and J. G. Andrews.  Modeling and analysis of $K$-tier downlink heterogeneous cellular networks.  IEEE Journal on Selected Areas in Communications, 30:550–560, 2012\.
  * [20] H. ElSawy, A. Sultan-Salem, M.-S. Alouini, and M. Z. Win.  Modeling and analysis of cellular networks using stochastic geometry: A tutorial.  IEEE Communications Surveys & Tutorials, 19:167–203, 2017.
  * [21] R. K. Ganti and M. Haenggi.  Interference and outage in clustered wireless ad hoc networks.  IEEE Transactions on Information Theory, 55:4067–4086, 2009.
  * [22] D. Gentner and G. Last.  Palm pairs and the general mass-transport principle.  Mathematische Zeitschrift, 267:695–716, 2011.
  * [23] M. Haenggi.  Stochastic Geometry for Wireless networks.  Cambridge University Press, 2013.
  * [24] M. Haenggi and R. K. Ganti.  Interference in large wireless networks.  Foundations and Trends in Networking, 3:127–248, 2009.
  * [25] E. B. Hamida, G. Chelius, A. Busson, and E. Fleury.  Neighbor discovery in multi-hop wireless networks: Evaluation and dimensioning with interference considerations.  Discrete Mathematics and Theoretical Computer Science, 10:87–114, 2008.
  * [26] Y. Hmamouche, M. Benjillali, S. Saoudi, H. Yanikomeroglu, and M. Di Renzo.  New trends in stochastic geometry for wireless networks: A tutorial and survey.  Proceedings of the IEEE, 109:1200–1252, 2021.
  * [27] T. Kwon and J.-W. Choi.  Spatial performance analysis and design principles for wireless peer discovery.  IEEE Transactions on Wireless Communications, 13:4507–4519, 2014\.
  * [28] T. Kwon, H. Ju, and H. Lee.  Performance study for random access-based wireless mutual broadcast networks with Ginibre point processes.  IEEE Communications Letters, 24:1581–1585, 2020.
  * [29] G. Last.  Modern random measures: Palm theory and related models.  In W. S. Kendall and I. Molchanov, editors, New Perspectives in Stochastic Geometry, pages 77–110. Oxford University Press, 2010.
  * [30] G. Last and M. Penrose.  Lectures on the Poisson Process.  Cambridge University Press, 2017.
  * [31] G. Last and H. Thorisson.  Invariant transports of stationary random measures and mass-stationarity.  The Annals of Probability, 37:790–813, 2009.
  * [32] X. Lu, M. Salehi, M. Haenggi, E. Hossain, and H. Jiang.  Stochastic geometry analysis of spatial-temporal performance in wireless networks: A tutorial.  IEEE Communications Surveys & Tutorials, 23:2753–2801, 2021.
  * [33] N. Miyoshi.  Downlink coverage probability in cellular networks with Poisson-Poisson cluster deployed base stations.  IEEE Wireless Communications Letters, 8:5–8, 2019.
  * [34] N. Miyoshi and T. Shirai.  A cellular network model with Ginibre configured base stations.  Advances in Applied Probability, 46:832–845, 2014.
  * [35] N. Miyoshi and T. Shirai.  Spatial modeling and analysis of cellular networks using the Ginibre point process: A tutorial.  IEICE Transactions on Communications, E99-B:2247–2255, 2016.
  * [36] S. Mukherjee.  Analytical Modeling of Heterogeneous Cellular Networks: Geometry, Coverage, and Capacity.  Cambridge University Press, 2014.
  * [37] K. Pandey, H. S. Dhillon, and A. K. Gupta.  On the contact and nearest-neighbor distance distributions for the $n$-dimensional Matérn cluster process.  IEEE Wireless Communications Letters, 9:394–397, 2020.
  * [38] C. Saha, M. Afshang, and H. S. Dhillon.  Enriched $K$-tier HetNet model to enable the analysis of user-centric small cell deployments.  IEEE Transactions on Wireless Communications, 16:1593–1608, 2017\.
  * [39] C. Saha, M. Afshang, and H. S. Dhillon.  3GPP-inspired HetNet model using Poisson cluster process: Sum-product functionals and downlink coverage.  IEEE Transactions on Communications, 66:2219–2234, 2018.
  * [40] C. Saha, H. S. Dhillon, N. Miyoshi, and J. G. Andrews.  Unified analysis of HetNets using Poisson cluster process under max-power association.  IEEE Transactions on Wireless Communications, 18:3797–3812, 2019\.
  * [41] V. Suryaprakash, J. Møller, and G. Fettweis.  On the modeling and analysis of heterogeneous radio access networks using a Poisson cluster process.  IEEE Transactions on Wireless Communications, 14:1035–1047, 2015\.
  * [42] U. Tanaka, Y. Ogata, and D. Stoyan.  Parameter estimation and model selection for Neyman-Scott point processes.  Biometrical Journal, 50:43–57, 2008.
  * [43] E. Turgut and M. C. Gursoy.  Downlink analysis in unmanned aerial vehicle (UAV) assisted cellular networks with clustered users.  IEEE Access, 6:36313–36324, 2018.
  * [44] L. Yang, T. J. Lim, J. Zhao, and M. Motani.  Modeling and analysis of HetNets with interference management using Poisson cluster process.  IEEE Transactions on Vehicular Technology, 70:12039–12054, 2021\.


# Neural-adjoint method for the inverse design of all-dielectric metasurfaces

Yang Deng  1 Simiao Ren  1 Kebin Fan  1 Jordan M. Malof  1 and Willie J.
Padilla1  1Department of Electrical and Computer Engineering, Duke University,
Durham, North Carolina 27708, USA<EMAIL_ADDRESS>
###### Abstract

All-dielectric metasurfaces exhibit exotic electromagnetic responses, similar
to those obtained with metal-based metamaterials. Research in all-dielectric
metasurfaces currently uses relatively simple unit-cell designs, but increased
geometrical complexity may yield even greater scattering states. Although
machine learning has recently been applied to the design of metasurfaces with
impressive results, the much more challenging task of finding a geometry that
yields the desired spectra remains largely unsolved. We explore and adapt a
recent deep learning approach – termed neural-adjoint – and find it is capable
of accurately and efficiently estimating complex geometry needed to yield a
targeted frequency-dependent scattering. We also show how the neural-adjoint
method can intelligently grow the design search space to include designs that
increasingly and accurately approximate the desired scattering response. The
neural-adjoint method is not restricted to the case demonstrated and may be
applied to plasmonics, photonic bandgap, and other structured material
systems.

††journal: oe††articletype: Research Article

##  1 Introduction

Electromagnetic metamaterials have been used to realize many exotic scattering
responses over the last two decades. Effects including negative refractive
index and cloaking have generated significant interest and have served to
drive the community [1, 2, 3, 4, 5]. A more applied, but still relevant
metamaterials achievement, is that of graded designs [6, 7]. It was shown
early on that unit-cells which form metamaterials may be designed with a
spatial dependence across a surface and / or volume, and various lensing
effects were shown that utilize spatial degrees of freedom. A principle of
gradient metasurfaces is that their scattering properties change slowly as a
function of spatial coordinate. Broadband metamaterial absorbers [8, 9], and
metamaterial spatial light modulators [10], also make use of dissimilar
neighboring unit-cells, however there is no such requirement to make neighbors
as alike as possible, and designs are simply cobbled together to achieve the
desired response. This fact highlights a design feature of sub-wavelength
metal-based metamaterials, i.e. their scattering response is primarily
connected to their geometry and – due to minimal neighbor interaction – the
unit-cell largely governs the electromagnetic properties of the array.
Conventional optimization techniques built-in to modern day electromagnetic
mode solvers are sufficient to achieve a desired response, and designs of
differing units may be assembled with little change to the overall response.

In principle, all-dielectric metasurface (ADM) unit-cells may also be used to
tesselate a surface in an arbitrary fashion similar to that achieved with
metal-based metamaterials. However, ADMs resonators are only marginally sub-
wavelength, and modes utilized are often not confined within their physical
bound – with an evanescent tail lying just outside their surface – resulting
in significant neighbor interaction. Further, surface modes related to the
periodicity are not restricted to nearest neighbors, and may persist over
several spatial periods. Despite the success of ADMs [11, 12], it is
conceivable that still-richer electromagnetic scattering can be achieved if
more complex geometries are employed. However, physical understanding for such
metasurfaces is poor – simple functional relationships, or even heuristic
guidance regarding super unit-cell geometry and final electromagnetic
properties is unavailable. The only contemporary means to estimate such
metasurface properties, given a candidate geometry, is simulation or
fabrication. However, given the vast space of potential designs and the speed
of conventional simulation and fabrication, it is completely infeasible to
iterate over designs in order to achieve a desired response.

A viable design alternative to numerical simulation for structured materials,
including metamaterials [13, 14, 15, 16, 17], photonic band gap [18, 19, 20,
21], and plasmonics studies [22, 23, 24], is deep learning. Deep neural
networks (DNNs) have successfully learned a forward mapping $s=f(g)$ between a
metasurface geometry $g$ and the resulting frequency dependent electromagnetic
scattering $s$, where $f$ is an unknown complex (e.g., highly nonlinear)
function. [25] A DNN – once trained – can effectively act as a high-speed
simulator that may be used to find the electromagnetic scattering of a
candidate geometry substantially faster (e.g., by six orders of magnitude
[16]) than conventional simulation.

##  2 Neural-Adjoint Method

While deep learning enables substantially faster evaluation of candidates,
given the vast number of possibilities for many problems (e.g., approx.
$10^{12}$ in our case), it is nonetheless difficult or impossible to iterate
over all or most candidate designs. From a design perspective, what would be
of greatest utility would be to instead specify a desired frequency dependent
electromagnetic response $s$, and have a model or solver compute a specific
approximate geometry $\hat{g}$, which yields $s$ (here $\hat{}$ denotes
approximate). The solution used to search for such an ideal geometry may be
cast as an inverse problem written as $\hat{g}=\hat{f}^{-1}(s)$. Hadamard
described a solution method as well-posed if it met three criteria – $H_{1}$:
existence (at least one solution), $H_{2}$: uniqueness (only one solution),
and $H_{3}$: stability (solution must depend continuously on $s$). A solution
approach is then ill-posed if one of the $H_{1}-H_{3}$ conditions fails [26].
Past work has shown that finding a specific ADM consisting of distinct
neighboring resonators in a super unit-cell geometry for a desired scattering
state is an ill-posed non-linear inverse problem and, in particular,
conditions $H_{1}$ and $H_{2}$ are not met [13].

Figure 1: (a) 3D illustration of an all-dielectric metasurface absorber. A
schematic of the all-dielectric metasurfaces super-cell is shown – top right.
(b) Training phase of the neural-adjoint method which uses fixed geometry
inputs and corresponding simulated spectrum. (c) The inference phase of the
neural-adjoint method, in which the target spectrum drives the search of the
optimal geometry. The geometry is shown in a green box here to indicate it is
the only trainable parameter.

Here we explore a recently-proposed deep inverse modeling approach, called the
neural-adjoint (NA) method [27], which outperformed many other methods for
solving ill-posed inverse problems. Here we adapt the NA method to solve a
challenging problem involving the design of ADMs geometries with 14 free
geometrical parameters; a much larger space than most recent work (e.g.,
$5-10$ free parameters in [18, 14, 16]). Furthermore, existing inverse methods
often involve complex training procedures, and ultimately produce sub-optimal
solutions. In contrast, the NA only requires that we train a single
conventional feed-forward neural network, and – as we show – appears to find
close approximations to the globally optimal solution within just one minute,
even for our complex ADM design problem. Furthermore, we demonstrate how the
NA method can be utilized to expand the design search space, essentially
providing a form of active learning that is specifically tailored to solve
inverse problems. We explore an example inverse problem, where a frequency
dependent infrared absorptivity ($A(\omega)=s$) is desired, which we would
like to achieve using an ADM consisting of a square array of 2$\times$2
resonators, each with an elliptical generatrix lying perpendicular to its
directrix (line of length $h$) – depicted in Fig. 1. Although past experience
and single unit-cell simulations guide us to choose approximate metasurface
dimensions which yield resonances in the desired spectral range, our inverse
problem is ill-posed, since we do not know if a solution exists, i.e.
condition $H_{1}$ is not met. Further, we do not meet the condition of
uniqueness, $H_{2}$, since many metasurface solutions (disparate geometries)
may yield the same spectra – at least within the accuracy of our chosen loss
metric, and numerical precision. The proposed ADMs consists of a geometry
space of fourteen parameters: [$h$, $p$, $r_{x_{1}}$–$r_{x_{4}}$,
$r_{y_{1}}$–$r_{y_{4}}$, $\theta_{1}$–$\theta_{4}$]. As shown in Fig. 1, $h$
is the height of all four elliptical resonators, $p$ is the periodicity of the
super cell, $r_{x_{1}}$ – $r_{x_{4}}$ and $r_{y_{1}}$ – $r_{y_{4}}$ are the
x-axis and y-axis radii of four elliptical resonators respectively, and
$\theta_{1}$ – $\theta_{4}$ are the rotational angles with respect to the
center axis of each elliptical resonator. All geometry values are randomly
sampled from the data grid shown in Table S1, which is included in the
Supplemental. The numerical simulations use SiC for the proposed ADMs, and the
Supplemental contains more details.

Figure 2: Forward model mean-squared-error (MSE) between DNN predicted spectra
and numerical simulations. The DNN result with (a) the best, (b) average, and
(c) below average performance. (d) Histogram detailing the distribution of MSE
for the entire validation set.

##  3 Results

The DNN forward model is trained by pair of 14 geometrical parameters and 2000
frequency points from 100-500 THz. Results of the forward model trained on 60k
$\\{g,s\\}$ pairs (approximately $5.75\times 10^{-8}$ of the entire geometry
space) are shown in Fig. 2 (a)-(c), where red curves are numerical
simulations, and blue curves are the DNN predictions. The absolute error
between simulation and prediction is plotted as the gray curve on the right
vertical axis of Fig. 2 (a)-(c). In Fig. 2 (d) we show the histogram of the
mean-squared-error (MSE) for the entire validation set, and find that 95% have
an MSE$\leq 2.49\times 10^{-3}$ (dashed gray line), and an average MSE of
$1.2\times 10^{-3}$. Having an accurate trained forward model we next turn
toward the inverse solution. The NA method [27] finds the optimal inverse
solution by fixing all the weights and biases of $\hat{f}$, and computing the
forward model’s gradient solely with respect to the input to the network
(i.e., the geometry), starting from randomly chosen values, denoted
$\hat{g}_{0}$. It is important to highlight that $\hat{f}$ is a closed-form
differentiable expression, and thus calculation of $\partial\hat{f}/\partial
g$ is trivial. Further, we can estimate the gradient of the input geometry
with respect to a loss function $\mathcal{L}$ that we are free to specify
(e.g., mean squared error). Therefore, if $s$ is our desired spectrum, and
$\hat{g}_{i}$ is our current best estimate of the metasurface geometry, we can
iteratively move along the loss surface to find a better solution using,

$\left.\hat{g}_{i+1}=\hat{g}_{i}+\alpha\frac{\partial\mathcal{L}(\hat{f}(\hat{g}_{i}),s)}{\partial
g}\right|_{g=\hat{g}_{i}}\\\ $ (1)

where $\alpha$ is the learning rate. The $\hat{g}_{i}$ can then be evaluated
iteratively until some convergence criteria is satisfied (e.g., $\mathcal{L}$
changes very little after each iteration). Because this is a gradient-based
procedure, it will only converge to a locally optimal solution. As a result,
the NA method prescribes that this search process be repeated $T$ times, each
time starting from a different randomly chosen value of $\hat{g}_{0}$. In
practice we find that we may run greater than $10^{4}$ $T$ initializations in
parallel with no cost in speed – only limited by available memory. Therefore
the NA method produces $T$ candidate designs, and we can choose the best
design (or several designs, if desired) by passing each design back into
$\hat{f}$ and evaluating their similarity to the target scattering properties.

Figure 3: Neural adjoint inverse results for prediction of target spectra that
exist within the geometry space. Example spectra from geometry predictions
that have (a) the best, (b) average, and (c)-(d) below average performance
compared to the target spectrum.

Notably, while the user can specify many different choices for $\mathcal{L}$,
the NA method prescribes that a so-called "boundary loss", $\mathcal{L}_{bnd}$
should be added to any user-chosen loss [27] and is given by:
$\mathcal{L}_{bnd}=ReLU(|\hat{g}-\mu_{g}|-\frac{1}{2}R_{g})$, where rectified
linear unit (ReLU) is the activation function, $\mu_{g}$ is the mean of the
geometry training data, and $R_{g}$ is its range. This boundary loss punishes
the inference process with increasing loss if the geometry search process
steps out of the space of the training data, where the forward model may
produce inaccurate estimates of scattering parameters. In our experiments we
use the following total loss:
$\mathcal{L}=(s-\hat{f}(\hat{g}))^{2}+\mathcal{L}_{bnd}$. As an initial test
of the NA inverse method, we feed in frequency dependent absorptivities
$A(\omega)$ where we know apriori that a solution $g$ exists, i.e. $s$ is a
numerical simulation from which $\hat{f}$ originates from. In Fig. 3 we show
results of the NA inverse method and each sub-plot shows characteristic
results of (a) the best results, (b) average results, and (c) and (d) below
average results, all based on MSE. In all of the these examples, the NA method
identifies a close approximation to the correct solution. This is impressive
given the complexity of the spectra present in Fig. 3. We suspect the small
remaining errors in the predicted design are due largely to the limited
precision of gradient descent as it nears solutions (i.e., minima) in the
error space; due to the non-zero learning rate in Eq. 1, it cannot converge to
the exact minimum point. We note however that learning rate can be gradually
reduced during the search process, at the cost of additional computation time,
until a solution of desired precision is obtained.

We next turn to the significantly more challenging task of applying the NA
method to a spectra where we are unaware if a solution exists within the
chosen geometrical parameter space, i.e. criterion $H_{1}$ for inverse
problems is violated. We chose the frequency dependent external quantum
efficiency (EQE) of gallium antimonide (GaSb) as $s$, shown as the gray curves
of Fig. 4. The metasurface will operate at elevated temperature, and thus we
consider the so-called graybody spectra – also termed the spectral exitance
$M_{e,\nu}(T)$ – which is given by the blackbody radiation curve times the
absorptivity. We keep the top 16,000 neural adjoint solutions (spectra) and
determine $M_{e,\nu}(T)$ for each of these at 100 temperatures between 1500
and 2500 k – a total of $1.6\times 10^{6}$ candidates. The shape of the EQE
curve differs significantly from typical spectra we see in our geometry space
(Fig. 3). None-the-less we find a best solution resulting from the NA method
at a temperature of $T=2100k$ that consists of a geometry of [$h=0.566$,
$p=1.440$, $r_{x_{1}}=0.180$, $r_{x_{2}}=0.155$, $r_{x_{3}}=0.214$,
$r_{x_{4}}=0.278$, $r_{y_{1}}=0.285$, $r_{y_{2}}=0.253$, $r_{y_{3}}=0.146$,
$r_{y_{4}}=0.256$, $\theta_{1}=-0.901^{\circ}$, $\theta_{2}=-20.677^{\circ}$,
$\theta_{3}=-37.982^{\circ}$, and $\theta_{4}=39.046^{\circ}$]. The spectral
exitance resulting from this geometry – calculated from $\hat{f}$ – is shown
as the blue curve of Fig. 4 (b), and we find an MSE, compared to the EQE of
GaSb, of 1.06$\times$10-2. We also apply a weighting function
$W(\nu)=1\chi_{[100,275]}+0\chi_{(275,300]}$ on the MSE forcing the NA method
to focus on the region of interest for energy harvesting purposes. To verify
our neural adjoint results, we numerically simulate the predicted geometry and
plot the resulting $M_{e,\nu}(T)$ in Fig. 5 (a) as the red curve – again
compared to the EQE of GaSb (gray curve). As can be seen, the simulated curve
has many relatively sharp peaks that are not present in Fig. 4(b). This is
because, as noted earlier, $\hat{f}$ does not perfectly match the numerical
simulator, and this will introduce errors in the design process. Thus since
the NA method relies on $\hat{f}$ to search for designs, it is also limited by
the accuracy of the forward model estimate. We also found that since $\hat{f}$
is trained from geometries constrained to a grid, the discrepancy between NA
solutions and numerical simulation arises because NA solutions are not
confined to the grid, where our model is most accurate. None-the-less we find
that our simulated $M_{e,\nu}(T)$ achieves an MSE of 1.65$\times 10^{-2}$, as
shown in Fig. 5 (a).

Figure 4: Neural adjoint inverse results for matching the EQE of GaSb. The
target spectrum (gray), DNN prediction for $M_{e,\nu}(T=2100K)$ (blue), for
the expanded geometry space (a), and the original geometry space (b) explored.
Neural adjoint predictions versus $p$ and $h$ for (c) the expanded geometry
space, and (d) the original space. Symbol colors indicate MSE, with values
given in the colorbar.

Another major obstacle is that our design search space does not contain a
geometry that can realize our target spectrum (i.e., Hadamard’s criteria
$H_{1}$). This is suggested by the fact that our best NA solution, shown in
Fig. 4(a), still does not match our target spectrum. However, we can use the
NA output to identify where we should expand our search space so that it will
include better designs. We can do this by visualizing the error of all inverse
solutions returned by NA, and looking for trends e.g., we may find that all
the best solutions are bunched up against some edge of our initial search
space, suggesting that expanding along that dimension may yield better
results. However, since we have a 14 dimensional design space, we are unable
to easily visualize these data. To address this problem, we use the Uniform
Manifold Approximation and Projection (UMAP) [28], which is a type of
dimensionality reduction method permitting us to visualize the distribution of
our inverse solutions performance in 2D, so that we may more easily identify
patterns. From this initial investigation with UMAP in Fig. S1 – shown in
Supplemental – we find that our best NA inverse solutions are limited by
height. Shown in Fig. 4 (d) are NA solutions as a function of height and
periodicity color mapped by corresponding MSE values. It is evident that not
only are our best solutions grouped at the maximum height allowed in our
geometrical space, $h=0.6$ $\mu$m, but also that the solutions improve as a
function of height. Encouraged by these results, we expanded our original
geometry space to include increased height values from 0.6 to 0.75 $\mu$m, by
simulating an additional 24k $\\{g,s\\}$ pairs. The NA model now trained on
the expanded geometrical space indeed finds an improved solution, shown as the
blue curve in Fig. 4 (a), where we realize an MSE that is reduced by a factor
of 2.7. The simulated red curve in Fig. 5 (b) further validates the result
that the MSE of numerical simulations is also reduced – here by a factor of
4.8. A plot of the 1000 best NA solutions in the expanded geometrical space
shown in Fig. 4 (c), however, indicate that we may be able to make continued
improvements, since we still have a gradient pushing for greater heights –
although the periodicity seems to be honing in on a value of 1.2 $\mu$m.

Figure 5: Numerically simulated $M_{e,\nu}(T=2100K)$ (red) of the optimal
geometry predicted by the neural-adjoint method to match the EQE of GaSb
(gray), for the original geometry space (a), and the expanded geometry space
(b).

##  4 Conclusion

We have adapted the neural-adjoint inverse design method [27] to accurately
predict the high-dimensional all-dielectric metasurface geometry needed to
produce a targeted infrared absorptivity spectra. When the geometry needed to
produce a desired spectrum lies outside of the bounds, the NA method appears
to find the best possible solution within the permissible search space. Unlike
other adjoint inverse approaches [29], the NA method does not require any
domain knowledge specific to the problem. In the event that the inverse
solution does not exist in the parameter space explored, NA may be used to
guide one to a better solution, through expansion of the design parameter
search space. This may help to reduce the initial required number of numerical
simulations, and to instead use NA guided simulation exploration. With its
exceptional computational speed, high accuracy, and potential use in active
learning that is explored here, the neural-adjoint method has an impressive
future in not only ADMs thermal emitter but also any ADMs inverse problems.
The NA method is not restricted to the case presented, but may be applied to
many other systems including photonic band gap and plasmonics.

## References

  * [1] D. R. Smith, W. J. Padilla, D. Vier, S. C. Nemat-Nasser, and S. Schultz, “Composite medium with simultaneously negative permeability and permittivity,” Physical review letters 84, 4184 (2000).
  * [2] J. B. Pendry, “Negative refraction makes a perfect lens,” Physical review letters 85, 3966 (2000).
  * [3] R. A. Shelby, D. R. Smith, and S. Schultz, “Experimental verification of a negative index of refraction,” science 292, 77–79 (2001).
  * [4] D. R. Smith, J. B. Pendry, and M. C. Wiltshire, “Metamaterials and negative refractive index,” Science 305, 788–792 (2004).
  * [5] D. Schurig, J. J. Mock, B. Justice, S. A. Cummer, J. B. Pendry, A. F. Starr, and D. R. Smith, “Metamaterial electromagnetic cloak at microwave frequencies,” Science 314, 977–980 (2006).
  * [6] R. Greegor, C. Parazzoli, J. Nielsen, M. Thompson, M. Tanielian, and D. Smith, “Simulation and testing of a graded negative index of refraction lens,” Applied Physics Letters 87, 091114 (2005).
  * [7] D. R. Smith, J. J. Mock, A. Starr, and D. Schurig, “Gradient index metamaterials,” Physical Review E 71, 036609 (2005).
  * [8] N. I. Landy, S. Sajuyigbe, J. J. Mock, D. R. Smith, and W. J. Padilla, “Perfect metamaterial absorber,” Physical review letters 100, 207402 (2008).
  * [9] Y. Avitzour, Y. A. Urzhumov, and G. Shvets, “Wide-angle infrared absorber based on a negative-index plasmonic metamaterial,” Physical Review B 79, 045131 (2009).
  * [10] D. Shrekenhamer, J. Montoya, S. Krishna, and W. J. Padilla, “Four-color metamaterial absorber thz spatial light modulator,” Advanced Optical Materials 1, 905–909 (2013).
  * [11] J. Sautter, I. Staude, M. Decker, E. Rusak, D. N. Neshev, I. Brener, and Y. S. Kivshar, “Active tuning of all-dielectric metasurfaces,” ACS nano 9, 4308–4315 (2015).
  * [12] D. Headland, S. Nirantar, W. Withayachumnankul, P. Gutruf, D. Abbott, M. Bhaskaran, C. Fumeaux, and S. Sriram, “Terahertz magnetic mirror realized with dielectric resonator antennas,” Advanced Materials 27, 7137–7144 (2015).
  * [13] Z. Liu, D. Zhu, S. P. Rodrigues, K.-T. Lee, and W. Cai, “Generative model for the inverse design of metasurfaces,” Nano letters 18, 6570–6576 (2018).
  * [14] W. Ma, F. Cheng, and Y. Liu, “Deep-learning-enabled on-demand design of chiral metamaterials,” ACS nano 12, 6326–6334 (2018).
  * [15] J. Jiang and J. A. Fan, “Global optimization of dielectric metasurfaces using a physics-driven neural network,” Nano letters 19, 5366–5372 (2019).
  * [16] C. C. Nadell, B. Huang, J. M. Malof, and W. J. Padilla, “Deep learning for accelerated all-dielectric metasurface design,” Optics express 27, 27523–27535 (2019).
  * [17] S. D. Campbell, D. Sell, R. P. Jenkins, E. B. Whiting, J. A. Fan, and D. H. Werner, “Review of numerical optimization techniques for meta-device design,” Optical Materials Express 9, 1842–1863 (2019).
  * [18] J. Peurifoy, Y. Shen, L. Jing, Y. Yang, F. Cano-Renteria, B. G. DeLacy, J. D. Joannopoulos, M. Tegmark, and M. Soljačić, “Nanophotonic particle simulation and inverse design using artificial neural networks,” Science advances 4, eaar4206 (2018).
  * [19] D. Liu, Y. Tan, E. Khoram, and Z. Yu, “Training deep neural networks for the inverse design of nanophotonic structures,” ACS Photonics 5, 1365–1369 (2018).
  * [20] Y. Kiarashinejad, M. Zandehshahvar, S. Abdollahramezani, O. Hemmatyar, R. Pourabolghasem, and A. Adibi, “Knowledge discovery in nanophotonics using geometric deep learning,” Advanced Intelligent Systems 2, 1900132 (2020).
  * [21] L. Huang, L. Xu, and A. E. Miroshnichenko, “Deep learning enabled nanophotonics,” in _Advances in Deep Learning,_ (IntechOpen, 2020).
  * [22] S. Inampudi and H. Mosallaei, “Neural network based design of metagratings,” Applied Physics Letters 112, 241102 (2018).
  * [23] J. Jiang, D. Sell, S. Hoyer, J. Hickey, J. Yang, and J. A. Fan, “Free-form diffractive metagrating design based on generative adversarial networks,” ACS nano 13, 8872–8878 (2019).
  * [24] P. R. Wiecha and O. L. Muskens, “Deep learning meets nanophotonics: A generalized accurate predictor for near fields and far fields of arbitrary 3d nanostructures,” Nano Letters 20, 329–338 (2019).
  * [25] J. Gareth, W. Daniela, H. Trevor, and T. Robert, _An introduction to statistical learning: with applications in R_ (Spinger, 2013).
  * [26] J. Mueller, _Linear and nonlinear inverse problems with practical applications_ (Society for Industrial and Applied Mathematics, Philadelphia, 2012).
  * [27] S. Ren, W. J. Padilla, and J. Malof, “Benchmarking deep inverse models over time, and the neural-adjoint method,” (2020).
  * [28] L. McInnes, J. Healy, and J. Melville, “Umap: Uniform manifold approximation and projection for dimension reduction,” arXiv preprint arXiv:1802.03426 (2018).
  * [29] C. M. Lalau-Keraly, S. Bhargava, O. D. Miller, and E. Yablonovitch, “Adjoint shape optimization applied to electromagnetic design,” Optics express 21, 21693–21701 (2013).