# Sorting Lists with Equal Keys Using Mergesort in Linear Time

Albert Tedja

###### Abstract

This article introduces a new optimization method to improve mergesort’s
runtime complexity, when sorting sequences that have equal keys to
$O(nlog_{2}k)$, where $k$ is the number of distinct keys in the sequence. When
$k$ is constant, it is evident that mergesort is capable of achieving linear
time by utilizing linked lists as its underlying data structure. Mergesort
linked list implementations can be optimized by introducing a new mechanism to
group elements with equal keys together, thus allowing merge algorithm to
achieve linear time.

##  1 Introduction

The incremental nature of linked list allows only linear search and traversal,
thus proving less useful compared to random-access arrays where they can
benefit more from the locality of reference and $O(log_{2}{n})$ search
algorithms such as binary or exponential search. Nevertheless, linked lists
$O(1)$ insertion and deletion can sometimes prove to be useful in some
situations, in particular, mergesort.

Mergesort algorithm is usually implemented in two different ways. The most
notable one is the top-down approach, in which the input sequence is split in
half and each half is sorted independently and recursively, and finally merged
to form the final sorted sequence. Top-down mergesort requires the length of
the input sequence to be known in advance to optimally divide the sequence
equally in half. The second method is the bottom-up approach, where it
iteratively creates sublists of length one and inserts them into a stack
implementation. Sublists in stack are merged depending on certain heuristics.
This removes the requirement of knowing the length of the list in advance.
Panny and Prodinger [PP95] provide a detailed analysis of top-down vs bottom-
up mergesort and conclude that top-down mergesort provides a better
performance although only by a slim margin. Many linked list implementations
keep track of an internal number of items in the list that it takes $O(1)$ to
retrieve its length. However, due to the iterative nature of linked lists, the
bottom-up approach of mergesort is chosen as the baseline of the optimization
described here.

##  2 Existing Optimizations and Implementations

Natural mergesort is a variant of the bottom-up approach with optimizations to
detect existing monotonic sublists in the original input sequence, also
colloquially known as runs. A run is not split but rather kept and merged with
its adjacent sublists. Brady [Bra05] introduces Runsort, that takes a step
further by also identifying decreasing monotonic sublists and reverses their
order before merging them.

Part of implementing an optimal natural mergesort algorithm is to determine
the heuristic of when to merge sublists to maximize optimality in both space
and time because each run has varying lengths. Shivers [Shi02] proposes a
solution by tracking the length of sublists and prioritizing merging shorter
sublists. Timsort [Pet02] approaches the problem in a similar fashion, but
adds an additional layer by keeping a minimum threshold value that if a
sublist’s length is below a certain threshold, insertion sort is used before
it is pushed into the stack. Knuth [Knu98, p.163] noted such optimization is
worthwhile to reduce administrative overhead of short merges without affecting
the asymptotic running time. These hybrid methods of using multiple sorting
algorithms are quite common, and can be found in industrial-grade sorting
libraries such as C++ STL’s std::sort or some implementations of the GNU C’s
qsort based on the works of Bentley and Douglas [BM93].

Wegner’s trisort [Weg82] attempts to sort linked lists with equal keys using
the quicksort method with the average cost in the order of
$O((m+n)log_{2}(\frac{n}{m}))$, where $n$ is the number of items occuring $m$
times. Similar to quicksort, however, it has the possibility to suffer from a
bad pivot. It is especially even more likely to happen in linked lists as only
either the $head$ or $tail$ can be chosen as the pivot.

##  3 Baseline Algorithm

Algorithm 1 shown in Appendix A describes a straight-forward bottom-up
mergesort algorithm implemented using linked lists without any optimization,
such as detecting existing runs or comparing lengths of sublists. Unoptimized
mergesort is chosen to demonstrate that the optimization described here is
isolated and performs without the compounding effect of other optimizations.
The algorithm hereafter shall serve as the baseline of the optimization
presented in this article.

To summarize, the baseline algorithm treats each element as a sublist of
length 1, and inserts sublists into a stack implementation. The algorithm does
not maintain additional administrative computational complexities such as
minimum runs in Timsort [Pet02] nor identifies shorter sublists in Shivers’
[Shi02] implementation, but rather resorts to a simpler approach similar to
O’Keefe’s samsort [O’K82]

The process in which the stack determines if it needs to merge sublists is
simply decided by the binary pattern of the stack counter. Let $c$ be the
stack counter, and it counts the number of sublists inserted into the stack.
The binary pattern of $c$, starting from the least significant digit,
indicates the number of merges required before inserting the $(c+1)$-th
sublist. For example, if $c=7$, then $c=0111_{2}$, then 3 merges are performed
across three levels 1, 2, and 3, when inserting the 8th sublist. Merging is
terminated once a 0 is encountered. This solution is therefore opportunistic
and performs merges as soon as possible, and allows minimal constant space
requirement for up to $2^{\alpha}-1$ number of sublists for a stack with
maximum depth of $\alpha$.

##  4 Merging Using Linked Lists

Using arrays as the underlying structure and buffer, mergesort’s merging
operation is known to be proportional to $m+n$ [Knu98, p.159], where $m$ and
$n$ are the lengths of the sublists to be merged. The disadvantage of this
approach is the $O(n)$ space requirement for the buffer. Attempts to reduce
the $O(n)$ space requirement, such as by Katajainen et al. [KPT96] or Huang et
al. [HL88], increase the merging time complexity to $O(nlog_{2}n)$, which is
not favorable.

Linked lists, however, naturally allow merging without an additional $O(n)$
space requirement. This comes at a cost that linked list implementations do
not normally benefit from the principle of locality, and searching for
insertion points must be done linearly. Carlsson et al. [CLP93] introduce
sublinear merging using a combination of arrays and linked lists by grouping
blocks from input sequences $X$ and $Y$, where each block is computed by using
binary and exponential searches alternatingly between $X$ and $Y$. The
computed ranges of blocks are organized inside a linked list $Z$, such that
traversing $Z$ outputs the merged sorted sequence, thereby preventing the need
to copy or swap elements.

It is possible to merge two sublists with cost function proportional to $k$,
where $k$ is the total number of distinct keys in the sublists, by using
linked list as the underlying data structure. Instead of identifying monotonic
ranges in input sequences, it identifies elements with equal keys. This
reduces the overall mergesort time complexity to linear when sorting such
particular datasets.

##  5 Optimizing with Hop Pointers

In addition to the ubiquitous linked list’s next pointer, and its contrast
previous pointer in doubly linked lists, this article introduces a new concept
called the hop pointer. The hop pointer is remarkably simple as by default it
simply points to itself, its host node. Given a sequence, and there exist
elements with equal keys adjacent to each other, we define these contiguous
set of elements with equal keys as a segment. When such a segment exists in a
sequence, then the hop pointer of the first node in the segment shall point to
the last node in the segment.

Figure 1: Shows two lists $L_{1}$ and $L_{2}$ with hop pointers. $L_{3}=$
MERGE$(L_{1},L_{2})$, from Algorithm 2. The second $3$ in $L_{3}$ is the first
$3$ from $L_{2}$. Its hop pointer can be optionally set to itself after the
merge, but not required for the optimization discussed here.

As observed in Figure 1, the hop pointer of the first element in a segment
identifies the last node in the segment. The purpose of hop pointers is to
compress the list down to its distinct keys, allowing an iterative algorithm
to skip forward to the end of the segment, essentially treating an entire
segment as one step.

######  Lemma 5.1.

Let $p$ and $q$ be the number of distinct keys in sorted lists $L_{1}$ and
$L_{2}$, respectively. The cost of merging $L_{1}$ and $L_{2}$ using hop
pointers is $O(p+q)$.

###### Proof.

Hop pointers group elements with equal keys as one element. Consequently, the
cost of iterating $L_{1}$ and $L_{2}$ is $O(p)$ and $O(q)$, respectively. An
efficient merge algorithm requires iterating $L_{1}$ and $L_{2}$ each exactly
once [Knu98, p.159]. The cost of merging is therefore $O(p+q)$. ∎

Merging two sequences no longer requires $O(|L_{1}|+|L_{2}|)$ but $O(p+q)$.
Since $L_{3}=$ MERGE$(L_{1},L_{2})$, we know that $k$, the number of distinct
keys in $L_{3}$, cannot be greater than $p+q$, therefore $k\leq p+q$.

This optimization extends to subsequent merges thereafter, such as merging
$L_{3}$ with $L_{4}$, $L_{5}$, …, $L_{n}$. The complexity of merging is no
longer dependent on the length of the lists but by the number of their
distinct keys, thus allowing merge algorithm to perform more efficiently.

Implementing hop pointers on top of the baseline algorithm is quite straight-
forward, and does not involve complex restructuring of the code or logic
rearrangement. Algorithm 2 in Appendix A shows the optimized merge using hop
pointers. A single hop pointer in every node is mostly sufficient for singly
linked lists. An additional hop-back pointer can be optionally added if
identifying a segment’s head node is necessary, for example, when appending
one linked list to another, or in doubly linked lists where reverse traversals
are permitted.

##  6 Linear Time Complexity

When $n=k$, every element is distinct, the conventional recurrence relation
expression for mergesort $T(n)=n+nlog_{2}n$ is assumed. When $n>k$, i.e. there
exist elements with equal keys, we define the recurrence relation to be
$T(n)=2T(\frac{n}{2})+k$, because the cost of the merge function can be no
more than $O(k)$. Therefore, the recurrence relation notation is:

$T(n)=\begin{cases}n+nlog_{2}n&n=k\\\
2T\left(\frac{n}{2}\right)+k&n>k\end{cases}$

The base case $T(n)=1$ for $n=1$ is omitted because it is already solved and
implied in the $n=k$ relation, which becomes the new base case.

Solving the recurrence relation for $n>k$ by expanding
$T\left(\frac{n}{2}\right)$ yields:

$\begin{split}T(n)&=2\left[2T\left(\frac{n}{4}\right)+k\right]+k\\\
&=4T\left(\frac{n}{4}\right)+3k\\\
&=4\left[2T\left(\frac{n}{8}\right)+k\right]+3k\\\
&=8T\left(\frac{n}{8}\right)+7k\\\ &=\dots\\\
&=pT\left(\frac{n}{p}\right)+(p-1)k\end{split}$

When $\frac{n}{p}=k$, $p=\frac{n}{k}$:

$\begin{split}T(n)&=\frac{n}{k}T(k)+(\frac{n}{k}-1)k\\\
&=\frac{n}{k}(k+klog_{2}k)+n-k\\\ &=n+nlog_{2}k+n-k\\\ &=2n+nlog_{2}k-k\\\
\end{split}$

When $k=n$, i.e. every element is distinct, the relation above reduces back to
the conventional mergesort recurrence relation $T(n)=n+nlog_{2}n$. However, if
$k$ remains constant, the mergesort runtime complexity becomes $O(nlog_{2}k)$
or simply $O(n)$.

##  7 Empirical Result

The following tests are conducted to measure the number of comparisons between
the baseline implementation and the hoptimized111In the spirit of
optimization, given a string ”hop-optimized”, omitting the redundancy of the
dash delimiter from the English language, finding the longest repeated non-
overlapping substring yields ”op”. The string can then be reduced to
”hoptimized”. implementation. There are 3 types of sample sets tested:

  * •

Shuffled \- Shuffled sequences where every element is distinct.

  * •

Sawtooth \- Repeated sawtooth pattern with values up to $k$.

  * •

K-Distinct \- Shuffled version of Sawtooth.

Each dataset is executed over 100 random permutations with input length $n$
exponentially increasing from $2^{7}$ to $2^{22}$. The total numbers of
comparison are recorded and presented in Table 1.

$n$ | Shuffled | Sawtooth | K-Distinct
---|---|---|---
|  Baseline |  Hoptimized |  Baseline |  Hoptimized |  Baseline |  Hoptimized
$2^{7}$ |  735 |  735 |  448 |  448 |  736 |  736
$2^{8}$ |  1 725 |  1 725 |  1 024 |  1 024 |  1 727 |  1 727
$2^{9}$ |  3 961 |  3 961 |  2 304 |  2 304 |  3 964 |  3 964
$2^{10}$ |  8 946 |  8 946 |  5 120 |  5 120 |  8 948 |  8 948
$2^{11}$ |  19 942 |  19 942 |  12 287 |  11 275 |  19 936 |  17 918
$2^{12}$ |  43 974 |  43 974 |  28 668 |  23 556 |  43 969 |  36 068
$2^{13}$ |  96 131 |  96 131 |  65 524 |  48 139 |  96 135 |  72 415
$2^{14}$ |  208 673 |  208 673 |  147 424 |  97 306 |  208 640 |  145 153
$2^{15}$ |  450 094 |  450 094 |  327 600 |  195 641 |  450 028 |  290 687
$2^{16}$ |  720 704 |  720 704 |  720 704 |  392 312 |  965 586 |  581 725
$2^{17}$ |  2 062 483 |  2 062 483 |  1 572 416 |  785 655 |  2 062 221 |  1 163 812
$2^{18}$ |  4 387 094 |  4 387 094 |  3 406 848 |  1 572 342 |  4 386 381 |  2 328 143
$2^{19}$ |  9 298 502 |  9 298 502 |  7 337 728 |  3 145 717 |  9 296 875 |  4 656 799
$2^{20}$ |  19 645 532 |  19 645 532 |  15 723 520 |  6 292 468 |  19 641 712 |  9 314 026
$2^{21}$ |  41 388 301 |  41 388 301 |  33 543 168 |  12 585 971 |  41 379 540 |  18 628 520
$2^{22}$ |  86 971 029 |  86 971 029 |  71 278 592 |  25 172 978 |  86 950 900 |  37 257 365
Table 1: Total number of comparisons in Bottom-Up Mergesort: Baseline vs
Hoptimized. $k=1024$

Shuffled sample sets do not show any improvement as shown in Table 1. This is
expected as hop pointers do not offer any benefit when sorting sequences where
each element is distinct. In Sawtooth and K-Distinct sample sets, there are
noticeable improvements in that the number of comparisons slowly turns linear
once $n$ grows larger than $k$. Dividing the total comparisons by $n$ yields a
more interesting observation shown in Table 2. The average number of
comparisons per element slowly plateaus even as $n$ grows exponentially.

$n$ | Sawtooth | K-Distinct
---|---|---
|  Baseline |  Hoptimized |  Baseline |  Hoptimized
$2^{7}$ |  3.50000 |  3.50000 |  5.75000 |  5.75000
$2^{8}$ |  4.00000 |  4.00000 |  6.74609 |  6.74609
$2^{9}$ |  4.50000 |  4.50000 |  7.74219 |  7.74219
$2^{10}$ |  5.00000 |  5.00000 |  8.73828 |  8.73828
$2^{11}$ |  5.99951 |  5.50049 |  9.73438 |  8.74902
$2^{12}$ |  6.99902 |  5.75098 |  10.73462 |  8.80566
$2^{13}$ |  7.99854 |  5.87634 |  11.73523 |  8.83972
$2^{14}$ |  8.99805 |  5.93909 |  12.73438 |  8.85944
$2^{15}$ |  9.99756 |  5.97049 |  13.73376 |  8.87106
$2^{16}$ |  10.99707 |  5.98621 |  14.73367 |  8.87918
$2^{17}$ |  11.99658 |  5.99407 |  15.73350 |  8.87918
$2^{18}$ |  12.99609 |  5.99801 |  16.73272 |  8.88116
$2^{19}$ |  13.99561 |  5.99998 |  17.73238 |  8.88214
$2^{20}$ |  14.99512 |  6.00097 |  18.73180 |  8.88255
$2^{21}$ |  15.99463 |  6.00146 |  19.73130 |  8.88277
$2^{22}$ |  16.99414 |  6.00171 |  20.73031 |  8.88285
Table 2: Average number of comparisons per element. Baseline vs Hoptimized.
$k=1024$

##  8 Conclusion

Even though the implementation described here does not take advantage of the
locality of reference or other optimizations, it is demonstrated that linear
time mergesort is achievable. By simply adding a $hop$ pointer in every node,
mergesort can achieve linear time when sorting permutations of sequences with
equal keys.

Many real-world datasets have equal keys, such as sorting cars by their makes
or models, merchandise by their brands, or citizens by cities or countries.
The stability of mergesort allows this optimization to reach an even greater
purpose because it keeps the original ordering intact.

Hop pointers can be implemented in tandem with other sorting optimizations,
such as detecting existing runs and better merging strategies, thus creating
an even faster mergesort. Hop pointers, additionally, can be used to optimize
other sorting algorithms such as insertion sort, or quickly determine the
number of unique keys after sorting. They can also increase the speed of
binary and exponential searches by further reducing the number of elements to
be searched.

## References

  * [BM93] Jon L. Bentley and M. Douglas McIlroy.  Engineering a sort function.  Software: Practice and Experience, 23(11):1249–1265, 1993.
  * [Bra05] Mike Brady.  ’runsort’ - an adaptive mergesort for prolog.  Technical Report TCD-CS-2005-34, TCD Computer Science Department, May 2005\.
  * [CLP93] Svante Carlsson, Christos Levcopoulos, and Ola Petersson.  Sublinear merging and natural mergesort.  Algorithmica, 9(6):629–648, Jun 1993.
  * [HL88] Bing-Chao Huang and Michael A. Langston.  Practical in-place merging.  Commun. ACM, 31(3):348–352, Mar 1988.
  * [Knu98] Donald E. Knuth.  Art of Computer Programming, Volume 3: Sorting and Searching.  Addison Wesley, 1998.
  * [KPT96] Jyrki Katajainen, Tomi Pasanen, and Jukka Teuhola.  Practical in-place mergesort.  Nordic Journal of Computing, 3(1):27–40, March 1996.
  * [O’K82] Richard O’Keefe.  A smooth applicative merge sort.  Dept. of Artificial Intelligence, University of Edinburgh, 1982.
  * [Pet02] Tim Peters.  https://svn.python.org/projects/python/trunk/Objects/listsort.txt, 2002\.  Online; accessed August 2019.
  * [PP95] W. Panny and H. Prodinger.  Bottom-up mergesort — a detailed analysis.  Algorithmica, 14(4):340–354, Oct 1995.
  * [Shi02] Olin Shivers.  A simple and efficient natural merge sort.  Technical report, Georgia Institute of Technology, Mar 2002.  Online; accessed on August 2019.
  * [Weg82] Lutz M. Wegner.  Sorting a linked list with equal keys.  Information Processing Letters, 15:205–208, 1982.

##  Appendix A Algorithms

Algorithm 1 Baseline Bottom-Up Mergesort

1:procedure merge($a,b$) $\triangleright$ Merge two lists $a$ and $b$

2: if $a=\oslash$ then return $b$

3: else if $b=\oslash$ then return $a$

4: $head\leftarrow\oslash$

5: if $a.value\leq b.value$ then

6: $head\leftarrow a$

7: $a\leftarrow a.next$

8: else

9: $head\leftarrow b$

10: $b\leftarrow b.next$

11:

12: $p\leftarrow head$

13: while $a\neq\oslash\land b\neq\oslash$ do

14: if $a.value\leq b.value$ then

15: $p.next\leftarrow a$

16: $p\leftarrow a$

17: $a\leftarrow a.next$

18: else

19: $p.next\leftarrow b$

20: $p\leftarrow b$

21: $b\leftarrow b.next$

22:

23: if $a=\oslash$ then $p.next\leftarrow b$

24: else if $b=\oslash$ then $p.next\leftarrow a$

25: return $head$

26:procedure mergesort($node$)

27: if $node=\oslash\lor node.next=\oslash$ then return $node$

28: $S\leftarrow Stack$

29: $c\leftarrow 0$ $\triangleright$ Stack counter

30: while $node\neq\oslash$ do

31: $next\leftarrow node.next$

32: $node.next\leftarrow\oslash$

33: $b\leftarrow c$

34: while $b\bmod 2=1$ do

35: $a\leftarrow\textsc{pop}(S)$

36: $node\leftarrow\textsc{merge}(a,node)$

37: $b\leftarrow\lfloor b/2\rfloor$

38:

39: push(S, node)

40: $c\leftarrow c+1$

41: $node\leftarrow next$

42:

43: if $node=\oslash$ then $node\leftarrow\textsc{pop}(S)$

44: while $S$ not empty do

45: $a\leftarrow\textsc{pop}(S)$

46: $node\leftarrow\textsc{merge}(a,node)$

47:

48: return $node$

Algorithm 2 Optimized Bottom-Up Mergesort with Hop Pointers

1:procedure merge($a,b$)

2: if $a=\oslash$ then return $b$

3: else if $b=\oslash$ then return $a$

4: $head\leftarrow\oslash$

5: if $a.value\leq b.value$ then

6: $head\leftarrow a$

7: $a\leftarrow a.hop.next$

8: else

9: $head\leftarrow b$

10: $b\leftarrow b.hop.next$

11:

12: $p\leftarrow head.hop$

13: while $a\neq\oslash\land b\neq\oslash$ do

14: if $a.value<b.value$ then

15: $p.next\leftarrow a$

16: $p\leftarrow a.hop$

17: $a\leftarrow a.hop.next$

18: else if $a.value>b.value$ then

19: $p.next\leftarrow b$

20: $p\leftarrow b.hop$

21: $b\leftarrow b.hop.next$

22: else

23: $p.next\leftarrow a$

24: $p\leftarrow b.hop$

25: $tmp\leftarrow a.hop.next$

26: $a.hop.next\leftarrow b$

27: $a.hop\leftarrow b.hop$

28: $a\leftarrow tmp$

29: $b\leftarrow b.hop.next$

30:

31: if $a=\oslash$ then $p.next\leftarrow b$

32: else if $b=\oslash$ then $p.next\leftarrow a$

33: return $head$


#  Capacity Bounds for Vertically-Drifted First Arrival Position Channels
under a Covariance Constraint

Yun-Feng Lo  Georgia Institute of Technology
Atlanta, Georgia
<EMAIL_ADDRESS>Yen-Chi Lee  Hon Hai (Foxconn) Research Institute
Taipei, Taiwan
<EMAIL_ADDRESS>Min-Hsiu Hsieh  Hon Hai (Foxconn) Research Institute
Taipei, Taiwan
<EMAIL_ADDRESS>
###### Abstract

In this paper, we delve into the capacity problem of additive vertically-
drifted first arrival position noise channel, which models a communication
system where the position of molecules is harnessed to convey information.
Drawing inspiration from the principles governing vector Gaussian interference
channels, we examine this capacity problem within the context of a covariance
constraint on input distributions. We offer analytical upper and lower bounds
on this capacity for a three-dimensional spatial setting. This is achieved
through a meticulous analysis of the characteristic function coupled with an
investigation into the stability properties. The results of this study
contribute to the ongoing effort to understand the fundamental limits of
molecular communication systems.

###### Index Terms:

molecular communication (MC), vertical drift (VD), first arrival position
(FAP), channel capacity, covariance constraint, characteristic function (CF).

##  I Introduction

Molecular communication (MC) is a communication paradigm that employs message
molecules (MM) as information carriers [1, 2]. Due to its biocompatibility and
feasibility at the nanoscale, MC is considered a promising approach for nano-
networks [3, 4]. In MC systems, MMs act as the information carriers, and a
propagation mechanism is necessary to transport them to the receiver. This
mechanism can be diffusion-based [5], flow-based [6], or an engineered
transport system such as molecular motors [7]. Among these mechanisms,
_diffusion-based MC_ , often in combination with advection and chemical
reaction networks (CRNs), has been the prevalent approach in the literature
[4, Table 4].

The channel characteristics of a diffusive MC system are dependent not only on
the physical properties of the propagation medium but also on the reception
mechanism. For our study, we adopt the fully-absorbing receiver, which is the
common type of active reception [5]. First arrival models for absorbing
receivers can be roughly categorized into three types: time modulation,
position modulation, and joint position-time modulation. In the following, we
briefly review the progress made in channel capacity research for first
arrival time (FAT) and first arrival position (FAP) channels.

FAT channels can be represented as time-invariant additive channels for which
$T_{\text{out}}=T_{\text{in}}+T_{\text{n}},$ where the arrival time
$T_{\text{out}}$ is given by the releasing time $T_{\text{in}}$ plus a random
time delay $T_{\text{n}}$ caused by the propagation mechanisms [8]. The
additive inverse Gaussian noise (AIGN) channel is FAT channel with an inverse
Gaussian distribution for the random delay $T_{\text{n}}$. The capacity of the
AIGN channel has been studied in the literature, with both upper and lower
bounds on the capacity developed [8, 9, 10]. Moreover, [11] explores the
capacity-achieving input of the AIGN channel, considering both average and
peak release time constraints.

For FAP channels in D-dimensional (dim) spaces, the one-shot channel model can
also be written in additive vector form [12, 13] as

$\boldsymbol{X}_{\text{out}}=\boldsymbol{X}_{\text{in}}+\boldsymbol{X}_{\text{n}},$
(1)

where $\boldsymbol{X}_{\text{in}}$ is the emission position,
$\boldsymbol{X}_{\text{out}}$ is the arrival position, and
$\boldsymbol{X}_{\text{n}}$ is the random position bias incurred by the
propagation mechanisms. Note that $\boldsymbol{X}_{\text{in}}$,
$\boldsymbol{X}_{\text{out}}$, and $\boldsymbol{X}_{\text{n}}$ are all
Euclidean vectors in $\mathbb{R}^{d}$. Throughout this paper, we define
$d:=\text{D}-1$.

Although the density function of $\boldsymbol{X}_{\text{n}}$ has been derived
in previous works [14, 12, 13], an exact characterization of the capacity of
FAP channels remains elusive, except for a special case discussed (via
numerical simulation) in [12, Section IV-B], where equally spaced $M$-ary
modulation is assumed for a 2D FAP channel, and the transition probabilities
are discretized. This paper aims to fill this gap by providing analytic
capacity bounds for 3D FAP channels without assuming any specific modulation.

We study the capacity of additive Vertically-Drifted First Arrival Position
(VDFAP) noise channels

$\displaystyle\begin{split}C=\sup_{f(\boldsymbol{X}_{\text{in}}):\leavevmode\nobreak\
\mathbb{E}\left[\boldsymbol{X}_{\text{in}}^{\text{
}}\boldsymbol{X}_{\text{in}}^{\mathsf{\scriptstyle{T}}}\right]\preceq\Sigma}I(\boldsymbol{X}_{\text{in}};\boldsymbol{X}_{\text{out}}),\end{split}$
(2)

where the objective function is the mutual information $I$ between
$\boldsymbol{X}_{\text{in}}$ and $\boldsymbol{X}_{\text{out}}$, and the
supremum is taken over all input distributions $f(\boldsymbol{X}_{\text{in}})$
satisfying a _covariance matrix constraint_ [15] expressed as
$\mathbb{E}\left[\boldsymbol{X}_{\text{in}}^{\text{
}}\boldsymbol{X}_{\text{in}}^{\mathsf{\scriptstyle{T}}}\right]\preceq\Sigma$,
where $\Sigma\succ 0$.111We define the notation $A\succeq B$ as $A-B$ being
positive semi-definite, and $A\succ B$ as $A-B$ being positive definite.
Similar definitions apply for $\preceq$ and $\prec$. To obtain analytic bounds
for this capacity, we first derive the characteristic function (CF) of the
VDFAP noise distribution. Using the CF, we not only derive formulas for the
first two moments but also prove a weak stability property of VDFAP
distributions. We then apply the moments and weak stability (see Section
III-B) to derive lower and upper bounds for the capacity.

The structure of this paper is as follows. In Section II, we present the
system model used in our analysis. In Section III, we examine the
characteristic function of the VDFAP distribution. We then use the derived
moments and weak stability property to provide lower and upper bounds for the
capacity of VDFAP channels in Section IV. Finally, we summarize our results
and provide concluding remarks in Section V.

##  II System Model

Figure 1:  This figure depicts a 3D MC system where both the transmitter and
receiver have plane shapes. The communication takes place in a fluid medium
with a vertical drift $v_{3}<0$, and the transmission distance is denoted by
$\lambda>0$. The orthonormal basis vectors in 3D space are denoted by
$\\{\mathbf{\hat{e}}_{i}\\}_{i=1}^{3}$.

In this paper, we investigate a diffusion-based molecular communication
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 used.
We consider a scenario where there is a constant drift [6] in the ambient
space. The diffusion effects of MMs are captured by the diffusion coefficient
$D=\sigma^{2}/2$ [4, 13].

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

In this study, we assume an ideal MC system that satisfies the following
assumptions:

  * •

The transmitter has perfect control over the emission position of the
molecules.

  * •

The transmitter plane is transparent, allowing MMs to move through it without
experiencing any force after they are released.

  * •

The receiver can perfectly measure the first arrival positions of the
molecules.

  * •

Upon first arrival at the receiver plane, molecules are captured and removed
from the system.

  * •

The movement of every molecule is independent.

The positional information is $d$-dim, where $d=\text{D}-1$ for brevity and
consistency throughout the paper. The first arrival position $\boldsymbol{Y}$
of a molecule released at position $\boldsymbol{X}$ can be expressed as
$\boldsymbol{Y}=\boldsymbol{X}+\boldsymbol{N}$, where $\boldsymbol{N}$ denotes
the deviation of the first arrival position. The probability density function
of $\boldsymbol{N}$ is known as the _FAP density_ [13, 12].

In order to express the FAP density, we denote the D-dim drift vector
$\mathbf{v}=[v_{1},\ldots,v_{\text{D}}]^{\mathsf{\scriptstyle{T}}}$ as
$\mathbf{v}=[\boldsymbol{v}_{\text{par}}^{\mathsf{\scriptstyle{T}}},v_{\text{D}}]^{\mathsf{\scriptstyle{T}}}$,
where
$\boldsymbol{v}_{\text{par}}:=[v_{1},\ldots,v_{d}]^{\mathsf{\scriptstyle{T}}}$
contains the drift components parallel to the transmitter and receiver planes,
and $v_{\text{D}}$ is the drift component perpendicular to the transmitter and
receiver planes. As a rule of thumb, we use bold fonts to stand for column
vectors; slanted vectors are $d$-dim, while non-slanted vectors are D-dim. The
operator $(\cdot)^{\mathsf{\scriptstyle{T}}}$ represents transposition. We
also introduce the notation $\mathbf{u}:=\frac{\mathbf{v}}{\sigma^{2}}$ and
similarly for $\boldsymbol{u}_{\text{par}}$ and $u_{1},\ldots,u_{\text{D}}$.
These $\\{u_{j}\\}_{j=1}^{\text{D}}$ can be interpreted as normalized drift.

In [12, 14], the 2D and 3D FAP densities were separately derived based on the
above listed system assumptions. By the methodology proposed in [13], we can
unify the two FAP densities (i.e., $\text{D}=2,3$) into a single expression:

$\displaystyle\begin{split}&f^{(d)}_{\boldsymbol{N}}(\boldsymbol{n})\\\
=\leavevmode\nobreak\
&2\lambda\left(\frac{\|\mathbf{u}\|}{\sqrt{2\pi}}\right)^{d+1}e^{\boldsymbol{u}_{\text{par}}^{\intercal}\boldsymbol{n}-u_{\text{D}}\lambda}\frac{K_{\tfrac{d+1}{2}}\left(\|\mathbf{u}\|\sqrt{\norm{\boldsymbol{n}}^{2}+\lambda^{2}}\right)}{\left(\|\mathbf{u}\|\sqrt{\norm{\boldsymbol{n}}^{2}+\lambda^{2}}\right)^{\tfrac{d+1}{2}}},\end{split}$
(3)

where $K_{\nu}(\cdot)$ denotes the order-$\nu$ modified Bessel function of the
second kind [16], and $\|\cdot\|$ denotes the Euclidean norm.

This paper focuses on a sub-family of FAP distributions with two specific
properties:

  1. 1.

The _parallel drift_ components are zero, i.e.,
$\boldsymbol{v}_{\text{par}}=\boldsymbol{0}$.

  2. 2.

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.)

This sub-family is referred to as Vertically-Drifted First Arrival Position
(VDFAP) distributions because only the vertical drift component is non-zero.
We can set $\boldsymbol{u}_{\text{par}}=\boldsymbol{0}$ so that
$\|\mathbf{u}\|$ in (3) 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)}_{\boldsymbol{N}}(\boldsymbol{n})=2\lambda\left(\frac{|u|}{\sqrt{2\pi}}\right)^{d+1}e^{\lambda|u|}\frac{K_{\tfrac{d+1}{2}}\left(|u|\sqrt{\norm{\boldsymbol{n}}^{2}+\lambda^{2}}\right)}{\left(|u|\sqrt{\norm{\boldsymbol{n}}^{2}+\lambda^{2}}\right)^{\tfrac{d+1}{2}}}.\end{split}$
(4)

We can thus view $u$ as a parameter of VDFAP distribution. As a shorthand, we
denote $\boldsymbol{X}\sim\text{VDFAP}^{(d)}(u,\lambda)$ if a $d$-dim random
vector $\boldsymbol{X}$ follows the VDFAP distribution defined by (4) 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 The Characteristic Function of VDFAP Distribution

From probability theory, distributions can be characterized by its
_characteristic function_ (CF). The CF viewpoint facilitates the computation
of moments and the analysis of stability properties for VDFAP distributions.

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

$\displaystyle\Phi^{(d)}_{\boldsymbol{N}}(\boldsymbol{\omega}):=\mathbb{E}[\exp(i\boldsymbol{\omega}^{\mathsf{\scriptstyle{T}}}\boldsymbol{N})],\text{\
for\ }\boldsymbol{\omega}\in\mathbb{R}^{d}.$ (5)

We have derived a novel closed-form expression for (5) when $d=1,2$ as
detailed in Appendix A. The resulting closed-form expression for the CF (5) is

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

where we have included the parameters $u,\lambda$ to emphasize the dependency.
It is important to note that letting $|u|\to 0$ in the CF formula (6) recovers
the $d$-dim Cauchy CF [17, 18]:

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

The parameters $u,\lambda$ will be suppressed hereafter when the context
allows for clear understanding. The observation in (7) is consistent with
[19], which indicates that as the drift approaches zero, the density of FAP
distribution given by (3) reduces to a Cauchy density.

Via the closed-form CF expression (6), we can obtain the first two moments, as
well as a weak stability property of VDFAP distributions.

###  III-A Mean Vector and Covariance Matrix

Based on (5), the mean vector of
$\boldsymbol{N}\sim\text{VDFAP}^{(d)}(u,\lambda)$ can be obtained as

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

and the correlation matrix of $\boldsymbol{N}$ as

$\displaystyle\mathbb{E}\left[\boldsymbol{N}\boldsymbol{N}^{\mathsf{\scriptstyle{T}}}\right]=-\nabla^{2}\Phi^{(d)}_{\boldsymbol{N}}(\boldsymbol{\omega})\bigg{|}_{\boldsymbol{\omega}=\boldsymbol{0}}.$
(9)

Due to the radial symmetry of $\boldsymbol{N}$, its mean is zero, and thus the
correlation matrix coincides with the covariance matrix. We first calculate
the gradient of the CF:

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

Then we can proceed from (10) to calculate the Hessian of the CF:

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

where $\mathbb{I}_{d}$ denotes the $d\times d$ identity matrix. Hence, the
covariance matrix of $\boldsymbol{N}$ can be expressed as:

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

It is interesting to note that the variance of each component of a random
vector $\boldsymbol{N}\sim\text{VDFAP}^{(d)}(u,\lambda)$ is
$\frac{\lambda}{|u|}$. Also, the components of $\boldsymbol{N}$ are pairwise
uncorrelated, but not independent, as the PDF (4) cannot be decomposed as the
product of its marginal PDFs.

###  III-B A Weak Stability Property

Cauchy distributions, along with Gaussian and Lévy distributions, are known to
be examples of _stable_ distributions. These distributions possess a specific
form of characteristic function (CF) [20], which is not satisfied by (6).
Therefore, we can conclude that VDFAP distributions are _not_ stable. However,
despite not being strictly stable, VDFAP distributions still exhibit certain
“weaker” stability properties.

One of the weaker stability properties of VDFAP distributions is that if
$\boldsymbol{N}_{1}\sim\text{VDFAP}^{(d)}(u,\lambda_{1})$ and
$\boldsymbol{N}_{2}\sim\text{VDFAP}^{(d)}(u,\lambda_{2})$ are two independent
VDFAP random vectors with the same (normalized) drift $u$, then their sum
$\boldsymbol{N}_{1}+\boldsymbol{N}_{2}$ follows the
$\text{VDFAP}^{(d)}(u,\lambda_{1}+\lambda_{2})$ distribution.

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

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

Since the CF uniquely characterizes the distribution, we can compare this
result to the CF formula (6) again to conclude that
$\boldsymbol{N}_{1}+\boldsymbol{N}_{2}\sim\text{VDFAP}^{(d)}(u,\lambda_{1}+\lambda_{2})$.

This _weak stability property_ will be useful in deriving a lower bound for
the capacity of VDFAP channels.

##  IV Bounds on the Capacity of 3D VDFAP Channel

Borrowing the wisdom from vector Gaussian interference channels we investigate
the VDFAP channel capacity under the covariance matrix constraint [15, 21] as
follows:

$\displaystyle\begin{split}C=\sup_{f(\boldsymbol{X}):\leavevmode\nobreak\
\mathbb{E}\left[\boldsymbol{X}\boldsymbol{X}^{\mathsf{\scriptstyle{T}}}\right]\preceq\Sigma}I(\boldsymbol{X};\boldsymbol{Y}),\end{split}$
(14)

where $\Sigma$ is a prescribed _positive definite_ matrix. Note that the
channel output $\boldsymbol{Y}$ equals to $\boldsymbol{X}+\boldsymbol{N}$,
where the noise $\boldsymbol{N}$ follows $\text{VDFAP}^{(2)}(u,\lambda)$
distribution, and is _independent_ of the input $\boldsymbol{X}$. The capacity
appeared in (14) depends on the parameter triplet $(u,\lambda,\Sigma)$.

Applying the additive channel structure
$\boldsymbol{Y}=\boldsymbol{X}+\boldsymbol{N}$, equation (14) can be
simplified as:

$\displaystyle\begin{split}C=\sup_{f(\boldsymbol{X}):\leavevmode\nobreak\
\mathbb{E}\left[\boldsymbol{X}\boldsymbol{X}^{\mathsf{\scriptstyle{T}}}\right]\preceq\Sigma}h(\boldsymbol{Y})-h(\boldsymbol{N}),\end{split}$
(15)

where $h(\cdot)$ denotes the differential entropy. We have calculated the
differential entropy of a random vector $\boldsymbol{N}$ following
$\text{VDFAP}^{(2)}(u,\lambda)$ in closed-form:

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

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

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

The calculation of (16) is given in Appendix B.

###  IV-A Lower Bound

Using the weak stability property of VDFAP distributions proved in III-B, we
obtain a lower bound on the capacity (14) of the VDFAP channel:

$\displaystyle C(u,\lambda,\Sigma)$ $\displaystyle\geq
2\log\Big(1+\tfrac{|u|}{\lambda}\sigma_{\text{min}}\Big{missing})-\log\left(1+\tfrac{|u|^{2}}{1+\lambda|u|}\sigma_{\text{min}}\right)$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
+\lambda|u|e^{\lambda|u|}\big{(}e\cdot\text{Ei}(-1-\lambda|u|)-3\cdot\text{Ei}(-\lambda|u|)\big{)}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
-\Big{(}\lambda+|u|\sigma_{\text{min}}\Big{)}|u|\exp\left\\{\left(\lambda+|u|\sigma_{\text{min}}\right)|u|\right\\}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
\times\bigg{\\{}e\cdot\text{Ei}\Big{(}-1-\big{(}\lambda+|u|\sigma_{\text{min}}\big{)}|u|\Big{)}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
-3\cdot\text{Ei}\Big{(}-\big{(}\lambda+|u|\sigma_{\text{min}}\big{)}|u|\Big{)}\bigg{\\}}$
(18) $\displaystyle>0,$

where $\sigma_{\text{min}}$ denotes the minimum eigenvalue of $\Sigma$.

###### Proof.

To obtain a lower bound on the capacity, we pick particular input
distributions $f(\boldsymbol{X})$ satisfying the second-moment constraint
$\mathbb{E}\left[\boldsymbol{X}\boldsymbol{X}^{\mathsf{\scriptstyle{T}}}\right]\preceq\Sigma$,
calculate $h(\boldsymbol{Y})-h(\boldsymbol{N})$, and then optimize the result
via taking supremum.

Take $\boldsymbol{X}\sim\text{VDFAP}^{(2)}(u,\lambda^{\prime})$ for some
positive $\lambda^{\prime}\leq|u|\sigma_{\text{min}}$. By (12), we have that
$\mathbb{E}\left[\boldsymbol{X}\boldsymbol{X}^{\mathsf{\scriptstyle{T}}}\right]=\frac{\lambda^{\prime}}{|u|}\mathbb{I}_{2}\preceq\sigma_{\text{min}}\mathbb{I}_{2}\preceq\Sigma,$
so $\boldsymbol{X}$ satisfies the covariance constraint. Also, by weak
stability property in III-B, $\boldsymbol{Y}=\boldsymbol{X}+\boldsymbol{N}$
follows $\text{VDFAP}^{(2)}(u,\lambda+\lambda^{\prime})$. Then, a lower bound
$C\geq h(\boldsymbol{Y})-h(\boldsymbol{N})$ follows from (15).

To analyze this lower bound, we introduce an ancillary function
$h_{0}:\mathbb{R}^{+}\to\mathbb{R}$ defined by

$\displaystyle\begin{split}h_{0}(s)&:=2\log(s)-\log(1+s)\\\ &\quad-
se^{s}\left(e\cdot\text{Ei}(-1-s)-3\cdot\text{Ei}(-s)\right),\text{\ for\
}s>0,\end{split}$ (19)

so that the differential entropy (16) of a random vector $\boldsymbol{Z}$,
which follows $\text{VDFAP}^{(2)}(u,\lambda^{\prime\prime})$, can be written
as

$\displaystyle h(\mathbf{Z})=h_{0}(|u|\lambda^{\prime\prime})+\log(2\pi
e^{3})-2\log(|u|).$ (20)

Therefore, the lower bound can be expressed via the function $h_{0}(\cdot)$ as
$C\geq h_{0}(|u|(\lambda+\lambda^{\prime}))-h_{0}(|u|\lambda),$ which holds
for any $0<\lambda^{\prime}\leq|u|\sigma_{\text{min}}$. Hence, we can maximize
over $\lambda^{\prime}$ to get that

$\displaystyle
C\geq\sup_{0<\lambda^{\prime}\leq|u|\sigma_{\text{min}}}h_{0}(|u|(\lambda+\lambda^{\prime}))-h_{0}(|u|\lambda).$
(21)

It is shown in Appendix C that $h_{0}(\cdot)$ is strictly increasing.
Therefore, (21) is equivalent to

$\displaystyle C\geq
h_{0}\bigg{(}|u|\Big{(}\lambda+|u|\sigma_{\text{min}}\Big{)}\bigg{)}-h_{0}(|u|\lambda).$
(22)

Plugging (19) into (22) yields (18). The lower bound in (22), and hence (18),
is strictly positive due to the two facts that $h_{0}(\cdot)$ is strictly
increasing, and that $\sigma_{\text{min}}$ is strictly positive as
$\Sigma\succ 0$. ∎

###  IV-B Upper Bound

It is well-established in the literature that the differential entropy is
maximized by multivariate Gaussian under a prescribed covariance matrix [22,
Chapter 12],[23]. This fact can be applied to establish an upper bound on the
capacity. Under the constraint in (14), the bound can be expressed as:

$\displaystyle C(u,\lambda,\Sigma)$
$\displaystyle\leq\frac{1}{2}\log\left(\det\Big(\tfrac{1}{\lambda^{2}}\Sigma+\tfrac{1}{\lambda|u|}\mathbb{I}_{2}\Big{missing})\right)+\log(1+\lambda|u|)$
$\displaystyle\quad+\lambda|u|e^{\lambda|u|}\left(e\cdot\text{Ei}(-1-\lambda|u|)-3\cdot\text{Ei}(-\lambda|u|)\right)$
$\displaystyle\quad-2,$ (23)

where $\det(\cdot)$ denotes the determinant.

###### Proof.

Because $\boldsymbol{X}$ and $\boldsymbol{N}$ are independent and that
$\mathbb{E}[\boldsymbol{N}]=\boldsymbol{0}$,

$\displaystyle\mathbb{E}\left[\boldsymbol{Y}\boldsymbol{Y}^{\mathsf{\scriptstyle{T}}}\right]$
$\displaystyle=\leavevmode\nobreak\ $
$\displaystyle\mathbb{E}\left[\left(\boldsymbol{X}+\boldsymbol{N}\right)\left(\boldsymbol{X}+\boldsymbol{N}\right)^{\mathsf{\scriptstyle{T}}}\right]$
$\displaystyle=\leavevmode\nobreak\ $
$\displaystyle\mathbb{E}\left[\boldsymbol{X}\boldsymbol{X}^{\mathsf{\scriptstyle{T}}}+\boldsymbol{X}\boldsymbol{N}^{\mathsf{\scriptstyle{T}}}+\boldsymbol{N}\boldsymbol{X}^{\mathsf{\scriptstyle{T}}}+\boldsymbol{N}\boldsymbol{N}^{\mathsf{\scriptstyle{T}}}\right]$
$\displaystyle\overset{(a)}{=}\leavevmode\nobreak\ $
$\displaystyle\mathbb{E}\left[\boldsymbol{X}\boldsymbol{X}^{\mathsf{\scriptstyle{T}}}\right]+\mathbb{E}\left[\boldsymbol{X}\right]\mathbb{E}\left[\boldsymbol{N}\right]^{\mathsf{\scriptstyle{T}}}+\mathbb{E}\left[\boldsymbol{N}\right]\mathbb{E}\left[\boldsymbol{X}\right]^{\mathsf{\scriptstyle{T}}}+\mathbb{E}\left[\boldsymbol{N}\boldsymbol{N}^{\mathsf{\scriptstyle{T}}}\right]$
$\displaystyle\overset{(b)}{=}\leavevmode\nobreak\ $
$\displaystyle\mathbb{E}\left[\boldsymbol{X}\boldsymbol{X}^{\mathsf{\scriptstyle{T}}}\right]+\mathbb{E}\left[\boldsymbol{N}\boldsymbol{N}^{\mathsf{\scriptstyle{T}}}\right]$
$\displaystyle\overset{(c)}{\preceq}\leavevmode\nobreak\ $
$\displaystyle\Sigma+\tfrac{\lambda}{|u|}\mathbb{I}_{2},$ (24)

where (a) is due to independence, (b) is due to
$\mathbb{E}[\boldsymbol{N}]=\boldsymbol{0}$, and (c) is due to the constraint
$\mathbb{E}\left[\boldsymbol{X}\boldsymbol{X}^{\mathsf{\scriptstyle{T}}}\right]\preceq\Sigma$
and (12).

Under the covariance constraint (24) on the output $\boldsymbol{Y}$, its
differential entropy $h(\boldsymbol{Y})$ is maximized when $\boldsymbol{Y}$
follows the bivariate Gaussian distribution
$\mathcal{N}\Big{(}\mathbf{0},\Sigma+\tfrac{\lambda}{|u|}\mathbb{I}_{2}\Big{)}$.
Hence, from (15) we have

$\displaystyle C\leq
h\bigg{(}\mathcal{N}\Big{(}\mathbf{0},\Sigma+\tfrac{\lambda}{|u|}\mathbb{I}_{2}\Big{)}\bigg{)}-h\left(\text{VDFAP}^{(2)}(u,\lambda)\right).$
(25)

Since the differential entropy of the bivariate Gaussian is

$\displaystyle
h\bigg{(}\mathcal{N}\Big{(}\mathbf{0},\Sigma+\tfrac{\lambda}{|u|}\mathbb{I}_{2}\Big{)}\bigg{)}\bigg{)}$
$\displaystyle=\frac{1}{2}\log\Big((2\pi
e)^{2}\det\big(\Sigma+\tfrac{\lambda}{|u|}\mathbb{I}_{2}\big{missing})\Big{missing}),$
(26)

together with (16) we can simplify (25) into (23). ∎

##  V Conclusion

The concept of First Arrival Position (FAP) channels has been recently
introduced in the field of diffusive molecular communication. This paper has
advanced the understanding of this novel concept by addressing the capacity
problem of additive vertically-drifted (VD) FAP noise channels, specifically
under the covariance matrix constraint.

The primary contribution of our research lies in providing explicit
expressions for both the upper and lower bounds of this capacity as the
spatial dimension equals $3$. To arrive at these expressions, we carried out a
meticulous analysis of the characteristic function of VDFAP distributions.
Furthermore, our research has uncovered a novel stability property of VDFAP.

Our findings contribute to the ongoing efforts to comprehend the fundamental
limits of molecular communication systems. They open new avenues for future
research, including extending the analysis to other spatial dimensions and
exploring the implications of the discovered stability property.

## References

  * [1] T. Nakano, A. W. Eckford, and T. Haraguchi, _Molecular Communication_. Cambridge University Press, 2013.
  * [2] 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.
  * [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] 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.
  * [5] 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.
  * [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. 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.
  * [8] 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.
  * [9] A. W. Eckford, K. Srinivas, and R. S. Adve, “The peak constrained additive inverse Gaussian noise channel,” in _Proc. IEEE Int. Symp. Information Theory_ , Cambridge, MA, USA, July 2012, pp. 2983–2987.
  * [10] H.-T. Chang and S. M. Moser, “Bounds on the capacity of the additive inverse Gaussian noise channel,” in _Proc. IEEE Int. Symp. Information Theory_ , Cambridge, MA, USA, July 2012, pp. 299–303.
  * [11] 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.
  * [12] 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.
  * [13] Y.-C. Lee, J.-M. Wu, and M.-H. Hsieh, “A unified framework for calculating first arrival position density in molecular communication,” arXiv:2201.04476v3 [cs.IT], Apr. 2023.
  * [14] 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_ , Barcelona, Spain, July 2016, pp. 1033–1037.
  * [15] X. Shang, B. Chen, G. Kramer, and H. V. Poor, “Capacity regions and sum-rate capacities of vector Gaussian interference channels,” _IEEE Transactions on Information Theory_ , vol. 56, no. 10, pp. 5030–5044, 2010.
  * [16] I. S. Gradshteyn and I. M. Ryzhik, _Table of Integrals, Series, and Products_. Academic Press, 2014.
  * [17] 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.
  * [18] S. Kotz and S. Nadarajah, _Multivariate t-distributions and Their Applications_. Cambridge University Press, 2004.
  * [19] Y.-C. Lee and M.-H. Hsieh, “On the capacity of zero-drift first arrival position channels in diffusive molecular communication,” arXiv:2201.11383v2 [cs.IT], May 2023.
  * [20] 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.
  * [21] M. Vu, “MISO capacity with per-antenna power constraint,” _IEEE Transactions on Communications_ , vol. 59, no. 5, pp. 1268–1274, 2011.
  * [22] T. M. Cover, _Elements of Information Theory_. John Wiley & Sons, 1999.
  * [23] N. Ebrahimi, E. S. Soofi, and R. Soyer, “Multivariate maximum entropy identification, transformation, and dependence,” _Journal of Multivariate Analysis_ , vol. 99, no. 6, pp. 1217–1231, 2008.
  * [24] 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.
  * [25] Y. L. Luke, _The Special Functions and Their Approximations, Vol. 2_. Elsevier, 1969.

##  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 (4).
Specifically, a random variable $N$ following the Student’s $t$-distribution
of $\nu$ degrees of freedom ($\nu>0$) has the following PDF [24]:

$\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},$ (27)

where $\Gamma(\cdot)$ is the gamma function, and the corresponding CF:

$\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}.$ (28)

This implies that for $\omega\in\mathbb{R}\setminus\\{0\\}$,

$\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.$ (29)

Fix an arbitrary $s>0$. Applying the substitutions $\tilde{n}=(s/\sqrt{\nu})n$
and $\tilde{\omega}=(\sqrt{\nu}/s)\omega$ to (29), 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}$ (30)

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

$\displaystyle\begin{split}f^{(d)}_{\boldsymbol{N}}(\boldsymbol{n})&=\frac{\Gamma\left(\frac{d+2}{2}\right)}{\pi^{\frac{d+2}{2}}}\lambda
e^{\lambda|u|}\\\
&\quad\cdot\int_{\mathbb{R}}\left(\norm{\boldsymbol{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}$ (31)

Thus for the $d$-dim VDFAP distribution, its CF

$\displaystyle\Phi^{(d)}_{\boldsymbol{N}}(\boldsymbol{\omega})$
$\displaystyle=\int_{\mathbb{R}^{d}}f^{(d)}_{\boldsymbol{N}}(\boldsymbol{n})\leavevmode\nobreak\
e^{i\boldsymbol{\omega}^{\intercal}\boldsymbol{n}}\leavevmode\nobreak\
\mathrm{d}\boldsymbol{n}$ (32)

can be obtained by plugging (31) into (32), resulting in

$\displaystyle\begin{split}\Phi^{(d)}_{\boldsymbol{N}}(\boldsymbol{\omega})=&\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}$ (33)

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

We recognize that the integral in (33) 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 $\boldsymbol{\mu}=\boldsymbol{0}$ and scale
$\boldsymbol{\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}}},$
(34)

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

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

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

##  Appendix B Calculation of the 3D Differential Entropy

From (4), we can write $f^{(2)}_{\boldsymbol{N}}(\boldsymbol{n})=c\cdot
f(\norm{\boldsymbol{n}})$ where

$\displaystyle\begin{split}c&:=\dfrac{\lambda}{\sqrt{2\pi^{3}}}|u|^{3}e^{\lambda|u|};\leavevmode\nobreak\
f(r):=\dfrac{K_{3/2}\left(|u|\sqrt{r^{2}+\lambda^{2}}\right)}{\left(|u|\sqrt{r^{2}+\lambda^{2}}\right)^{3/2}}.\end{split}$
(36)

Then, the differential entropy of a random vector $\boldsymbol{N}$ following
$\text{VDFAP}^{(2)}(u,\lambda)$, by definition, is

$\displaystyle
h(\boldsymbol{N}):=-\int_{\mathbb{R}^{2}}f^{(2)}_{\boldsymbol{N}}(\boldsymbol{n})\log\left(f^{(2)}_{\boldsymbol{N}}(\boldsymbol{n})\right)\mathrm{d}\boldsymbol{n}.$
(37)

Since VDFAP distributions are radially symmetric, converting to polar
coordinates yields:

$\displaystyle h(\boldsymbol{N})=-\int_{0}^{2\pi}\int_{0}^{\infty}c\cdot
f(r)\cdot(\log c+\log f(r))\cdot r\leavevmode\nobreak\
\mathrm{d}r\leavevmode\nobreak\ \mathrm{d}\theta,$ (38)

where we have applied (36). Using the fact that
$f^{(2)}_{\boldsymbol{N}}(\boldsymbol{n})$ is a PDF on $\mathbb{R}^{2}$ and
thus integrate to one, we can further express:

$\displaystyle h(\boldsymbol{N})=-\log c-2\pi c\int_{0}^{\infty}rf(r)\log
f(r)\leavevmode\nobreak\ \mathrm{d}r.$ (39)

Thus, it remains to evaluate the improper integral. Applying the change of
variable

$\displaystyle\begin{split}\rho:=|u|\sqrt{r^{2}+\lambda^{2}};\leavevmode\nobreak\
\leavevmode\nobreak\
\mathrm{d}\rho=|u|\frac{r}{\sqrt{r^{2}+\lambda^{2}}}\mathrm{d}r=|u|^{2}\frac{r}{\rho}\mathrm{d}r\end{split}$
(40)

to this improper integral results in

$\displaystyle\begin{split}&\int_{0}^{\infty}rf(r)\log
f(r)\leavevmode\nobreak\
\mathrm{d}r=\int_{|u|\lambda}^{\infty}\frac{K_{3/2}(\rho)}{|u|^{2}\rho^{1/2}}\log\left(\frac{K_{3/2}(\rho)}{\rho^{3/2}}\right)\mathrm{d}\rho,\end{split}$
(41)

which can be evaluated using the following integral formula: 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}$ (42)

To further express (39) solely in terms of parameters $u$ and $\lambda$, we
calculate: $-\log
c=-\log\lambda-3\log|u|-|u|\lambda+\frac{1}{2}\log(2\pi^{3})$ and $2\pi
c=\sqrt{\frac{2}{\pi}}\lambda|u|^{3}e^{|u|\lambda}.$ Therefore, taking
$a=|u|\lambda$ in (42), plugging the result into (41) and then into (39)
yields

$\displaystyle h(\boldsymbol{N})$
$\displaystyle=-\log\lambda-3\log|u|-|u|\lambda+\tfrac{1}{2}\log(2\pi^{3})$
$\displaystyle\quad-\lambda|u|e^{|u|\lambda}\left(e\cdot\text{Ei}(-1-|u|\lambda)-3\cdot\text{Ei}(-|u|\lambda)\right)$
$\displaystyle\quad+\lambda|u|+\tfrac{1}{2}\left(6\log(|u|\lambda)-2\log(1+|u|\lambda)\right)$
$\displaystyle\quad+\tfrac{1}{2}\left(6+\log\tfrac{2}{\pi}\right).$ (43)

After collecting terms and simplifying, we arrived at (16).

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

Observe that $h_{0}(s)=2\log(s)-\log(1+s)-g(s)$, where we introduce 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.$ (44)

Taking the derivative of (44) and using the definition of $\text{Ei}(\cdot)$
given in (17), we obtain the formula

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

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)).$ (46)

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

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

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}$
(48)

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