# Characterizing First Arrival Position Channels: Noise Distribution and
Capacity Analysis

Yen-Chi Lee∗, Yun-Feng Lo, Jen-Ming Wu, and Min-Hsiu Hsieh,  This manuscript
refines and extends the concepts from our unpublished preprints [1, 2],
presenting a more comprehensive and integrated narrative of the research. Yen-
Chi Lee, Jen-Ming Wu and Min-Hsiu Hsieh are all affiliated with Hon Hai
(Foxconn) Research Institute, Taipei, Taiwan. Emails<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>Lo is affiliated
with Georgia Institute of Technology, Atlanta, Georgia, United States. Email:
<EMAIL_ADDRESS>
###### Abstract

This paper introduces a novel mathematical model for Molecular Communication
(MC) systems, utilizing First Arrival Position (FAP) as a fundamental mode of
information transmission. We address two critical challenges: the
characterization of FAP density and the establishment of capacity bounds for
channels with vertically-drifted FAP. Our method relate macroscopic Partial
Differential Equation (PDE) models to microscopic Stochastic Differential
Equation (SDE) models, resulting in a precise expression that links FAP
density with elliptic-type Green’s function. This formula is distinguished by
its wide applicability across any spatial dimensions, any drift directions,
and various receiver geometries. We demonstrate the practicality of our model
through case studies: 2D and 3D planar receivers. The accuracy of our formula
is also validated by particle-based simulations. Advancing further, the
explicit FAP density forms enable us to establish closed-form upper and lower
bounds for the capacity of vertically-drifted FAP channels under a second-
moment constraint, significantly advancing the understanding of FAP channels
in MC systems.

###### Index Terms:

molecular communication (MC), partial differential equation (PDE), stochastic
analysis, first arrival position (FAP), channel capacity

##  I Introduction

The challenge of transmitting information over a distance while maintaining
the fidelity of the information has been present throughout human history,
from ancient times to the present day. While modern communication systems have
addressed this problem using electromagnetic signals, such techniques
encounter limitations in nano-scale applications due to restrictions on
wavelength, antenna size, and energy [3, 4]. To circumvent these restrictions,
molecular signals have been proposed as a viable alternative, particularly for
nanonetwork applications [5, 6]. This has given birth to a new field of
communication called molecular communication (MC).

###  I-A Background on Molecular Communication

In MC systems, information is conveyed through tiny molecules, referred to as
message molecules (MMs), that act as information carriers [7]. To transport
these MMs through the physical channel toward a receiver, a propagation
mechanism is necessary. Various forms of propagation mechanisms exist, such as
diffusion-based [8], flow-based [9], or an engineered transport system like
molecular motors [10]. Among these mechanisms, diffusion-based propagation
combined with advection (or somtimes chemical reaction [11]), has been the
prevalent approach considered in the MC literature thus far. In this study, we
will focus on _diffusion-based_ (also known as _diffusive_) MC systems. Note
that we use the terms “molecules” and “particles” interchangeably to refer to
the MMs, as their shape is not relevant to our discussion.

In diffusive MC systems, information can be transmitted by modulating various
physical properties of the MMs [12]. Once these signaling molecules reach the
vicinity of the receiver, they can be detected and processed to extract the
required information for detection and decoding [13]. The reception mechanism
of an MC receiver can be classified into two categories: i) passive reception,
and ii) active reception [14]. A widely used simple active reception model is
the fully-absorbing receiver [15]. For our following discussion, we will focus
on the fully-absorbing receiver, which often referred to as an “absorbing
receiver” for brevity. Note that this assumption is also commonly employed in
recent research on molecular MIMO [16] and molecular index modulation [17,
18].

Absorbing receivers harness the information extracted during the initial
contact between MMs and the receiver boundary. The stochastic properties of
MMs can be captured by first-hitting process (see [16]). The modulation
schemes [12] applicable to absorbing receivers can be roughly categorized into
three types: timing-based [19], position-based [20, 21], and joint timing-
position-based schemes [22]. In this study, we delve deep into position-based
information type and explore a novel channel family known as First Arrival
Position (FAP) channels.

###  I-B Main Results of this Work

In this paper, we extend upon the materials introduced in our previously
unpublished preprints [1, 2]. We have dedicated substantial effort to refine
and integrate the fundamental ideas and discoveries from these earlier
versions. In addition, we have enriched the content to bolster the motivation
behind this study and have proposed potential applications (Section I-B).
Furthermore, in this journal version, we have performed particle-based
simulations to verify the validity of our derived FAP density formula, and
performed graphical illustrations of our derived capacity bounds with respect
to key parameters. Last but not least, we have organized our key results about
capacity in a theorem-style format (see Section V) to facilitate easy
comprehension for readers.

The main contributions of this paper are as follows:

  * •

Modeling FAP-type MC systems: By imposing some ideal assumptions, we model
FAP-type MC systems as additive vector channels (see Eq. (4)), and the
particle’s evolution is captured by an Itô diffusion. Please see Fig. 1 and
Fig. 2.

  * •

FAP Channel Characterization: We introduce an unified framework (Section III)
capable of calculating the FAP density in diffusive MC systems. This framework
is applicable to MC systems of arbitrary spatial dimension, arbitrary drift
direction, and various receiver shapes including the commonly used sphere and
plane shapes.

  * •

Upper and Lower Bounds on Capacity: Under a second-moment constraint, we
derive the upper and lower bounds for VDFAP channel capacity. To the best of
our knowledge, this is the first work in MC to analytically explore the
fundamental data transmission rate (i.e., channel capacity) for FAP-type MC
channels.

###  I-C Organization of this Paper

The remainder of this paper is structured as follows. In Section II, we
present the physical model and system assumptions used in our analysis.
Section III introduces our novel unified approach for finding FAP density. In
Section IV, we examine the characteristic function of the VDFAP noise
distribution, as well as its weak stability property. We then use the derived
formulas for moments and the weak stability property to establish lower and
upper bounds for the capacity of VDFAP channels in Section V. Numerical
results are provided in Section VI. Finally, we summarize our results and
provide concluding remarks in Section VII.

##  II System Model

###  II-A Physical Model

There are two different perspectives for modeling diffusion phenomenon in MC
systems. The macroscopic viewpoint utilizes the diffusion equation, also known
as the heat equation (both are partial differential equation (PDE) models), to
depict the evolution of the concentration $c(\mathbf{x},t)$ of message
molecules. Based on Fick’s law of diffusion [13], the temporal evolution of
$c(\mathbf{x},t)$, assuming no chemical reactions, can be described as [13,
eq. (17)]:

$\partial_{t}c\left(\mathbf{x},t;\mathbf{x}_{0}\right)+\nabla\cdot\left(\mathbf{v}(\mathbf{x},t)c\left(\mathbf{x},t;\mathbf{x}_{0}\right)\right)=D\Delta
c\left(\mathbf{x},t;\mathbf{x}_{0}\right),$ (1)

where $\mathbf{x}_{0}$ denotes the emission point, $\mathbf{v}$ denotes the
background velocity field, $\nabla\cdot$ represents the divergence operator,
and $\Delta\equiv\nabla^{2}$ is the Laplacian operator[23, Section 12.11].
Here, $D$ is the diffusion coefficient [13], which we assume to be constant
with respect to space and time. The value of $D$ is dependent on factors such
as temperature, fluid viscosity, and the molecule’s Stokes radius, as outlined
in [7]. Note that the shorthand notations $\partial_{t}$ stands for
$\frac{\partial}{\partial t}$ and $\partial_{x_{j}}$ stands for
$\frac{\partial}{\partial{x_{j}}}$.

In contrast, from a microscopic perspective, a suitable model for an
individual trajectory $\mathbf{X}_{t}$ of a message molecule is an Itô
diffusion process[24, Section 7.1]. We will assume that the parameters in the
diffusion channel vary slowly over time, allowing the MC channel to be viewed
as an approximately time-invariant system [13]. This assumption corresponds to
the time-homogeneous Itô diffusion model. Namely, a _time-homogeneous Itô
diffusion_ process is a stochastic process that satisfies a stochastic
differential equation (SDE) of the form [25, Section 6.2.1]:

$\ \mathrm{d}\mathbf{X}_{t}=\mathbf{b}(\mathbf{X}_{t})\
\mathrm{d}t+\sigma(\mathbf{X}_{t})\ \mathrm{d}\mathbf{B}_{t},$ (2)

where $\mathbf{B}_{t}$ is a D-dimensional standard Brownian motion. The first
term, $\mathbf{b}(\mathbf{X}_{t})\ \mathrm{d}t$, is the deterministic
component, which describes the drift effect induced by an external potential
field. The second term, $\sigma(\mathbf{X}_{t})\ \mathrm{d}\mathbf{B}_{t}$,
represents the random fluctuations caused by the constant bombardment of
background molecules. Notice that $\mathbf{b}(\cdot)$ is a D-dimensional
vector, and $\sigma(\cdot)$ is a $\text{D}\times\text{D}$ matrix. Both are
functions of the spatial variable $\mathbf{x}$.

These two models, PDE and SDE, possess an important relation. The microscopic
perspective via the SDE (2) can be connected to the macroscopic perspective
via the PDE (1) through the Fokker-Planck equation, also known as the
Kolmogorov forward equation[24, Chapter 8]. Specifically, for an Itô diffusion
process $\mathbf{X}_{t}$ satisfying (2) for $t>0$ and the initial condition
$\mathbf{X}_{0}=\mathbf{x}_{0}$, the probability density of $\mathbf{X}_{t}$,
denoted by $p(\mathbf{x},t;\mathbf{x}_{0})$, satisfies the Kolmogorov forward
equation:

$\displaystyle\begin{split}\partial_{t}p(\mathbf{x},t;\mathbf{x}_{0})&=-\sum_{j=1}^{\text{D}}\partial_{x_{j}}\left[\mathbf{b}(\mathbf{x})p(\mathbf{x},t;\mathbf{x}_{0})\right]\\\
&\quad+\sum_{j=1}^{\text{D}}\sum_{k=1}^{\text{D}}\partial_{x_{j}}\partial_{x_{k}}\left[D_{jk}(\mathbf{x})p(\mathbf{x},t;\mathbf{x}_{0})\right]\end{split}$
(3)

where $b_{j}(\mathbf{x})$ denotes the $j$-th component of
$\mathbf{b}(\mathbf{x})$ and $D_{jk}(\mathbf{x})$ is the $(j,k)$-th entry of
the matrix
$D(\mathbf{x}):=\tfrac{1}{2}\sigma(\mathbf{x})\sigma(\mathbf{x})^{\mathsf{\scriptstyle{T}}}$.
To establish a connection between (3) and (1), we make two assumptions that
are common in the MC literature [7, 13]: i) the background velocity field is
constant with respect to time, i.e., assuming
$\mathbf{v}(\mathbf{x},t)=\mathbf{v}(\mathbf{x})$ in the macroscopic PDE
model; ii) the diffusion coefficient $D$ is constant throughout space and
time. That is, we assume the matrix $\sigma(\mathbf{x})$ to be a scalar
multiple of the identity matrix, where the scalar is a constant. To simplify
notations, we also denote this scalar constant by $\sigma$. Under these two
assumptions, we can establish equivalence between (1) and (3) by identifying
the following quantities: $c=p$, $\mathbf{v}(\cdot)=\mathbf{b}(\cdot)$ and
$D=\sigma^{2}/2$.111Note that since (1) and (3) are _linear_ PDEs, any scalar
multiple of a solution is still a solution. In concentration-encoded MC where
MMs move independently, the relation $c=N_{\text{total}}\cdot p$ holds
($N_{\text{total}}$ denotes the total number of MMs released) [13].

###  II-B Communication Model

Figure 1:  This figure illustrates a 2D FAP channel with line-shaped Tx and
Rx, where the Tx is assumed to be transparent, allowing particles to move
through it without experiencing any force. The emission point is located on
the Tx line. Note that the two red particles are highlighted as exemplary
particles to indicate the emission position $(x_{1},x_{2})$ and the hitting
position $(\xi,\eta)$, respectively.

Figure 2:  This figure illustrates a 3D FAP channel with plane-shaped Tx and
Rx, where the Tx is assumed to be transparent, allowing particles to move
through it without experiencing any force. The emission point is located on
the Tx plane. Note that the two red particles are highlighted as exemplary
particles to indicate the current position $(x_{1},x_{2},x_{3})$ and the
hitting position $(\xi,\eta,\zeta)$, respectively.

In this paper, we investigate a diffusive MC system, which employs the
emission positions of MMs in a D-dimensional fluid medium to convey
information. In practice, $\text{D}=2,3$ are commonly selected. See Fig. 1 and
Fig. 2. We consider a scenario where there is a constant velocity field [9] in
the ambient space. This corresponds to further assuming
$\mathbf{v}(\cdot)=\mathbf{b}(\cdot)=\mathbf{v}$, where $\mathbf{v}$ is a
constant vector. We call $\mathbf{v}$ the drift (vector) hereafter.
Additionally, we assume all the MMs are of the same type [26]. The diffusion
effects of MMs are captured by a diffusion coefficient $D=\sigma^{2}/2$, which
is assumed to be a constant throughout space and time [7]. We call this the
isotropic assumption on the fluid medium.

The MC system we considered comprises a transmitter (Tx) and a receiver (Rx),
which are modeled as parallel $(\text{D}-1)$-dimensional hyperplanes
(hereafter referred to as _planes_) separated by a distance $\lambda>0$.
Without loss of generality, we can set the Tx plane at coordinate
$x_{\text{D}}=\lambda$ and the Rx plane at coordinate $x_{\text{D}}=0$.

In this study, we assume some simplified assumptions to our FAP MC system.

  * (A1)

The movement of every MM is independent. (Most MC works adopt this
assumption.)

  * (A2)

Upon first arrival at the receiver plane, MMs are captured and removed from
the system. (Most MC works using an fully-absorbing receiver adopt this
assumption.)

  * (A3)

The receiver can perfectly measure the first arrival positions of the MMs.
(Note: We add this assumption in order to have a tractable analysis for
channel capacity. However, the main results in Section III and Section IV do
not depend on this assumption.)

Although the ambient space is D-dimensional, the position information is
$d$-dimensional, where $d:=\text{D}-1$ for brevity and consistency throughout
this paper. The first arrival position $\bm{Y}$ of a molecule released at
position $\bm{X}$ can be expressed as

$\bm{Y}=\bm{X}+\bm{N},$ (4)

where $\bm{N}$ denotes the deviation (i.e., noise) of the first arrival
position. The probability density function of $\bm{N}$ is defined to be the
_FAP density_ [20, 21], and the additive channel (4) is defined to be the _FAP
channel_.

Please note that the realization of the random vector $\bm{N}$ is represented
by the lowercase letter $\bm{n}$. In instances where $\bm{n}$ resides in a
$\mathbb{R}^{d}$ space, we express it as $\bm{n}=(n_{1},n_{2},\ldots,n_{d})$,
which can also be equivalently denoted as $[n_{1},n_{2},\ldots,n_{d}]$. The
same convention will also be applied to other random vectors, such as $\bm{X}$
and $\bm{Y}$.

###  II-C The FAP Noise Distribution

In [21, 20], the 2D and 3D FAP densities were separately derived based on
different methods. By the framework proposed in this paper, we can unify the
derivations in different scenarios into a standard procedure. Our new
derivation, which recovers the FAP formulas in [21, 20], is presented in
Appendix D. Notice that the special function $K_{\frac{3}{2}}$ has an
elementary form[27, Section 8.432]:

$K_{\frac{3}{2}}(x)=\sqrt{\dfrac{\pi}{2}}e^{-x}\left(\dfrac{1+x}{x^{\frac{3}{2}}}\right),$
(5)

where $K_{\nu}(\cdot)$ denotes the order-$\nu$ modified Bessel function of the
second kind [28]. By denoting $\bm{n}=\bm{y}-\bm{x}$, we can merge the two FAP
densities, for $\text{D}=2$ and $\text{D}=3$, into a single expression:

$\displaystyle\begin{split}f^{(d)}_{\bm{N}}(\bm{n};\mathbf{u},\lambda)=\
&2\lambda\frac{\|\mathbf{u}\|^{\tfrac{d+1}{2}}}{(2\pi)^{\tfrac{d+1}{2}}}\cdot
e^{\bm{u}_{\text{par}}^{\intercal}\bm{n}-u_{\text{D}}\lambda}\\\
\cdot&\frac{K_{\tfrac{d+1}{2}}\left(\|\mathbf{u}\|\sqrt{\norm{\bm{n}}^{2}+\lambda^{2}}\right)}{\left(\sqrt{\norm{\bm{n}}^{2}+\lambda^{2}}\right)^{\tfrac{d+1}{2}}},\text{\
for }\bm{n}\in\mathbb{R}^{d}\end{split}$ (6)

where $\|\cdot\|$ denotes the Euclidean norm. We present here some special
notational conventions used throughout this paper. We write the D-dimensional
drift vector
$\mathbf{v}=[v_{1},\ldots,v_{\text{D}}]^{\mathsf{\scriptstyle{T}}}$ as
$\mathbf{v}=[\bm{v}_{\text{par}}^{\mathsf{\scriptstyle{T}}},v_{\text{D}}]^{\mathsf{\scriptstyle{T}}}$,
where $\bm{v}_{\text{par}}=[v_{1},\ldots,v_{d}]^{\mathsf{\scriptstyle{T}}}$
contains the drift components parallel to the Tx and Rx planes, and
$v_{\text{D}}$ is the drift component perpendicular to the Tx and Rx planes.
As a rule of thumb, we use bold fonts to stand for (column) vectors; slanted
vectors are $d$-dimensional, while non-slanted vectors are D-dimensional. The
operator $(\cdot)^{\mathsf{\scriptstyle{T}}}$ represents transposition. We
also introduce the notation $\mathbf{u}=\frac{\mathbf{v}}{\sigma^{2}}$ ,
$\bm{u}_{\text{par}}=\frac{\bm{v}_{\text{par}}}{\sigma^{2}}$ and
$u_{j}=\frac{v_{j}}{\sigma^{2}}$ for each $j\in\\{1,\ldots,\text{D}\\}$. These
$\\{u_{j}\\}_{j=1}^{\text{D}}$ can be interpreted as _normalized drift_ : they
are the ratio of the drift “strength” $v_{j}$ over the diffusion “strength”
$\sigma^{2}$. Since $v_{j}$ has the unit of $\mu\mathrm{m}/s$ and
$\sigma^{2}=D/2$ has the unit of $\mu\mathrm{m}^{2}/s$, the normalized drift
$u_{j}$ has the unit of $\mu\mathrm{m}^{-1}$.

For the channel capacity investigation, we focus on a sub-family of FAP
(noise) distributions with two specific properties:

  * i)

The parallel drift components are zero, i.e., $\bm{v}_{\text{par}}=\bm{0}$.

  * ii)

The _vertical drift_ (VD) component points from the Tx to the Rx, i.e.,
$v_{\text{D}}<0$. (Intuitively, this vertical drift helps the transmission of
information.)

We name this sub-family as vertically-drifted first arrival position (VDFAP)
distributions because only the vertical drift component is non-zero. We can
set $\bm{u}_{\text{par}}=\bm{0}$ so that $\|\mathbf{u}\|$ in (6) becomes
$|u_{\text{D}}|=-u_{\text{D}}$. Abbreviating $u_{\text{D}}$ as $u$, the VDFAP
densities can be expressed as:

$\displaystyle\begin{split}f^{(d)}_{\bm{N}}(\bm{n};u,\lambda)=\
&2\lambda\left(\frac{|u|}{\sqrt{2\pi}}\right)^{d+1}\cdot e^{\lambda|u|}\\\
\cdot&\frac{K_{\tfrac{d+1}{2}}\left(|u|\sqrt{\norm{\bm{n}}^{2}+\lambda^{2}}\right)}{\left(|u|\sqrt{\norm{\bm{n}}^{2}+\lambda^{2}}\right)^{\tfrac{d+1}{2}}},\text{\
for }\bm{n}\in\mathbb{R}^{d}.\end{split}$ (7)

We also denote the VDFAP densities by $f^{(d)}_{\bm{N}}(\bm{n})$ when the
parameters $u,\lambda$ are clear from the context. As a shorthand, we denote
$\bm{N}\sim\text{VDFAP}^{(d)}(u,\lambda)$ if a $d$-dimensional random vector
$\bm{N}$ follows the VDFAP distribution defined by (7) with parameters $u<0$
and $\lambda>0$. Note that we refer to a random variable $X\in\mathbb{R}$ as a
“one-dimensional random vector” to maintain consistency in terminology.

##  III A Novel Framework for Efficiently Determining FAP Density

In this section, we propose a novel and unified framework for determining the
FAP density. To achieve this goal, we introduce several stochastic analytical
tools, such as the generator of diffusion semigroup and Dynkin’s formula[25,
Section 9.3].

Our new method for calculating the FAP density can be divided into three steps
as summarized in Table I. Specifically, the FAP density can be computed as
$f_{\mathbf{Y}|\mathbf{X}}(\mathbf{y}|\mathbf{x})=\big{|}\frac{\partial
G(\mathbf{x},\mathbf{y})}{\partial\mathbf{n}_{\mathbf{y}}}\big{|}$ once the
elliptic Green’s function $G(\mathbf{x},\mathbf{y})$ is known. Note that the
“previous method” in Table I refers to those methods that tackle the diffusion
equation (1) directly, such as [21].

For the previous method, a parabolic-type PDE (which is time-dependent) is
used, while the new method removes the time-dependency by considering the
generator of the diffusion semigroup. As a consequence, using the previous
method often requires doing integration with respect to time, which is
unnecessary while applying our new method.

Table I: Comparison Between Previous and New Method of Calculating FAP Density |  Previous Method |  New Method (we propose in this work)
---|---|---
Step 1 |  Calculates the free space Green’s function for parabolic-type[29, Chapter 3] PDE. |  Removes time-dependency via considering the generator of the diffusion semigroup.
Step 2 |  Accommodates the free space solution to absorbing boundary conditions (e.g., using the image method). |  Solves elliptic-type[29, Chapter 8] boundary value problem with appropriate boundary conditions.
Step 3 |  Calculates the flux density (see [21]) and does integration with respect to time. |  Obtains FAP density directly from the Green’s function $G(\mathbf{x},\mathbf{y})$. No time-integration is required.

In the remainder of this section, we will establish the validity of the
proposed FAP density computing framework.

###  III-A Infinitestimal Generator of Itô Diffusion

Recall that a D-dimensional Itô diffusion is a stochastic process that
satisfies (2): $\ \mathrm{d}\mathbf{X}_{t}=\mathbf{b}(\mathbf{X}_{t})\
\mathrm{d}t+\sigma(\mathbf{X}_{t})\ \mathrm{d}\mathbf{B}_{t}.$ In this
subsection, we will first assume that $\mathbf{b}(\cdot)$ and $\sigma(\cdot)$
are constant only with respect to time, but not to space, allowing the
corresponding channel to be regarded as a time-invariant system. We will make
further assumptions on $\sigma(\cdot)$ as the deductions progress.
Nevertheless, for the purpose to calculate FAP density in specific scenarios
(as presented in Appendix D), we only need the stronger assumptions as stated
in Section II.

For each Itô diffusion $\mathbf{X}_{t}$, a corresponding operator called the
infinitesimal generator, or simply generator, can be associated. Let
$\mathcal{D}(A)$ denote the domain of the generator $A$ [30]. We define the
notation $\mathbb{E}^{\mathbf{x}}[\cdot]$ as taking expectation conditioned on
$\mathbf{X}_{0}=\mathbf{x}$, so that

$\mathbb{E}^{\mathbf{x}}[f(\mathbf{X}_{t})]:=\mathbb{E}[f(\mathbf{X}_{t})|\mathbf{X}_{0}=\mathbf{x}].$
(8)

The generator $A$ (which operates on function $f$) of a process
$\mathbf{X}_{t}$ can be defined as

$Af=A\\{f(\mathbf{x})\\}:=\lim\limits_{t\searrow
0}\dfrac{\mathbb{E}^{\mathbf{x}}[f(\mathbf{X}_{t})]-f(\mathbf{x})}{t}\text{\
for any\ }f\in\mathcal{D}(A).$ (9)

For time-homogeneous Itô processes $\mathbf{X}_{t}$, the time evolution is a
Markov process. By letting
$T_{t}f:=\mathbb{E}^{\mathbf{x}}[f(\mathbf{X}_{t})],$ we can obtain a
semigroup of operators $\mathcal{T}=\\{T_{t}\\}_{t\geq 0}$. Notice that the
term

$Af=\lim\limits_{t\searrow 0}\dfrac{T_{t}f-f}{t}$ (10)

can be interpreted as the rate of linear increment with respect to $t$ through
the semigroup evolution.

To explicitly calculate the generator $A$ based on (2), we can assume that $f$
is of class $C^{2}$, which means that it has both a continuous first
derivative and a continuous second derivative. Using Taylor expansion [25,
Formula (9.2.2)] and Itô’s formula[24, Chapter 4], we have:

$\displaystyle\begin{split}\ \mathrm{d}f\left(\mathbf{X}_{t}\right)=\ &(\nabla
f)^{\mathsf{\scriptstyle{T}}}(\mathrm{d}\mathbf{X}_{t})+\frac{1}{2}(\mathrm{d}\mathbf{X}_{t})^{\mathsf{\scriptstyle{T}}}\text{Hess}(f)(\mathrm{d}\mathbf{X}_{t})\\\
=\ &\left[(\nabla
f)^{\mathsf{\scriptstyle{T}}}(\mathbf{b}\left(\mathbf{X}_{t}\right))+\tr(\text{Hess}(f)\frac{\sigma\sigma^{\mathsf{\scriptstyle{T}}}\left(\mathbf{X}_{t}\right)}{2})\right]\mathrm{d}t\\\
&+(\nabla
f)^{\mathsf{\scriptstyle{T}}}\sigma\left(\mathbf{X}_{t}\right)\mathrm{d}\mathbf{B}_{t},\end{split}$
(11)

where $\tr(\cdot)$ denotes the trace operation, $\text{Hess}(f)$ denotes the
Hessian of $f$ and the matrix notation
$\sigma\sigma^{\mathsf{\scriptstyle{T}}}\left(\mathbf{X}_{t}\right)$ is short
for
$\sigma\left(\mathbf{X}_{t}\right)\left(\sigma\left(\mathbf{X}_{t}\right)\right)^{\mathsf{\scriptstyle{T}}}$.
Note that in (11), all the gradient and Hessian terms are evaluated at
$\bm{X}_{t}$. After taking expectation to both sides of (11), and then
plugging into (9), we use the fact that
$\mathbb{E}[\mathrm{d}\mathbf{B}_{t}]=0$ to obtain an explicit expression for
the generator $A$ as:

$A\\{f(\mathbf{x})\\}=\mathbf{b}(\mathbf{x})^{\mathsf{\scriptstyle{T}}}(\nabla
f(\mathbf{x}))+\dfrac{\sigma^{2}(\mathbf{x})}{2}\Delta f(\mathbf{x}),$ (12)

where $\mathbf{x}$ can be viewed as the starting point of a diffusion process.
From (12), we can see that $A$ is a second-order differential operator of
elliptic-type. Notice that in (12) we have assumed that the matrix
$\sigma(\cdot)$ is a scalar multiple of the identity matrix. To simplify the
notation, we also denoted this scalar function by $\sigma(\cdot)$ in (12).

We proceed to consider a boundary value problem (BVP) with an unknown function
$\phi\in C^{2}$, given by:

$\left\\{\begin{array}[]{ll}A\phi=0&\text{ in }\Omega\\\ \phi=g&\text{ on
}\partial\Omega\end{array}\right.,$ (13)

where $\Omega$ is the domain in which the MMs can diffuse, and
$\partial\Omega$ is the boundary of the domain. In (13), $A$ denotes a partial
differential operator and $g$ is the prescribed boundary data, which is also a
$C^{2}$ function. Note that for elliptic-type partial differential operators,
representation formulas for commonly encountered boundary conditions can be
directly found in PDE handbooks such as [29].

###  III-B A Novel Relation Between FAP Density and Elliptic Green’s Function

Recall the infinitesimal generator defined in (9). By rearranging the terms in
(9), we can obtain the following integral expression:

$\mathbb{E}^{\mathbf{x}}\left[f\left(\mathbf{X}_{t}\right)\right]=f(\mathbf{x})+\mathbb{E}^{\mathbf{x}}\left[\int_{0}^{t}A\\{f\left(\mathbf{X}_{s}\right)\\}\
\mathrm{d}s\right].$ (14)

It is important to note that the variable $t$ appearing in (14) is
deterministic, not a random variable.

In the field of stochastic analysis, Dynkin’s formula relates the expected
value of a function $f$ of a Markov process at a stopping time $\tau$ [31,
Definition 8.2.1] to the initial value $f(\mathbf{x})$ of the function and an
integral involving the generator $A$ of the process. Technically, we have to
make some assumptions on $f$. As we are considering a second-order
differential operator (12), it is common to require $f$ to be twice
continuously differentiable (i.e., $f\in C^{2}$). Letting $\tau$ be a stopping
time such that $\mathbb{E}^{\mathbf{x}}[\tau]<+\infty$, the Dynkin’s formula
can be expressed as:

$\mathbb{E}^{\mathbf{x}}\left[f\left(\mathbf{X}_{\tau}\right)\right]=f(\mathbf{x})+\mathbb{E}^{\mathbf{x}}\left[\int_{0}^{\tau}A\\{f\left(\mathbf{X}_{s}\right)\\}\
\mathrm{d}s\right].$ (15)

Notably, $\tau$ in (15) is now a random variable. Additional details regarding
this formula can be found in [25, eq. (9.3.10)].

Let $g$ denote a smooth function defined on the boundary $\partial\Omega$. We
can express the expected value of $g(\mathbf{X}_{\tau})$ as follows:

$\displaystyle\begin{split}\mathbb{E}^{\mathbf{x}}\big{[}g(\mathbf{X}_{\tau})\big{]}&=\mathbb{E}\Big{[}g(\mathbf{X}_{\tau})\Big{|}\mathbf{X}_{0}=\mathbf{x}\Big{]}=\int_{\partial\Omega}f_{\mathbf{Y}|\mathbf{X}}(\mathbf{y}|\mathbf{x})g(\mathbf{y})\
\mathrm{d}S_{\mathbf{y}},\end{split}$ (16)

where $\mathrm{d}S_{\mathbf{y}}$ is the magnitude of surface element (see[23,
Section 10.6]) at $\mathbf{y}$. Here, $\mathbf{x}\in\Omega$ represents the
starting point (corresponding to the emission point) of the diffusion, while
the hitting position $\mathbf{y}$ belongs to the boundary, i.e.,
$\mathbf{y}\in\partial\Omega$.

The conditional probability density function $f_{\mathbf{Y}|\mathbf{X}}$ in
(16) can be interpreted as the distribution of the position of
$\mathbf{X}_{t}$ at the stopping time $\tau$. In MC applications,
$\mathbf{X}_{t}$ is the trajectory of a MM, and the stopping event corresponds
to the capture by Rx. Thus, $f_{\mathbf{Y}|\mathbf{X}}$ is the desired FAP
density on the receiver boundary $\partial\Omega$. Our remaining task is to
express $f_{\mathbf{Y}|\mathbf{X}}$ using the Green’s function.

Suppose we have a solution $\phi(\mathbf{x})$ for the BVP in (13) with
prescribed boundary data $g$ such that $A\phi=0$ inside $\Omega$. By setting
$f(\mathbf{x})=\phi(\mathbf{x})$ in (15), we obtain
$\int_{0}^{\tau}A\\{\phi(\mathbf{X}_{s})\\}\ \mathrm{d}s=0.$ Here, $\tau$
represents the first hitting time. The physical interpretation of this
condition is that $\mathbf{X}_{s}$ is located inside $\Omega$ prior to the
hitting event, which leads to $A\\{\phi(\mathbf{X}_{s})\\}=0$ for $0<s<\tau$.
Furthermore, since $\phi(\mathbf{x})$ coincides with $g(\mathbf{x})$ on the
boundary, we can express the left-hand side of (15) as
$\mathbb{E}^{\mathbf{x}}[g(\mathbf{X}_{\tau})]$. Combining these two facts
yields:

$\mathbb{E}^{\mathbf{x}}[g(\mathbf{X}_{\tau})]=\phi(\mathbf{x})+0=\int_{\partial\Omega}f_{\mathbf{Y}|\mathbf{X}}(\mathbf{y}|\mathbf{x})g(\mathbf{y})\
\mathrm{d}S_{\mathbf{y}}.$ (17)

This holds true for any $\mathbf{x}\in\Omega$.

For a general elliptic-type BVP defined in a domain $\Omega$ with boundary
$\partial\Omega$, we have the following formula[29, Section 7.2.1]:

$\phi(\mathbf{x})=\int_{\Omega}\Xi(\mathbf{y})H_{1}(\mathbf{x},\mathbf{y})\
\mathrm{d}V_{\mathbf{y}}+\int_{\partial\Omega}g(\mathbf{y})H_{2}(\mathbf{x},\mathbf{y})\
\mathrm{d}S_{\mathbf{y}}.$ (18)

Here, $\Xi$ represents the source term, $\mathrm{d}V_{\mathbf{y}}$ is the
volume element[23, Section 10.7], $\mathrm{d}S_{\mathbf{y}}$ is the magnitude
of surface element[23, Section 10.6], and $H_{1}(\mathbf{x},\mathbf{y})$,
$H_{2}(\mathbf{x},\mathbf{y})$ depend on the type of boundary conditions under
consideration. For MC systems without molecule reproduction or annihilation
inside the channel, we can set $\Xi(\mathbf{y})=0$, and hence
$H_{1}(\mathbf{x},\mathbf{y})$ is irrelevant. Notice that BVP (13) has the
Dirichlet-type boundary condition [28, Chapter 13]. According to[29, Section
7.2.1], we have

$H_{2}(\mathbf{x},\mathbf{y})=-\dfrac{\partial
G(\mathbf{x},\mathbf{y})}{\partial\mathbf{n}_{\mathbf{y}}}$ (19)

in the solution formula (18), where $G$ represents the Green’s function of BVP
(13). The notation $\frac{\partial}{\partial\mathbf{n}_{\mathbf{y}}}$ means
the directional derivative along the outward normal at the point $\mathbf{y}$
on $\partial\Omega$. Letting $\Xi=0$ and plugging (19) into (18), we obtain a
simpler form:

$\phi(\mathbf{x})=\int_{\mathbf{y}\in\partial\Omega}\left|\dfrac{\partial
G(\mathbf{x},\mathbf{y})}{\partial\mathbf{n}_{\mathbf{y}}}\right|g(\mathbf{y})\
\mathrm{d}S_{\mathbf{y}}.$ (20)

Note that we have added an absolute value since the term
$\big{|}\frac{\partial
G(\mathbf{x},\mathbf{y})}{\partial\mathbf{n}_{\mathbf{y}}}\big{|}$ has the
physical meaning of conditional probability density, so it must be non-
negative. Finally, by comparing (17) and (20), we can establish a fundamental
relation between the FAP density and the elliptic-type Green’s function:

$f_{\mathbf{Y}|\mathbf{X}}(\mathbf{y}|\mathbf{x})=\left|\dfrac{\partial
G(\mathbf{x},\mathbf{y})}{\partial\mathbf{n}_{\mathbf{y}}}\right|\Bigg{|}_{\partial\Omega},$
(21)

where $\big{|}_{\partial\Omega}$ means $\mathbf{y}$ is evaluated on
$\partial\Omega$. Notice that (21) gives a concise relation between the FAP
density function and the Green’s function. How to explicitly calculate the FAP
density using (21) in specific scenarios is demonstrated in Appendix D.

##  IV Properties of VDFAP Distributions via Frequency Domain Analysis

To bound the VDFAP channel capacity under a second-moment constraint in Sec.
V, we conduct a frequency-domain analysis of the VDFAP distribution via its
_characteristic function_ (CF). From the CF we investigate the first and
second moments of VDFAP noise for analyzing the output second moment of VDFAP
channels. Also, we discover a property we call _weak stability_ exhibited by
VDFAP distributions. This property turns out to be essential for proving a
lower bound for the VDFAP channel capacity.

Based on probability theory [32, Section 3.3], a distribution can be
characterized by its CF, which is the Fourier transform of its PDF, i.e., a
frequency-domain representation of the distribution. Interestingly, via the
frequency-domain viewpoint, it is easier to compute the moments as well as
analyze the (weak) stability properties for VDFAP distributions.

For a random vector $\bm{N}$ following $\text{VDFAP}^{(d)}(u,\lambda)$, we
denote its CF as

$\displaystyle\Phi^{(d)}_{\bm{N}}(\bm{\omega};u,\lambda)=\mathbb{E}[\exp(i\bm{\omega}^{\mathsf{\scriptstyle{T}}}\bm{N})],\text{\
for\ }\bm{\omega}\in\mathbb{R}^{d},$ (22)

where $\mathbb{E}$ is the expectation operator and $i$ the imaginary unit. We
also write the CF simply as $\Phi^{(d)}_{\bm{N}}(\bm{\omega})$ when the
parameters are clear from the context. To our best knowledge, we have derived
a novel closed-form expression for the CF (22) of VDFAP when $d=1,2$, detailed
in Appendix A, as:

$\displaystyle\Phi^{(d)}_{\bm{N}}(\bm{\omega};u,\lambda)=\exp\left(-\lambda\left(\sqrt{\|\bm{\omega}\|^{2}+|u|^{2}}-|u|\right)\right).$
(23)

Note that letting $|u|\to 0$ in (23) reduces to the $d$-dimensional Cauchy CF
form:

$\displaystyle\lim_{|u|\to
0}\Phi^{(d)}_{\bm{N}}(\bm{\omega};u,\lambda)=\exp\left(-\lambda\|\bm{\omega}\|\right),$
(24)

which is discussed more in [33].

Via the closed-form CF expression (23) of a VDFAP distribution, we are able to
obtain, relatively easily, formulas for its moments in Subsec. IV-A, as well
as study its (weak) stability properties in Subsec. IV-B.

###  IV-A The Moments

In this subsection, we aim to find the first two moments of VDFAP
distributions. In particular, we derive the first moment $\mathbb{E}[\bm{N}]$
and the second moment222Whenever $\bm{N}$ has component random variables
$N_{1},\ldots,N_{d}$, we call
$\mathbb{E}\left[\|\bm{N}\|^{2}\right]=\sum_{i=1}^{d}\mathbb{E}\left[|N_{i}|^{2}\right]$
the “second moment” of $\bm{N}$. $\mathbb{E}\left[\|\bm{N}\|^{2}\right]$ for
their further use in analyzing the VDFAP channel capacity under a second-
moment constraint.

In principle, moments can be calculated by integrating the product of
monomials with the PDF. However, due to the complicated PDF (7), it is not
clear how to find closed-form expressions for these integrals. Due to its
similarity to the moment-generating function, the CF is useful for the
calculation of moments. In particular, from (22), the first moment of
$\bm{N}\sim\text{VDFAP}^{(d)}(u,\lambda)$ can be calculated as

$\displaystyle\mathbb{E}[\bm{N}]=-i\nabla\Phi^{(d)}_{\bm{N}}(\bm{\omega})\bigg{|}_{\bm{\omega}=\bm{0}}$
(25)

(where $\nabla$ is the gradient operator, and $\big{|}_{\bm{\omega}=\bm{0}}$
denotes evaluation at $\bm{\omega}=\bm{0}$), and its correlation matrix333We
are interested in the correlation matrix because its trace is exactly the
second moment:
$\Tr{\mathbb{E}\left[\bm{N}\bm{N}^{\mathsf{\scriptstyle{T}}}\right]}=\mathbb{E}\left[\Tr{\bm{N}\bm{N}^{\mathsf{\scriptstyle{T}}}}\right]=\mathbb{E}\left[\Tr{\bm{N}^{\mathsf{\scriptstyle{T}}}\bm{N}}\right]=\mathbb{E}\left[\Tr{\|\bm{N}\|^{2}}\right]=\mathbb{E}\left[\|\bm{N}\|^{2}\right]$,
where we applied the cyclic property of trace.
$\mathbb{E}\left[\bm{N}\bm{N}^{\mathsf{\scriptstyle{T}}}\right]$ as

$\displaystyle\mathbb{E}\left[\bm{N}\bm{N}^{\mathsf{\scriptstyle{T}}}\right]=-\text{Hess}\Big{(}\Phi^{(d)}_{\bm{N}}\Big{)}(\bm{\omega})\bigg{|}_{\bm{\omega}=\bm{0}}.$
(26)

To this end, we calculate the gradient of CF (23) as

$\displaystyle\nabla\Phi^{(d)}_{\bm{N}}(\bm{\omega})=\frac{-\lambda\Phi^{(d)}_{\bm{N}}(\bm{\omega})}{\sqrt{\|\bm{\omega}\|^{2}+|u|^{2}}}\bm{\omega}$
(27)

and the Hessian matrix of CF as

$\displaystyle\begin{split}\text{Hess}\Big{(}\Phi^{(d)}_{\bm{N}}\Big{)}(\bm{\omega})=\leavevmode\nobreak\
&\frac{-\lambda\Phi^{(d)}_{\bm{N}}(\bm{\omega})}{\sqrt{\|\bm{\omega}\|^{2}+|u|^{2}}}\mathbb{I}_{d}\\\
&+\frac{\lambda\Phi^{(d)}_{\bm{N}}(\bm{\omega})\left(1+\lambda\sqrt{\|\bm{\omega}\|^{2}+|u|^{2}}\right)}{(\|\bm{\omega}\|^{2}+|u|^{2})^{\frac{3}{2}}}\bm{\omega}\bm{\omega}^{\mathsf{\scriptstyle{T}}},\end{split}$
(28)

where $\mathbb{I}_{d}$ denotes the $d$-by-$d$ identity matrix. From (25) and
(27), we can deduce that

$\displaystyle\mathbb{E}[\bm{N}]=\bm{0}.$ (29)

Alternatively, due to the radial symmetry of $\bm{N}$, we can also deduce that
its first moment is zero. Also, from (26) and (28), the correlation matrix of
$\bm{N}$ is

$\displaystyle\mathbb{E}\left[\bm{N}\bm{N}^{\mathsf{\scriptstyle{T}}}\right]=\frac{\lambda}{|u|}\mathbb{I}_{d},$
(30)

which, upon taking the trace, yields

$\displaystyle\mathbb{E}[\|\bm{N}\|^{2}]=\frac{\lambda d}{|u|}.$ (31)

Another consequence of (30) is that components of
$\bm{N}\sim\text{VDFAP}^{(d)}(u,\lambda)$ are pairwise uncorrelated. However,
they are _not_ independent, since the PDF (7) cannot be decomposed as the
product of its marginals.

###  IV-B (Weak) Stability

In this subsection, we aim to prove a weak stability property for the VDFAP
distributions, which will be used later to establish a lower bound on the
VDFAP channel capacity.

Gaussian distributions, along with Cauchy distributions and Lévy
distributions, are known to be examples of _stable_ distributions [26].
Roughly speaking, a distribution is stable if any linear combination of two
independent random vectors from this distribution still follows the same
distribution. That is, it is akin to the _closure_ property of vector spaces.
By definition, a distribution is called stable if and only if it possesses a
specific form of CF [34, Eq. (2)], which is not satisfied by (23). Therefore,
VDFAP distributions are _not_ stable under this definition. However, despite
not being (strictly) stable, VDFAP distributions still exhibit certain
“weaker” stability properties.

One such _weak stability property_ that is instrumental in deriving a lower
bound (refer to the proof of Theorem 1) for the capacity of VDFAP channels is
encapsulated in the following lemma.

######  Lemma (Weak Stability Property for VDFAP).

If $\bm{N}_{1}\sim\text{VDFAP}^{(d)}(u,\lambda_{1})$ and
$\bm{N}_{2}\sim\text{VDFAP}^{(d)}(u,\lambda_{2})$ are two independent VDFAP
random vectors with the same normalized drift $u$, then their sum
$\bm{N}_{1}+\bm{N}_{2}\sim\text{VDFAP}^{(d)}(u,\lambda_{1}+\lambda_{2}).$

###### Proof.

To prove this property, note that $\bm{N}_{1}$ and $\bm{N}_{2}$ are
independent. Therefore, the CF of their sum is the product of their respective
CFs [35, Theorem 3.3.1], i.e.,
$\mathbb{E}[\exp(i\bm{\omega}^{\mathsf{\scriptstyle{T}}}(\bm{N}_{1}+\bm{N}_{2}))]=\Phi^{(d)}_{\bm{N}_{1}}(\bm{\omega};u,\lambda_{1})\Phi^{(d)}_{\bm{N}_{2}}(\bm{\omega};u,\lambda_{2}).$
Using the CF formula (23) for $\bm{N}_{1}$ and $\bm{N}_{2}$, we obtain

$\displaystyle\begin{split}&\mathbb{E}[\exp(i\bm{\omega}^{\mathsf{\scriptstyle{T}}}(\bm{N}_{1}+\bm{N}_{2}))]\\\
=&\exp\left(-\lambda_{1}\left(\sqrt{\|\bm{\omega}\|^{2}+|u|^{2}}-|u|\right)\right)\\\
&\cdot\exp\left(-\lambda_{2}\left(\sqrt{\|\bm{\omega}\|^{2}+|u|^{2}}-|u|\right)\right)\\\
=&\exp\left(-\left(\lambda_{1}+\lambda_{2}\right)\left(\sqrt{\|\bm{\omega}\|^{2}+|u|^{2}}-|u|\right)\right).\end{split}$
(32)

Since the CF uniquely characterizes the distribution [35, Theorem 3.1], we can
compare this result to the CF formula (23) to conclude that
$\bm{N}_{1}+\bm{N}_{2}\sim\text{VDFAP}^{(d)}(u,\lambda_{1}+\lambda_{2})$. ∎

##  V Bounds on the Capacity of VDFAP Channel

Borrowing the wisdom from parallel Gaussian channels with a common power
constraint444Although the noise components are assumed to be independent in
the setting of parallel Gaussian channels, this is not the case for VDFAP
channels when $d\geq 2$, as we argued at the end of Sec. IV-A. (see [36,
Section 10.4]), we investigate the VDFAP channel capacity under the second-
moment constraint as follows:

$\displaystyle\begin{split}C=\sup_{p(\bm{X}):\leavevmode\nobreak\
\mathbb{E}\left[\|\bm{X}\|^{2}\right]\leq P}I(\bm{X};\bm{Y}),\end{split}$ (33)

where the objective function $I(\bm{X};\bm{Y})$ is the _mutual information_
[36] between two $d$-dim random vectors $\bm{X}$ and $\bm{Y}$. The supremum is
taken over all input distributions $p(\bm{X})$ satisfying a constraint
$\mathbb{E}\left[\|\bm{X}\|^{2}\right]\leq P$ we call the _second-moment
constraint_ 555Mathematically, this is identical to the _common power
constraint_ for parallel Gaussian channels. But for FAP channels,
$\mathbb{E}\left[\|\bm{X}\|^{2}\right]$, which has the unit of squared length,
cannot be interpreted as power in the physical sense., where $P>0$ is an upper
bound on the second moment of the input random vector $\bm{X}$. Note that the
VDFAP channel output $\bm{Y}$ equals to $\bm{X}+\bm{N}$, where the noise
$\bm{N}$ follows $\text{VDFAP}^{(d)}(u,\lambda)$ distribution, and is
_independent_ of the input $\bm{X}$. The capacity appeared in (33) depends on
the dimension $d$ and the parameter triplet $(u,\lambda,P)$, so we sometimes
write it as $C^{(d)}(u,\lambda,P)$ to stress these dependencies.

Applying the additive channel structure $\bm{Y}=\bm{X}+\bm{N}$ and the
independence of $\bm{X}$ and $\bm{N}$, we have (assuming all differential
entropies $h(\cdot)$ involved are well-defined and finite)
$I(\bm{X};\bm{Y})=h(\bm{Y})-h(\bm{Y}|\bm{X})=h(\bm{Y})-h(\bm{N}|\bm{X})=h(\bm{Y})-h(\bm{N})$,
and thus (33) can be simplified as:

$\displaystyle\begin{split}C=\sup_{p(\bm{X}):\leavevmode\nobreak\
\mathbb{E}\left[\|\bm{X}\|^{2}\right]\leq
P}h(\bm{X}+\bm{N})-h(\bm{N}),\end{split}$ (34)

where $h(\bm{W})$ denotes the differential entropy of a $d$-dim random vector
$\bm{W}$ with support $S$ and PDF $f(\bm{w})$. Namely,
$h(\bm{W}):=-\int_{S}f(\bm{w})\log f(\bm{w})\leavevmode\nobreak\
\mathrm{d}\bm{w},$ where $\log(\cdot)$ denotes the _natural logarithm_
function throughout this paper. Two things to note here:

  * i)

In this Sec. V and Appendix B, the differential entropy $h(\cdot)$, the
capacity $C$ and its bounds are all measured in the unit of _nats_ for ease of
expression. In Sec. VI, when plotting bounds on $C$ we use the unit of _bits_
for ease of interpretation. $1\leavevmode\nobreak\
\text{nats}=\log_{2}e\leavevmode\nobreak\ \text{bits}\approx
1.44\leavevmode\nobreak\ \text{bits}$.

  * ii)

Abusing the notation, we allow the argument of $h(\cdot)$ to be a random
vector or the name of a distribution. For example, we write
$h(\text{VDFAP}^{(d)}(u,\lambda))$ to denote the differential entropy of the
$\text{VDFAP}^{(d)}(u,\lambda)$ distribution.

For $d=1$, to our best knowledge, no closed-form expression for the
differential entropy of $\text{VDFAP}^{(1)}(u,\lambda)$ is known. However, for
$d=2$, we have calculated a closed-form expression for differential entropy of
a random vector $\bm{N}$ following $\text{VDFAP}^{(2)}(u,\lambda)$. The
formula is

$\displaystyle\begin{split}&h(\text{VDFAP}^{(2)}(u,\lambda))\\\
=\leavevmode\nobreak\ &\log\left(2\pi
e^{3}\right)+2\log(\lambda)-\log(1+\lambda|u|)\\\
&-\lambda|u|e^{\lambda|u|}\big{(}e\cdot\text{Ei}(-1-\lambda|u|)-3\cdot\text{Ei}(-\lambda|u|)\big{)},\end{split}$
(35)

where the exponential integral function $\text{Ei}(\cdot)$ is defined as [27,
Section 8.21]:

$\displaystyle\text{Ei}(x)=-\int_{-x}^{\infty}\frac{e^{-t}}{t}\
\mathrm{d}t,\text{\ for\ }x<0.$ (36)

The detailed calculation of (35) is presented in Appendix B.

Next, we aim to establish two theorems characterizing the lower and upper
bounds for the VDFAP channel capacity $C$. By using the derived closed-form
expression (35), we further derive closed-form lower and upper bounds for $C$
in 3D space ($d=2$).

###  V-A Lower Bound on Capacity

Our strategy is to apply the weak stability property of VDFAP distributions
proved in Section IV-B to obtain a lower bound of the VDFAP channel capacity
$C=C^{(d)}(u,\lambda,P)$. The result is stated in the following Theorem.

######  Theorem 1 (Lower Bound on $C$).

The capacity $C^{(d)}(u,\lambda,P)$ is lower bounded as follows:

$\displaystyle\begin{split}C^{(d)}(u,\lambda,P)\geq\
&\sup_{0<\lambda^{\prime}\leq|u|P/d}h\left(\text{VDFAP}^{(d)}\left(u,\lambda+\lambda^{\prime}\right)\right)\\\
&-h\left(\text{VDFAP}^{(d)}\left(u,\lambda\right)\right).\end{split}$ (37)

###### Proof.

A general strategy to obtain lower bounds on the capacity (34) of additive
channels is to pick a family of parametrized input distributions
$p_{\theta}(\bm{X})$ satisfying the specified constraint(s), calculate the
corresponding $h(\bm{X}+\bm{N})-h(\bm{N})$ for each parameter $\theta$, and
then optimize the result via taking supremum over $\theta$.

Applying this general strategy to VDFAP channels under the second-moment
constraint while having the _weak stability_ Lemma (proved in Sec. Lemma) in
mind, we select the family of input distributions to be VDFAP with the _same_
$u$ as the noise but a generally different $\lambda^{\prime}$ such that the
second-moment constraint is satisfied. That is,
$\bm{X}\sim\text{VDFAP}^{(d)}(u,\lambda^{\prime})$ for some strictly positive
$\lambda^{\prime}\leq|u|P/d$. By (31), we have that
$\mathbb{E}\left[\|\bm{X}\|^{2}\right]=\frac{\lambda^{\prime}d}{|u|}\leq\frac{(|u|P/d)d}{|u|}=P,$
so $\bm{X}$ satisfies the second-moment constraint. Also, by _weak stability_
Lemma, $\bm{Y}=\bm{X}+\bm{N}$ follows
$\text{VDFAP}^{(d)}(u,\lambda+\lambda^{\prime})$. Then, a lower bound that
holds for any $0<\lambda^{\prime}\leq|u|P/d$ is

$\displaystyle\begin{split}&C^{(d)}(u,\lambda,\Sigma)\\\
\geq\leavevmode\nobreak\ &h(\bm{X}+\bm{N})-h(\bm{N})\\\ =\leavevmode\nobreak\
&h\left(\text{VDFAP}^{(d)}(u,\lambda+\lambda^{\prime})\right)-h\left(\text{VDFAP}^{(d)}(u,\lambda)\right),\end{split}$
(38)

where the inequality follows directly from (34). Taking the supremum over all
admissible $\lambda^{\prime}$ proves (37). ∎

######  Corollary 1 (Lower Bound on $C$ for $d=2$).

A closed-form expression of the lower bound given in Theorem 1 for $d=2$ is

$\displaystyle\begin{split}C^{(2)}(u,\lambda,P)&\geq
2\log\Big(1+\tfrac{P|u|}{2\lambda}\Big{missing})-\log\left(1+\tfrac{P|u|^{2}}{2(1+\lambda|u|)}\right)\\\
&\leavevmode\nobreak\ \leavevmode\nobreak\
+\lambda|u|e^{\lambda|u|}\big{(}e\cdot\text{Ei}(-1-\lambda|u|)-3\cdot\text{Ei}(-\lambda|u|)\big{)}\\\
&\leavevmode\nobreak\ \leavevmode\nobreak\
-\Big{(}\lambda+\tfrac{P|u|}{2}\Big{)}|u|\exp\left\\{\left(\lambda+\tfrac{P|u|}{2}\right)|u|\right\\}\\\
&\leavevmode\nobreak\ \leavevmode\nobreak\
\times\bigg{\\{}e\cdot\text{Ei}\Big{(}-1-\big{(}\lambda+\tfrac{P|u|}{2}\big{)}|u|\Big{)}\\\
&\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
-3\cdot\text{Ei}\Big{(}-\big{(}\lambda+\tfrac{P|u|}{2}\big{)}|u|\Big{)}\bigg{\\}}\\\
&>0.\end{split}$ (39)

###### Proof.

To obtain the closed-form expression for the lower bound in (37) for $d=2$, we
introduce an ancillary function $h_{0}:\mathbb{R}^{+}\to\mathbb{R}$ (where
$\mathbb{R}^{+}=\\{x\in\mathbb{R}:x>0\\}$) defined by

$\displaystyle\begin{split}h_{0}(s):=\leavevmode\nobreak\
&2\log(s)-\log(1+s)\\\
&-se^{s}\big{(}e\cdot\text{Ei}(-1-s)-3\cdot\text{Ei}(-s)\big{)},\end{split}$
(40)

so that the differential entropy (35) of a random vector
$\bm{Z}\sim\text{VDFAP}^{(2)}(u,\lambda^{\prime\prime})$ can be expressed as
$h(\mathbf{Z})=h_{0}(|u|\lambda^{\prime\prime})+\log(2\pi e^{3})-2\log(|u|).$
Therefore, the lower bound (37) can be expressed via the function
$h_{0}(\cdot)$ as

$\displaystyle
C^{(2)}(u,\lambda,P)\geq\sup_{0<\lambda^{\prime}\leq|u|P/2}h_{0}(|u|(\lambda+\lambda^{\prime}))-h_{0}(|u|\lambda).$
(41)

It is shown in Appendix C that $h_{0}(\cdot)$ is strictly monotonically
increasing. Therefore, (41) boils down to plugging the maximum
$\lambda^{\prime}$, i.e., $|u|P/2$, into (41):

$\displaystyle\begin{split}C^{(2)}(u,\lambda,P)&\geq
h_{0}\bigg{(}|u|\Big{(}\lambda+|u|P/2\Big{)}\bigg{)}-h_{0}(|u|\lambda)>0.\end{split}$
(42)

The strict positivity of the lower bound in (42) also follows from the
strictly monotonically increasing property of $h_{0}(\cdot)$ and the strict
positivity of $|u|$ and $P$. Also note that the input distribution achieving
this lower bound is exactly $\text{VDFAP}(u,|u|P/2)$, according to the proof
of Theorem 1.

Last but not least, the closed-form lower bound in (39) can be obtained by
directly plugging the expression (40) of $h_{0}$ into the lower bound in (42).
The strict positivity of the lower bound in (39) follows directly. ∎

###  V-B Upper Bound on Capacity

It is well known [37, Theorem 2.7, (2.20)] in information theory that under a
prescribed upper bound $Q>0$ on the second moment $\mathbb{E}[\|\bm{Y}\|^{2}]$
of a continuous random vector $\bm{Y}\in\mathbb{R}^{d}$, the differential
entropy $h(\bm{Y})$ is maximized by the multivariate Gaussian
$\mathcal{N}\left(\bm{0},\frac{Q}{d}\mathbb{I}_{d}\right)$. That is,

$\displaystyle\max_{q(\bm{Y}):\mathbb{E}[\|\bm{Y}\|^{2}]\leq
Q}h(\bm{Y})=h\left(\mathcal{N}\left(\bm{0},\tfrac{Q}{d}\mathbb{I}_{d}\right)\right)=\tfrac{d}{2}\log\left(2\pi
e\tfrac{Q}{d}\right),$ (43)

where $q(\bm{Y})$ denotes the distribution of $\bm{Y}$. Our strategy is to
apply this well-known fact to establish an upper bound for the VDFAP channel
capacity $C=C^{(d)}(u,\lambda,P)$. The result is stated in the following
theorem.

######  Theorem 2 (Upper Bound on $C$).

The capacity $C^{(d)}(u,\lambda,P)$ is upper bounded as follows:

$\displaystyle\begin{split}C^{(d)}(u,\lambda,P)\leq\leavevmode\nobreak\
&\tfrac{d}{2}\log\left(2\pi
e\left(\tfrac{P}{d}+\tfrac{\lambda}{|u|}\right)\right)\\\
&-h\left(\text{VDFAP}^{(d)}(u,\lambda)\right).\end{split}$ (44)

###### Proof.

Consider (34) as the maximization problem

$\displaystyle C=\sup_{p(\bm{X})\in\mathcal{F}_{P}}h(\bm{X}+\bm{N})-h(\bm{N})$
(45)

where the constraint set $\mathcal{F}_{P}$ is the family of input
distributions $\mathcal{F}_{P}=\left\\{p(\bm{X}):\leavevmode\nobreak\
\mathbb{E}\left[\|\bm{X}\|^{2}\right]\leq P\right\\}.$ For the VDFAP channel,
we can calculate the second moment of $\bm{Y}=\bm{X}+\bm{N}$ as

$\displaystyle\begin{split}\mathbb{E}[\|\bm{Y}\|^{2}]&=\mathbb{E}[\|\bm{X}+\bm{N}\|^{2}]\\\
&=\mathbb{E}[\|\bm{X}\|^{2}+2\bm{X}^{\mathsf{\scriptstyle{T}}}\bm{N}+\|\bm{N}\|^{2}]\\\
&=\mathbb{E}[\|\bm{X}\|^{2}]+2\mathbb{E}[\bm{X}^{\mathsf{\scriptstyle{T}}}\bm{N}]+\mathbb{E}[\|\bm{N}\|^{2}]\\\
&\overset{(a)}{=}\mathbb{E}[\|\bm{X}\|^{2}]+2\mathbb{E}[\bm{X}]^{\mathsf{\scriptstyle{T}}}\mathbb{E}[\bm{N}]+\mathbb{E}[\|\bm{N}\|^{2}]\\\
&\overset{(b)}{=}\mathbb{E}[\|\bm{X}\|^{2}]+\mathbb{E}[\|\bm{N}\|^{2}]\\\
&\overset{(c)}{=}\mathbb{E}[\|\bm{X}\|^{2}]+\tfrac{\lambda
d}{|u|},\end{split}$ (46)

where (a) is due to the independence of $\bm{X}$ and $\bm{N}$, (b) is due to
$\mathbb{E}[\bm{N}]=\bm{0}$ and (c) is due to (31). Defining $Q=P+\lambda
d/|u|$ and the corresponding family of second-moment constrained output
distributions as $\mathcal{G}_{Q}=\left\\{q(\bm{Y}):\leavevmode\nobreak\
\mathbb{E}\left[\|\bm{Y}\|^{2}\right]\leq Q\right\\},$ we have that (46)
implies: for any $P>0$,

$\displaystyle\begin{split}&p(\bm{X})\in\mathcal{F}_{P}\text{ and
}\bm{Y}=\bm{X}+\bm{N}\leavevmode\nobreak\ (\bm{X},\bm{N}\text{
independent})\\\ \Longrightarrow\leavevmode\nobreak\
&q(\bm{Y})\in\mathcal{G}_{Q}.\end{split}$ (47)

In other words, if there is an input distribution $p(\bm{X})$ satisfying the
second-moment constraint $P$, then the corresponding output
$\bm{Y}=\bm{X}+\bm{N}$ (with $\bm{X},\bm{N}$ independent) has a distribution
$q(\bm{Y})$ satisfying the second-moment constraint $Q$. However, the converse
statement does _not_ hold: even if a desired output distribution $q(\bm{Y})$
satisfies a second-moment constraint, there is _no_ guarantee that a
corresponding input distribution $p(\bm{X})$ should exist such that
$\bm{Y}=\bm{X}+\bm{N}$. Hence, the following maximization problem

$\displaystyle
C^{\prime}=\sup_{q(\bm{Y})\in\mathcal{G}_{Q}}h(\bm{Y})-h(\bm{N})$ (48)

is a _relaxation_ of the original problem (45), and because of this
relaxation, we have $C\leq C^{\prime}$. To conclude, from (43) and (48) we
have

$\displaystyle C\leq C^{\prime}=\tfrac{d}{2}\log\left(2\pi
e\tfrac{Q}{d}\right)-h\left(\text{VDFAP}^{(d)}(u,\lambda)\right),$ (49)

and (44) follows upon replacing $Q$ by its definition. ∎

######  Corollary 2 (Upper Bound on $C$ for $d=2$).

$\displaystyle\begin{split}C^{(2)}(u,\lambda,P)&\leq\log\left(\tfrac{P}{2\lambda^{2}}+\tfrac{1}{\lambda|u|}\right)+\log(1+\lambda|u|)-2\\\
&+\lambda|u|e^{\lambda|u|}\left(e\cdot\text{Ei}(-1-\lambda|u|)-3\cdot\text{Ei}(-\lambda|u|)\right).\end{split}$
(50)

###### Proof.

To get a closed-form upper bound when $d=2$, we plug the differential entropy
expression (35) into (44) and simplify to get (50). ∎

(a) Particle-based simulated FAP density on a planar receiver when drift is
zero. The vertical height represents the relative frequency, estimated via the
Monte Carlo method.

(b) Theoretical FAP density on a planar receiver when drift is zero.

Figure 3:  Comparison of particle-based simulated and theoretical FAP
densities on a planar receiver (i.e., $d=2$) for zero drift
$\mathbf{u}=[0,0,0]$.

(a) Particle-based simulated FAP density on a planar receiver when
$\mathbf{u}$ is $[2,-3,-1]$ $(\mu\text{m}^{-1})$. The vertical height
represents the relative frequency, estimated via the Monte Carlo method.

(b) Theoretical FAP density on a planar receiver when $\mathbf{u}$ is
$[2,-3,-1]$ $(\mu\text{m}^{-1})$.

Figure 4:  Comparison of particle-based simulated and theoretical FAP
densities on a planar receiver (i.e., $d=2$) for a case of non-zero normalized
drift $\mathbf{u}=[2,-3,-1]\leavevmode\nobreak\ (\mu\mathrm{m}^{-1})$.

##  VI Numerical Results

###  VI-A Particle-based Verification of the Derived FAP Density Formulas

In the 3D scenario ($d=2$) depicted in Fig. 2, where molecules are released
from a point source at $[0,0,1]\leavevmode\nobreak\ (\mu\mathrm{m})$ (i.e.,
the Tx-Rx distance $\lambda=1\leavevmode\nobreak\ (\mu\mathrm{m})$), we
carried out particle-based simulations (PBS) and numerical evaluations of the
theoretical FAP densities as per (6). We examined two different normalized
drift vectors, one being zero $\mathbf{u}=[0,0,0]$ and the other non-zero
$\mathbf{u}=[2,-3,-1]\leavevmode\nobreak\ (\mu\mathrm{m}^{-1})$. The resulting
PDFs of position deviation, referred to as FAP noise $[n_{1},n_{2}]$, are
displayed in Fig. 3 and Fig. 4, respectively.

For the PBS, we utilized the Monte Carlo method as detailed in [7, Eq.
(9)-(11)] to simulate $M=10^{6}$ independent trajectories of dissolved Ethanol
molecules diffusing in water with a diffusion coefficient $D\approx
840\leavevmode\nobreak\ (\mu\mathrm{m}^{2}/\mathrm{sec})$ and a time step of
$10^{-6}\leavevmode\nobreak\ (\mathrm{sec})$.666To avoid possible confusion
with the noise $N$, we use $M$ to denote the number of particles. The first
hitting position of each molecule was approximated by recording the first two
coordinates when the third coordinate first becomes negative, indicating
crossing the receiver plane at $x_{3}=0$. The derived FAP density estimates
are presented in Fig. 3(a) and Fig. 4(a), corresponding to each drift
scenario. The theoretical predictions were computed using the derived FAP
density formula (6). In producing these figures, we truncated values below
$10^{-3}$ to emphasize the approximate support of the distribution, although
theoretically, the FAP distributions could be supported on the whole
$\mathbb{R}^{2}$ plane.

The PBS results in the no-drift scenario (Fig. 3(a)) show a symmetrical
distribution of particles around the origin, aligning with our theoretical
predictions (Fig. 3(b)). In the scenario with drift, the particle distribution
skews in the direction of the parallel normalized drift
$\bm{u}_{\text{par}}=[2,-3]\leavevmode\nobreak\ (\mu\mathrm{m}^{-1})$, as seen
in Fig. 4(a), which is consistent with our theoretical expectations (Fig.
4(b)).

Figures 3 and 4 demonstrate a strong alignment between our theoretical
predictions and the PBS results, with a simulation error in density estimate
on the order of $1/\sqrt{M}=10^{-3}$. This agreement provides strong evidence
for the accuracy and reliability of our derived FAP density function.

###  VI-B Graphical Illustration of the Trends of Capacity Bounds

This subsection presents numerical evaluations illustrating the impact of key
parameters on the capacity bounds derived in (39) and (50) for a 3D scenario
($d=2$). In Fig. 5, we choose the operating point $[P=1\leavevmode\nobreak\
(\mu\mathrm{m}^{2}),\lambda=1\leavevmode\nobreak\
(\mu\mathrm{m}),u=-1\leavevmode\nobreak\ (\mu\mathrm{m}^{-1})]$ and vary one
parameter while keeping the other two fixed, to investigate the influence from
each key parameter alone while the others are held constant.

In Fig. 5(a), we demonstrate how capacity increases with the input second-
moment constraint $P$ , when $\lambda=1\leavevmode\nobreak\ (\mu\mathrm{m})$
and $u=-1\leavevmode\nobreak\ (\mu\mathrm{m}^{-1})$ are fixed. From the
figure, the trend is logarithmic (also can be seen from the upper bound
formula (50)), similar to the case for parallel Gaussian channels (see [36,
Section 10.4]). Also, it can be seen that the gap between upper and lower
bounds closes as $P$ grows, and the maximum gap is roughly $0.2$ (bits/channel
use) for the range $P\in(0,10]\leavevmode\nobreak\ (\mu\mathrm{m}^{2})$.

Fig. 5(b) illustrates the inverse relationship between capacity and the Tx-Rx
distance $\lambda$, when $P=1\leavevmode\nobreak\ (\mu\mathrm{m}^{2})$ and
$u=-1\leavevmode\nobreak\ (\mu\mathrm{m}^{-1})$ are fixed. This inverse
relationship is intuitive: the farther the Tx-Rx distance, the “worse” the FAP
channel is. It is also evident that the gap between upper and lower bounds
closes as $\lambda$ grows, and the largest gap is less than $0.5$
(bits/channel use) for the range $\lambda\in(0,10]\leavevmode\nobreak\
(\mu\mathrm{m})$.

Finally, Fig. 5(c) highlights the effect of normalized drift $u$ on capacity,
when $P=1\leavevmode\nobreak\ (\mu\mathrm{m}^{2})$ and
$\lambda=1\leavevmode\nobreak\ (\mu\mathrm{m})$ are fixed. It makes intuitive
sense that when $|u|$ is large, i.e., when the relative strength of drift over
diffusion is large, the VDFAP channel is more capable, as it takes a shorter
duration to travel from the Tx to the Rx, and the diffusion effect, which
introduces positional uncertainty during the same duration, is less
pronounced. From the figure, it is apparent that for $|u|$ larger than roughly
$0.5\leavevmode\nobreak\ (\mu\mathrm{m}^{-1})$, the gap between upper and
lower bounds closes as $|u|$ grows, and the maximum gap is roughly $0.3$
(bits/channel use) for this range of $|u|$. However, for $|u|$ smaller than
$0.5\leavevmode\nobreak\ (\mu\mathrm{m}^{-1})$, it can be observed that the
effectiveness of the upper bound deteriorates as $|u|\to 0$. In fact, from
(50) we see directly that the upper bound goes to $+\infty$ as $|u|\to 0$
while keeping $P$ and $\lambda$ fixed. This makes sense since the upper bound
is obtained by “matching” the output second moment to a multivariate Gaussian.
However, as $|u|\to 0$, the second moment of VDFAP noise $\bm{N}$ approaches
$+\infty$, and by (46) so does the second moment of the multivariate Gaussian.
This raises the interesting question of whether the second-moment constraint
is suitable when the VDFAP noise distribution approaches a multivariate Cauchy
distribution (known for having an infinite second moment) as $|u|\to 0$. This
topic is addressed in a separate work [33] by some of the authors.

These figures elucidate the comprehensiveness and tightness of the derived
bounds, underscoring how the upper and lower bounds diverge or converge under
various key parameters. This visual representation aids in understanding the
nuanced effects of these parameters on the derived capacity bounds for VDFAP
channels.

(a) Capacity bounds versus second-moment constraint $P$ with
$\lambda=1\leavevmode\nobreak\ (\mu\mathrm{m})$ and $u=-1\leavevmode\nobreak\
(\mu\mathrm{m}^{-1})$ fixed.

(b) Capacity bounds versus transmitter-receiver distance $\lambda$ with
$P=1\leavevmode\nobreak\ (\mu\mathrm{m}^{2})$ and $u=-1\leavevmode\nobreak\
(\mu\mathrm{m}^{-1})$ fixed.

(c) Capacity bounds versus absolute normalized drift $|u|$ with
$P=1\leavevmode\nobreak\ (\mu\mathrm{m}^{2})$ and
$\lambda=1\leavevmode\nobreak\ (\mu\mathrm{m})$ fixed.

Figure 5:  Numerical evaluations of derived lower bound (39) and upper bound
(50) on the capacity of VDFAP channels when $d=2$, illustrating the influence
of $P$, $\lambda$, and $u$ on the trends of capacity.

##  VII Conclusion and Future Work

###  VII-A Summary and Conclusion

In this paper, we have developed a mathematical rigorous model for MC systems
that harnesses FAP as a mode of information transmission. We tackled two
fundamental challenges: i) characterizing the FAP density and ii) establishing
capacity bounds for vertically-drifted FAP channels.

Our approach to characterizing FAP density in MC systems with a fully-
absorbing receiver is based on bridging macroscopic PDE models with
microscopic SDE models, resulting in a concise and versatile expression, as
presented in Eq. (21). This formula, which connects FAP density to elliptic-
type Green’s function, is notable for its broad applicability across any
spatial dimensions, any drift directions, and various receiver geometries.

The practicality and accuracy of our approach were demonstrated through
specific case studies involving various receiver geometries, including 2D
linear, 3D planar, and 3D spherical receivers (as presented in Appendix D).
These cases exemplify the versatility of our model. Further validation of our
model’s accuracy was achieved through particle-based simulations (in Section
VI-A), affirming the precision of our derived formulas.

On the other hand, armed with explicit expressions for the FAP density, we
have successfully established closed-form upper and lower bounds (Theorem 1
and Theorem 2) on the capacity of vertically-drifted FAP channels under a
second-moment constraint. This represents a significant advancement in our
understanding of the fundamental information transmission rate for FAP
channels in MC systems.

###  VII-B Potential Applications and Future Work

The development of FAP channels in molecular communication systems is not only
a theoretical advance but also holds significant practical implications. Our
research has particularly vital implications for Molecular MIMO [16] and
Molecular Index Modulation (Molecular-IM) systems.

In the field of Molecular-IM systems, as detailed in [17] and [18], there has
been a widespread adoption of “Rx box” models and Uniform Circular Array (UCA)
topologies, which are illustrated in [18, Fig. 1]. These models typically
utilize spherical-shaped Rx antennas. Nonetheless, there are considerable
practical difficulties in implementing perfectly spherical antennas,
particularly when integrating them into flat surfaces, owing to the limited
contact area they provide. This preference for spherical shapes stems from the
foundational research by Yilmaz et al. [15] in 2014, which introduced a
_point-to-sphere_ density formula. This formula has become a cornerstone in
most contemporary Molecular MIMO related research, influencing studies like
those conducted by Sabu et al. [38] and Kwak et al. [39]. These investigations
are instrumental in determining the proportion of MMs received by each antenna
in both Molecular MIMO and Molecular-IM systems. However, this approach seems
to be constrained to scenarios where the receivers or antennas have a
spherical configuration.

Our significant advancement in this area is the development of a _point-to-
plane_ position transition model, as depicted in Fig. 2. This model enables
the accommodation of various antenna geometries777 Recently, Brand et al. also
discussed the scenario of point Tx and cylindrical-shaped Rx in [40]. , such
as an ellipsoid-shaped antenna (which has not yet been explored in the
literature). By utilizing our point-to-plane formula and applying integration,
we can approximate a curved surface with numerous small planes. This approach
offers a promising solution for different antenna shapes. Additionally, our
model innovatively incorporates a non-zero background drift, a feature not
addressed in the point-to-sphere formula proposed by Yilmaz et al. [15]. This
inclusion marks a notable improvement over the existing methodologies. We
believe our new model represents a crucial step towards resolving the
complexities associated with more general antenna geometries in Molecular MIMO
and IM systems. Finally, note that a very recent work [41] proposed a
heuristic analytical model (utilizing negative particles) to describe the
impulse response of a Molecular SIMO (i.e., Single-Input Multiple-Output)
system with multiple Fully-Absorbing (FA) receivers. An interesting future
direction is to integrate our theoretical model with this heuristic model, and
compare their predictive performances with PBS.

##  Appendix A  Derivation of the Characteristic Function

We apply the Fourier transform (FT) pair relationship between the PDF and CF
of Student’s $t$-distribution to obtain an integral representation for (7).
Specifically, a random variable $N$ following the Student’s $t$-distribution
of $\nu$ degrees of freedom ($\nu>0$) has the PDF [42]:

$\displaystyle f^{(\nu)}_{N}(n)$
$\displaystyle=\frac{\Gamma\left(\frac{\nu+1}{2}\right)}{\sqrt{\nu\pi}\Gamma\left(\frac{\nu}{2}\right)}\left(1+\tfrac{n^{2}}{\nu}\right)^{-\tfrac{\nu+1}{2}},\text{\
for\ }n\in\mathbb{R},$ (51)

(where $\Gamma(\cdot)$ is the gamma function) and the CF [42]:

$\displaystyle\mathbb{E}[e^{i\omega
N}]=\frac{K_{\nu/2}(\sqrt{\nu}|\omega|)\cdot(\sqrt{\nu}|\omega|)^{\tfrac{\nu}{2}}}{\Gamma\left(\frac{\nu}{2}\right)2^{\frac{\nu-2}{2}}},\text{\
for\ }\omega\in\mathbb{R}.$ (52)

The FT pair relationship yields that

$\displaystyle K_{\nu/2}(\sqrt{\nu}|\omega|)$
$\displaystyle=\frac{2^{\frac{\nu-2}{2}}\Gamma\left(\frac{\nu+1}{2}\right)}{(\sqrt{\nu}|\omega|)^{\frac{\nu}{2}}\sqrt{\nu\pi}}\int_{\mathbb{R}}\left(1+\tfrac{n^{2}}{\nu}\right)^{-\tfrac{\nu+1}{2}}e^{i\omega
n}\leavevmode\nobreak\ \mathrm{d}n,$ (53)

for $\omega\in\mathbb{R}\setminus\\{0\\}$. Fix any $s>0$. Applying the
substitutions $\tilde{n}=(s/\sqrt{\nu})n$ and
$\tilde{\omega}=(\sqrt{\nu}/s)\omega$ to (53), we have

$\displaystyle\frac{K_{\nu/2}(s|\tilde{\omega}|)}{(s|\tilde{\omega}|)^{\frac{\nu}{2}}}$
$\displaystyle=\frac{2^{\frac{\nu-2}{2}}\Gamma\left(\frac{\nu+1}{2}\right)}{\sqrt{\pi}|\tilde{\omega}|^{\nu}}\int_{\mathbb{R}}\left(s^{2}+\tilde{n}^{2}\right)^{-\tfrac{\nu+1}{2}}e^{i\tilde{\omega}\tilde{n}}\leavevmode\nobreak\
\mathrm{d}\tilde{n},$ (54)

for $\tilde{\omega}\in\mathbb{R}\setminus\\{0\\}$. Setting $\nu=d+1$,
$s=(\norm{\bm{n}}^{2}+\lambda^{2})^{1/2}$ and $\tilde{\omega}=|u|$ in (54),
and then plugging into (7) yields

$\displaystyle\begin{split}f^{(d)}_{\bm{N}}(\bm{n};u,\lambda)&=\frac{\Gamma\left(\frac{d+2}{2}\right)}{\pi^{\frac{d+2}{2}}}\lambda
e^{\lambda|u|}\\\
&\quad\cdot\int_{\mathbb{R}}\left(\norm{\bm{n}}^{2}+\tilde{n}^{2}+\lambda^{2}\right)^{-\tfrac{d+2}{2}}e^{i\tilde{\omega}\tilde{n}}\leavevmode\nobreak\
\mathrm{d}\tilde{n}.\end{split}$ (55)

Thus for the $d$-dim VDFAP distribution, its CF (22):

$\displaystyle\Phi^{(d)}_{\bm{N}}(\bm{\omega};u,\lambda)$
$\displaystyle=\int_{\mathbb{R}^{d}}f^{(d)}_{\bm{N}}(\bm{n};u,\lambda)\leavevmode\nobreak\
e^{i\bm{\omega}^{\intercal}\bm{n}}\leavevmode\nobreak\ \mathrm{d}\bm{n}$ (56)

can be obtained by plugging (55) into (56), resulting in

$\displaystyle\begin{split}\Phi^{(d)}_{\bm{N}}(\bm{\omega};u,\lambda)=&\frac{\Gamma\left(\frac{d+2}{2}\right)}{\pi^{\frac{d+2}{2}}}\lambda
e^{\lambda|u|}\\\
&\cdot\int_{\mathbb{R}^{d+1}}\big{(}\norm{\mathbf{n}}^{2}+\lambda^{2}\big{)}^{-\frac{d+2}{2}}e^{i\mathbf{w}^{\intercal}\mathbf{n}}\leavevmode\nobreak\
\mathrm{d}\mathbf{n},\end{split}$ (57)

where
$\mathbf{n}=\left[\bm{n}^{\mathsf{\scriptstyle{T}}},\tilde{n}\right]^{\mathsf{\scriptstyle{T}}}$
and
$\mathbf{w}=\left[\bm{\omega}^{\mathsf{\scriptstyle{T}}},\tilde{\omega}\right]^{\mathsf{\scriptstyle{T}}}$
are vectors in $\mathbb{R}^{d+1}$.

We recognize that the integral in (57) is proportional to the CF of a
$(d+1)$-variate Cauchy distribution, with the last frequency variable
$\tilde{\omega}$ fixed at $|u|$. The PDF and CF of a $(d+1)$-variate Cauchy
vector $\mathbf{X}$ with location $\bm{\mu}=\bm{0}$ and scale
$\bm{\Sigma}=\lambda^{2}\mathbb{I}_{d+1}$ can be expressed as:

$\displaystyle
f_{\mathbf{X}}(\mathbf{x})=\frac{\Gamma\left(\frac{d+2}{2}\right)}{\pi^{\frac{d+2}{2}}}\frac{\lambda}{\left(\|\mathbf{x}\|^{2}+\lambda^{2}\right)^{\frac{d+2}{2}}},$
(58)

and
$\mathbb{E}\left[e^{i\mathbf{w}^{\intercal}\mathbf{X}}\right]=\exp(-\lambda\|\mathbf{w}\|),$
as shown in [43, 44]. Hence, using the FT pair relationship between the PDF
and CF of $\mathbf{X}$, we can evaluate (57) and get

$\displaystyle\begin{split}\Phi^{(d)}_{\bm{N}}(\bm{\omega})&=e^{\lambda|u|}\exp(-\lambda\|\mathbf{w}\|)\\\
&=\exp\bigg(-\lambda\Big{(}\sqrt{\norm{\bm{\omega}}^{2}+|u|^{2}}-|u|\Big{)}\bigg{missing}),\end{split}$
(59)

where we made use of the substitution $\tilde{\omega}=|u|$.

##  Appendix B  Closed-form Calculation of the Differential Entropy for $d=2$
VDFAP Distribution

We directly calculate the differential entropy as the integral in the
definition. After converting to polar coordinates (suggested by radial
symmetry) and making a substitution, the integral reduces to the form of the
following integral formula[27, Section 8.432]: for $a>0$,

$\displaystyle\begin{split}&\int_{a}^{\infty}\frac{K_{3/2}(\rho)}{\rho^{1/2}}\log\left(\frac{K_{3/2}(\rho)}{\rho^{3/2}}\right)\mathrm{d}\rho\\\
=\leavevmode\nobreak\
&\sqrt{\tfrac{\pi}{2}}\Big{(}e\cdot\text{Ei}(-1-a)-3\cdot\text{Ei}(-a)-e^{-a}\\\
&-(2a)^{-1}e^{-a}\left(6\log
a-2\log(1+a)+6+\log\tfrac{2}{\pi}\right)\Big{)}.\end{split}$ (60)

Using this formula, we are able to write the differential entropy in a
simplified closed-form as (35).

##  Appendix C  Strictly Monotonically Increasing Property of the Ancillary
Function $h_{0}(\cdot)$

Observe that $h_{0}(s)=2\log(s)-\log(1+s)-g(s)$, where we introduced a
function $g:\mathbb{R}^{+}\to\mathbb{R}$ defined by

$\displaystyle g(s)=se^{s+1}\text{Ei}(-(s+1))-3se^{s}\text{Ei}(-s),\text{ for
}s>0.$ (61)

Taking the derivative of (61) and using the definition of $\text{Ei}(\cdot)$
given in (36), we obtain the following formula by direct differentiation. That
is,

$\displaystyle
g^{\prime}(s)=\derivative{s}g(s)=\frac{s+1}{s}g(s)+\frac{s}{s+1}-3.$ (62)

As a consequence, the derivative of $h_{0}(\cdot)$ can be expressed as

$\displaystyle
h_{0}^{\prime}(s)=\derivative{s}h_{0}(s)=\frac{s+1}{s}(2-g(s)).$ (63)

To show that $h_{0}^{\prime}(s)>0$ for any $s>0$, and hence $h_{0}(\cdot)$ is
strictly monotonically increasing, it suffices to show that $g(s)<2$ for any
$s>0$. Using the inequalities [45, p.201]:

$\displaystyle\frac{s}{s+1}<-se^{s}\text{Ei}(-s)<\frac{s+1}{s+2}\quad(\forall
s>0),$ (64)

we can deduce that

$\displaystyle\begin{split}g(s)&<3\cdot\frac{s+1}{s+2}-\frac{s}{s+1}\cdot\frac{s+1}{(s+1)+1}=\frac{2s+3}{s+2}<2,\end{split}$
(65)

which holds for any $s>0$. Therefore, we have established the previously
mentioned sufficient condition for $h_{0}(\cdot)$ to be strictly monotonically
increasing.

##  Appendix D Calculating FAP Density in Specific Scenarios

This section demonstrates the application of formula (21) and the methodology
presented in Section III to compute the FAP density in some practical
scenarios that are relevant to MC.

###  D-A 2D Linear Receiver and 3D Planar Receiver

In this subsection, we unify the FAP density calculation for 2D linear
receivers and 3D planar receivers by employing the innovative approach
introduced in Section III.

To calculate the $d$-dimensional FAP density, we consider an Itô diffusion
$\mathbf{X}_{t}\in\mathbb{R}^{\text{D}}$ with a semigroup generator
$A:=\sum_{i=1}^{\text{D}}v_{i}\dfrac{\partial}{\partial
x_{i}}+\dfrac{\sigma^{2}}{2}\sum_{i=1}^{\text{D}}\dfrac{\partial^{2}}{\partial
x_{i}^{2}}$ (recall $d:=\text{D}-1$). We denote the Laplacian operator[29,
Section 8.5.1] in D dimensions by $\Delta_{\text{D}}$ and consider the
following drift-diffusion BVP in Cartesian coordinates:

$\left\\{\begin{array}[]{ll}A\phi=\sum_{i=1}^{\text{D}}v_{i}\frac{\partial\phi}{\partial
x_{i}}+\frac{\sigma^{2}}{2}\Delta_{\text{D}}\phi=0&\text{ in }\Omega\\\
\phi=g&\text{ on }\partial\Omega\end{array}\right.,$ (66)

where $\phi$ is a (dummy) unknown function,
$\Omega=\mathbb{R}^{\text{D}}\cap\left\\{x_{\text{D}}>0\right\\}$ is the
domain of (66), and
$\partial\Omega=\mathbb{R}^{\text{D}}\cap\left\\{x_{\text{D}}=0\right\\}$ is
the boundary of the domain. Here, we use the notation $x_{j}$ to indicate the
$j$-th component of the position vector
$\mathbf{x}=[x_{1},\ldots,x_{\text{D}}]^{\mathsf{\scriptstyle{T}}}$ in an
D-dimensional space. We temporarily set $\sigma^{2}=1$, and then treat the
general case later.

To simplify the problem further, we introduce a change of variable
$\psi(\mathbf{x})=\gamma(\mathbf{x})\phi(\mathbf{x}).$ The _drift factor_
$\gamma(\cdot)$ is defined as
$\gamma(\mathbf{x}):=\exp{\mathbf{v}^{\mathsf{\scriptstyle{T}}}\mathbf{x}},$
which depends only on the drift components. After changing the variables, we
obtain a new BVP with the same $\Omega$ as

$\left\\{\begin{array}[]{ll}\tilde{A}(\psi)=\Delta_{\text{D}}\psi-s^{2}\psi=0&\text{
in }\Omega\\\ \psi=\tilde{g}&\text{ on }\partial\Omega\end{array}\right.,$
(67)

where $s=\|\mathbf{v}\|$ denotes the magnitude of the drift vector. This new
operator $\tilde{A}$ is in the form of a Helmholtz operator (see[29, Section
7.3 and Section 8.3]), so we call this new BVP (67) the Helmholtz BVP.

According to [29], for both $\text{D}=2$ and $\text{D}=3$, the solution to the
Helmholtz BVP (67) can be compactly written as

$\psi(\mathbf{x})=\int_{\mathbf{y}\in\partial\Omega}\Big{|}\frac{\partial\tilde{G}(\mathbf{x},\mathbf{y})}{\partial\mathbf{n}_{\mathbf{y}}}\Big{|}\tilde{g}(\mathbf{y})\leavevmode\nobreak\
\mathrm{d}\mathbf{y},$ (68)

where the notation
$\big{|}\frac{\partial\tilde{G}(\mathbf{x},\mathbf{y})}{\partial\mathbf{n}_{\mathbf{y}}}\big{|}$
denotes the absolute value of the directional derivative of $\tilde{G}$ with
respect to the outward normal on the boundary $\partial\Omega$. The Green’s
function888We use $\tilde{G}$ to emphasize that this is the Green’s function
corresponding to the Helmholtz BVP, not to the original BVP. $\tilde{G}$
corresponding to Helmholtz BVP (67) is known and can be separately written
down for $\text{D}=2$ and $\text{D}=3$.

In order to obtain a representation formula for the solution $\phi$ of the
original BVP (66), we substitute
$\psi(\mathbf{x})=\gamma(\mathbf{x})\phi(\mathbf{x})$ back into (68), yielding

$\begin{aligned}
\phi(\mathbf{x})&=\frac{\psi(\mathbf{x})}{\gamma(\mathbf{x})}=\exp\\{-\mathbf{v}^{\mathsf{\scriptstyle{T}}}\mathbf{x}\\}\int_{\partial\Omega}\Big{|}\frac{\partial\tilde{G}(\mathbf{x},\mathbf{y})}{\partial\mathbf{n}_{\mathbf{y}}}\Big{|}\tilde{g}(\mathbf{y})\leavevmode\nobreak\
\mathrm{d}\mathbf{y}\\\
&\overset{(a)}{=}\int_{\partial\Omega}\exp\left\\{\mathbf{v}^{\mathsf{\scriptstyle{T}}}\mathbf{y}-\mathbf{v}^{\mathsf{\scriptstyle{T}}}\mathbf{x}\right\\}\Big{|}\frac{\partial\tilde{G}(\mathbf{x},\mathbf{y})}{\partial\mathbf{n}_{\mathbf{y}}}\Big{|}g(\mathbf{y})\leavevmode\nobreak\
\mathrm{d}\mathbf{y}\\\
&\overset{(b)}{=}\int_{\partial\Omega}\mathcal{K}_{\mathbf{v}}(\mathbf{x},\mathbf{y})g(\mathbf{y})\leavevmode\nobreak\
\mathrm{d}\mathbf{y}\end{aligned}.$ (69)

In $(a)$, we have used a relation between the boundary data before and after
the change of variable:
$\tilde{g}(\mathbf{y})=e^{\mathbf{v}^{\mathsf{T}}\mathbf{y}}g(\mathbf{y}).$ In
$(b)$, we use the notation $\mathcal{K}_{\mathbf{v}}(\mathbf{x},\mathbf{y})$
to represent the overall “integral kernel” (see [23, Section 6.1] for this
term). Note that this integral kernel can be also expressed by
$\big{|}\frac{\partial
G(\mathbf{x},\mathbf{y})}{\partial\mathbf{n}_{\mathbf{y}}}\big{|}$, where $G$
corresponds to the Green’s function of original BVP (66). From (69), we see
that

$\displaystyle\Big{|}\frac{\partial
G(\mathbf{x},\mathbf{y})}{\partial\mathbf{n}_{\mathbf{y}}}\Big{|}=\mathcal{K}_{\mathbf{v}}(\mathbf{x},\mathbf{y})=\exp\left\\{\mathbf{v}^{\mathsf{\scriptstyle{T}}}\left(\mathbf{y}-\mathbf{x}\right)\right\\}\Big{|}\frac{\partial\tilde{G}(\mathbf{x},\mathbf{y})}{\partial\mathbf{n}_{\mathbf{y}}}\Big{|}.$
(70)

In the following, we separately write down, for $\text{D}=2$ and $\text{D}=3$,
the Green’s function $\tilde{G}$ corresponding to Helmholtz BVP (67) from [29]
and explicitly carry out the calculation in (70).

####  D-A1 2D Linear Receiver

According to[29, Section 7.3], when $\text{D}=2$, the Green’s function
corresponding to Helmholtz BVP (67) has the form

$\tilde{G}(\mathbf{x},\mathbf{y})=\tilde{G}(x_{1},x_{2},\xi,\eta)=\frac{1}{2\pi}[K_{0}(s\rho_{1})-K_{0}(s\rho_{2})],$
(71)

where $\mathbf{x}=[x_{1},x_{2}]^{\mathsf{\scriptstyle{T}}}$ and
$\mathbf{y}=[\xi,\eta]^{\mathsf{\scriptstyle{T}}}$ are general positions in
2D, and we have defined $\rho_{1}:=\sqrt{(x_{1}-\xi)^{2}+(x_{2}-\eta)^{2}}$
and $\rho_{2}:=\sqrt{(x_{2}-\xi)^{2}+(x_{2}+\eta)^{2}}.$ From (71), we can
directly calculate

$\displaystyle\begin{split}\Big{|}\frac{\partial\tilde{G}(\mathbf{x},\mathbf{y})}{\partial\mathbf{n}_{\mathbf{y}}}\Big{|}\Bigg{|}_{\partial\Omega}&=\left[\dfrac{\partial\tilde{G}}{\partial\rho_{1}}\dfrac{\partial\rho_{1}}{\partial\eta}+\dfrac{\partial\tilde{G}}{\partial\rho_{2}}\dfrac{\partial\rho_{2}}{\partial\eta}\right]_{\eta=0}\\\
&=\dfrac{s\lambda}{\pi}\dfrac{K_{1}\left(s\sqrt{(x_{1}-\xi)^{2}+\lambda^{2}}\right)}{\sqrt{(x_{1}-\xi)^{2}+\lambda^{2}}},\end{split}$
(72)

where we have specified the emission point
$\mathbf{x}=[x_{1},\lambda]^{\mathsf{\scriptstyle{T}}}$ as well as the arrival
point $\mathbf{y}=[\xi,0]^{\mathsf{\scriptstyle{T}}}$ hereafter for the 2D
linear receiver. Notice that $\mathbf{y}$ now lies on the boundary
$\partial\Omega$. The parameter $\lambda>0$ is the distance from the emission
point to the linear receiver, as illustrated in Fig. 1. Then according to (70)
and $s=\|\mathbf{v}\|$, we obtain

$\displaystyle\begin{split}\Big{|}\frac{\partial
G(\mathbf{x},\mathbf{y})}{\partial\mathbf{n}_{\mathbf{y}}}\Big{|}\Bigg{|}_{\partial\Omega}&=\mathcal{K}_{\mathbf{v}}(\mathbf{x},\mathbf{y})\Big{|}_{\partial\Omega}\\\
&=\exp\left\\{\mathbf{v}^{\mathsf{\scriptstyle{T}}}\left(\mathbf{y}-\mathbf{x}\right)\right\\}\dfrac{\|\mathbf{v}\|\lambda}{\pi}\dfrac{K_{1}\left(\|\mathbf{v}\|\|\mathbf{y}-\mathbf{x}\|\right)}{\|\mathbf{y}-\mathbf{x}\|},\end{split}$
(73)

where we have substituted
$\|\mathbf{y}-\mathbf{x}\|=\sqrt{(x_{1}-\xi)^{2}+\lambda^{2}}$ into (72). This
$\mathcal{K}_{\mathbf{v}}(\mathbf{x},\mathbf{y})\Big{|}_{\partial\Omega}$ is
just the desired FAP density according to relation (21) presented in Section
III. We will denote this integral kernel by
$f_{\mathbf{Y}\mid\mathbf{X}}(\mathbf{y}|\mathbf{x})$ hereafter.

For the general case $\sigma^{2}\neq 1$, we only need to replace $v_{i}$ with
$\frac{v_{i}}{\sigma^{2}}=u_{i}$ (recall
$\mathbf{u}:=\frac{\mathbf{v}}{\sigma^{2}}$), yielding the desired FAP density
that allows arbitrary drift directions:

$\displaystyle\begin{split}f_{\mathbf{Y}\mid\mathbf{X}}(\mathbf{y}|\mathbf{x})=\dfrac{\|\mathbf{u}\|\lambda}{\pi}\exp\left\\{\mathbf{u}^{\mathsf{\scriptstyle{T}}}\left(\mathbf{y}-\mathbf{x}\right)\right\\}\dfrac{K_{1}\left(\|\mathbf{u}\|\|\mathbf{y}-\mathbf{x}\|\right)}{\|\mathbf{y}-\mathbf{x}\|}.\end{split}$
(74)

Equation (74) is a new formula in the MC field. Setting $v_{1}=0$ in (74)
recovers the restrictive version in [21]:

$\displaystyle\begin{split}&f_{\mathbf{Y}\mid\mathbf{X}}(\mathbf{y}|\mathbf{x})\\\
=\leavevmode\nobreak\ &\dfrac{|v_{2}|\lambda}{2\pi
D}\exp\left\\{\dfrac{-v_{2}\lambda}{2D}\right\\}\dfrac{K_{1}\left(\dfrac{|v_{2}|}{2D}\sqrt{(x_{1}-\xi)^{2}+\lambda^{2}}\right)}{\sqrt{(x_{1}-\xi)^{2}+\lambda^{2}}},\end{split}$
(75)

where we have used the relation $\sigma^{2}=2D$.999 In the literature, the
diffusion coefficient $D$ may have different definitions, which can vary by a
constant factor. However, the unit of $D$ consistently corresponds to
$\sigma^{2}$.  To conclude, our formula (74) improves the formula [21, eq.
(19)] by allowing arbitrary drift directions.

####  D-A2 3D Planar Receiver

According to[29, Section 8.3], when $\text{D}=3$, the Green’s function
corresponding to Helmholtz BVP (67) has the form:

$\displaystyle\begin{split}\tilde{G}(\mathbf{x},\mathbf{y})&=\tilde{G}(x_{1},x_{2},x_{3},\xi,\eta,\zeta)\\\
&=\dfrac{\exp\left(-s\mathcal{R}_{1}\right)}{4\pi\mathcal{R}_{1}}-\dfrac{\exp\left(-s\mathcal{R}_{2}\right)}{4\pi\mathcal{R}_{2}},\end{split}$
(76)

where $\mathbf{x}=[x_{1},x_{2},x_{3}]^{\mathsf{\scriptstyle{T}}}$ and
$\mathbf{y}=[\xi,\eta,\zeta]^{\mathsf{\scriptstyle{T}}}$ are general positions
in 3D, and we have defined
$\mathcal{R}_{1}:=\sqrt{(x_{1}-\xi)^{2}+(x_{2}-\eta)^{2}+(x_{3}-\zeta)^{2}}$
and
$\mathcal{R}_{2}:=\sqrt{(x_{1}-\xi)^{2}+(x_{2}-\eta)^{2}+(x_{3}+\zeta)^{2}}$.
From (76), we can directly differentiate to get

$\displaystyle\begin{split}\Big{|}\frac{\partial\tilde{G}(\mathbf{x},\mathbf{y})}{\partial\mathbf{n}_{\mathbf{y}}}\Big{|}\Bigg{|}_{\partial\Omega}&=\left[\frac{\partial}{\partial\zeta}\tilde{G}(x_{1},x_{2},x_{2},\xi,\eta,\zeta)\right]_{\zeta=0}\\\
&=\frac{\lambda}{2\pi}\exp\\{-s\|\mathbf{y}-\mathbf{x}\|\\}\left[\frac{1+s\|\mathbf{y}-\mathbf{x}\|}{\|\mathbf{y}-\mathbf{x}\|^{3}}\right]\end{split}$
(77)

where we have specified the emission point
$\mathbf{x}=[x_{1},x_{2},\lambda]^{\mathsf{\scriptstyle{T}}}$ and the arrival
point $\mathbf{y}=[\xi,\eta,0]^{\mathsf{\scriptstyle{T}}}$ hereafter for the
3D planar receiver. Notice that $\mathbf{y}$ now lies on the boundary
$\partial\Omega$. The parameter $\lambda>0$ is the distance from the emission
point to the planar receiver, as illustrated in Fig. 2. Then according to (70)
and $s=\|\mathbf{v}\|$, we obtain

$\displaystyle\begin{split}&\Big{|}\frac{\partial
G(\mathbf{x},\mathbf{y})}{\partial\mathbf{n}_{\mathbf{y}}}\Big{|}\Bigg{|}_{\partial\Omega}=\mathcal{K}_{\mathbf{v}}(\mathbf{x},\mathbf{y})\Big{|}_{\partial\Omega}\\\
=\leavevmode\nobreak\
&\exp\left\\{\mathbf{v}^{\mathsf{\scriptstyle{T}}}\left(\mathbf{y}-\mathbf{x}\right)\right\\}\\\
&\frac{\lambda}{2\pi}\exp\\{-\|\mathbf{v}\|\|\mathbf{y}-\mathbf{x}\|\\}\left[\frac{1+\|\mathbf{v}\|\|\mathbf{y}-\mathbf{x}\|}{\|\mathbf{y}-\mathbf{x}\|^{3}}\right].\end{split}$
(78)

This integral kernel
$\mathcal{K}_{\mathbf{v}}(\mathbf{x},\mathbf{y})\Big{|}_{\partial\Omega}$ can
be interpreted as the probability density function of the FAP (see (21),
Section III), so we denote it by
$f_{\mathbf{Y}\mid\mathbf{X}}(\mathbf{y}|\mathbf{x})$ hereafter.

Lastly, to treat the general case $\sigma^{2}\neq 1$, we can simply replace
$v_{i}$ with $\frac{v_{i}}{\sigma^{2}}=u_{i}$ (recall
$\mathbf{u}:=\frac{\mathbf{v}}{\sigma^{2}}$). This yields:

$\displaystyle\begin{split}f_{\mathbf{Y}\mid\mathbf{X}}(\mathbf{y}|\mathbf{x})=\leavevmode\nobreak\
&\frac{\lambda}{2\pi}\exp\left\\{\mathbf{u}^{\mathsf{\scriptstyle{T}}}\left(\mathbf{y}-\mathbf{x}\right)\right\\}\\\
&\exp\left\\{-\|\mathbf{u}\|\norm{\mathbf{y}-\mathbf{x}}\right\\}\left[\frac{1+\|\mathbf{u}\|\norm{\mathbf{y}-\mathbf{x}}}{\norm{\mathbf{y}-\mathbf{x}}^{3}}\right].\end{split}$
(79)

Equation (79) is the desired 3D FAP density function, allowing arbitrary drift
directions.

###  D-B 3D Spherical Receiver

In this subsection, we tackle the 3D spherical receiver case with the
assumption of zero drift velocity, the same setting as in [22]. Note that the
result we are going to show is already known in the literature, so we present
a very brief version of how to obtain the angular density using our proposed
methodology. The main goal of this subsection is to demonstrate the capability
of our method to handle receivers with various shapes.

We consider an Itô diffusion $\mathbf{X}_{t}\in\mathbb{R}^{3}$ with zero
drift, i.e., all $v_{i}=0$. The corresponding generator $A$ is simply the 3D
Laplacian operator. It is well known that the 3D Laplace equation in the
spherical coordinate [23, Section 12.11] can be expressed as

$\displaystyle\begin{split}\frac{1}{r^{2}}\frac{\partial}{\partial
r}\left(r^{2}\frac{\partial\phi}{\partial
r}\right)&+\frac{1}{r^{2}\sin\theta}\frac{\partial}{\partial\theta}\left(\sin\theta\frac{\partial\phi}{\partial\theta}\right)\\\
&+\frac{1}{r^{2}\sin^{2}\theta}\frac{\partial^{2}\phi}{\partial\varphi^{2}}=0\end{split}$
(80)

where $r=\sqrt{x^{2}+y^{2}+z^{2}}$. Assuming a prescribed boundary data
$g(\theta,\varphi)$ at the surface ($r=R$) of the spherical receiver with
radius $R>0$, i.e., $\phi=g(\theta,\varphi)\text{\ at\ }r=R,$ the solution of
the outer problem for $r\geq R$ can be obtained as [29, Section 8.1.3]:

$\displaystyle\begin{split}&\phi(r,\theta,\varphi)\\\ =&\
\frac{R}{4\pi}\int_{0}^{2\pi}\int_{0}^{\pi}g\left(\theta_{0},\varphi_{0}\right)\cdot\\\
&\frac{r^{2}-R^{2}}{\left(r^{2}-2rR\cdot\kappa(\theta,\varphi;\theta_{0},\varphi_{0})+R^{2}\right)^{3/2}}\sin\theta_{0}\
\mathrm{d}\theta_{0}\ \mathrm{d}\varphi_{0},\end{split}$ (81)

where

$\displaystyle\kappa(\theta,\varphi;\theta_{0},\varphi_{0})=\cos\theta\cos\theta_{0}+\sin\theta\sin\theta_{0}\cos(\varphi-\varphi_{0}).$
(82)

A comparison between formula (81) and equation (20) allows us to determine
$\big{|}\frac{\partial
G}{\partial\mathbf{n}_{\mathbf{y}}}(\mathbf{x},\mathbf{y})\big{|}$, which is
exactly the FAP density $f_{\mathbf{Y}|\mathbf{X}}(\mathbf{y}|\mathbf{x})$ due
to our proposed formula (21). That is, we have already obtained the (marginal)
_angular density_ of the hitting position of molecules, which was expressed in
[22, eq. (3)-(4)] and [46, eq. (6.3.3a)]. This demonstrates the versatility of
our FAP calculating methodology in handling various receiver shapes, including
spherical receivers.

## References

  * [1] Y.-C. Lee, J.-M. Wu, and M.-H. Hsieh, “A unified framework for calculating first arrival position density in molecular communication,” _arXiv preprint arXiv:2201.04476v3_ , 2022.
  * [2] Y.-F. Lo, Y.-C. Lee, and M.-H. Hsieh, “Capacity bounds for vertically-drifted first arrival position channels under a covariance constraint,” _arXiv preprint arXiv:2305.02706v4_ , 2023.
  * [3] T. Nakano, A. W. Eckford, and T. Haraguchi, _Molecular Communication_. Cambridge University Press, 2013.
  * [4] W. Guo, C. Mias, N. Farsad, and J.-L. Wu, “Molecular versus electromagnetic wave propagation loss in macro-scale environments,” _IEEE Transactions on Molecular, Biological and Multi-Scale Communications_ , vol. 1, no. 1, pp. 18–25, 2015.
  * [5] I. F. Akyildiz, F. Brunetti, and C. Blázquez, “Nanonetworks: A new communication paradigm,” _Computer Networks_ , vol. 52, no. 12, pp. 2260–2279, 2008.
  * [6] M. Pierobon and I. F. Akyildiz, “A physical end-to-end model for molecular communication in nanonetworks,” _IEEE Journal on Selected Areas in Communications_ , vol. 28, no. 4, pp. 602–611, 2010.
  * [7] N. Farsad, H. B. Yilmaz, A. Eckford, C.-B. Chae, and W. Guo, “A comprehensive survey of recent advancements in molecular communication,” _IEEE Communications Surveys & Tutorials_, vol. 18, no. 3, pp. 1887–1919, 2016.
  * [8] M. Pierobon and I. F. Akyildiz, “Capacity of a diffusion-based molecular communication system with channel memory and molecular noise,” _IEEE Transactions on Information Theory_ , vol. 59, no. 2, pp. 942–954, 2012.
  * [9] S. Kadloor, R. S. Adve, and A. W. Eckford, “Molecular communication using Brownian motion with drift,” _IEEE Transactions on NanoBioscience_ , vol. 11, no. 2, pp. 89–99, 2012.
  * [10] M. Moore, A. Enomoto, T. Nakano, R. Egashira, T. Suda, A. Kayasuga, H. Kojima, H. Sakakibara, and K. Oiwa, “A design of a molecular communication system for nanomachines using molecular motors,” in _Proc. 4th Annu. IEEE Int. Conf. Pervasive Comput. Commun. Workshops (PERCOMW)_ , Pisa, Italy, Mar. 2006, pp. 554–559.
  * [11] D. Bi, Y. Deng, M. Pierobon, and A. Nallanathan, “Chemical reactions-based microfluidic transmitter and receiver design for molecular communication,” _IEEE Transactions on Communications_ , vol. 68, no. 9, pp. 5590–5605, 2020.
  * [12] M. Ş. Kuran, H. B. Yilmaz, I. Demirkol, N. Farsad, and A. Goldsmith, “A survey on modulation techniques in molecular communication via diffusion,” _IEEE Communications Surveys & Tutorials_, vol. 23, no. 1, pp. 7–28, 2020.
  * [13] V. Jamali, A. Ahmadzadeh, W. Wicke, A. Noel, and R. Schober, “Channel modeling for diffusive molecular communication—a tutorial review,” _Proceedings of the IEEE_ , vol. 107, no. 7, pp. 1256–1301, 2019.
  * [14] A. Noel, Y. Deng, D. Makrakis, and A. Hafid, “Active versus passive: receiver model transforms for diffusive molecular communication,” in _2016 IEEE Global Communications Conference (GLOBECOM)_ , Washington, D.C., USA, Dec. 2016, pp. 1–6.
  * [15] H. B. Yilmaz, A. C. Heren, T. Tuğcu, and C.-B. Chae, “Three-dimensional channel characteristics for molecular communications with an absorbing receiver,” _IEEE Communications Letters_ , vol. 18, no. 6, pp. 929–932, 2014.
  * [16] B.-H. Koo, C. Lee, H. B. Yilmaz, N. Farsad, A. Eckford, and C.-B. Chae, “Molecular MIMO: From theory to prototype,” _IEEE Journal on Selected Areas in Communications_ , vol. 34, no. 3, pp. 600–614, 2016.
  * [17] M. C. Gursoy, E. Basar, A. E. Pusane, and T. Tuğcu, “Pulse position-based spatial modulation for molecular communications,” _IEEE Communications Letters_ , vol. 23, no. 4, pp. 596–599, 2019.
  * [18] M. C. Gursoy, A. Celik, E. Basar, A. E. Pusane, and T. Tuğcu, “Molecular index modulation with space-time equalization,” _IEEE Wireless Communications Letters_ , vol. 9, no. 5, pp. 702–705, 2020.
  * [19] K. V. Srinivas, A. W. Eckford, and R. S. Adve, “Molecular communication in fluid media: The additive inverse Gaussian noise channel,” _IEEE Transactions on Information Theory_ , vol. 58, no. 7, pp. 4678–4692, 2012.
  * [20] Y.-C. Lee, C.-C. Chen, P.-C. Yeh, and C.-H. Lee, “Distribution of first arrival position in molecular communication,” in _Proc. IEEE Int. Symp. Information Theory (ISIT)_ , Barcelona, Spain, July 2016, pp. 1033–1037.
  * [21] N. Pandey, R. K. Mallik, and B. Lall, “Molecular communication: The first arrival position channel,” _IEEE Wireless Communications Letters_ , vol. 8, no. 2, pp. 508–511, 2018.
  * [22] B. C. Akdeniz, N. A. Turgut, H. B. Yilmaz, C.-B. Chae, T. Tuğcu, and A. E. Pusane, “Molecular signal modeling of a partially counting absorbing spherical receiver,” _IEEE Transactions on Communications_ , vol. 66, no. 12, pp. 6237–6246, 2018.
  * [23] E. Kreyszig, _Advanced Engineering Mathematics_ , 10th ed. Hoboken, NJ, USA: Wiley, 2011.
  * [24] B. Øksendal, _Stochastic Differential Equations: An Introduction with Applications_ , 6th ed. Springer Science & Business Media, 2002.
  * [25] O. Calin, _An Informal Introduction to Stochastic Calculus with Applications_. World Scientific, 2015.
  * [26] N. Farsad, W. Guo, C.-B. Chae, and A. Eckford, “Stable distributions as noise models for molecular communication,” in _Proc. IEEE Global Communications Conference (GLOBECOM)_ , San Diego, CA, USA, Dec. 2015, pp. 1–6.
  * [27] I. S. Gradshteyn and I. M. Ryzhik, _Table of Integrals, Series, and Products_ , 8th ed. Academic Press, 2014.
  * [28] D. G. Zill, _Advanced Engineering Mathematics_. Jones & Bartlett Learning, 2020.
  * [29] A. D. Polyanin, _Handbook of Linear Partial Differential Equations for Engineers and Scientists_ , 1st ed. Chapman and hall/crc, 2001.
  * [30] M. H. Stone, “On one-parameter unitary groups in Hilbert space,” _Annals of Mathematics_ , pp. 643–648, 1932.
  * [31] S. E. Shreve, _Stochastic Calculus for Finance II: Continuous-Time Models_. Springer, 2004.
  * [32] R. Durrett, _Probability: Theory and Examples_. Cambridge University Press, 2019.
  * [33] Y.-C. Lee and M.-H. Hsieh, “On the capacity of zero-drift first arrival position channels in diffusive molecular communication,” to appear in _IEEE International Conference on Communications (ICC)_ , 2024.
  * [34] H. Fallahgoul, S. Hashemiparast, F. J. Fabozzi, and Y. S. Kim, “Multivariate stable distributions and generating densities,” _Applied Mathematics Letters_ , vol. 26, no. 3, pp. 324–329, 2013.
  * [35] E. Lukacs, _Characteristic Function_. Charles Griffin & Company Limited, 1970.
  * [36] T. M. Cover, _Elements of Information Theory_. John Wiley & Sons, 1999.
  * [37] Y. Polyanskiy and Y. Wu, _Information Theory: From Coding to Learning_. Cambridge University Press, 2022+ (book draft).
  * [38] N. V. Sabu, N. Varshney, and A. K. Gupta, “3-D diffusive molecular communication with two fully-absorbing receivers: Hitting probability and performance analysis,” _IEEE Transactions on Molecular, Biological and Multi-Scale Communications_ , vol. 6, no. 3, pp. 244–249, 2020.
  * [39] J. W. Kwak, H. B. Yilmaz, N. Farsad, C.-B. Chae, and A. J. Goldsmith, “Two-way molecular communications,” _IEEE Transactions on Communications_ , vol. 68, no. 6, pp. 3550–3563, 2020.
  * [40] L. Brand, S. Lotter, V. Jamali, R. Schober, and M. Schäfer, “Area rate efficiency in multi-link molecular communications,” _IEEE Transactions on Molecular, Biological and Multi-Scale Communications_ , vol. 9, no. 4, pp. 391–407, 2023.
  * [41] F. Vakilipoor, A. N. Ansari, and M. Magarini, “Heuristic barycenter modeling of fully absorbing receivers in diffusive molecular communication channels,” _IEEE Transactions on Communications_ , vol. 72, no. 1, pp. 133–145, 2024.
  * [42] R. E. Gaunt, “A simple proof of the characteristic function of student’s $t$-distribution,” _Communications in Statistics-Theory and Methods_ , vol. 50, no. 14, pp. 3380–3383, 2021.
  * [43] H.-Y. Lee, H.-J. Park, and H.-M. Kim, “A clarification of the Cauchy distribution,” _Communications for Statistical Applications and Methods_ , vol. 21, no. 2, pp. 183–191, 2014.
  * [44] S. Kotz and S. Nadarajah, _Multivariate t-distributions and Their Applications_. Cambridge University Press, 2004.
  * [45] Y. L. Luke, _The Special Functions and Their Approximations, Vol. 2_. Elsevier, 1969.
  * [46] S. Redner, _A Guide to First-Passage Processes_. Cambridge University Press, 2001.