# First Arrival Position in Molecular Communication Via Generator of Diffusion
Semigroup

Yen-Chi Lee  Hon Hai (Foxconn) Research Institute
Taipei, Taiwan
<EMAIL_ADDRESS>Jen-Ming Wu  Hon Hai (Foxconn) Research Institute
Taipei, Taiwan
<EMAIL_ADDRESS>Min-Hsiu Hsieh  Hon Hai (Foxconn) Research Institute
Taipei, Taiwan
<EMAIL_ADDRESS>4th Given Name Surname  dept. name of organization
(of Aff.)
name of organization (of Aff.)
City, Country
email address or ORCID  5th Given Name Surname  dept. name of organization (of
Aff.)
name of organization (of Aff.)
City, Country
email address or ORCID  6th Given Name Surname  dept. name of organization (of
Aff.)
name of organization (of Aff.)
City, Country
email address or ORCID

###### Abstract

We consider the problem of characterizing the first arrival position (FAP)
density in molecular communication (MC) with a diffusion-advection channel
that permits a constant drift velocity pointed to arbitrary direction. The
advantage of FAP modulation lies in the fact that it could encode more
information into higher dimensional spatial variables, compared to other
modulation techniques using time or molecule numbers. However, effective
methods to characterize the FAP density in a general framework do not exist.
In this paper, we devise a methodology that fully resolves the FAP density
with planar absorbing receivers in arbitrary dimensions. Our work recovers
existing results of FAP in 2D and 3D as special cases. The key insight of our
approach is to remove the time dependence of the MC system evolution based on
the generator of diffusion semigroups.

###### Index Terms:

Molecular communication, diffusion, first arrival time (FAT), first arrival
position (FAP), semigroup, generator.

##  I Introduction

Since ancient times, the problem of conveying information over a distance has
been an important issue in human history. Modern communication systems resolve
this problem with electromagnetic (EM) signals. However, EM-based
communication techniques are challenging for tiny (for instance, nano-scale)
devices; if not feasible, due to the constraints such as the wavelength,
antenna size, or energy issues [1, 2].

In molecular communication (MC) systems, small molecules111We omit the shape
of a message molecule and think of it as a mass point as most works in
molecular communication do. Thus we use both terms “particle” and “molecular”
interchangeably. called message molecules (MMs) act as chemical signals
conveying the information [3, 4]. After the information-carrying particles are
released in the channel, a propagation mechanism is necessary for transporting
them to the receiver. This mechanism can be diffusion-based [5], flow-
based[6], or an engineered transport system such as molecular motors [7, 8].
Among these different propagation mechanisms, diffusion-based MC, sometimes in
combination with advection and chemical reaction networks, has been the most
prevalent approach for MC systems, and will be the main focus of study in this
paper.

In order to transmit information in a diffusion-based MC system, we can
modulate different physical properties of the MMs [9, 10, 11, 12]. Signaling
molecules that reach the vicinity of the receiver can be observed and
processed by the receiver to extract the information that is necessary for
performing detection and decoding [13]. The reception mechanism of a MC
receiver can be categorized into two classes: i) passive reception, and ii)
active reception. The simplest active reception is the fully-absorbing
receiver [14] which has the ability to measure the time and position of
arrival of each molecule, and to remove the MM right after it is received. For
MC systems with a fully-absorbing receiver, how to completely characterize the
arrival time or the arrival position is no doubt a very important issue.

In Section I-A, we briefly review some important results concerning first
arrival time (FAT) in MC literature. Until very recently, the first arrival
position (FAP) information was also introduced into the MC realm, see [15,
16]. We will give a short introduction to the current status of FAP results in
Section I-B.

###  I-A First Arrival Time

When information is encoded on the time of release of MM, we refer to this
subclass of MC channels as molecular timing channels (MTCs). The main
difference between timing-based modulation and other modulations, such as the
number or type modulation, is that: the channel input is continuous instead of
discrete.

The MTC model was first proposed in [17]. In its simplest form, a MTC can be
realized through a single MM released by the transmitter (Tx) at time
$t_{\text{in}}$ with information encoded on this release time. The MM goes
through some random propagation and arrive at the receiver (Rx) at time
$t_{\text{out}}$. We have

$t_{\text{out}}=t_{\text{in}}+t_{\text{n}},$ (1)

where $t_{\text{n}}$ is some random delay due to MM propagation mechanisms.

In [18], an MTC with additive inverse Gaussian noise (AIGN) was introduced
into MC realm, where $t_{\text{n}}$ in Eq. (1) is inverse Gaussian
distributed. The authors of [18] derived upper and lower bounds on capacity
per channel use. Later in [19], tighter bounds on capacity of the same MTC
with AIGN noise were derived. The capacity-achieving MM input distribution was
also characterized. In [20, 21, 22], single particle release is extended to
multiple, say $M$ MM releases. In this scenario, information is encoded in a
vector t, in which each component is the time of release of one of the $M$
MMs. For multiple particle case, the channel model changes to:

$\textbf{t}_{\text{out}}=\text{sort}(\textbf{t}_{\text{in}}+\textbf{t}_{\text{n}}),$
(2)

where $\textbf{t}_{\text{n}}$ is an $\mathbb{R}^{M}$ vector of random delays
associated with each MM. Notice that $\text{sort}(\mathbf{t})$ is the sort
operator that permutes the input vector into an ascending order. The sort
operation is necessary when time information of each MM is considered since
the MMs may arrive out of order.

###  I-B First Arrival Position

Apart from the first arrival time, there is another important degree of
freedom to carry information: the first arrival position. The simplest FAP
channel can be written as

$\mathbf{x}_{\text{out}}=\mathbf{x}_{\text{in}}+\mathbf{x}_{\text{n}},$ (3)

where $\mathbf{x}_{\text{in}}$, $\mathbf{x}_{\text{out}}$ and
$\mathbf{x}_{\text{n}}$ are all Euclidean vectors in $\mathbb{R}^{n}$ with $n$
being the dimension of the underlying space.

To the best of our knowledge, the first paper in MC society that mentions FAP
as an information carrying property is [15], the density function of FAP is
derived in closed-form for both pure diffusion and diffusion with a constant
drift in 3D spaces. However, this approach is limited in scope when higher
dimensional FAP is considered. Later, the authors in [16] use Green’s function
and the method of images to derive the FAP density function for 2D MC systems,
where the fluid medium has a constant drift restricted to the longitudinal
direction from the transmitter towards the receiver. However, the method of
images in [16] cannot be generalized to the setting of arbitrary drift
velocity field, which is commonly encountered in fluid environments, because
it uses symmetry heavily, but an arbitrary direction of drift will break down
the required symmetry. Finally, the capacity of molecular position channels
remains completely unknown except some very specific cases, see [16] for
capacity of $M$-ary modulation scheme.

###  I-C Our Contributions

Although most works in MC consider FAT modulation for absorbing receivers,
there are at least two reasons why FAP is preferable.

  * •

For each message particle, the FAT information is only one-dimensional, while
the FAP could have higher dimensions to carry information. Hence, the capacity
of FAP could be larger then FAT per single message particle.

  * •

Practically, “guard intervals” between two consecutive transmissions are
required to avoid cross-over effects, as described in Eq. (2). Due to these
guard intervals, the total transmission time for FAT-type modulation will
increase roughly proportional to the number of MMs to be transmitted.
Consequently, for applications in which the time efficiency plays an important
role, the FAP-type modulation is arguably a better solution.

The main contributions of this paper are two-fold. 1) We relax the constraint
imposed on the drift direction in previous works. 2) We introduce the
semigroup approach to facilitate characterization of the FAP density for
planar absorbing receivers in higher dimensional diffusion channels. Our
methodology unifies the process toward finding FAP in 2D and 3D with a solid
theoretical backup (see Section III for details). Technically, the semigroup
approach avoids complicated time integration in finding FAP. We believe that
our work opens the possibility toward general resolution of the capacity of
molecular position channels.

###  I-D Structure of this Paper

The remainder of this paper is structured as follows. Section II describes the
system model we choose, and some related papers working on FAP. Section III
presents the new methodology we proposed to find FAP. Finally, concluding
remarks are made in Section IV.

##  II System Model and Related Works

###  II-A System Model

We consider a MC system located in an $n$-dimensional Euclidean space,
consists of a point Tx and a planar absorbing Rx. The distance between Tx and
Rx is $d$. Notice that we can always arrange the basis of $\mathbb{R}^{n}$ so
that the last vector in the basis set is parallel to the transmission
direction. Without loss of generality, we consider that the transmitter is a
point located at the Cartesian coordinate $(0,\cdots,0,d)$, and the receiver
is an infinite large absorbing plane located at
$\\{x\in\mathbb{R}^{n}:x_{n}=0\\}$. Abstract MC system figures for 2D and 3D
spaces are illustrated in Figure 1 and Figure 2, respectively. The physical
channel of our MC model is composed of a diffusive fluid medium with a
constant drift.

In theoretical MC literature, there are two different viewpoints to model
diffusion channels. The macroscopic viewpoint uses diffusion equation (aka
heat equation) to capture the evolution of the continuous distribution of
concentration of MMs. While in the microscopic viewpoint, individual message
molecule is monitored using its own trajectory. For the equivalence of these
two viewpoints, readers can refer to classical potential theory as in [23], or
a good review paper in MC such as [9] or [13].

In the macroscopic viewpoint, the physical channel is modeled by Eq. (4)
below. The most important quantity in macroscopic MC system analysis is the
concentration field222Here we use the word “field” to mean a function of space
and time. $c(\mathbf{r},t)$ of message molecules at spatial position
$\mathbf{r}$ and time $t$. By Fick’s law of diffusion [13], the evolution of
$c(\mathbf{r},t)$ can be captured by the diffusion-advection equation:

$\partial_{t}c\left(\mathbf{r},t;\mathbf{r}_{0}\right)+\mathbf{v}(\mathbf{r},t)\cdot\nabla
c\left(\mathbf{r},t;\mathbf{r}_{0}\right)=D\nabla^{2}c\left(\mathbf{r},t;\mathbf{r}_{0}\right),$
(4)

where $\mathbf{r}_{0}$ is the point where the diffusion starts, $\mathbf{v}$
is the velocity field of the fluid medium which is assumed to be
incompressible, $\nabla$ and $\nabla^{2}$ are the gradient and the Laplace
operators, respectively, and $D$ is the diffusion coefficient. The value of
$D$ is determined by the temperature, the fluid viscosity, and the molecule’s
Stokes radius, see [9].

On the other hand in the microscopic viewpoint, a common assumption is that
the trajectory of each molecule can be well distinguished. A mainstream model
in MC for the trajectory $X_{t}$ of a MM is the Itô diffusion process. An Itô
diffusion in Euclidean space $\mathbb{R}^{n}$ is a stochastic process
satisfying a stochastic differential equation (SDE) of the form

$dX_{t}=b(X_{t})dt+\sigma(X_{t})dB_{t},$ (5)

where $B_{t}$ is an $n$ dimensional standard Brownian motion. Throughout this
paper, we assume that $b(\cdot)$ and $\sigma(\cdot)$ are both constants and
can be determined by the properties of fluid and MMs. We further assume that
the receiver for our MC system is perfectly absorbing and has the ability to
correctly measure the time and position at first-arrival for each individual
molecule [16].

In the following, we denote the initial position of a MM by x, the final (or
received) position by y, and the FAP distribution by a conditional density
function $f(\mathbf{y}|\mathbf{x})$. The main problem is to characterize the
function $f(\mathbf{y}|\mathbf{x})$.

###  II-B Related Works

####  II-B1 2D FAP density

The authors of [16] solved this problem in 2D partially333Using the image
method in [16], we can only tackle with drifts in longitudinal direction
parallel to the direction of transmission, but not in arbitrary direction..

Figure 1: System Model of 2D First Arrival Position Channels.

In [16], a 2D MC system with a constant drift velocity $\mathbf{v}$ pointing
precisely to the receiver is considered. The system model is illustrated in
Figure 1. The authors of [16] adopt the method of images to construct Green’s
function for an absorbing boundary. The resulting FAP density function can be
found as in [24, Eq. (42)] of Appendix A of the online version of this paper.
Notice that we cannot recover the FAP expression as appeared in [16, Eq.
(19)]. Instead, we have provided a self-contained proof in Appendix A of [24].

####  II-B2 3D FAP density

In [15], a 3D MC system with a constant drift velocity $\mathbf{v}$ pointing
to arbitrary direction is considered. The system model is illustrated in
Figure 2.

Figure 2: System Model of 3D First Arrival Position Channels.

The authors of [15] use Dynkin’s formula to link the microscopic properties at
first arrival time to the solution of macroscopic diffusion equation. A
closed-form FAP density function in 3D is obtained in [15], which has the
following form:

$\displaystyle f_{Y|X}(\mathbf{y}|\mathbf{x})$ (6) $\displaystyle=$
$\displaystyle\frac{\lambda}{2\pi}\exp\left\\{-\dfrac{v_{3}\lambda}{\sigma^{2}}\right\\}\exp\left\\{\dfrac{v_{1}}{\sigma^{2}}\left(\xi-
x_{1}\right)+\dfrac{v_{2}}{\sigma^{2}}\left(\eta-x_{2}\right)\right\\}$
$\displaystyle\cdot\exp\left\\{-\dfrac{|v|}{\sigma^{2}}\left\lVert\mathbf{y}-\mathbf{x}\right\rVert\right\\}\left[\frac{1+\frac{|v|}{\sigma^{2}}\left\lVert\mathbf{y}-\mathbf{x}\right\rVert}{\left\lVert\mathbf{y}-\mathbf{x}\right\rVert^{3}}\right],$

where $\sigma^{2}=2D$ is the microscopic diffusion coefficient, $d=\lambda$ is
the distance between transmitter plane and receiver plane, $v_{3}$ is the
component of drift velocity parallel (or longitudinal) to the information
transmission direction, and $v_{1}$ and $v_{2}$ are the drift components
perpendicular (or transverse) to the transmission. Here we use the symbol
$\left\lVert\cdot\right\rVert$ to represent the Euclidean norm. Namely,
$\left\lVert\mathbf{y}-\mathbf{x}\right\rVert=\sqrt{(\xi-x_{1})^{2}+(\eta-
x_{2})^{2}+\lambda^{2}},$ where $\textbf{x}=(x_{1},x_{2},\lambda)$,
$\textbf{y}=(\xi,\eta,0)$ are position vectors in $\mathbb{R}^{3}$.

##  III Finding FAP Via Generator of Diffusion Semigroup

By suitably examine macroscopic and microscopic relations for the diffusion
process, we come up with new ideas which can simplify the procedure towards
finding FAP density. After a careful examination, one can realize that FAP
density itself has no time dependency, so there must exist some redundant
calculations in old methods playing with heat equations, for instance the
method used in [16]. The key insight of our approach is to remove the time
dependence of the MC system evolution based on the generator of diffusion
semigroups.

Our proposed semigroup approach towards finding the FAP density can be
summarized in three main steps:

  * 1.

By utilizing the generator of diffusion semigroup, we can remove the time
variable from the very beginning. This can be done because of the Markov
property of diffusion semigroup for our MC system.

  * 2.

Look up representation formulas from partial differential equation (PDE)
literature for the corresponding elliptic boundary value problems (BVP) we
obtained in Step 1. Since there are plentiful well known results for solutions
of elliptic PDEs (even in high dimensional cases, see [25, 26]), we do not
need to derive representation formulas from scratch.

  * 3.

Via properly interpretating the Dynkin’s formula, we claim that the FAP
density can be obtained directly by the relation
$f_{Y|X}(\textbf{y}|\textbf{x})=\left|\dfrac{\partial
G(\textbf{y},\textbf{x})}{\partial n_{\textbf{y}}}\right|,$ once the elliptic
Green’s function $G(\textbf{y},\textbf{x})$ is known. Through this process, we
can bypass the complicated integration step as appeared in [16, Eq. (17)].

The comparisons between the old and new methods are listed in Table I. For the
sake of completeness, we present the old method concisely in Appendix B of
[24] in a self-contained fashion.

Table I: Comparison Between Old and New Methods
for Finding FAP Density |  Old Method (with $t$) |  New Method (without $t$)
---|---|---
Step 1 |  Finding the free space Green’s function for parabolic PDE. |  Removing time $t$ from the beginning by considering the generator of diffusion.
Step 2 |  Solving for absorbing Green’s function based on image method. |  Looking up the solution forms for certain types of elliptic BVPs.
Step 3 |  Calculating FAP density flux and doing time integration. |  Interpretating the Dynkin’s formula and obtaining the FAP density directly.

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

In physics, an Itô diffusion is used to describe the Brownian motion of a
particle subjected to a potential in a fluid medium. In mathematics, an $n$
dimensional Itô diffusion444We merely consider time-homogeneous Itô diffusion
throughout this paper. That is, $b(X_{t})$ and $\sigma(X_{t})$ does not depend
on $t$ explicitly. is a stochastic process satisfying a specific type of SDE
of the form $dX_{t}=b(X_{t})dt+\sigma(X_{t})dB_{t},$ as described in Eq. (5).
Throughout this paper, we assume that $b(\cdot)$ and $\sigma(\cdot)$ are both
constants.

To each Itô diffusion $X_{t}$, we can associate a corresponding operator
called infinitestimal generator, or simply generator, which is defined as
follows. Let $\mathcal{D}(A)$ be the domain of the generator $A$. We define
$\mathbb{E}^{\mathbf{x}}[f(X_{t})]:=\mathbb{E}[f(X_{t})|X_{0}=\mathbf{x}],$
where the notation $\mathbb{E}^{\mathbf{x}}[\cdot]$ stands for taking
expectation conditioned on $X_{0}=\mathbf{x}$. The generator $A$ of a process
$X_{t}$ can be then defined as

$Af(\mathbf{x})=\lim\limits_{t\searrow
0}\dfrac{\mathbb{E}^{\mathbf{x}}[f(X_{t})]-f(\mathbf{x})}{t}\quad\text{for}\quad
f\in\mathcal{D}(A).$ (7)

For time-homogeneous Itô process $X_{t}$, the time evolution is a Markov
process. By letting $T_{t}:=\mathbb{E}^{\mathbf{x}}[f(X_{t})],$ we obtain a
semigroup of operators: $T=(T_{t})_{t\geq 0}$. Thinking from another
perspective, $Af=\lim\limits_{t\searrow 0}\dfrac{T_{t}f-f}{t}$ can be regarded
as the linear increment term (omitting all higher order terms) of the
semigroup evolution.

The next step is to calculate the generator $A$ explicitly from the Itô
diffusion SDE: $dX_{t}=b(X_{t})dt+\sigma(X_{t})dB_{t}$. Supposing $f$ is of
class $C^{2}$, by Taylor expansion and Itô’s formula555Itô’s formula and
Taylor expansion of an Itô process are standard results in stochastic
analysis, see [27, 28]., we have:

$\displaystyle\begin{split}df\left(X_{t}\right)=\
&f^{\prime}\left(X_{t}\right)dX_{t}+\frac{1}{2}f^{\prime\prime}\left(X_{t}\right)d\langle
X\rangle_{t}\\\ =\
&\left[b\left(X_{t}\right)f^{\prime}\left(X_{t}\right)+\frac{\sigma^{2}\left(X_{t}\right)}{2}f^{\prime\prime}\left(X_{t}\right)\right]dt\\\
&+f^{\prime}\left(X_{t}\right)\sigma\left(X_{t}\right)dB_{t}.\end{split}$ (8)

The notation $d\langle X\rangle_{t}$ represents the quadrature variation of
random process $X_{t}$, see [28]. Directly plugging Eq. (8) into Eq. (7)
yields

$Af(\mathbf{x})=b(\mathbf{x})f^{\prime}(\mathbf{x})+\dfrac{\sigma^{2}(\mathbf{x})}{2}f^{\prime\prime}(\mathbf{x}).$
(9)

Eq. (9) is the generator of Itô diffusion corresponding to our MC system.

###  III-B Boundary Value Problem for Elliptic PDE

For MC systems located in a time-invariant and spatial-homogeneous fluid
environment, we can restrict ourselves to the case that both $b(\cdot)$ and
$\sigma(\cdot)$ are constants. Consider a PDE BVP with unknown function $u$:

$\left\\{\begin{array}[]{ll}Au=0&\text{ in }\Omega\\\ u=g&\text{ on
}\partial\Omega\end{array}\right.$ (10)

where $g$ is some prescribed boundary data. Note that the right hand side of
the first equation in (10) is zero. This means that there are no particle
production and annihilation process in the fluid channel. For elliptic
differential operators, BVPs are well studied, and the solution of (10) under
common boundary conditions can be found in literature of PDE, for instance
[25].

Let us discuss a simple example in more detail. Consider a half-space 3D
domain with $A=\nabla^{2}$ and a boundary condition:

$u(x,y,z)=g(x,y)\quad\text{at}\quad z=0.$ (11)

The solution of BVP (10) with (11) as boundary condition is well known and can
be found in standard handbooks of PDE:

$u(x,y,z)=\frac{1}{2\pi}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\frac{zg(\xi,\eta)d\xi
d\eta}{\left[(x-\xi)^{2}+(y-\eta)^{2}+z^{2}\right]^{3/2}}.$ (12)

Eq. (12) is called an integral representation formula in PDE language. Using
Green’s function, we can rewrite (12) as

$u(x,y,z)=\int_{\partial\Omega}\left|\dfrac{\partial
G(\xi,\eta,0;x,y,z)}{\partial\mathbf{n}}\right|g(\xi,\eta)d\xi d\eta,$ (13)

where $\mathbf{n}$ is the normal unit vector of the boundary. We claim that
the term $\left|\dfrac{\partial
G(\xi,\eta,0;x,y,z)}{\partial\mathbf{n}}\right|$ is already the FAP density on
the receiving plane. In this example, the BVP corresponds to no-drift pure
diffusion. We will prove the above claim in Section III-C. The calculations
for diffusion-with-drift case in 2D space are presented in Appendix A of [24].

###  III-C Dynkin’s Formula and FAP Density

In this subsection, we will use the Dynkin’s formula to clarify the relations
between macroscopic and microscopic viewpoints of diffusion phenomena. During
this clarification process, the FAP density is obtained incidentally.

Recall that the definition of an infinitestimal generator is
$Af(\mathbf{x})=\lim\limits_{t\searrow
0}\dfrac{\mathbb{E}^{\mathbf{x}}[f(X_{t})]-f(\mathbf{x})}{t}.$ Rewriting the
above equation into an integral form, we have

$\mathbb{E}^{\mathbf{x}}\left[f\left(X_{t}\right)\right]=f(\mathbf{x})+\mathbb{E}^{\mathbf{x}}\left[\int_{0}^{t}Af\left(X_{s}\right)ds\right].$
(14)

Note that the time $t$ appeared in the above formula is a deterministic
variable, not a random variable.

In stochastic analysis, Dynkin’s formula is a theorem giving information about
a diffusion process at a stopping time $\tau$ (which is a random variable
[28]). Let $f$ be of class $C^{2}$ with compact support and $\tau$ be a
stopping time with $\mathbb{E}^{\mathbf{x}}[\tau]<+\infty$. The Dynkin’s
formula can be stated as

$\mathbb{E}^{\mathbf{x}}\left[f\left(X_{\tau}\right)\right]=f(\mathbf{x})+\mathbb{E}^{\mathbf{x}}\left[\int_{0}^{\tau}Af\left(X_{s}\right)ds\right].$
(15)

With all these background materials at hand, we now show that how arbitrary
FAP density is obtained by our newly proposed approach in principle. Let $g$
be a smooth data defined on the boundary (i.e. the receiver plane) of the
domain. We can write

$\mathbb{E}^{\mathbf{x}}[g(X_{\tau})]=\mathbb{E}[g(X_{\tau})\mid
X_{0}=\mathbf{x}]=\int_{\partial\Omega}f_{Y|X}(\mathbf{y}|\mathbf{x})g(\mathbf{y})d\mathbf{y}$
(16)

by the definition of conditional probability density function (PDF). In Eq.
(16), $\mathbf{x}\in\Omega$ is the starting point of the diffusion; and since
$\tau$ is the hitting time, the hitting position
$\mathbf{y}\in\partial\Omega$. Notice that the conditional PDF $f_{Y|X}$ is
exactly the desired FAP density distributed on the receiver boundary
$\partial\Omega$.

Consider an Itô diffusion $X_{t}$ and its corresponding generator $A$. Suppose
we have a solution $u(\mathbf{x})$ of the BVP problem (10) for some prescribed
$g$. Then we plug $f(\mathbf{x})=u(\mathbf{x})$ into Eq. (15). The last term
of Eq. (15) becomes

$\int_{0}^{\tau}Au(X_{s})ds=0.$ (17)

This is because $0<s<\tau$ and $\tau$ is the first hitting time, so that
$X_{s}$ lies in $\Omega$, making $u(X_{s})=0$. On the other hand, the term on
the left hand side of Eq. (15) becomes $\mathbb{E}^{\mathbf{x}}[g(X_{\tau})]$
because $u(\mathbf{x})$ coincides with $g(\mathbf{x})$ on the boundary.
Combining these two facts, we can get

$\mathbb{E}^{\mathbf{x}}[g(X_{\tau})]=u(\mathbf{x})=\int_{\partial\Omega}f_{Y|X}(\mathbf{y}|\mathbf{x})g(\mathbf{y})d\mathbf{y}$
(18)

for any $\mathbf{x}\in\Omega$.

Generically, for elliptic-type BVP defined in a domain $V$ with boundary
$S=\partial V$, we have

$u(\mathbf{x})=\int_{V}\Phi(\mathbf{y})G(\mathbf{x},\mathbf{y})dV_{\mathbf{y}}+\int_{S}g(\mathbf{y})H(\mathbf{x},\mathbf{y})dS_{\mathbf{y}}$
(19)

where $\Phi$ represents the source term, $dV_{\mathbf{y}}$ is the volume
element, $dS_{\mathbf{y}}$ is the surface element, and the term
$H(\mathbf{x},\mathbf{y})$ depends on which type of boundary conditions we are
considering. If there is no molecule source reproducing MMs in the domain, we
can set $\Phi(\mathbf{y})=0$. For the purpose of deriving FAP density, we
consider the Dirichlet type boundary, which corresponds to

$H(\mathbf{x},\mathbf{y})=-\dfrac{\partial
G}{\partial\mathbf{n}_{\mathbf{y}}}(\mathbf{x},\mathbf{y}).$ (20)

As a result, formula (19) reduces to

$u(\mathbf{x})=\int_{S}-\dfrac{\partial
G}{\partial\mathbf{n}_{\mathbf{y}}}(\mathbf{x},\mathbf{y})g(\mathbf{y})dS_{\mathbf{y}}.$
(21)

Finally, by comparing Eq. (18) and Eq. (21), we reveals an important relation
between the FAP density and elliptic Green’s function:

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

Note that the left hand side of Eq. (22) always has a positive sign, so we add
an absolute value on the right hand side for emphasis.

To conclude, Eq. (22) can be used to solve for FAP density as long as we know
the Green’s function. This method in principle works for higher $n$
dimensional diffusion. By using generator of diffusion semigroup, we transform
the original heat-type BVP (with time variable) into elliptic-type BVP
(without time variable). Then we can look up solutions for the corresponding
elliptic Green’s functions in literature, such as [25, 26].

##  IV Conclusions

In this paper, we fully resolve the FAP density problem for planar absorbing
receivers in MC. Our approach is based on the theory of diffusion semigroups.
The comparisons between the existing and our methods are listed in Table I.
Through this semigroup approach, we integrate the existing FAP results in MC,
and unified the process toward finding FAP in 2D and 3D.

Secondly, this is the first paper in MC to clarify in a rigorous way the
mathematical link between the macroscopic and microscopic viewpoints for the
diffusion mechanism in diffusion-based MC systems. Previous works in MC treat
these two viewpoints separately. We hope this clarification could bring new
insights into future theoretical MC system analysis.

Finally, we provided a self-contained proof for the formula proposed in [16]
for 2D FAP density problem. We also make some improvements via relieving the
restriction that the drift direction must point exactly toward the receiver.
The calculation details are presented in Appendix A of [24].

## References

  * [1] T. Nakano, A. W. Eckford, and T. Haraguchi, _Molecular communication_. Cambridge University Press, 2013.
  * [2] 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.
  * [3] I. F. Akyildiz, F. Brunetti, and C. Blázquez, “Nanonetworks: A new communication paradigm,” _Computer Networks_ , vol. 52, no. 12, pp. 2260–2279, 2008.
  * [4] P.-C. Yeh, K.-C. Chen, Y.-C. Lee, L.-S. Meng, P.-J. Shih, P.-Y. Ko, W.-A. Lin, and C.-H. Lee, “A new frontier of wireless communication theory: diffusion-based molecular communications,” _IEEE Wireless Communications_ , vol. 19, no. 5, pp. 28–35, 2012.
  * [5] H. C. Berg, _Random walks in biology_. Princeton University Press, 2018.
  * [6] 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.
  * [7] 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 _Fourth Annual IEEE International Conference on Pervasive Computing and Communications Workshops (PERCOMW’06)_. IEEE, 2006, pp. 6–pp.
  * [8] M. Gregori and I. F. Akyildiz, “A new nanonetwork architecture using flagellated bacteria and catalytic nanomotors,” _IEEE Journal on selected areas in communications_ , vol. 28, no. 4, pp. 612–619, 2010.
  * [9] 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.
  * [10] M. S. Kuran, H. B. Yilmaz, T. Tugcu, and I. F. Akyildiz, “Modulation techniques for communication via diffusion in nanonetworks,” in _2011 IEEE international conference on communications (ICC)_. IEEE, 2011, pp. 1–5.
  * [11] 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\.
  * [12] Y.-P. Hsieh, Y.-C. Lee, P.-J. Shih, P.-C. Yeh, and K.-C. Chen, “On the asynchronous information embedding for event-driven systems in molecular communications,” _Nano Communication Networks_ , vol. 4, no. 1, pp. 2–13, 2013.
  * [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] H. B. Yilmaz, A. C. Heren, T. Tugcu, 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\.
  * [15] Y.-C. Lee, C.-C. Chen, P.-C. Yeh, and C.-H. Lee, “Distribution of first arrival position in molecular communication,” in _2016 IEEE International Symposium on Information Theory (ISIT)_. IEEE, 2016, pp. 1033–1037.
  * [16] 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.
  * [17] A. W. Eckford, “Nanoscale communication with brownian motion,” in _2007 41st Annual Conference on Information Sciences and Systems_. IEEE, 2007, pp. 160–165.
  * [18] 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.
  * [19] H. Li, S. M. Moser, and D. Guo, “Capacity of the memoryless additive inverse gaussian noise channel,” _IEEE Journal on Selected Areas in Communications_ , vol. 32, no. 12, pp. 2315–2329, 2014.
  * [20] A. W. Eckford, “Molecular communication: Physically realistic models and achievable information rates,” _arXiv preprint arXiv:0812.1554_ , 2008.
  * [21] C. Rose and I. S. Mian, “Inscribed matter communication: Part i,” _IEEE Transactions on Molecular, Biological and Multi-Scale Communications_ , vol. 2, no. 2, pp. 209–227, 2016.
  * [22] ——, “Inscribed matter communication: Part ii,” _IEEE Transactions on Molecular, Biological and Multi-Scale Communications_ , vol. 2, no. 2, pp. 228–239, 2016.
  * [23] J. L. Doob and J. Doob, _Classical potential theory and its probabilistic counterpart_. Springer, 1984, vol. 549\.
  * [24] Y.-C. Lee, J.-M. Wu, and M.-H. Hsieh, “First arrival position in molecular communication via generator of diffusion semigroup,” _arXiv preprint arXiv:2201.04476_ , 2022.
  * [25] A. D. Polyanin, _Handbook of linear partial differential equations for engineers and scientists_. Chapman and hall/crc, 2001.
  * [26] P. Grohs and L. Herrmann, “Deep neural network approximation for high-dimensional elliptic pdes with boundary conditions,” _arXiv preprint arXiv:2007.05384_ , 2020.
  * [27] B. Øksendal, _Stochastic differential equations: an introduction with applications_. Springer Science & Business Media, 2013.
  * [28] O. Calin, _An informal introduction to stochastic calculus with applications_. World Scientific, 2015.
  * [29] A. Tikhonov and A. Samarskii, “Methods of mathematical modeling and computational diagnostics,” _Methods of Mathematical Modeling and Computational Diagnostics_ , 1990.
  * [30] J. D. Jackson, “Classical electrodynamics,” 1999.
  * [31] I. S. Gradshteyn, I. M. Ryzhik, and R. H. Romer, “Tables of integrals, series, and products,” 1988.

##  Appendix A Calculation of FAP Density in 2D Space Using Semigroup Approach

In this appendix section, we demonstrate how to calculate FAP density using
semigroup approach proposed in Section III.

For two dimensional (2D) FAP density, we consider a 2D Itô diffusion $X_{t}$
with its semigroup generator $A$ as

$A=\sum_{i=1}^{2}v_{i}\dfrac{\partial}{\partial
x_{i}}+\dfrac{\sigma^{2}}{2}\sum_{i=1}^{2}\dfrac{\partial^{2}}{\partial
x_{i}^{2}}.$ (23)

In the above equation, the term $v_{i}$ stands for the $i$-th component of the
drift velocity $\mathbf{v}$. Let us denote the Laplacian operator in 2D by
$\Delta_{2}$ and consider the following BVP in Cartesian coordinate:

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

where $u$ is a (dummy) unknown function. Set the domain of BVP (24) to be

$\Omega=\mathbb{R}^{2}\cap\left\\{x_{2}>0\right\\},$ (25)

and the boundary to be

$\partial\Omega=\mathbb{R}^{2}\cap\left\\{x_{2}=0\right\\}.$ (26)

Here we use the notation $x_{j}$ to denote the $j$-th component of the
position vector $\mathbf{x}=(x_{1},\cdots,x_{n})$. For the purpose of
compatibility to higher dimensional cases, we arrange the basis order of the
Cartesian coordinate so that the direction of the last component is parallel
to the transmission direction.

The BVP (24) is not in a common resolved form in standard PDE solution tables
such as [25]. However, we can use a change of variables to transform (24) into
a Helmholtz equation. Before proceeding further, we temporarily set
$\sigma^{2}=1$ to avoid cumbersome calculations. After obtaining the final
result, we will then replace $v_{i}$ with $\frac{v_{i}}{\sigma^{2}}$ to
recover the solution for general $\sigma^{2}$ which is not necessarily equal
to $1$.

Suppose $u$ is the solution of BVP (24) with $\sigma^{2}=1$. That is, $u$
satisfies $\sum_{i=1}^{2}v_{i}\frac{\partial u}{\partial
x_{i}}+\frac{1}{2}\Delta_{2}u=0$ in the domain $\Omega$, and coincides with
the function $g$ on the boundary $\partial\Omega$. We define a drift factor
$\gamma$ to facilitate later calculations; the drift factor is defined as
follows:

$\gamma(\mathbf{x}):=\exp\\{\mathbf{v}\cdot\mathbf{x}\\}=\exp\left\\{v_{1}x_{1}+v_{2}x_{2}\right\\},$
(27)

where $\mathbf{v}\cdot\mathbf{x}$ denotes the standard inner product in
$\mathbb{R}^{2}$. We consider an alternative function
$w(\mathbf{x})=\gamma(\mathbf{x})u(\mathbf{x})$. Plugging
$u(\mathbf{x})=\frac{w(\mathbf{x})}{\gamma(\mathbf{x})}$ into BVP (24) yields

$\frac{1}{\gamma(\mathbf{x})}\left[\Delta_{2}w(\mathbf{x})-\left(v_{1}^{2}+v_{2}^{2}\right)w(\mathbf{x})\right]=0.$
(28)

Next, by letting

$s=|\mathbf{v}|=\sqrt{v_{1}^{2}+v_{2}^{2}},$ (29)

we can deduce that $w$ satisfies the Helmholtz equation:

$\Delta_{2}w-s^{2}w=0\quad(s\geq 0).$ (30)

The Helmholtz equation can be regarded as the eigenvalue problem of Laplacian
operator; its solution form can be found in standard PDE books [25]. To solve
$w$, we consider another BVP described by the Helmholtz equation:

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

where $H(\cdot)$ is the Helmholtz operator. Note that new boundary data
$\tilde{g}$ satisfies

$\tilde{g}(\mathbf{y})=e^{v_{1}y_{1}}g(\mathbf{y})\quad\text{for}\quad\mathbf{y}\in\partial\Omega$
(32)

since $y_{2}=0$ on the boundary.

The next step is to look up a representation formula of solution of BVP (31)
in the literature. The following equations (33, 34, 35, 36) are all cited from
[25]. Let the domain be $-\infty<x_{1}<\infty$, $0\leq x_{2}<\infty$ and
consider a first-type boundary condition666First-type boundary condition, also
known as Dirichlet boundary condition, means that we prescribe the values of
the unknown function on the boundary. as:

$w(x_{1},0)=g(x_{1})\quad\text{ at }\quad x_{2}=0.$ (33)

The solution of BVP (31) can be written down using the following
representation formula:

$w(x_{1},x_{2})=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f(\xi)\left[\frac{\partial}{\partial\eta}G(x_{1},x_{2},\xi,\eta)\right]_{\eta=0}d\xi.$
(34)

where the Green’s function appeared in Eq. (34) has the form

$\begin{array}[]{c}G(x_{1},x_{2},\xi,\eta)=\dfrac{1}{2\pi}[K_{0}(s\rho_{1})-K_{0}(s\rho_{2})]\end{array}$
(35)

with

$\displaystyle\begin{split}\rho_{1}&=\sqrt{(x_{1}-\xi)^{2}+(x_{2}-\eta)^{2}},\\\
\rho_{2}&=\sqrt{(x_{2}-\xi)^{2}+(x_{2}+\eta)^{2}}.\end{split}$ (36)

Interested readers can refer to [29] for more details about the derivation of
this Green’s function.

Now we let $\mathbf{x}=(x_{1},x_{2})$, $\mathbf{y}=(\xi,\eta)$. Based on the
methodology presented in Section III, we know that

$w(\mathbf{x})=\int_{\partial\Omega}\left|\frac{\partial
G(\mathbf{x},\mathbf{y})}{\partial\mathbf{n}_{\mathbf{y}}}\right|\tilde{g}(\mathbf{y})d\mathbf{y}.$
(37)

Here, the notation $\left|\frac{\partial
G(\mathbf{x},\mathbf{y})}{\partial\mathbf{n}_{\mathbf{y}}}\right|$ is the
partial derivative of $G$ with respect to the unit vector normal to the
boundary $\partial\Omega$. Using the fact that $K_{0}^{\prime}(x)=-K_{1}(x)$,
we calculate

$\displaystyle\begin{split}&\left.\dfrac{\partial
G}{\partial\eta}\right|_{\eta=0}\\\ =&\left[\dfrac{\partial
G}{\partial\rho_{1}}\dfrac{\partial\rho_{1}}{\partial\eta}+\dfrac{\partial
G}{\partial\rho_{2}}\dfrac{\partial\rho_{2}}{\partial\eta}\right]_{\eta=0}\\\
=&\
\dfrac{-s}{2\pi}K_{1}\left(s\sqrt{(x_{1}-\xi)^{2}+x_{2}^{2}}\right)\dfrac{-x_{2}}{\sqrt{(x_{1}-\xi)^{2}+x_{2}^{2}}}\\\
&+\dfrac{s}{2\pi}K_{1}\left(s\sqrt{(x_{1}-\xi)^{2}+x_{2}^{2}}\right)\dfrac{x_{2}}{\sqrt{(x_{1}-\xi)^{2}+x_{2}^{2}}}\\\
=&\
\dfrac{|\mathbf{v}|x_{2}}{\pi}\dfrac{K_{1}\left(|\mathbf{v}|\sqrt{(x_{1}-\xi)^{2}+x_{2}^{2}}\right)}{\sqrt{(x_{1}-\xi)^{2}+x_{2}^{2}}}.\end{split}$
(38)

To obtain a representation formula for the solution $u$ of the original BVP
(24), we substitute $w(\mathbf{x})=\gamma(\mathbf{x})u(\mathbf{x})$ into Eq.
(37), yielding

$\begin{aligned}
u(\mathbf{x})=\frac{w(\mathbf{x})}{\gamma(\mathbf{x})}&=\exp\\{-\mathbf{v}\cdot\mathbf{x}\\}\int_{\partial\Omega}\left|\frac{\partial
G(\mathbf{x},\mathbf{y})}{\partial\mathbf{n}_{\mathbf{y}}}\right|\tilde{g}(\mathbf{y})d\mathbf{y}\\\
&=\int_{\partial\Omega}\exp\left\\{v_{1}y_{1}-v_{1}x_{1}-v_{2}x_{2}\right\\}\\\
&\quad\left[\dfrac{\partial
G}{\partial\rho_{1}}\dfrac{\partial\rho_{1}}{\partial\eta}+\dfrac{\partial
G}{\partial\rho_{2}}\dfrac{\partial\rho_{2}}{\partial\eta}\right]_{\eta=0}g(y_{1})dy_{1}\\\
&=\int_{\partial\Omega}K_{\mathbf{v}}(\mathbf{x},\xi)g(\xi)d\xi\end{aligned},$
(39)

where $K_{\mathbf{v}}(\mathbf{x},\xi)$ is the integral kernel which we want to
determine. Note that in the second equality of (39), we have used the relation
$\tilde{g}(\mathbf{y})=e^{\mathbf{v}\cdot\mathbf{x}}g(\mathbf{y})$.

Now we let $x_{2}=d>0$ to be the transmission distance between Tx and Rx. By
examining Eq. (38) and Eq. (39), we obtain an exact formula for the integral
kernel $K_{\mathbf{v}}(\mathbf{x},\xi)$:

$\displaystyle\begin{split}&K_{\mathbf{v}}(\mathbf{x},\xi)\\\ =&\
\exp\left\\{v_{1}\xi-v_{1}x_{1}-v_{2}d\right\\}\left[\dfrac{\partial
G}{\partial\rho_{1}}\dfrac{\partial\rho_{1}}{\partial\eta}+\dfrac{\partial
G}{\partial\rho_{2}}\dfrac{\partial\rho_{2}}{\partial\eta}\right]_{\eta=0}\\\
=&\ \dfrac{|v|d}{\pi}\exp\\{-v_{2}d\\}\exp\\{-v_{1}(x_{1}-\xi)\\}\\\ &\
\cdot\dfrac{K_{1}\left(|\mathbf{v}|\sqrt{(x_{1}-\xi)^{2}+d^{2}}\right)}{\sqrt{(x_{1}-\xi)^{2}+d^{2}}}.\end{split}$
(40)

To tackle with the general case $\sigma^{2}\neq 1$, we can simply replace
$v_{i}$ with $\frac{v_{i}}{\sigma^{2}}$, yielding

$\displaystyle\begin{split}f_{Y\mid X}(\mathbf{y}|\mathbf{x})=&\
\dfrac{|\mathbf{v}|d}{\sigma^{2}\pi}\exp\left\\{\dfrac{-v_{2}d}{\sigma^{2}}\right\\}\exp\left\\{\dfrac{-v_{1}(x_{1}-\xi)}{\sigma^{2}}\right\\}\\\
&\cdot\dfrac{K_{1}\left(\dfrac{|\mathbf{v}|}{\sigma^{2}}\sqrt{(x_{1}-\xi)^{2}+d^{2}}\right)}{\sqrt{(x_{1}-\xi)^{2}+d^{2}}},\end{split}$
(41)

where $\mathbf{x}=(x_{1},d)$, $\mathbf{y}=(\xi,0)$, and $d$ is the distance
between Tx and Rx. Eq. (41) is the desired conditional density function,
allowing arbitrary drift directions, for 2D FAP problem.

In order to make a comparison to the restrictive version in [16, Eq. (19)], we
can set $v_{1}=0$ in Eq. (41), yielding

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

where we use the relation $\sigma^{2}=2D$.

##  Appendix B The Old Method: Separation of Variables and Method of Image

For the sake of completeness, we review in a self-contained fashion the
original approach to finding FAP density function as proposed in [16].
Throughout this appendix section, we set the $x$-direction to be the direction
of transmission, i.e. the direction pointing from Tx to Rx.

###  B-A Finding the Green’s Function for Free Space

For a 2-dimensional (2D) diffusion channel with a constant drift velocity
$\mathbf{v}=(v,0)$ pointing from the transmitter to the receiver along the
$x$-axis, equation (4) can be written explicitly as

$\frac{\partial c(x,y,t)}{\partial t}+v\frac{\partial c(x,y,t)}{\partial
x}=D\left(\frac{\partial^{2}c(x,y,t)}{\partial
x^{2}}+\frac{\partial^{2}c(x,y,t)}{\partial y^{2}}\right)$ (43)

where $c(x,y,0)=\delta\left(x-x_{0}\right)\delta\left(y-y_{0}\right)$ is the
2D Dirac delta function. The boundary condition is set to be:

$c(x,y,t)=0\quad\text{for}\quad x,y\in\partial\Omega,$ (44)

representing an absorbing receiver.

Using separation of variables, the concentration field $c(x,y,t)$ can be
written as

$c(x,y,t)=c_{1}(x,t)c_{2}(y,t).$ (45)

By plugging Eq. (45) back into Eq. (43), we get two separated equations

$\displaystyle\frac{\partial c_{1}}{\partial t}+v\frac{\partial
c_{1}}{\partial x}-D\frac{\partial^{2}c_{1}}{\partial x^{2}}$
$\displaystyle=0,$ (46) $\displaystyle\frac{\partial c_{2}}{\partial
t}-D\frac{\partial^{2}c_{2}}{\partial y^{2}}$ $\displaystyle=0.$ (47)

The fundamental solutions (aka free space Green’s function) of Eq. (46) and
Eq. (47) are

$c_{1}(x,t)=\frac{1}{\sqrt{4\pi
Dt}}\exp\left(-\frac{\left(x-x_{0}-vt\right)^{2}}{4Dt}\right)$ (48)

and

$c_{2}(y,t)=\frac{1}{\sqrt{4\pi
Dt}}\exp\left(-\frac{\left(y-y_{0}\right)^{2}}{4Dt}\right)$ (49)

respectively. Hence, the free space Green’s function of (43) can be written as
the product of Eq. (48) and Eq. (49):

$\displaystyle\begin{split}&G\left(x,y,x_{0},y_{0},t\right)\\\ =\
&\frac{1}{4\pi
Dt}\exp\left(\frac{-\left(x-x_{0}-vt\right)^{2}-\left(y-y_{0}\right)^{2}}{4Dt}\right).\end{split}$
(50)

###  B-B Finding the Green’s Function for Absorbing Boundary

For Eq. (50), we have not yet taken the boundary conditions into account.
Since the receiver is assume to be perfectly absorbing, we require the
concentration to be zero on the boundary, as Eq. (44) shows.

There is one well-known method, originated in classical electrodynamics, that
works well in problems of finding static distributions for flat boundaries,
called the method of image, see [30]. Since we only care about the
concentration function inside the domain, we can effectively put an “mirror
image” of negative mass at the reflection point of the original releasing
point. Using this method, the domain of the original unknown function is not
extended, but the function is made to satisfy given boundary conditions, such
as absorbing boundary in our case.

To use the method of image, we write

$\displaystyle\begin{split}&G_{\text{abs}}(x,y,x_{0},y_{0},t)\\\ =\
&G(x,y,x_{0},y_{0},t)-a(x_{0})G(x,y,-x_{0},y_{0},t)\end{split}$ (51)

and then solve for $a(x_{0})$. Note that the subscript “abs” stands for
absorbing. The absorbing boundary condition requires that

$a(x_{0})=\dfrac{G(0,y,x_{0},y_{0},t)}{G(0,y,-x_{0},y_{0},t)}=\exp\left(-\dfrac{x_{0}v}{D}\right).$
(52)

We can now write down the Green’s function for absorbing boundary as

$\begin{aligned} &G_{\text{abs}}\left(x,y,x_{0},y_{0},t\right)\\\ =\
&\frac{1}{4\pi
Dt}\left[\exp\left(-\frac{\left(x-x_{0}-vt\right)^{2}+\left(y-y_{0}\right)^{2}}{4Dt}\right)\right]\\\
&\left.-\exp\left(-\frac{x_{0}v}{D}-\frac{\left(x+x_{0}-vt\right)^{2}+\left(y-y_{0}\right)^{2}}{4Dt}\right)\right]\end{aligned}.$
(53)

###  B-C From Green’s function to FAP Density

Consider a quantity called diffusive flux $J(0,y,t)$ at $x=0$ plane, which
possesses a physical unit of particle density per unit length per unit time.
By Fick’s law, we know that

$J(0,y,t)=-D\dfrac{\partial G_{\text{abs}}(x,y,x_{0},y_{0},t)}{\partial x}.$
(54)

Using the absorbing Green’s function formula as shown in Eq. (53), we can
calculate the flux explicitly, yielding

$J(0,y,t)=-\frac{x_{0}}{4\pi
Dt^{2}}\exp\left(-\frac{\left(x_{0}+vt\right)^{2}+\left(y-y_{0}\right)^{2}}{4Dt}\right).$
(55)

Note that the original formula [16, Eq. (14)] has some calculation error, so
we correct the formula herein.

Our goal is to find the FAP density $f(\mathbf{y}|\mathbf{x})$. Let
$\mathbf{x}=\left(x_{0},y_{0}\right)$ denote the point of release of the
information molecule and $\mathbf{y}=(x,y)$ its arrival position at the
receiver boundary. The conditional PDF of $Y$ at a receiver located at $x=0$
is given as

$f_{Y\mid X}(0,y\mid
x_{0},y_{0})=\left.\int_{0}^{\infty}J(0,y,t)\mathrm{d}t\right..$ (56)

Plugging the diffusion flux we just obtained in Eq. (55) into Eq. (56), we
have

$\displaystyle\begin{split}&f_{Y\mid X}(0,y\mid x_{0},y_{0})\\\ =\
&-\int_{0}^{\infty}\frac{x_{0}}{4\pi
Dt^{2}}\exp\left(-\frac{\left(x_{0}+vt\right)^{2}+\left(y-y_{0}\right)^{2}}{4Dt}\right)\mathrm{d}t.\end{split}$
(57)

In order to calculate this integration explicitly, we rearrange the terms to
get

$\displaystyle\begin{split}&f_{Y\mid X}(0,y\mid x_{0},y_{0})\\\ =\ &\
\frac{-x_{0}}{4\pi D}\exp\left(-\frac{x_{0}v}{2D}\right)\\\
&\times\int_{0}^{\infty}\frac{1}{t^{2}}\exp\left(-\frac{x_{0}^{2}+\left(y-y_{0}\right)^{2}}{4Dt}-\frac{v^{2}t}{4D}\right)\mathrm{d}t.\end{split}$
(58)

Letting $K_{s}(z)$ be the modified Bessel function of the second kind of order
$s$, for positive constants $a$ and $b$, we have the relation [31]:

$\int_{0}t^{-s-1}\exp\left(-\left(\frac{a}{t}+bt\right)\right)\mathrm{d}t=2\left(\frac{b}{a}\right)^{s/2}K_{s}(2\sqrt{ab}),$
(59)

where $s$ and $z$ are real and complex numbers respectively. Combining Eq.
(58) and Eq. (59) with $x_{0}=d$, we finally get

$\displaystyle\begin{split}&f(\mathbf{y}|\mathbf{x})\\\ =\ &f_{Y\mid
X}(0,y\mid d,y_{0})\\\ =\ &\frac{|v|d}{2\pi
D}\exp\left(\frac{-dv}{2D}\right)\dfrac{1}{\sqrt{\left(d^{2}+\left(y-y_{0}\right)^{2}\right)}}\\\
&\cdot
K_{1}\left(\frac{|v|\sqrt{\left(d^{2}+\left(y-y_{0}\right)^{2}\right)}}{2D}\right).\end{split}$
(60)

It is worth noting that the conditional FAP density as shown in Eq. (60)
depends only on $y-y_{0}$, but not on time variable $t$. We can conclude that
the FAP channel is an additive time-invariant channel.