# A Scalable Track-Before-Detect Method
With Poisson/Multi-Bernoulli Model

Thomas Kropfreiter1, Jason L. Williams2, and Florian Meyer1 This research was
supported by the Office of Naval Research under Grant N00014-21-1-2267.
1Scripps Institution of Oceanography and the Department of Electrical and
Computer Engineering,
University of California San Diego, La Jolla, CA, USA (e-mail:
<EMAIL_ADDRESS><EMAIL_ADDRESS> 2Data61, Commonwealth Scientific
and Industrial Research Organisation, Australia
<EMAIL_ADDRESS>
###### Abstract

We propose a scalable track-before-detect (TBD) tracking method based on a
Poisson/multi-Bernoulli model. To limit computational complexity, we
approximate the exact multi-Bernoulli mixture posterior probability density
function (pdf) by a multi-Bernoulli pdf. Data association based on the sum-
product algorithm and recycling of Bernoulli components enable the detection
and tracking of low-observable objects with limited computational resources.
Our simulation results demonstrate a significantly improved tracking
performance compared to a state-of-the-art TBD method.

###### Index Terms:

Multiobject tracking, track-before-detect, message passing, random finite
sets.

##  I Introduction

Multiobject tracking aims to estimate the time-dependent number and states of
multiple objects based on data provided by one or more sensors. The
measurements of conventional multiobject tracking methods are produced by a
detector that performs thresholding and, possibly, some further preparatory
processing of the raw sensor data [1, 2, 3, 4, 5, 6, 7, 8, 9]. This
preparatory processing step reduces data flow and computational complexity,
but results in a loss of tracking information. Thus, the conventional _detect-
than-track_ approach can lead to a significantly reduced tracking performance,
especially in scenarios with low signal-to-noise (SNR) ratio [10, 11, 12, 13,
14, 15, 16, 17].

In such scenarios, tracking methods that use raw sensor data as measurements,
can potentially achieve an improved tracking performance. Many existing
_track-before-detect_ (TBD) methods are based on batch processing using
maximum likelihood estimation [10], the Hough transform [11], and dynamic
programming techniques [12]. However, due to their high computational
complexity, they are often unsuitable for real-time operation. A well-
established real-time TBD method is the Histogram Probabilistic Multi-
Hypothesis Tracker (H-PMHT) [13], which is based on the expectation-
maximization algorithm. However, tuning the parameters of the H-PMHT is known
to be difficult [17]. Another class of real-time TBD tracking methods is those
derived using random finite sets (RFSs). These methods include the (single)
Bernoulli filter [14, 15] for single-object tracking and multi-Bernoulli (MB)
filters [16, 17] for multiobject tracking. Existing MB filters [16, 17] for
multiobject TBD rely on certain heuristics to introduce Bernoulli components
for objects that appear in the scene for the first time.

A state-of-the-art conventional tracking method is the track-oriented marginal
multi-Bernoulli/Poisson (TOMB/P) filter [5, 6, 8], which models an unknown
number of objects and their states by the union of a Poisson RFS and a MB RFS.
The TOMB/P filter relies on a high-dimensional marginalization operation
performed by means of message passing [6, 18]. A modified version of the
TOMB/P filter [19] “recycles” Bernoulli components with low existence
probability by transferring them to the Poisson RFS instead of pruning them.
This recycling step can significantly improve tracking performance in
scenarios with reduced sensor performance [19]. The TOMB/P filter has also
been derived and extended to multiple sensors using the framework of factor
graphs and the sum-product algorithm (SPA) [6, 7]. It has been demonstrated
that it can outperform existing detect-than-track methods and that it is
highly scalable in relevant system parameters [5, 6, 7, 8].

In this paper, we propose a scalable TBD method that is based on a Poisson/MB
model similar to the one of the TOMB/P filter. We first review the general TBD
measurement model and discuss a simplification that is suitable for the case
where at any discrete time step, the influence of an object is limited to a
single “bin” or “pixel” of the raw sensor data. Based on this simplified TBD
measurement model, we derive an exact update step. Since the resulting exact
posterior probability density function (pdf) is a MB mixture pdf, we
approximate it by an MB pdf to limit computational complexity. To guarantee
reliable multiobject tracking performance, we recycle Bernoulli components
with a low existence probability instead of pruning them. Finally, we verify
the performance of the proposed method numerically and compare it to a state-
of-the-art TBD tracking method [16].

The main contributions of this paper can be summarized as follows:

  * •

We establish a Bayesian Poisson/MB model for the development of scalable TBD
methods.

  * •

We derive a scalable Poisson/MB filter for TBD that relies on recycling of
Bernoulli components.

  * •

We demonstrate performance improvements of our TBD method compared to a state-
of-the-art technique.

_Notation:_ We will use the following basic notation. Random variables are
displayed in sans serif, upright fonts; their realizations in serif, italic
fonts. Vectors and matrices are denoted by bold lowercase and uppercase
letters, respectively. For example, $\mathsfbr{x}$ is a random variable, and
$x$ is its realization and $\bm{\mathsfbr{x}}$ is a random vector and $\bm{x}$
is its realization. Random sets and their realizations are denoted by upright
sans serif and calligraphic font, respectively. For example, $\mathsfbr{X}$ is
a random set and ${\cal{X}}$ is its realization. We denote pdfs by $f(\cdot)$
and probability mass functions (pmfs) by $p(\cdot)$.
${\cal{N}}(\bm{x};\bm{\mu},\bm{\varSigma})$ denotes the Gaussian pdf (of
random vector $\bm{\mathsfbr{x}}$) with mean $\bm{\mu}$ and covariance matrix
$\bm{\varSigma}$ and ${\cal{R}}(x;\sigma)$ denotes the Rayleigh pdf (of scalar
random variable $\mathsfbr{x}$) with scale parameter $\sigma$. The probability
of an event is denoted $\mathrm{Pr}\\{\cdot\\}$. The symbol $\propto$
indicates equality up to a normalization factor.

##  II RFS Fundamentals

A RFS $\mathsfbr{X}$ is a random variable whose realizations ${\cal{X}}$ are
finite sets $\big{\\{}\bm{x}^{(1)},\ldots,\bm{x}^{(n)}\big{\\}}$ of vectors
$\bm{x}^{(i)}\\!\in\mathbb{R}^{n_{x}}$. Both the vectors
$\bm{\mathsfbr{x}}^{(i)}$ and their number ${\sf n}=|\mathsfbr{X}|$ (the
cardinality of $\mathsfbr{X}$) are random. Thus, $\mathsfbr{X}$ consists of a
random number ${\sf n}$ of random vectors
$\bm{\mathsfbr{x}}^{(1)},\ldots,\bm{\mathsfbr{x}}^{({\sf n})}$. While the
conventional Riemann integral is not defined for sets, one can define the set
integral of a real-valued set function $g({\cal{X}})$ as [4]

$\displaystyle\int g({\cal{X}})\hskip 0.85358pt\hskip
0.85358pt\delta{\cal{X}}$ $\displaystyle\hskip
8.53581pt\triangleq\sum_{n=0}^{\infty}\hskip 0.85358pt\frac{1}{n!}\hskip
0.85358pt\int_{\mathbb{R}^{nn_{x}}}\\!g(\\{\bm{x}^{(1)},\ldots,\bm{x}^{(n)}\\})\hskip
0.85358pt\hskip 0.85358pt\text{d}\bm{x}^{(1)}\cdots\text{d}\bm{x}^{(n)}.$ (1)

Note that each term of the sum corresponds to one value of the cardinality
$n=|{\cal{X}}|$. The statistics of an RFS $\mathsfbr{X}$ can be described by
its multiobject pdf $f({\cal{X}})$. For any realization
${\cal{X}}=\big{\\{}\bm{x}^{(1)},\ldots,\bm{x}^{(n)}\big{\\}}$, the
multiobject pdf $f({\cal{X}})$ is given by

$f({\cal{X}})\hskip 0.85358pt=\hskip 0.85358ptn!\hskip 0.85358pt\rho(n)\hskip
0.85358ptf_{n}(\bm{x}^{(1)},\ldots,\bm{x}^{(n)})\hskip 0.85358pt.\vspace{1mm}$
(2)

Here, $\rho(n)\triangleq\mathrm{Pr}\\{|\mathsfbr{X}|=n\\}$,
$n\in\mathbb{N}_{0}$ is the pmf of ${\sf n}=|\mathsfbr{X}|$, and
$f_{n}(\bm{x}^{(1)},\ldots,\bm{x}^{(n)})$ is a joint pdf of the random vectors
$\bm{\mathsfbr{x}}^{(1)},\ldots,\bm{\mathsfbr{x}}^{(n)}$ that is invariant to
a permutation of its arguments. Note that based on the set integral (1),
$f({\cal{X}})$ integrates to one.

Next, we will review five classes of RFSs [4, 5, 6] that are relevant for the
derivation of the proposed method. A Poisson RFS $\mathsfbr{X}$ is
characterized by a cardinality that is Poisson distributed with mean $\mu$,
i.e., $\rho(n)=e^{-\mu}\mu^{n}/n!\hskip 0.85358pt$, $n\\!\in\\!\mathbb{N}_{0}$
and by elements $\bm{\mathsfbr{x}}^{(1)},\ldots,\bm{\mathsfbr{x}}^{(n)}$ that
are independent and identically distributed (iid) according to the spatial pdf
$f(\bm{x})$, i.e,
$f_{n}(\bm{x}^{(1)},\ldots,\bm{x}^{(n)})=\prod_{i=1}^{n}f(\bm{x}^{(i)})$.
Following (2) the multiobject pdf now reads

$f^{\text{P}}({\cal{X}})\hskip 0.85358pt=\hskip
0.85358pte^{-\int\lambda(\bm{x}^{\prime})\hskip
0.85358pt\text{d}\bm{x}^{\prime}}\prod_{\bm{x}\hskip
0.85358pt\in{\cal{X}}}\lambda(\bm{x})\vspace{-.5mm}$

where $\lambda(\bm{x})=\mu f(\bm{x})$ is called the probability hypothesis
density (PHD) or intensity function.

A Bernoulli RFS $\mathsfbr{X}$ is represented by an existence probability $r$
and a spatial pdf $f(\bm{x})$. It consists of either none or one element with
probability $1-r$ and $r$, respectively. According to (2), the multiobject pdf
thus reads

$f({\cal{X}})\hskip
0.85358pt=\begin{cases}1\\!-r,&{\cal{X}}\\!=\emptyset,\\\\[-1.70717pt] r\hskip
0.85358ptf(\bm{x})\hskip
0.85358pt,&{\cal{X}}\\!=\\!\\{\bm{x}\\},\\\\[-1.70717pt]
0,&\text{otherwise}.\end{cases}\vspace{-1.5mm}$ (3)

A MB RFS $\mathsfbr{X}$ is the union of a fixed number $J$ of statistically
independent Bernoulli RFSs $\mathsfbr{X}^{(j)}$, $j\in\\{1,\ldots,J\\}$ with
multiobject pdfs $f^{(j)}({\cal{X}})$ (cf. (3)) described by the existence
probabilities $r^{(j)}$ and the spatial pdfs $f^{(j)}(\bm{x})$. The
multiobject pdf $f^{\text{MB}}({\cal{X}})$ of the MB RFS can be obtained by
applying the set convolution [4] to the individual Bernoulli pdfs
$f^{(j)}({\cal{X}})$. For any realization
${\cal{X}}=\\{\bm{x}^{(1)},\ldots,\bm{x}^{(n)}\\}$ with $n\leqslant J$, the
multiobject pdf $f^{\text{MB}}({\cal{X}})$ can be evaluated as

$f^{\text{MB}}({\cal{X}})=\sum_{{\cal{X}}^{(1)}\uplus\ldots\uplus\hskip
0.85358pt{\cal{X}}^{(J)}={\cal{X}}}\prod_{j\hskip 0.85358pt=\hskip
0.85358pt1}^{J}\hskip 0.85358pt\hskip 0.85358ptf^{(j)}({\cal{X}}^{(j)})\hskip
0.85358pt.$ (4)

Here, $\sum_{{\cal{X}}^{(1)}\uplus\ldots\uplus{\cal{X}}^{(J)}={\cal{X}}}$
denotes the sum over all disjoint decompositions of ${\cal{X}}$ into sets
${\cal{X}}^{(j)}$, $j\in\\{1\ldots,J\\}$ such that
${\cal{X}}^{(1)}\cup\ldots\cup{\cal{X}}^{(J)}={\cal{X}}$. For example, for
${\cal{X}}=\big{\\{}\bm{x}^{(1)},\bm{x}^{(2)}\big{\\}}$ and $J=2$, the
multiobject pdf $f^{\text{MB}}({\cal{X}})$ can be evaluated as

$\displaystyle f^{\text{MB}}\big{(}\\{\bm{x}^{(1)},\bm{x}^{(2)}\\}\big{)}$
$\displaystyle=r^{(1)}\hskip 0.85358pt\hskip 0.85358ptr^{(2)}\hskip
0.85358pt\hskip 0.85358pt\Big{(}f^{(1)}\big{(}\bm{x}^{(1)}\big{)}\hskip
0.85358ptf^{(2)}\big{(}\bm{x}^{(2)}\big{)}+f^{(1)}\big{(}\bm{x}^{(2)}\big{)}\hskip
0.85358ptf^{(2)}\big{(}\bm{x}^{(1)}\big{)}\Big{)}.$

Similarly, for ${\cal{X}}=\big{\\{}\bm{x}^{(1)}\big{\\}}$ and $J=3$, the
multiobject pdf $f^{\text{MB}}({\cal{X}})$ reads

$\displaystyle f^{\text{MB}}\big{(}\\{\bm{x}^{(1)}\\}\big{)}$
$\displaystyle=r^{(1)}\hskip 0.85358pt\hskip
0.85358pt(1-r^{(2)})(1-r^{(3)})\hskip 0.85358pt\hskip
0.85358ptf^{(1)}\big{(}\bm{x}^{(1)}\big{)}$ $\displaystyle+(1-r^{(1)})\hskip
0.85358pt\hskip 0.85358ptr^{(2)}\hskip 0.85358pt\hskip
0.85358pt(1-r^{(3)})\hskip 0.85358pt\hskip
0.85358ptf^{(2)}\big{(}\bm{x}^{(1)}\big{)}$ $\displaystyle+(1-r^{(1)})\hskip
0.85358pt\hskip 0.85358pt(1-r^{(2)})\hskip 0.85358pt\hskip
0.85358ptr^{(3)}\hskip 0.85358pt\hskip
0.85358ptf^{(3)}\big{(}\bm{x}^{(1)}\big{)}.$

Note that for $n>J$, we have $f^{\text{MB}}({\cal{X}})=0$.

A Poisson/MB RFS $\mathsfbr{X}$ is the union of a Poisson RFS and a MB RFS.
The pdf of a Poisson/MB RFS can be obtained by applying the set convolution to
the pdf of the Poisson RFS and the pdf of the MB RFS. Let $J$ be the number of
components of the MB RFS. For any realization
${\cal{X}}=\\{\bm{x}^{(1)},\ldots,\bm{x}^{(n)}\\}$ with $n\in\mathbb{N}_{0}$,
the multiobject pdf $f^{\text{PMB}}({\cal{X}})$ can be evaluated as

$\displaystyle f^{\text{PMB}}({\cal{X}})\hskip
5.69054pt=\sum_{{\cal{X}}^{(0)}\uplus\hskip 0.85358pt{\cal{X}}^{(1)}=\hskip
0.85358pt{\cal{X}}}f^{\text{P}}\big{(}{\cal{X}}^{(0)}\big{)}\hskip
0.85358ptf^{\text{MB}}\big{(}{\cal{X}}^{(1)}\big{)}.$ (5)

Here, $\sum_{{\cal{X}}^{(0)}\uplus\hskip 0.85358pt{\cal{X}}^{(1)}=\hskip
0.85358pt{\cal{X}}}$ denotes the sum over all disjoint decompositions of
${\cal{X}}$ into two sets ${\cal{X}}^{(0)}$ and ${\cal{X}}^{(1)}$ such that
${\cal{X}}^{(0)}\cup\hskip 0.85358pt{\cal{X}}^{(1)}={\cal{X}}$.

Finally, a MB mixture RFS is a weighted sum of MB RFSs where without loss of
generality, we assume that all MB RFSs have the same number of Bernoulli
components $J$. For any realization
${\cal{X}}=\\{\bm{x}^{(1)},\ldots,\bm{x}^{(n)}\\}$ with $n\leqslant J$, the
multiobject pdf $f^{\text{MBM}}({\cal{X}})$ can be evaluated as

$f^{\text{MBM}}({\cal{X}})=\sum^{I}_{i\hskip 0.85358pt=\hskip
0.85358pt1}w_{i}\hskip 0.85358ptf_{i}^{\text{MB}}({\cal{X}})\vspace{-1mm}$ (6)

where $I$ is the number of different MB pdfs and $\sum^{I}_{i\hskip
0.85358pt=\hskip 0.85358pt1}w_{i}=1$. Note that in the implementation of our
proposed tracking method, the sums in (4), (5) and (6) are never explicitly
evaluated.

##  III System Model

In this section, we describe the system model underlying the proposed
algorithm. The multiobject state at time $k$ is represented by an RFS
$\mathsfbr{X}_{k}=\\{\bm{\mathsfbr{x}}_{k}^{(1)},\ldots,\bm{\mathsfbr{x}}_{k}^{(n)}\\}$.
The single-object state $\bm{\mathsfbr{x}}_{k}^{(i)}$, $i\in\\{1,\ldots,n\\}$
consists of the object’s intensity $\sgamma_{k}^{(i)}$, position
$\bm{\mathsfbr{p}}_{k}^{(i)}$, and possibly motion related parameters. The
measurement
$\bm{\mathsfbr{z}}_{k}\triangleq[\mathsfbr{z}_{k}^{(1)}\ldots\mathsfbr{z}_{k}^{(M)}]^{\text{T}}$
at time $k$, consists of $M$ non-negative scalars
$\mathsfbr{z}_{k}^{(m)}\geqslant 0$, $m\in\\{1,\ldots,M\\}$, that represent
the measured intensity of the $m$’th “bin” or “pixel”.

###  III-A State-Transition Model

We use the well-established conventional RFS state-transition model [4, 5, 6].
At time $k-\\!1$, an object with state
$\bm{\mathsfbr{x}}_{k-1}\\!\in\mathsfbr{X}_{k-1}$ either survives or dies with
probabilities $p_{\text{S}}$ and $1\\!-p_{\text{S}}$, respectively. If it
survives, its new state $\bm{\mathsfbr{x}}_{k}$ is distributed according to
the single-object state transition pdf $f(\bm{x}_{k}|\bm{x}_{k-1})$.

We assume that the states of different objects survive/die and evolve in time
independently, i.e., given $\bm{x}_{k-1}$, the single-object state
$\bm{\mathsfbr{x}}_{k}$ is conditionally independent of all the other single-
object states $\bm{\mathsfbr{x}}_{k}^{\prime}$. Thus, conditioned on the
multiobject state ${\cal{X}}_{k-1}$, the multiobject state of the survived
objects $\mathsfbr{X}_{k}^{\text{S}}$, can be modeled by an MB RFS, i.e.,
$\mathsfbr{X}_{k}^{\text{S}}=\bigcup_{\bm{x}_{k-1}\in{\cal{X}}_{k-1}}\\!\mathsfbr{S}_{k}(\bm{x}_{k-1})$.
The components $\mathsfbr{S}_{k}(\bm{x}_{k-1})$ of this MB RFS are Bernoulli
RFSs (cf. (3)) with existence probabilities $p_{\text{S}}$ and spatial pdfs
$f(\bm{x}_{k}|\bm{x}_{k-1})$.

Newborn objects are modeled by a Poisson RFS $\mathsfbr{X}_{k}^{\text{B}}$
with mean $\mu_{\text{B}}$, spatial pdf $f_{\text{B}}(\bm{x}_{k})$ and, hence,
PHD $\lambda_{\text{B}}(\bm{x}_{k})=\mu_{\text{B}}\hskip
0.85358ptf_{\text{B}}(\bm{x}_{k})$. Conditioned on ${\cal{X}}_{k-1}$,
surviving objects $\mathsfbr{X}_{k}^{\text{S}}$ are assumed independent of the
newborn objects $\mathsfbr{X}_{k}^{\text{B}}$. Thus, for ${\cal{X}}_{k-1}$
fixed, the overall multiobject state at time $k$, is given by

$\mathsfbr{X}_{k}=\mathsfbr{X}_{k}^{\text{S}}\cup\mathsfbr{X}_{k}^{\text{B}}=\Bigg{(}\bigcup_{\bm{x}_{k-1}\in\hskip
0.85358pt{\cal{X}}_{k-1}}\\!\\!\\!\mathsfbr{S}_{k}\big{(}\bm{x}_{k-1}\big{)}\Bigg{)}\cup\mathsfbr{X}_{k}^{\text{B}}\hskip
0.85358pt.$

This model defines the state-transition pdf
$f({\cal{X}}_{k}|{\cal{X}}_{k-1})$, which can be calculated explicitly via set
convolution [4].

###  III-B General Measurement Model

We consider a general superpositional111The considered superpositional
intensity model typically applies to sonar and radar tracking applications.
However, it can be easily altered to an occlusion model used in image tracking
applications. intensity model [17] where the influence of object
$\bm{\mathsfbr{x}}_{k}^{(i)}$, $i\in\\{1,\ldots,n\\}$ to the intensity
measurement for pixel $\mathsfbr{z}_{k}^{(m)}$, $m\in\\{1,\ldots,M\\}$ is
described by an arbitrary point spread function (PSF)
$d^{(m)}(\bm{x}_{k}^{(i)})$. The total influence of the multiobject state
$\mathsfbr{X}_{k}$ on measurement $\mathsfbr{z}_{k}^{(m)}$,
$m\in\\{1,\ldots,M\\}$ is modeled as the sum of all individual PSFs, i.e.,

$D^{(m)}({\cal{X}}_{k})\hskip 0.85358pt=\hskip 0.85358pt\sum_{i\hskip
0.85358pt=\hskip 0.85358pt1}^{n}\hskip
2.84526ptd^{(m)}(\bm{x}_{k}^{(i)})\hskip 0.85358pt.$ (7)

This total influence $D^{(m)}({\cal{X}}_{k})$ is a sufficient statistic with
respect to measurement $\mathsfbr{z}_{k}^{(m)}$, $m\in\\{1,\ldots,M\\}$.

For example, the Swerling 1 model [15] often used in radar applications,
models all contributions of objects and background noise as statistically
independent and circularly symmetric Gaussian random variables. The
measurement $\mathsfbr{z}_{k}^{(m)}$, $m\in\\{1,\ldots,M\\}$ is the magnitude
of the sum of all contributions. Following the Swerling 1 model, the general
likelihood function of measurement $\mathsfbr{z}_{k}^{(m)}$,
$m\in\\{1,\ldots,M\\}$ is thus given by the Rayleigh distribution

$f_{\mathrm{g}}\big{(}z^{(m)}_{k}|{\cal{X}}_{k}\big{)}={\cal{R}}\Big{(}z_{k}^{(m)};\sqrt{D^{(m)}({\cal{X}}_{k})+\sigma^{2}_{\text{n}}}\hskip
0.85358pt\Big{)}\vspace{0mm}$

where $\sigma^{2}_{\text{n}}$ is the variance of the background noise.

Conditioned on ${\cal{X}}_{k}$, all measurements are assumed statistically
independent. The general joint likelihood function can thus be obtained as

$\displaystyle f_{\mathrm{g}}(\bm{z}_{k}|{\cal{X}}_{k})\hskip 0.85358pt$
$\displaystyle=\hskip 0.85358pt\prod_{m\hskip 0.85358pt=\hskip
0.85358pt1}^{M}f_{\mathrm{g}}\big{(}z^{(m)}_{k}|{\cal{X}}_{k}\big{)}\hskip
0.85358pt.$

Note that according to this general intensity model, every object can
potentially contribute to every measurement.

###  III-C Considered Simplified Measurement Model

As an approximation for simplified inference, we consider a model where (i)
every object contributes to exactly one distinct measurement
$\mathsfbr{z}_{k}^{(m)}$, $m\in\\{1,\ldots,M\\}$, and (ii) the probability
that the object with state $\bm{\mathsfbr{x}}_{k}^{(i)}$ contributes to
measurement $\mathsfbr{z}_{k}^{(m)}$ is proportional to the value of the PSF
$d^{(m)}(\bm{x}_{k}^{(i)})$. This simplified model is suitable for scenarios
where the PSF of each object is highly concentrated around a single
measurement, and objects are unlikely to fall in the same pixel.

To describe the associations of objects and measurements at time $k$, we
introduce the vector
$\bm{\mathsfbr{a}}_{k}\triangleq[\mathsfbr{a}_{k}^{(1)},\ldots,\mathsfbr{a}_{k}^{(n)}]^{\text{T}}$
with elements $\mathsfbr{a}_{k}^{(i)}\in\\{1,\ldots,M\\}$. Here,
$\mathsfbr{a}_{k}^{(i)}=m\in\\{1,\ldots,M\\}$ indicates that the object with
state $\bm{\mathsfbr{x}}_{k}^{(i)}$, $i\in\\{1,\ldots,n\\}$ only contributes
to measurement $\mathsfbr{z}_{k}^{(m)}$. An object-measurement association
vector $\bm{a}_{k}$ is admissible if each measurement involves contributions
of at most one object. All admissible associations form the association
alphabet $\mathcal{A}_{n,M}$. For a fixed association vector
$\bm{a}_{k}\in\mathcal{A}_{n,M}$, the total influence of the multiobject state
$\mathsfbr{X}_{k}$ on measurement $\mathsfbr{z}_{k}^{(m)}$,
$m\in\\{1,\ldots,M\\}$ is thus given by (cf. (7))

$\tilde{D}^{(m)}({\cal{X}}_{k})\hskip 0.85358pt=\hskip
0.85358pt\begin{cases}d^{(m)}(\bm{x}_{k}^{(i)})\hskip 0.85358pt,&\exists
i\hskip 2.84526pt\text{such that}\hskip 2.84526pta_{k}^{(i)}=m\hskip
0.85358pt,\\\\[3.41432pt] 0\hskip 0.85358pt,&\text{otherwise}\hskip
0.85358pt.\end{cases}$

Hence, if an object $\bm{\mathsfbr{x}}_{k}^{(i)}$ is associated to measurement
$a_{k}^{(i)}=m$, then the total influence of the multiobject state
$\mathsfbr{X}_{k}$ on measurement $a_{k}^{(i)}=m$ is
$d^{(m)}(\bm{x}_{k}^{(i)})$. If no object is associated to measurement $m$,
then there is no influence of the multiobject state $\mathsfbr{X}_{k}$ on
measurement $m$.

As we show in Appendix A, the joint likelihood function of $\mathsfbr{z}_{k}$
given ${\cal{X}}_{k}=\\{\bm{x}_{k}^{(1)},\ldots,\bm{x}_{k}^{(n)}\\}$ is given
by

$\displaystyle f(\bm{z}_{k}|{\cal{X}}_{k})$
$\displaystyle\propto\sum_{\bm{a}_{k}\in\hskip
0.85358pt\mathcal{A}_{n,M}}\bigg{(}\prod_{i\hskip 0.85358pt=\hskip
0.85358pt1}^{n}\hskip 0.85358pt\hskip
0.85358ptd^{(a_{k}^{(i)})}(\bm{x}_{k}^{(i)})f_{1}\big{(}z_{k}^{(a_{k}^{(i)})}\big{|}\bm{x}_{k}^{(i)}\big{)}\bigg{)}$
$\displaystyle\hskip 42.67912pt\times\prod_{m\hskip 0.85358pt\in\hskip
0.85358pt\mathcal{M}_{\bm{a}_{k}}}f_{0}\big{(}z_{k}^{(m)}\big{)}$ (8)

for $n\leqslant M$ and according to $f(\bm{z}_{k}|{\cal{X}}_{k})=0$ for $n>M$,
respectively. Here,
$\mathcal{M}_{\bm{a}_{k}}=\big{\\{}1,\ldots,M\big{\\}}\setminus\big{\\{}a_{k}^{(1)},\dots,a_{k}^{(n)}\big{\\}}$
consists of the indexes of all measurements that are not associated to any
object state $\bm{\mathsfbr{x}}_{k}^{(i)}$. In addition, $f_{0}(z_{k}^{(m)})$
and $f_{1}(z_{k}^{(m)}|\bm{x}_{k}^{(i)})$ denote the pdfs of measurement
$z_{k}^{(m)},m\in\\{1,\ldots,M\\}$ for the cases where it is associated to
none or one object $i\in\\{1,\ldots,n\\}$, respectively. For example,
following the Swerling 1 model, the pdfs $f_{1}(z_{k}^{(m)}|\bm{x}_{k}^{(i)})$
and $f_{0}(z_{k}^{(m)})$ are given by

$\displaystyle f_{1}\big{(}z_{k}^{(m)}\big{|}\bm{x}_{k}^{(i)}\big{)}$
$\displaystyle\triangleq{\cal{R}}\Big{(}z^{(m)}_{k};\sqrt{d^{(m)}(\bm{x}^{(i)}_{k})+\sigma^{2}_{\text{n}}}\hskip
0.85358pt\Big{)}$ (9) $\displaystyle f_{0}\big{(}z_{k}^{(m)}\big{)}$
$\displaystyle\triangleq{\cal{R}}\big{(}z^{(m)}_{k};\sigma_{\text{n}}\big{)}\hskip
0.85358pt.$ (10)

##  IV Exact Estimation

In the Bayesian sequential estimation framework, the statistics of the state
$\mathsfbr{X}_{k}$ at time $k$, conditioned on all received measurements
$\bm{z}_{1:k}\triangleq[\bm{z}_{1}\ldots\bm{z}_{k}]^{\text{T}}$ up to time
$k$, can be described by the posterior pdf $f({\cal{X}}_{k}|\bm{z}_{1:k})$.
This pdf is calculated from the previous posterior pdf
$f({\cal{X}}_{k-1}|\bm{z}_{1:k-1})$ via a prediction and an update step. The
prediction step calculates the predicted posterior pdf
$f({\cal{X}}_{k}|\bm{z}_{1:k-1})$ from the previous posterior pdf
$f({\cal{X}}_{k-1}|\bm{z}_{1:k-1})$ based on the state-transition pdf
$f({\cal{X}}_{k}|{\cal{X}}_{k-1})$ discussed in Section III-A. It is identical
to the one performed by the TOMB/P filter [5, 6, 7, 8] and will be skipped.
Note that $f({\cal{X}}_{k}|\bm{z}_{1:k-1})$ is a Poisson/MB pdf with $J_{k-1}$
Bernoulli components and intensity function $\lambda_{k|k-1}(\bm{x}_{k})$.

The update step calculates the current posterior pdf
$f({\cal{X}}_{k}|\bm{z}_{1:k})$ from the predicted posterior pdf
$f({\cal{X}}_{k}|\bm{z}_{1:k-1})$ based on the likelihood function
$f(\bm{z}_{k}|{\cal{X}}_{k})$ in (8). As we show in Appendix B, the updated
posterior pdf $f({\cal{X}}_{k}|\bm{z}_{1:k})$ is no longer a Poisson/MB pdf
but a MB mixture pdf (cf. (6)). Contrary to the conventional TOMB/P filter,
the updated posterior pdf contains no Poisson part anymore. This is a direct
consequence of the fact that the TBD measurement model (cf. Section III-C)
does not consider missed detections. However, this is not restrictive for
practical scenarios, since objects with very low intensities are also admitted
by our model (and thus difficult to distinguish from pure noise). More
precisely, as derived in Appendix B, the exact updated posterior pdf can be
expressed as

$\displaystyle f({\cal{X}}_{k}|\bm{z}_{1:k})$ $\displaystyle\hskip
2.84526pt=\sum_{\bm{a}^{\prime}_{k}\hskip 0.85358pt\in\hskip
0.85358pt\mathcal{A}^{\prime}_{J_{k},M}}\hskip
0.85358ptp(\bm{a}^{\prime}_{k})\hskip 0.85358pt\hskip
0.85358ptf^{\text{MB}}_{\bm{a}^{\prime}_{k}}({\cal{X}}_{k})\hskip 0.85358pt,$
$\displaystyle\hskip
2.84526pt=\sum_{{\cal{X}}_{k}^{(1)}\uplus\ldots\uplus{\cal{X}}_{k}^{(J_{k})}={\cal{X}}_{k}}\hskip
0.85358pt\sum_{\bm{a}^{\prime}_{k}\hskip 0.85358pt\in\hskip
0.85358pt\mathcal{A}^{\prime}_{J_{k},M}}\hskip
0.85358ptp(\bm{a}^{\prime}_{k})\hskip 0.85358pt\prod_{j\hskip 0.85358pt=\hskip
0.85358pt1}^{J_{k}}\hskip
0.85358ptf^{(j,a_{k}^{\prime(j)})}({\cal{X}}_{k}^{(j)})$ (11)

where each MB mixture component (i) corresponds to one admissible object-
measurement association
$\bm{\mathsfbr{a}}^{\prime}_{k}\triangleq[\mathsfbr{a}_{k}^{\prime(1)}\ldots\mathsfbr{a}_{k}^{\prime(J_{k})}]^{\text{T}}$,
(ii) has $J_{k}=J_{k-1}+M$ Bernoulli components, and (iii) is weighted by the
probability $p(\bm{a}^{\prime}_{k})$. Note that each Bernoulli component
represents an object that potentially exists and that
$\bm{\mathsfbr{a}}^{\prime}_{k}$ associates each potential object
$j\in\\{1,\ldots,J_{k}\\}$ to a measurement $m\in\\{1,\ldots,M\\}$.

In contrast to the association vector $\bm{\mathsfbr{a}}_{k}$ used for the
derivation of the likelihood function in (8), here
$\bm{\mathsfbr{a}}^{\prime}_{k}$ accounts also for the possible non-existence
of objects. More precisely, $\bm{\mathsfbr{a}}^{\prime}_{k}$ has entries
$\mathsfbr{a}_{k}^{\prime(j)}\in\\{0,1,\ldots,M\\}$ for
$j\in\\{1,\ldots,J_{k-1}\\}$ and entries
$\mathsfbr{a}_{k}^{\prime(j)}\in\\{0,1\\}$ for
$j\in\\{J_{k-1}+1,\ldots,J_{k-1}+M\\}$. Here, for
$j\in\\{1,\ldots,J_{k-1}\\}$, $\mathsfbr{a}_{k}^{\prime(j)}=\hskip 0.85358pt0$
indicates that an object with state $\bm{\mathsfbr{x}}_{k}^{(j)}$ does not
exist and $\mathsfbr{a}_{k}^{\prime(j)}=m\in\\{1,\ldots,M\\}$ indicates that
it does exist and contribute to measurement $\mathsfbr{z}_{k}^{(m)}$.
Furthermore, for $j\in\\{J_{k-1}+1,\ldots,J_{k-1}+M\\}$,
$\mathsfbr{a}_{k}^{\prime(j)}=1$ indicates that a new object with state
$\bm{\mathsfbr{x}}_{k}^{(j)}$ contributes to measurement $m=j-J_{k-1}$ and
$\mathsfbr{a}_{k}^{\prime(j)}=0$ that no new object contributes to measurement
$m=j-J_{k-1}$. All admissible association vectors $\bm{a}^{\prime}_{k}$ form
the set $\mathcal{A}^{\prime}_{J_{k},M}$. The association pmf
$p(\bm{a}^{\prime}_{k})$ is given by

$p(\bm{a}^{\prime}_{k})\propto\prod_{j\hskip 0.85358pt=\hskip
0.85358pt1}^{J_{k}}\hskip
0.85358pt\beta_{k}^{(j,a_{k}^{\prime(j)})}\vspace{.8mm}$ (12)

for $\bm{a}^{\prime}_{k}\in\mathcal{A}^{\prime}_{J_{k},M}$ and by
$p(\bm{a}^{\prime}_{k})=0$ for
$\bm{a}^{\prime}_{k}\notin\mathcal{A}^{\prime}_{J_{k},M}$. In the following,
we provide expressions for the existence probabilities $r_{k}^{(j,m)}$, the
spatial pdfs $f^{(j,m)}(\bm{x}_{k})$, and the association weights
$\beta_{k}^{(j,m)}$ used in (11) and (12).

For each potential object $j\in\\{1,\ldots,J_{k-1}\\}$ and measurement
$m\in\\{1,\ldots,M\\}$, we have

$\displaystyle r_{k}^{(j,m)}$ $\displaystyle=1$ (13) $\displaystyle
f^{(j,m)}(\bm{x}_{k})$ $\displaystyle=\frac{d^{(m)}(\bm{x}_{k})\hskip
0.85358ptf_{1}(z_{k}^{(m)}|\bm{x}_{k})\hskip
0.85358ptf_{k|k-1}^{(j)}(\bm{x}_{k})}{c_{k}^{(j,m)}}$ (14)
$\displaystyle\beta_{k}^{(j,m)}$ $\displaystyle=r_{k|k-1}^{(j)}\hskip
0.85358ptc_{k}^{(j,m)}$ (15)

where we introduced the normalization constant $c_{k}^{(j,m)}\triangleq\int
d^{(m)}(\bm{x}_{k})\hskip 0.85358ptf_{1}(z_{k}^{(m)}|\bm{x}_{k})\hskip
0.85358ptf_{k|k-1}^{(j)}(\bm{x}_{k})\hskip 0.85358pt\text{d}\bm{x}_{k}\hskip
0.85358pt$. Here, (13) indicates that the object modeled by Bernoulli
component $j$ exists and contributes to measurement $m$. The state of this
object is distributed according to (14) and the likelihood of this event is
characterized by (15). Furthermore, for $m=0$ we have $r_{k}^{(j,0)}=0$,
$f^{(j,0)}(\bm{x}_{k})$ not defined, and

$\beta_{k}^{(j,0)}=1-r_{k|k-1}^{(j)}$ (16)

Here, $r_{k}^{(j,0)}=0$ indicates that the object modeled by Bernoulli
component $j$ does not exist. The likelihood of this event is given by (16).

For each new potential object $j=J_{k-1}+m$, $m\in\\{1,\dots,$ $M\\}$, we get

$\displaystyle r_{k}^{(j,1)}$
$\displaystyle=\frac{c_{k}^{(j)}}{f_{0}(z_{k}^{(m)})+c_{k}^{(j)}}$ (17)
$\displaystyle f^{(j,1)}(\bm{x}_{k})$
$\displaystyle=\frac{d^{(m)}(\bm{x}_{k})\hskip
0.85358ptf_{1}(z_{k}^{(m)}|\bm{x}_{k})\hskip
0.85358pt\lambda_{k|k-1}(\bm{x}_{k})}{c_{k}^{(j)}}$ (18)
$\displaystyle\beta_{k}^{(j,1)}$
$\displaystyle=f_{0}(z_{k}^{(m)})+c_{k}^{(j)}$ (19)

with normalization constant $c_{k}^{(j)}\triangleq\int
d^{(m)}(\bm{x}_{k})f_{1}(z_{k}^{(m)}|\bm{x}_{k})$
$\lambda_{k|k-1}(\bm{x}_{k})\hskip 0.85358pt\text{d}\bm{x}_{k}$. Here, (17) is
the probability that there is a new object that contributes to measurement $m$
(assuming that no other existing object contributed to measurement $m$). The
state of this object is distributed according to (18) and the likelihood of
this event is given by (19). Finally, we have $\beta_{k}^{(j,0)}=1$,
$f^{(j,0)}(\bm{x}_{k})$ not defined, and $r_{k}^{(j,0)}=0$.

##  V Approximate Update Step

In the following, we will describe the approximate update step used by our TBD
method. This update step approximates the exact MB mixture posterior pdf in
(11) by a Poisson/MB pdf. This limits computational complexity and, thus,
enables real-time tracking of multiple low-observable objects.

###  V-A MB Approximation

First, we approximate the exact MB mixture posterior pdf in (11) by an MB pdf.
This approximation is based on approximating the association pmf
$p(\bm{a}^{\prime}_{k})$ in (12) by the product of its marginals. We first
extend the association alphabet $\mathcal{A}^{\prime}_{J_{k},M}$ in (12) to
$\bar{\mathcal{A}}^{\prime}_{J_{k},M}\triangleq\\{0,1,\ldots,M\\}^{J_{k-1}}\times\\{0,1\\}^{M}$.
Note that $\bar{\mathcal{A}}^{\prime}_{J_{k},M}$ now also contains
inadmissible associations. This does not affect $p(\bm{a}^{\prime}_{k})$ since
by definition $p(\bm{a}^{\prime}_{k})=0$ for
$\bar{\mathcal{A}}^{\prime}_{J_{k},M}\setminus\mathcal{A}^{\prime}_{J_{k},M}$.
Next, we approximate $p(\bm{a}^{\prime}_{k})$ according to

$p(\bm{a}^{\prime}_{k})\approx\prod_{j=1}^{J_{k}}\hskip
0.85358ptp\big{(}a_{k}^{\prime(j)}\big{)}\hskip
0.85358pt,\quad\bm{a}^{\prime}_{k}\in\bar{\mathcal{A}}^{\prime}_{J_{k},M}$
(20)

with

$p\big{(}a_{k}^{\prime(j)}\big{)}=\sum_{\sim\bm{a}_{k}^{\prime(j)}}\hskip
0.85358ptp(\bm{a}^{\prime}_{k})\hskip 0.85358pt.\vspace{-1mm}$

Here, $\sim\bm{a}_{k}^{\prime(j)}$ denotes the vector of all
$a_{k}^{\prime(j^{\prime})}$ with $j^{\prime}\in\\{1,\ldots,$
$J_{k}\\}\setminus j$. Note that a fast and scalable calculation of
$p\big{(}a_{k}^{\prime(j)}\big{)}$ is enabled by the SPA [18, 5, 6].

Next, we insert (20) into (11), which yields the approximate posterior pdf
$\tilde{f}({\cal{X}}_{k}|\bm{z}_{1:k})$ given by

$\displaystyle\tilde{f}({\cal{X}}_{k}|\bm{z}_{1:k})$
$\displaystyle=\sum_{{\cal{X}}_{k}^{(1)}\uplus\ldots\uplus\hskip
0.85358pt{\cal{X}}_{k}^{(J_{k})}=\hskip 0.85358pt{\cal{X}}_{k}}\hskip
0.85358pt\sum_{\bm{a}^{\prime}_{k}\hskip 0.85358pt\in\hskip
0.85358pt\bar{\mathcal{A}}^{\prime}_{J_{k},M}}\hskip 0.85358pt\prod_{j\hskip
0.85358pt=\hskip 0.85358pt1}^{J_{k}}\hskip 0.85358pt\hskip
0.85358ptp\big{(}a_{k}^{\prime(j)}\big{)}\hskip
0.85358ptf^{(j,a_{k}^{(j)})}({\cal{X}}_{k}^{(j)})\hskip 0.85358pt.$

By further using the identity
$\sum_{\bm{a}^{\prime}_{k}\in\bar{\mathcal{A}}^{\prime}_{J_{k},M}}\\!\prod_{j\hskip
0.85358pt=\hskip 0.85358pt1}^{J_{k}}p(a_{k}^{\prime(j)})=\big{(}\prod_{j\hskip
0.85358pt=\hskip 0.85358pt1}^{J_{k-1}}\sum_{a_{k}^{\prime(j)}=\hskip
0.85358pt0}^{M}p(a_{k}^{\prime(j)})\big{)}\prod_{j\hskip 0.85358pt=\hskip
0.85358ptJ_{k-1}+1}^{J_{k}}\sum_{a_{k}^{\prime(j)}=\hskip
0.85358pt0}^{1}p(a_{k}^{\prime(j)})$, we can now rewrite (LABEL:eq:approx2) as

$\displaystyle\tilde{f}({\cal{X}}_{k}|\bm{z}_{1:k})\hskip 0.85358pt$
$\displaystyle=\sum_{{\cal{X}}_{k}^{(1)}\uplus\ldots\uplus\hskip
0.85358pt{\cal{X}}_{k}^{(J_{k})}=\hskip 0.85358pt{\cal{X}}_{k}}\prod_{j\hskip
0.85358pt=\hskip 0.85358pt1}^{J_{k-1}}\sum_{a_{k}^{\prime(j)}=\hskip
0.85358pt0}^{M}\hskip 0.85358ptp(a_{k}^{\prime(j)})\hskip
0.85358ptf^{(j,a_{k}^{(j)})}({\cal{X}}_{k}^{(j)})$
$\displaystyle\times\prod_{j\hskip 0.85358pt=\hskip
0.85358ptJ_{k-1}+1}^{J_{k}}\hskip 0.85358pt\sum_{a_{k}^{\prime(j)}=\hskip
0.85358pt0}^{1}\hskip 0.85358ptp(a_{k}^{\prime(j)})\hskip
0.85358ptf^{(j,a_{k}^{(j)})}({\cal{X}}_{k}^{(j)}).$ (22)

By comparing (22) with (4), it can be seen that this approximate posterior pdf
$\tilde{f}({\cal{X}}_{k}|\bm{z}_{1:k})$ is now a MB pdf, i.e., it can be
rewritten as

$\tilde{f}({\cal{X}}_{k}|\bm{z}_{1:k})\hskip 0.85358pt\hskip
0.85358pt=\sum_{{\cal{X}}_{k}^{(1)}\uplus\ldots\uplus\hskip
0.85358pt{\cal{X}}_{k}^{(J_{k})}=\hskip 0.85358pt{\cal{X}}_{k}}\prod_{j\hskip
0.85358pt=\hskip 0.85358pt1}^{J_{k}}\hskip
0.85358ptf^{(j)}({\cal{X}}_{k}^{(j)})$ (23)

where the existence probabilities and spatial pdfs of the Bernoulli pdfs
$f^{(j)}({\cal{X}}_{k}^{(j)})$ for $j\in\\{1,\ldots,J_{k-1}\\}$ are given by

$\displaystyle r_{k}^{(j)}$ $\displaystyle=\sum_{a_{k}^{\prime(j)}=\hskip
0.85358pt1}^{M}p(a_{k}^{\prime(j)})$ (24) $\displaystyle f^{(j)}(\bm{x}_{k})$
$\displaystyle=\frac{1}{r_{k}^{(j)}}\sum_{a_{k}^{\prime(j)}=\hskip
0.85358pt1}^{M}p(a_{k}^{\prime(j)})\hskip
0.85358ptf^{(j,a_{k}^{\prime(j)})}(\bm{x}_{k})$ (25)

and for $j\in\\{J_{k-1}+1,\ldots,J_{k}\\}$ by

$\displaystyle r_{k}^{(j)}$ $\displaystyle=p(a_{k}^{\prime(j)}=1)\hskip
2.84526ptr_{k}^{(j,1)}$ (26) $\displaystyle f^{(j)}(\bm{x}_{k})$
$\displaystyle=\hskip 0.85358ptf^{(j,1)}(\bm{x}_{k})\hskip 0.85358pt.$ (27)

Note that the idea of approximating a MB mixture pdf by an MB pdf was also
used in the derivation of the TOMB/P filter [5] and the SPA-based labeled MB
filter [20].

###  V-B Recycling of MB Components

The approximate posterior pdf in (23) consists of $J_{k}=J_{k-1}+M$ Bernoulli
components, i.e., the number of Bernoulli components increases by $M$ at each
time step $k$. This is because newborn objects, modeled by the intensity
function $\lambda_{\text{B}}(\bm{x}_{k})$, may appear in the scene and
potentially contribute to all of the $M$ measurements. However, many Bernoulli
components typically have a very low existence probability and are thus
unlikely to represent an existing object.

In order for real-time tracking to remain feasible, the number of Bernoulli
components has to be limited. Contrary to most multiobject tracking methods
(see e.g. [1, 2, 3, 4, 6, 7, 8] ), rather than pruning components with an
existence probability below a fixed threshold and potentially discarding
valuable tracking information, we employ the concept of recycling [19]. Here,
Bernoulli components
$j\in\mathcal{J}_{k}^{\text{R}}\subseteq\\{1,\ldots,J_{k}\\}$ with an
existence probability $r_{k}^{(j)}$ below a predefined threshold
$\eta_{\text{R}}$ are “transferred” to the Poisson part of the posterior pdf
by means of moment matching [19]. This yields the approximated posterior PHD

$\lambda(\bm{x}_{k})=\sum_{j\in\mathcal{J}_{k}^{\text{R}}}r_{k}^{(j)}\hskip
0.85358ptf^{(j)}(\bm{x}_{k})\vspace{-1.2mm}$

where $r_{k}^{(j)}$ is given by (24) or (26) and $f^{(j)}(\bm{x}_{k})$ by (25)
or (27), respectively. After applying this recycling step, the approximate
posterior pdf is again a Poisson/MB pdf. In particular, the Poisson part is
represented by the approximate posterior PHD $\lambda(\bm{x}_{k})$ in (V-B)
and the MB part by the existence probabilities $r_{k}^{(j)}$ and spatial pdfs
$f^{(j)}(\bm{x}_{k})$,
$j\in\\{1,\ldots,J_{k}\\}\setminus\mathcal{J}_{k}^{\text{R}}$ in (24)–(27).

##  VI Numerical Study

We consider a two-dimensional (2D) simulation scenario with a region of
interest (ROI) of
$[\text{0}\text{m},\text{64}\text{m}]\times[\text{0}\text{m},\text{64}\text{m}]$.
We simulated 10 objects during 200 time steps. The object states consist of 2D
position, 2D velocity, and the object’s intensity, i.e.,
$\bm{\mathsfbr{x}}_{k}\\!=[\bm{\mathsfbr{p}}_{k}\,\,\sgamma_{k}]^{\text{T}}$
with
$\bm{\mathsfbr{p}}_{k}=[\mathsfbr{p}_{k,1}\,\,\mathsfbr{p}_{k,2}\,\,\dot{\mathsfbr{p}}_{k,1}\,\,\dot{\mathsfbr{p}}_{k,2}]^{\text{T}}$.
The kinematic part of the object’s state $\bm{\mathsfbr{p}}_{k}$ evolves
according to the nearly constant velocity motion model [21] with iid driving
noise distributed according to
$\mathcal{N}(\bm{\epsilon}_{k,\text{p}};\bm{0}_{4},10^{-3}\bm{I}_{4})$, and
the object’s intensity according to a random walk model with iid driving noise
distributed according to $\mathcal{N}(\epsilon_{k,\text{I}};0,10^{-4})$,
respectively. The objects appear at various times before time step 30 and at
randomly chosen positions in the area
$[\text{17}\text{m},\text{48}\text{m}]\times[\text{17}\text{m},\text{48}\text{m}]$
at each simulation run, and they disappear at various times after time step
170 or when they leave the ROI. The object’s initial velocity is drawn from
$\mathcal{N}(\bm{v}_{k};0,\sigma^{2}_{\text{v}}\bm{I}_{2})$ with variance
$\sigma^{2}_{\text{v}}=10^{-2}$. We consider two scenarios. All objects appear
with an initial intensity of $\gamma_{\text{I}}=10$ and $\gamma_{\text{I}}=4$
in scenarios 1 and 2, respectively. While the scenario is not deliberately
constructed to cause objects to come into close proximity, this occurs
randomly, and the behavior of the compared methods is observed.

$1$$20$$40$$60$$1$$20$$40$$60$$p_{1}$$p_{2}$

Figure 1: Example of true trajectories for $\gamma_{\text{I}}=10$ (white
lines), as well as of trajectories estimated by the proposed T-TOMB/P filter
(red lines) and measurements acquired at time $k\\!=\\!100$.

The measurement is an image consisting of $\text{64}\times\text{64}$ cells
covering the ROI. Each cell is a square of 1m side length and has a scalar
intensity. We use the measurement model described in Section III-B and set the
PSF to

$d^{(a_{k}^{(i)})}(\bm{x}_{k}^{(i)})=\begin{cases}\gamma_{k}^{(i)},&\text{object
}\bm{x}_{k}^{(i)}\text{ is in cell }a_{k}^{(i)}\\\\[2.84526pt]
0,&\text{otherwise}\end{cases}\vspace{-2mm}$

Hence, $\mathsfbr{z}_{k}^{(m)}$ is distributed according to
$\mathcal{R}\big{(}z^{(m)}_{k};\sqrt{\gamma^{(i)}_{k}+\sigma^{2}_{\text{n}}}\big{)}$
(cf. (9)) if object $\bm{\mathsfbr{x}}_{k}^{(i)}$ is in cell $m$ and according
to $\mathcal{R}\big{(}z^{(m)}_{k};\sigma_{\text{n}}\big{)}$ (cf. (10)) if no
object is in cell $m$. Here, $\gamma^{(i)}_{k}$ is the intensity of object $i$
at time $k$ and $\sigma_{\text{n}}=1$ is the standard deviation of the
background noise. If two or more objects
$i\in\mathcal{I}^{\prime}_{k}\subseteq\\{1,\ldots,I_{k}\\}$ are in the same
cell at the same time, we select $i$ by drawing a sample from
$\mathsfbr{i}^{\prime}\sim\gamma_{k}^{(i)}$, $i\in\mathcal{I}^{\prime}_{k}$.
Hence, $\mathsfbr{z}_{k}^{(m)}$ can only have contributions of at most one
object (cf. Section III-B). Thus, the higher the intensity value
$\gamma_{k}^{(i)}$ of object $\bm{\mathsfbr{x}}_{k}^{(i)}$, the more likely it
is to contribute to $\mathsfbr{z}_{k}^{(m)}$.

$40$$80$$120$$160$$200$$0$$5$$10$$15$$20$0(a)$k$MOSPA errorT-MBT-TOMB/P
(proposed)

$40$$80$$120$$160$$200$$0$$5$$10$$15$$20$0(b)$k$MOSPA errorT-MBT-TOMB/P
(proposed)

Figure 2: MOSPA error of T-TOMB/P and T-MB versus time $k$ for (a)
$\gamma_{\text{I}}=10$ and (b) $\gamma_{\text{I}}=4$.

We employ a particle implementation of our proposed TBD algorithm, briefly
termed T-TOMB/P due to the close relation to the conventional TOMB/P filter
[5] for point measurements. We compare T-TOMB/P to the TBD-based MB filter
proposed in [16], briefly referred to as T-MB. Both T-TOMB/P and T-MB
represent the spatial pdf of each Bernoulli component by $3,000$ particles.
T-TOMB/P further represents the posterior PHD by $50,000$ particles, the birth
PHD by another $50,000$, where the resulting $100,000$ particles are again
reduced to $50,000$ after the update step. More precisely, the birth PHD
$\lambda_{\text{B}}(\bm{x}_{k})=\mu_{\text{B}}\hskip
0.85358ptf_{\text{B}}(\bm{x}_{k})$ is represented by drawing particles from
the pdf $f_{\text{B}}(\bm{x}_{k})=f(p_{k,1}\,\,p_{k,2})\hskip
0.85358ptf_{\text{v}}\big{(}\dot{p}_{k,1},\dot{p}_{k,2}\big{)}\hskip
0.85358ptf_{\text{I}}(\gamma_{k})$ and by setting $\mu_{\text{B}}=4/64^{2}$.
Here, $f(p_{k,1}\,\,p_{k,2})$ is uniform over the ROI,
$f_{\text{v}}\big{(}\dot{p}_{k,1},\dot{p}_{k,2}\big{)}$ is
$\mathcal{N}(\dot{p}_{k,1},\dot{p}_{k,2};\bm{0}_{2},\sigma^{2}_{\text{v}}\hskip
0.85358pt\textbf{I}_{2})$, and $f_{\text{I}}(\gamma_{k})$ is uniform from $0$
to $\eta_{\text{I}}=30$, respectively. For the generation of new Bernoulli
components of T-MB, we adapt the scheme of [17]. In fact, T-MB generates a new
Bernoulli component for each measurement $z_{k-1}^{(m)}$ whose intensity value
is above the threshold $\eta_{\text{new}}=1.5\hskip
0.85358pt\sqrt{\gamma_{\text{I}}+\sigma^{2}_{\text{n}}}$. More precisely, the
existence probability of each new Bernoulli component is set to $10^{-4}$ and
the spatial pdf is represented by drawing particles from the pdf
$f^{\prime}_{\text{B}}(\bm{x}_{k})\propto\int f(\bm{x}_{k}|\bm{x}_{k-1})\hskip
0.85358ptf\big{(}z^{(m)}_{k-1}\hskip 0.85358pt\big{|}\hskip
0.85358ptp_{k-1,1}\,\,p_{k-1,2}\big{)}\allowbreak\hskip
0.85358ptf_{\text{v}}\big{(}\dot{p}_{k-1,1},\dot{p}_{k-1,2}\big{)}f_{\text{I}}(\gamma_{k-1})\hskip
0.85358pt\text{d}\bm{x}_{k-1}$. Here, the function
$f\big{(}z_{k-1}^{(m)}\big{|}p_{k-1,1}\,\,p_{k-1,2}\big{)}$ is uniform over
the cell area of measurement $m$ and
$f_{\text{v}}\big{(}\dot{p}_{k-1,1},\dot{p}_{k-1,2}\big{)}$ and
$f_{\text{I}}(\gamma_{k-1})$ are chosen as in $f_{\text{B}}(\bm{x}_{k})$.
Furthermore, T-TOMB/P recycles Bernoulli components with an existence
probability below $\eta_{\text{R}}=10^{-1}$ and T-MB prunes Bernoulli
components with an existence probability below $\eta_{\text{T}}=10^{-4}$,
respectively. T-TOMB/P and T-MB set $p_{\text{S}}=0.999\hskip 0.85358pt$.

The example shown in Fig. 1 suggests excellent detection and estimation
performance of the proposed method. For a quantitative assessment and
comparison of the performance of both filters, we computed the mean Euclidean
distance based optimal subpattern assignment (MOSPA) metric [22] with cutoff
parameter $c=20$, order $p=\\!2$, and averaged over 1000 simulation Monte
Carlo runs. Fig. 2 shows the obtained results of T-TOMB/P and T-MB for
$\gamma_{\text{I}}=10$ and $\gamma_{\text{I}}=4$. It is seen that T-TOMB/P
consistently outperforms T-MB. This can be attributed to the excellent
behavior of the Bernoulli component generation/recycling scheme of T-TOMB/P
and the fact that T-TOMB/P considers data association, i.e., a measurement is
allowed to be associated to at most one object and an object to exactly one
measurement, respectively. The performance of T-MB depends greatly on the
choice of $\eta_{\text{new}}$. A smaller value of $\eta_{\text{new}}$ results
in a faster detection of newly appearing objects, but in a higher number of
false tracks, i.e., Bernoulli components modeling not existing objects, on the
other hand. Another cause for the performance difference can be understood by
considering the case where two objects come close to each other. Since T-MB
allows a measurement to be associated to more than one object, one of the two
objects can be tracked by both tracks and the other object by none after
object separation. This is a direct consequence of the fact that T-MB does not
consider data association.

##  VII Conclusion

We proposed a scalable track-before-detect (TBD) method for the tracking of
low-observable objects that relies on a Poisson/multi-Bernoulli model. To
limit the computational complexity, we approximated the exact posterior pdf by
a multi-Bernoulli pdf. For a reliable tracking in real time, a recycling of
Bernoulli components is performed. We demonstrated that with the proposed
method, a significant improvement in tracking performance can be achieved
compared to a state-of-the-art TBD tracking method. A possible direction for
future research is an extension of the proposed method to a more general
measurement model by using sum-product algorithms for data association with
extended objects [23, 24].

## Appendix A

In the following, we derive (8). Based on assumptions (i) and (ii) in Section
III-C, the joint likelihood function for
${\cal{X}}_{k}=\\{\bm{x}_{k}^{(1)},\ldots,\bm{x}_{k}^{(n)}\\}$ can be found as
[25]

$f(\bm{z}_{k}|{\cal{X}}_{k})\propto\sum_{{\cal{X}}^{(1)}\uplus\ldots\uplus{\cal{X}}^{(M)}=\hskip
0.85358pt{\cal{X}}}g(z_{k}^{(m)},{\cal{X}}_{k}^{(m)})$

where we have introduced

$g(z_{k}^{(m)},{\cal{X}}_{k}^{(i)})=\begin{cases}f_{0}\big{(}z_{k}^{(m)}\big{)},&{\cal{X}}_{k}^{(i)}=\emptyset\\\\[4.2679pt]
d^{(m)}(\bm{x}_{k}^{(i)})\hskip
0.85358ptf_{1}\big{(}z_{k}^{(m)}\big{|}\bm{x}_{k}^{(i)}\big{)},&{\cal{X}}_{k}^{(i)}=\\{\bm{x}_{k}^{(i)}\\}\\\\[4.2679pt]
\end{cases}\vspace{0.5mm}$ (28)

and $\sum_{{\cal{X}}^{(1)}\uplus\ldots\uplus{\cal{X}}^{(M)}={\cal{X}}}$
denotes the sum over all disjoint decompositions of ${\cal{X}}$ into sets
${\cal{X}}^{(m)}$, $m\in\\{1,\ldots,M\\}$ such that
${\cal{X}}^{(1)}\cup\ldots\cup{\cal{X}}^{(M)}={\cal{X}}$.

By further inserting (28), we get

$\displaystyle f(\bm{z}_{k}|{\cal{X}}_{k})$
$\displaystyle\propto\sum_{{\cal{X}}^{(1)}\uplus\ldots\uplus{\cal{X}}^{(M)}=\hskip
0.85358pt{\cal{X}}}\bigg{(}\prod_{m:\hskip 0.85358pt{\cal{X}}_{k}^{(m)}=\hskip
0.85358pt\\{\bm{x}_{k}^{(i)}\\}}d^{(m)}(\bm{x}_{k}^{(i)})\hskip
0.85358ptf_{1}\big{(}z_{k}^{(m)}\big{|}\bm{x}_{k}^{(i)}\big{)}\bigg{)}$
$\displaystyle\times\prod_{m^{\prime}:\hskip
0.85358pt{\cal{X}}_{k}^{(m^{\prime})}=\hskip
0.85358pt\emptyset}f_{0}\big{(}z_{k}^{(m^{\prime})}\big{)}\hskip 0.85358pt.$
(29)

Expression (29) can be reformulated by using the object-measurement
association vector $\bm{\mathsfbr{a}}_{k}$ introduced in Section III-C, which
yields the final expression for the joint likelihood function in (8).

## Appendix B

In this appendix, we derive the exact posterior pdf
$f({\cal{X}}_{k}|\bm{z}_{1:k})$ in (11). First, we obtain the posterior pdf by
applying Bayes theorem [4] according to

$f({\cal{X}}_{k}|\bm{z}_{1:k})\propto f(\bm{z}_{k}|{\cal{X}}_{k})\hskip
0.85358ptf({\cal{X}}_{k}|\bm{z}_{1:k-1})\hskip 0.85358pt.\vspace{1mm}$ (30)

Next, we reformulate the likelihood function in (8) as follows. We define the
measurement-object association vector
$\bm{\mathsfbr{b}}_{k}\triangleq[\mathsfbr{b}_{k}^{(1)},\ldots,\mathsfbr{b}_{k}^{(M)}]$
with entries $\mathsfbr{b}_{k}^{(m)}\in\\{0,1,\ldots,n\\}$, which expresses
the same information as $\bm{\mathsfbr{a}}_{k}$ but in a different form [6,
18]. In particular, $\mathsfbr{b}_{k}^{(m)}=i\in\\{1,\ldots,n\\}$ indicates
that object $i$ contributes to measurement $m$ and $\mathsfbr{b}_{k}=0$
indicates that no object contributes to measurement $m$. An association
$\bm{b}_{k}$ is admissible if any object contributes to exactly one
measurement $m$, and at most one object contributes to each measurement. All
admissible measurement-object associations form the set $\mathcal{B}_{M,n}$.

The combination of both $\bm{\mathsfbr{a}}_{k}$ and $\bm{\mathsfbr{b}}_{k}$ is
now used to reformulate (8) for each
${\cal{X}}_{k}=\\{\bm{x}_{k}^{(1)}\ldots,\bm{x}_{k}^{(n)}\\}$ according to

$\displaystyle f(\bm{z}_{k}|{\cal{X}}_{k})$
$\displaystyle\propto\sum_{{\cal{X}}_{k,0}\uplus{\cal{X}}_{k,1}={\cal{X}}_{k}}\sum_{\bm{a}_{k}\hskip
0.85358pt\in\hskip 0.85358pt\mathcal{A}_{n_{1},M}}\prod_{i\hskip
0.85358pt=\hskip 0.85358pt1}^{n_{1}}\hskip 0.85358pt\hskip
0.85358ptf_{1}\big{(}z_{k}^{(a_{k}^{(i)})}|\bm{x}_{k,1}^{(i)}\big{)}\hskip
0.85358pt\hskip 0.85358ptd^{(a_{k}^{(i)})}(\bm{x}_{k,1}^{(i)})$
$\displaystyle\times\hskip
1.42262pt\sum_{\bm{b}_{k}\in\mathcal{B}_{\mathcal{M}_{\bm{a}_{k}},n_{0}}}\hskip
1.70717pt\prod_{m\hskip 0.85358pt\in\hskip
0.85358pt\mathcal{M}_{\bm{a}_{k}}}g\big{(}z_{k}^{(m)},{\cal{X}}_{k,0}^{(b_{k}^{(m)})}\big{)}$
(31)

with $n\leqslant M$. Here, ${\cal{X}}_{k}$ is decomposed into arbitrary
subsets
${\cal{X}}_{k,0}=\\{\bm{x}_{k,0}^{(1)},\ldots,\bm{x}_{k,0}^{(n_{0})}\\}$ and
${\cal{X}}_{k,1}=\\{\bm{x}_{k,1}^{(1)},\ldots,\bm{x}_{k,1}^{(n_{1})}\\}$ with
${\cal{X}}_{k}={\cal{X}}_{k,0}\cup{\cal{X}}_{k,1}$ and $n=n_{0}+n_{1}$.
Further, $g\big{(}z_{k}^{(m)},{\cal{X}}_{k,0}^{(b_{k}^{(m)})}\big{)}$ is given
by (28) with ${\cal{X}}_{k,0}^{(b_{k}^{(m)})}=\emptyset$ for $b_{k}^{(m)}=0$
and ${\cal{X}}_{k,0}^{(b_{k}^{(m)})}=\\{\bm{x}_{k}^{(i)}\\}$ for
$b_{k}^{(m)}=i$ and we introduced
$\mathcal{B}_{\mathcal{M}_{\bm{a}_{k}},n_{0}}=\big{\\{}\bm{b}_{k}\in\mathcal{B}_{M,n_{0}}\hskip
0.85358pt|\hskip 0.85358ptb_{k}^{(m)}=0,\hskip 0.85358pt\forall
m\notin\mathcal{M}_{\bm{a}_{k}}\big{\\}}$. Recall that
$f(\bm{z}_{k}|{\cal{X}}_{k})$ in (31) is $0$ for $n>M$.

Next, we recap that the predicted posterior pdf is of Poisson/MB form and can
be represented by

$\displaystyle f({\cal{X}}_{k}|z_{1:k-1})$
$\displaystyle=\sum_{{\cal{X}}_{k,0}\uplus{\cal{X}}_{k}^{(1)}\uplus\ldots\uplus{\cal{X}}_{k}^{(J_{k-1})}\hskip
0.85358pt=\hskip 0.85358pt{\cal{X}}_{k}}f^{\text{P}}({\cal{X}}_{k,0})\hskip
0.85358pt\prod_{j\hskip 0.85358pt=\hskip
0.85358pt1}^{J_{k-1}}f^{(j)}({\cal{X}}_{k,1}^{(j)})$ (32)

where the Bernoulli pdfs $f^{(j)}({\cal{X}}_{k})$ are parametrized by the
existence probabilities $r_{k|k-1}^{(j)}$ and spatial pdfs
$f_{k|k-1}^{(j)}(\bm{x}_{k})$ and the Poisson pdf
$f^{\text{P}}({\cal{X}}_{k,0})$ is represented by the posterior PHD
$\lambda_{k|k-1}(\bm{x}_{k})$. By plugging (31) and (32) into (30) and
performing certain reformulations, we obtain [25]

$\displaystyle f({\cal{X}}_{k}|\bm{z}_{1:k})$ $\displaystyle\hskip
2.84526pt\propto\sum_{{\cal{X}}_{k,0}\uplus{\cal{X}}_{k,1}={\cal{X}}_{k}}\sum_{\bar{\bm{a}}_{k}\in\bar{\mathcal{A}}_{J_{k-1},M}}h_{1}({\cal{X}}_{k,1}^{(1)},\ldots,{\cal{X}}_{k,1}^{(J_{k-1})},\bar{\bm{a}}_{k})$
$\displaystyle\hskip 9.10486pt\times\hskip
0.85358pt\sum_{\mathcal{B}_{\mathcal{M}_{\bm{a}_{k}},n_{0}}}h_{0}({\cal{X}}_{k,0}^{(1)},\ldots,{\cal{X}}_{k,0}^{(M)},\bm{b}_{k},\bar{\bm{a}}_{k})\hskip
0.85358pt.$ (33)

Here the association vector
$\bar{\bm{a}}_{k}\triangleq[\bar{a}^{(1)}_{k},\ldots,\bar{a}^{(J_{k-1})}_{k}]$,
$j\in\\{1,\ldots,J_{k-1}\\}$ has entries
$\bar{a}^{(j)}_{k}\in\\{0,1,\ldots,M\\}$ where $\bar{a}^{(j)}_{k}=0$ indicates
that an object with state $\bm{\mathsfbr{x}}_{k}^{(j)}$ does not exist and
$\bar{\mathsfbr{a}}_{k}^{(j)}=m\in\\{1,\ldots,M\\}$ indicates that it does
exist and contribute to measurement $\mathsfbr{z}_{k}^{(m)}$. All admissible
association vectors $\bar{\bm{a}}_{k}$ form the set
$\bar{\mathcal{A}}_{J_{k},M}$.

In addition, the function $h_{1}({\cal{X}}_{k,1},\bar{\bm{a}}_{k})$ is given
by

$h_{1}({\cal{X}}_{k,1}^{(1)},\ldots,{\cal{X}}_{k,1}^{(J_{k-1})},\bar{\bm{a}}_{k})\triangleq\prod_{j\hskip
0.85358pt=\hskip
0.85358pt1}^{J_{k-1}}g_{1}\big{(}z_{k}^{(\bar{a}_{k}^{(j)})}|{\cal{X}}_{k,1}^{(j)}\big{)}\hskip
0.85358ptf^{(j)}\big{(}{\cal{X}}_{k,1}^{(j)}\big{)}.\vspace{0.5mm}$ (34)

Here, $f^{(j)}({\cal{X}}_{k,1}^{(j)})$ is the Bernoulli pdf in (32) and
$g_{1}\big{(}z_{k}^{(\bar{a}_{k}^{(j)})}|{\cal{X}}_{k,1}^{(j)}\big{)}$ is
obtained as

$\displaystyle
g_{1}\big{(}z_{k}^{(\bar{a}_{k}^{(j)})}|{\cal{X}}_{k,1}^{(j)}\big{)}$
$\displaystyle\hskip
2.84526pt=\begin{cases}d^{(m)}\big{(}\bm{x}_{k,1}^{(i)})\hskip
0.85358ptf_{1}(z_{k}^{(m)}|\bm{x}_{k,1}^{(i)}\big{)}\hskip
0.85358pt,&\bar{a}_{k}^{(j)}=m,\hskip
0.85358pt{\cal{X}}_{k,1}^{(j)}=\\{\bm{x}_{k,1}^{(i)}\\}\\\\[2.84526pt] 1\hskip
0.85358pt,&\bar{a}_{k}^{(j)}=0,\hskip
0.85358pt{\cal{X}}_{k,1}^{(j)}=\emptyset\\\\[2.84526pt] 0\hskip
0.85358pt,&\text{otherwise.}\end{cases}\vspace{0.5mm}$

Expression (34) is a product of weighted Bernoulli pdfs with parameters as in
(13)–(16). Furthermore, the function
$h_{0}({\cal{X}}_{k,0},\bm{b}_{k},\bar{\bm{a}}_{k})$ reads

$\displaystyle
h_{0}({\cal{X}}_{k,0}^{(1)},\ldots,{\cal{X}}_{k,0}^{(M)},\bm{b}_{k},\bar{\bm{a}}_{k})$
$\displaystyle\hskip 42.67912pt\triangleq\prod_{m\hskip 0.85358pt\in\hskip
0.85358pt\mathcal{M}_{\bar{\bm{a}}_{k}}}g\big{(}z_{k}^{(m)},{\cal{X}}_{k,0}^{(b_{k}^{(m)})}\big{)}\hskip
0.85358ptf^{\text{P}}\big{(}{\cal{X}}_{k,0}^{(b_{k}^{(m)})}\big{)}$ (35)

where $g\big{(}z_{k}^{(m)},{\cal{X}}_{k,0}^{(b_{k}^{(m)})}\big{)}$ is given by
(28) with ${\cal{X}}_{k,0}^{(b_{k}^{(m)})}=\emptyset$ for $b_{k}^{(m)}=0$ and
${\cal{X}}_{k,0}^{(b_{k}^{(m)})}=\\{\bm{x}_{k}^{(i)}\\}$ for $b_{k}^{(m)}=i$,
and $f^{\text{P}}\big{(}{\cal{X}}_{k,0}^{(b_{k}^{(m)})}\big{)}$ is the Poisson
pdf in (32). Expression (35) is a product of weighted Bernoulli pdfs with the
parameters as in (18)–(19). By using the association vector
$\bm{\mathsfbr{a}}^{\prime}_{k}$ as defined in Section IV, we can now
reformulate (33) according to

$\displaystyle f({\cal{X}}_{k}|\bm{z}_{1:k})$ $\displaystyle\hskip
8.53581pt\propto\sum_{\uplus_{i\hskip 0.85358pt=\hskip
0.85358pt1}^{J_{k}}{\cal{X}}_{k}^{(i)}={\cal{X}}_{k}}\sum_{\bm{a}^{\prime}_{k}\in\mathcal{A}_{J_{k},M}}\prod_{j\hskip
0.85358pt=\hskip 0.85358pt1}^{J_{k}}\hskip
0.85358pt\beta_{k}^{(j,a_{k}^{\prime(j)})}\hskip
0.85358ptf^{(j,a_{k}^{\prime(j)})}({\cal{X}}_{k}^{(j)})\hskip 0.85358pt.$

Finally, we introduce the association pmf $p(\bm{a}^{\prime}_{k})$ (cf. (12))
and obtain the final expression for the posterior pdf in (11).

## References

  * [1] Y. Bar-Shalom, P. K. Willett, and X. Tian, _Tracking and Data Fusion: A Handbook of Algorithms_. Storrs, CT, USA: Yaakov Bar-Shalom, 2011.
  * [2] S. Challa, M. R. Morelande, D. Musicki, and R. Evans, _Fundamentals of Object Tracking_. Cambridge, UK: Cambridge University Press, 2011.
  * [3] W. Koch, _Tracking and Sensor Data Fusion: Methodological Framework and Selected Applications_. Berlin, Germany: Springer, 2014.
  * [4] R. P. S. Mahler, _Statistical Multisource-Multitarget Information Fusion_. Boston, MA, USA: Artech House, 2007.
  * [5] J. L. Williams, “Marginal multi-Bernoulli filters: RFS derivation of MHT, JIPDA, and association-based MeMBer,” _IEEE Trans. Aerosp. Electron. Syst._ , vol. 51, no. 3, pp. 1664–1687, Jul. 2015.
  * [6] F. Meyer, T. Kropfreiter, J. L. Williams, R. A. Lau, F. Hlawatsch, P. Braca, and M. Z. Win, “Message passing algorithms for scalable multitarget tracking,” _Proc. IEEE_ , vol. 106, no. 2, pp. 221–259, Feb. 2018.
  * [7] F. Meyer, P. Braca, P. Willett, and F. Hlawatsch, “A scalable algorithm for tracking an unknown number of targets using multiple sensors,” _IEEE Trans. Signal Process._ , vol. 65, no. 13, pp. 3478–3493, Jul. 2017.
  * [8] T. Kropfreiter, F. Meyer, and F. Hlawatsch, “Sequential Monte Carlo implementation of the track-oriented marginal multi-Bernoulli/Poisson filter,” in _Proc. FUSION-16_ , Heidelberg, Germany, Jul. 2016, pp. 972–979.
  * [9] T. Kropfreiter, F. Meyer, S. Coraluppi, C. Carthel, R. Mendrzik, and P. Willett, “Track coalescence and repulsion: MHT, JPDA, and BP,” in _Proc. FUSION-21_ , Sun City, South Africa, Nov. 2021.
  * [10] S. Tonissen and Y. Bar-Shalom, “Maximum likelihood track-before-detect with fluctuating target amplitude,” _IEEE Trans. Aerosp. Electron. Syst._ , vol. 34, no. 3, pp. 796–809, Jul. 1998.
  * [11] L. R. Moyer, J. Spak, and P. Lamanna, “A multi-dimensional Hough transform-based track-before-detect technique for detecting weak targets in strong clutter backgrounds,” _IEEE Trans. Aerosp. Electron. Syst._ , vol. 47, no. 4, pp. 3062–3068, Oct. 2011.
  * [12] Y. Barniv, “Dynamic programming solution for detecting dim moving targets,” _IEEE Trans. Aerosp. Electron. Syst._ , vol. AES-21, no. 1, pp. 144–156, Jan. 1985.
  * [13] S. J. Davey and H. X. Gaetjens, _Track-Before-Detect Using Expectation Maximisation_. Singapore: Springer, 2018\.
  * [14] B. Ristic, B.-T. Vo, B.-N. Vo, and A. Farina, “A Tutorial on Bernoulli Filters: Theory, Implementation and Applications,” _IEEE Trans. Signal Process._ , vol. 61, no. 13, pp. 3406–3430, Jul. 2013.
  * [15] B. Ristic, L. Rosenberg, D. Y. Kim, X. Wang, and J. Williams, “Bernoulli filter for track-before-detect: Swerling-1 target in K-distributed clutter,” in _Proc. RADAR-19_ , Toulon, France, Sep. 2019, pp. 1–5.
  * [16] B. Vo, B. Vo, N. Pham, and D. Suter, “Joint detection and estimation of multiple objects from image observations,” _IEEE Trans. Signal Process._ , vol. 58, no. 10, pp. 5129–5141, Oct. 2010.
  * [17] D. Y. Kim, B. Ristic, R. Guan, and L. Rosenberg, “A Bernoulli track-before-detect filter for interacting targets in maritime radar,” _IEEE Trans. Aerosp. Electron. Syst._ , pp. 1–10, Jan. 2021.
  * [18] J. L. Williams and R. Lau, “Approximate evaluation of marginal association probabilities with belief propagation,” _IEEE Trans. Aerosp. Electron. Syst._ , vol. 50, no. 4, pp. 2942–2959, Oct. 2014.
  * [19] J. L. Williams, “Hybrid Poisson and multi-Bernoulli filters,” in _Proc. FUSION-12_ , Singapore, Jul. 2012, pp. 1103–1110.
  * [20] T. Kropfreiter, F. Meyer, and F. Hlawatsch, “A fast labeled multi-Bernoulli filter using belief propagation,” _IEEE Trans. Aerosp. Electron. Syst._ , vol. 56, no. 3, pp. 2478–2488, Jun. 2020.
  * [21] Y. Bar-Shalom, X.-R. Li, and T. Kirubarajan, _Estimation with Applications to Tracking and Navigation_. New York, NY, USA: Wiley, 2002.
  * [22] D. Schuhmacher, B.-T. Vo, and B.-N. Vo, “A consistent metric for performance evaluation of multi-object filters,” _IEEE Trans. Signal Process._ , vol. 56, no. 8, pp. 3447–3457, Aug. 2008.
  * [23] F. Meyer and M. Z. Win, “Scalable data association for extended object tracking,” _IEEE Trans. Signal Inf. Process. Netw._ , vol. 6, pp. 491–507, May 2020.
  * [24] F. Meyer and J. L. Williams, “Scalable detection and tracking of geometric extended objects,” 2021, arXiv:2103.11279.
  * [25] A. F. Garcia-Fernandez, J. L. Williams, K. Granström, and L. Svensson, “Poisson multi-Bernoulli mixture filter: Direct derivation and implementation,” _IEEE Trans. Aerosp. Electron. Syst._ , vol. 54, no. 4, pp. 1883–1901, Feb. 2018.


# Soft-wall induced structure and dynamics of partially confined supercritical
fluids

Kanka Ghosh<EMAIL_ADDRESS>C.V.Krishnamurthy<EMAIL_ADDRESS>Department of Physics, Indian Institute of Technology Madras, Chennai-600036,
India

###### Abstract

The interplay between structure and dynamics of partially confined Lennard
Jones (LJ) fluids, deep into the supercritical phase, are studied over a wide
range of densities in the context of Frenkel line (FL), which separates rigid
liquidlike and non-rigid gaslike regimes in the phase diagram of the
supercritical fluids. Extensive molecular dynamics simulations carried out at
the two ends of the FL (P = $5000$ bar, T = $300$ K and T = $1500$ K) reveal
intriguing features in supercritical fluids as a function of stiffness of the
partially confining atomistic walls. The liquidlike regime of a LJ fluid (P =
$5000$ bar, T = $300$ K), mimicking Argon, partially confined between walls
separated by $10$ Å along the z-axis, and otherwise unconstrained, reveals
amorphous and liquidlike structural signatures in the radial distribution
function parallel to the walls and enhanced self-diffusion as the wall
stiffness is decreased. In sharp contrast, in the gas-like regime (P = $5000$
bar, T = $1500$ K), soft walls lead to increasing structural order hindering
self-diffusion. Further, the correlations between structure and self-diffusion
are found to be well captured by excess entropy. The rich behaviour shown by
supercritical fluids under partial confinement, even with simple interatomic
potentials, is found to be fairly independent of hydrophilicity and
hydrophobicity. The study identifies persisting sub-diffusive features over
intermediate time scales, emerging from the strong interplay between density
and confinement, to dictate the evolution and stabilization of structures. It
is anticipated that these results may help gain a better understanding of the
behaviour of partially confined complex fluids found in nature.

††preprint: AIP/123-QED

Walls that are soft and flexible, play a crucial role in modifying the
dynamics of the fluids and its flow behavior. The experimental evidence of the
occurrence of flow instability due to fluid flow past a soft surface has been
reported Kumaran and Muralikrishnan (2000). Linear stability analyses have
been performed for flow behavior through soft-gel coated walls Dinesh and
Pushpavanam (2018). Further, soft wall generated turbulence in the flow has
also been studied thoroughly Srinivas and Kumaran (2017). These results
suggest an interesting physics happening under the soft wall confinements.
Furthermore, MD simulations had been carried out for water-like core-softened
fluids under both fixed and flexible boundaries Krott and Bordin (2013).
Diffusion coefficients have been found to vary non-monotonically for core-
softened fluids under fixed walls, while monotonic trend is observed for
fluctuating walls Krott and Bordin (2013). The structural properties of core-
softened fluid are also changed while making walls fluctuating. The
thermodynamic behavior of water-like anomalous fluids, modeled using core-
softened potential, show a very different dependencies on the rigidity or the
flexibility of the walls Bordin, Krott, and Barbosa (2014).

In sufficiently thin spacings, where separation between the two confining
walls are of the order of few atomic diameters of the fluid particles,
interesting and anomalous behavior of fluid has been found to exist Mittal,
Errington, and Truskett (2006); Mittal  _et al._ (2008); Vanderlick and Davis
(1987); Schoen, Diestler, and Cushman (1987); Schoen  _et al._ (1988);
Pagonabarraga  _et al._ (1999); Kumar  _et al._ (2005); Varadarajan, Dasgupta,
and Ayappa (2018); Liu, Harder, and Berne (2004); J.J.Magda, M.Tirrell, and
H.T.Davis (1985). Besides, these model systems mimic a plethora of real
systems, implying a wide range of applicability from physics, chemistry and
biology Sievers  _et al._ (2016); Minton (2001); Patra  _et al._ (2014).
Typically, dense fluids under confinement act quite abnormally compared to its
bulk counterpart as wall-fluid interactions in a confined domain produce layer
like structural patterns, which in turn significantly modify the dynamics of
fluid. Even, at room temperature, extremely dilute gas under partial
confinement shows appreciable changes in dynamics with a pronounced
suppression of density fluctuations with respect to the bulk gaseous systems
Ghosh and Krishnamurthy (2018a).

Our recent studyGhosh and Krishnamurthy (2018b) on the effect of confinement
on the structural behavior of supercritical fluids indicates the possibility
of unusual dynamic behavior of supercritical fluids under confinement. Though
the dynamics of supercritical fluids across the Frenkel line (FL) in bulk
phase have been studied for last few years Brazhkin  _et al._ (2012, 2013);
Prescher  _et al._ (2017), the domain of confined fluids in supercritical
phase is vastly unexplored.

For the confinement, we employ a model for an atomistic wall that is aimed at
capturing the thermal motion of the wall particles and ensure that the wall
presents a ’rough’ surface that is dynamic. The particles are arranged in a
lattice as in a solid and interact with each other through a well-defined
inter-atomic potential. The motion of wall atoms is coupled to a thermostat of
Nose-Hoover type to maintain the same temperature as that of the supercritical
fluid to avoid unnecessary heat flow through the fluid.
$10^{5}$ atoms interacting with truncated and shifted Lennard Jones (LJ)
potential mimicking argon ($\epsilon$/$k_{B}$ = $120$ K, $\sigma$ = $3.4$ Å,
cutoff = $20$ Å$\sim$ $6\sigma$) have been simulated considering a set of
thermodynamics states in the supercritical regime (P = $5000$ bar,$240$ K
$\leqslant$ T $\leqslant$ $1500$ K) of the bulk phase. The FL crossover has
been determined through velocity autocorrelation function (VACF) and radial
distribution function (RDF) calculations to be at around T $\sim$ $600$-$700$
K. We explore the changes in the dynamic and structural features across the FL
through a series of MD simulations of supercritical fluid (SCF) under partial
confinement using atomistic walls. To determine the FL and consistency checks
in bulk, constant pressure temperature (NPT) ensembles are used while NVT
ensembles are used to simulate fluids under atomistic boundaries. Walls are
simulated along z, while periodic boundary conditions are applied along other
two directions (x, y). Throughout the study, a fixed value of $10$ Å ($\sim$
$3$ $\sigma$) has been chosen as a

Figure 1:  A snapshot of supercritical argon confined between (a) rigid
atomistic walls and (b) soft, flexible atomistic walls with a spacing H
($10$Å) between the walls at $300$ K. Grey colored atoms denote argon and the
blue atoms represent wall atoms. OVITO software is used to visualize the
snapshot Ovi ; Stukowski (2010).

representative spacing to understand the interplay between the density and
confinement on the dynamics and strong packing structures. After an energy
minimization, a standard velocity-verlet algorithm with a time step
($\Delta$t) of 0.0001 ps has been used to equilibrate the system up to $50$ ps
followed by a $1000$ ps production run to calculate and analyze the properties
of interests.

The solid, atomistic walls are made of three layers of the face-centered cubic
(fcc) lattice. The number of wall atoms are $583444$ for T = $300$ K and
$1076404$ for T = $1500$ K. Each of these wall atoms are attached to the
lattice sites by harmonic springs. For our study, we chose a fixed spacing (H)
of $10$ Å($\sim$ $3\sigma$) between the walls at both side of FL (T = $300$ K
and T = $1500$ K) and vary the spring constant (k) for these springs from a
higher value $5000$ ev/Å2 to a very low value of $0.005$ ev/Å2 to investigate
the effect of wall softness on the structure and dynamics of supercritical
fluid. Higher spring constants (k=$5000$ ev/Å2) mimic rigid walls, by
restricting the MSD of the wall atoms with respect to their lattice sites to a
very low value, while lower spring constants (k=$0.005$ ev/Å2) mimic soft
walls, by increasing the MSD of the wall atoms with respect to their lattice
sites. To define rigid wall, we allow the root mean squared displacement
(RMSD) of wall particles to be $40$ times and $100$ times smaller than a
typical distance traversed by a fluid particle between two collisions Brazhkin
_et al._ (2012) ($\sim$ $1$ Å), at $1500$ K and $300$ K respectively. Further
details on parameters used for wall-wall, wall-fluid interactions, and
complete computational exercises to locate FL can be found in our previous
work Ghosh and Krishnamurthy (2018b). Figure 1 shows a typical snapshot of the
systems with rigid and soft walls. All MD simulations have been carried out
using LAMMPS software package Plimpton (1995); Lam .

Figure 2 presents the mean squared displacement (MSD) components, parallel and
perpendicular to the walls, for the liquid-like (T = $300$ K) and gas-like (T
= $1500$ K) regimes of the partially confined SCF. The main results, shown for
rigid (k = $5000$ eV/Å2) and soft (k = $0.005$ eV/Å2) atomistic wall
conditions, feature a ballistic regime ($\sim$ t2) and a diffusive ($\sim$ t)
or sub-diffusive ($\sim$ tα; $\alpha$ = $0$) regime separated by an
intermediate regime displaying a variety of cross-over characteristics. Figure
2 includes the bulk MSD for purposes of comparison.

On short time scales, the fluid particles appear to be largely unaffected by
the confinement. While in the bulk, particles undergo thermal collisions
resulting in diffusive motion, under partial confinement, particles undergo
wall-mediated collisions in addition to thermal collisions resulting in
diffusive or sub-diffusive motion. These wall-mediated collisions, being
governed by the fluid density, temperature, spacing between the walls and the
stiffness of the walls, are seen to produce significant changes in the MSD
over intermediate and long-time scales.

The wall-mediated collisions are seen to constrain the z-component of the
particle motion completely in the gas-like regime at T = $1500$ K (Fig.2.(d))
over long time scales. The constraint arises due to frequent velocity
reversals produced by wall-mediated collisions reinforced by the strong
confinement (see Velocity autocorrelation function (VACF) in Fig S1.(c) in the
supplement where such velocity reversals produce minima in the VACF) that lead
to progressively smaller changes in the z-component of the particle position.
The effect of the wall stiffness manifests on intermediate time scales
influencing the details of the cross-over from ballistic to diffusive or sub-
diffusive motion. Onset of corrections to ballistic motion due to wall-
mediated collisions is earlier for soft wall than for rigid wall (green and
red curves in Fig 2.(d)) due to the larger amplitude of thermal motion of wall
particles in the former.

The motion parallel to the walls in the gas-like regime (T = $1500$ K) is
diffusive over long time scales. Since the particles of the soft wall have a
much larger root mean squared displacement than those of the rigid wall at T =
$1500$ K, the wall-mediated collisions from the soft wall are significantly
more non-specular in nature. Consequently, the in-plane scattering due to
wall-mediated collisions induce faster decay in velocity autocorrelation
function when the wall is soft than when the wall is rigid (see VACF in Fig
S1.(d) in the supplement). It may indeed be noted from Fig S1.(d) of the
supplement, that partial confinement induces faster decay than what is
possible through purely thermal collisions. What is more interesting

Figure 2:  Variation of MSD∥ (along x, y) and MSD⟂ (along z) as a function of
time in a double logarithmic scale. Fig. 2.(a), (b) represent phase space
point at $300$ K (liquidlike regime) for two extreme k (stiffness coefficient
of walls) values and Fig. 2. (c), (d) represent phase space point at $1500$ K
(gaslike regime) for same two extreme k values. For all the cases,
corresponding bulk values are shown in blue dotted lines.

is that as the wall stiffness decreases, the non-specular nature of the wall-
mediated collisions increases and results in the VACF∥ developing almost a
shallow minimum and vanishing much earlier. The diffusion coefficients
parallel to the walls in the gas-like and liquid-like regimes evaluated over a
range of wall stiffness are shown in Table 1.

As can be seen from Table 1, at $1500$ K, the diffusion coefficients parallel
to the walls (D∥) are reduced by about $20$ $\%$ and $70$ $\%$ from the bulk
when the wall is rigid and soft respectively. While in-plane scattering could
contribute to a lowering of the diffusion coefficient, it will be seen later
that an increased structural order also plays a role in the net

Table 1:  Variation of Diffusion coefficients parallel to the walls (D∥) with different wall-stiffnesses (k values) for a fixed confined spacing (H =$10$Å) at $300$ K and $1500$ K. For reference corresponding bulk diffusion coefficients are also shown in column 4 for these two phase space points. | k (ev/Å2) |  $D_{\parallel}$(Å2/ps) |  $D_{bulk}$(Å2/ps)
---|---|---|---
T = 300 K | 5000 | 0.076 |
(liquidlike regime) | 1000 | 0.077 |
| 10 | 0.129 | 0.74
H = 10 Å | 0.5 | 0.234 |
| 0.05 | 0.358 |
| 0.005 | 0.509 |
T = 1500 K | 5000 | 7.191 |
(gaslike regime) | 1000 | 7.118 |
| 10 | 6.489 | 9.117
H = 10 Å | 0.5 | 5.861 |
| 0.05 | 4.366 |
| 0.005 | 2.656 |

reduction of the diffusion coefficient.

In the liquidlike regime (T = $300$ K), we note surprising trends in both
parallel and perpendicular components of the motion particularly over
intermediate time scales. From Table 1, we note that bulk diffusion
coefficient is significantly lowered with respect to that in the gas-like
regime (T = $1500$ K) due to the higher fluid density. The velocity reversals,
along the z-axis, produced by wall-mediated collisions acting as it does at
liquid-like densities lead to caging effects over intermediate time scales and
vanishes soon after for both rigid and soft walls. It is important to note
that while the caging effect exists even in the bulk SCF due to the high
densities, it is more pronounced when the wall is rigid than when the wall is
soft (see VACF Fig S1.(a) in the supplement). This pronounced caging effect
seems to arise from a nearly specular nature of wall-mediated collisions
induced by the rigid wall and compounded by the strong confinement at such
high densities. The surprising aspect of this pronounced caging effect is that
it manifests as a sub-diffusive feature in the MSD⟂ resembling that of a
glassy phase and persists over a fairly wide intermediate time scale
suggestive of the characteristic slowing down encountered in the glass forming
liquids Varadarajan, Dasgupta, and Ayappa (2018). The sub-diffusive feature
exhibited by MSD⟂ for the soft wall confinement which develops more gradually
over time is distinctly different from that for the rigid wall although
persisting over almost as wide an intermediate time scale. Over long time
scales, however, the z-component motion is constrained for both rigid and soft
wall conditions.

Motion parallel to the walls exhibits surprising features as well. At T =
$300$ K, while the particles of the soft wall and the rigid wall would both
have relatively smaller root mean squared displacements than at T = $1500$ K,
the relatively smaller root mean squared displacement associated with the
rigid wall would induce wall-mediated collisions that are almost specular in
nature. With strong confinement compounding the nearly specular scattering,
the VACF∥ exhibits deeper minima as wall stiffness is increased at such high
densities (see VACF Fig S1.(b) in the supplement). The existence of a well-
defined minimum, followed by the vanishing of the VACF∥ soon after, manifests
as a sub-diffusive feature in the MSD∥ when the wall is rigid. From Table 1,
we note that the diffusion coefficient is nearly $90$ $\%$ lower than that in
the bulk, under rigid wall confinement, while the diffusion coefficient is
only $30$ $\%$ lower than that in the bulk, under soft wall confinement. Such
a large reduction of the diffusion coefficient appears possible only if in-
plane scattering produced by wall-mediated collisions combines with an
increased structural order imposing restrictions to particle motion.

The constraints on MSD⟂ at T = $300$ K as well as at T = $1500$ K, have
important structural implications – at long time, the fluid particles have
root mean squared displacements of $2.9$ Å and $3.7$ Å at T = $300$ K, and
$2.8$ Å and $3.3$ Å at T = $1500$ K, for the rigid and soft wall confinements
respectively along z. Although the separation between the walls is $10$ Å,

Figure 3: Distributions of argon particles in supercritical state, averaged
over several timesteps at $300$K (Fig 3.(a)) and at $1500$ K (Fig 3.(b)) are
calculated for two extreme k values, normal to the walls for $10$Å confined
spacing. For comparison, the bulk density distribution is shown.

these constraints on the root mean squared displacements of fluid particles
force them to adopt non-uniform but symmetric spatial distributions with
respect to the z-axis. Fig 3.(a) and (b) show the normalized particle counts
along the z-axis at T = $300$ K and T = $1500$ K for rigid and soft wall
confinements. At T = $300$ K, the higher fluid density imposes further packing
restrictions resulting in a finer layering across the z-axis. The scale over
which local number density varies is about $3$ Å consistent with the average
displacement of fluid particles in both the gas-like regime as well as the
liquid-like regime of the SCF. As mentioned in our previous work Ghosh and
Krishnamurthy (2018b), the number of layers formed due to confinement scales
linearly with the ratio of spacing to the atomic diameter ($\frac{H}{\sigma}$
= $10$/$3.4$ $\sim$ $3$) in the liquidlike regime of SCF.

It is intriguing to see how motion parallel to the walls influences the
lateral structure. This structure-dynamics interrelationship is examined using
the radial distribution function parallel to the walls, g∥(r), and the excess
entropy associated with it.

The radial distribution function parallel to the walls (g∥(r)) is computed for
specified intervals along the z-axis from

$\displaystyle g_{\parallel}(r)\equiv$
$\displaystyle\frac{1}{\rho^{2}V}\sum_{i\neq
j}\delta\left(r-r_{ij}\right)\Bigg{[}\theta\left(|z_{i}-z_{j}|+\frac{\delta
z}{2}\right)$ $\displaystyle-\theta\left(|z_{i}-z_{j}|-\frac{\delta
z}{2}\right)\Bigg{]}$ (1)

where $V$ is the volume, $\rho$ is the density, rij is the distance parallel
to the walls between molecules i and j, zi is the $z$ coordinate of the
molecule $i$, and $\delta(x)$ is the Dirac $\delta$ function. The Heaviside
function $\theta(x)$ restricts the sum to a pair of particles located in the
same slab of thickness $\delta z$. We have considered $\delta z$ to be same as
the width of each layer. We use a uniform bin width and bin number of $80$ to
calculate g∥(r) for all the cases. Figure 4 and Figure 5 show g∥(r) for soft
and rigid wall confinements at the two temperatures, $1500$ K and $300$ K
respectively. For purposes of comparison, the g(r) for the bulk phases are
also included.

The two-body excess entropy, parallel to the walls, is defined as Nettleton
and Green (1958); Truskett, Torquato, and Debenedetti (2000)

$s^{(2)}=-2\pi\rho
k_{B}\int[g_{\parallel}(r)lng_{\parallel}(r)-g_{\parallel}(r)+1]r^{2}dr$ (2)

Table 2 shows the two-body excess entropy evaluated as a function of wall
stiffness (k) for T = $300$ K and T = $1500$ K for the same spatial regions
considered in the computation of g∥(r). The bulk values of s(2) are negative
since the reference is an ideal gas. The values for the confined systems are
more negative due to the more ordered structures formed. As the density is
kept constant representing a particular phase space point for both bulk and
confined systems, the quantity $\frac{s^{(2)}}{2\pi\rho k_{B}}$ is used for
the study.

Figure 4:  Figure 4.(a) shows g∥(r) of the central region in the number
distribution of SCF (argon) particles at $1500$ K, with H = $10$Å, for two
extreme stiffness coefficients (k). Figure 4.(b) shows g∥(r) of the region
close to the left wall (wall at z = $-5$ Å) for two extreme k values. For
comparison bulk g(r) has also been shown.

Although Fig 4.(a) and (b) show that g∥(r) appears not to be influenced
significantly by the cross-sectional density variation under rigid wall
confinement at $1500$ K, compared to the bulk, the excess entropy values from
Table 2 clearly indicate that there is some degree of ordering in the central
region and even more ordering in the higher density region of the cross-
section (see Fig. 3.(b)). It may be recalled that at T = $1500$ K, under rigid
wall confinement, the perpendicular component of

Figure 5:  Figure 5 shows g∥(r) of the layer close to the left wall (wall at z
= $-5$ Å) in the number distribution of SCF (argon) particles at $300$ K, with
H = $10$Å, for two extreme stiffness coefficients (k). For comparison bulk
g(r) has also been shown.

the wall-mediated collisions is seen to (a) bring about deviations from the
ballistic regime a little earlier than that in the bulk, and (b) limit the
values of MSD⟂ soon after, within a short time. This rapid change in the
nature of motion is seen to produce cross-sectional density variations early
on. In other words, the long-time scale density variations shown in Fig 3.(b)
appear to set in much earlier. The increase in structural order, captured by
the lowering of the excess entropy, leads to restricted particle motion
contributing to the overall reduction of $20$ $\%$ in the diffusion
coefficient with respect to the bulk. The enhanced non-specular wall mediated
collisions along with an increased structural order (seen in Fig 4 and Table
2) parallel to the soft walls at $1500$ K, play a crucial role in the
decrement of the diffusion coefficient to $70$ $\%$ compared to the bulk.

At T = $300$ K, the persistence of sub-diffusive features over a wide
intermediate time scale seen in the MSD⟂ appears to be the key to the
formation of spatial structures and signal the mutual influence between
particle motion and structure as time evolves. The ‘slowing down’ of the
dynamics implies the emergence of spatially ordered structures which in-turn
restricts the dynamics and this in turn reinforces the ordering and so on till
MSD⟂ no longer changes. As can be seen from Fig 3.(a), the cross-sectional
density variations are stronger when the wall is rigid due to the combined
effect of high fluid density and strong confinement.

The onset of the sub-diffusive features can be seen in the parallel component
of MSD (MSD∥) as well when the wall is rigid. The large reduction in the
excess entropy indicates significant structural ordering and is indeed evident
from g∥(r) shown in Fig 5. The appearance of several peaks in g∥(r) and the
distorted features of these successive maxima and minima suggest that the
structure is more amorphous than liquid-like. The $90$ $\%$ reduction in the
diffusion coefficient appears to be largely due to the higher degree of
structural order.

When the wall is soft, the features in g∥(r) appear relatively smoother with
fewer peaks resembling a more liquid-like phase. Table 2 indicates that there
is considerable structural ordering even when the wall is soft and appears to
play a major role in reducing the diffusion coefficient by $30$ $\%$ with
respect to the bulk. This again seems to be an outcome of the combined effect
of high fluid density and strong confinement.

We note that all the results described in the article have been obtained with
the ratio of fluid-wall ($\epsilon_{w-f}$) and fluid-fluid

Table 2:  Variation of Pair-excess entropy (s(2)/$2\pi\rho k_{B}$), parallel to the walls, with different wall-stiffnesses (k values) for a fixed confined spacing (H =$10$Å) at $300$ K and $1500$ K. For reference corresponding bulk values are also shown in column 6 for these two phase space points. | k (ev/Å2) | s(2)/$2\pi\rho k_{B}$ | s(2)/$2\pi\rho k_{B}$ | s(2)/$2\pi\rho k_{B}$ | s(2)/$2\pi\rho k_{B}$
---|---|---|---|---|---
|  | (central region) | (Region near left wall) | (Region near right wall) | (Bulk)
T = 300 K | 5000 | -63.82 | -53.48 | -52.54 |
(liquidlike regime) | 1000 | -65.31 | -53.41 | -54.31 |
| 10 | -61.94 | -57.64 | -56.86 |
H = 10 Å | 0.5 | -43.11 | -58.32 | -56.29 | -9.39
| 0.05 | -20.16 | -23.11 | -23.03 |
| 0.005 | -13.68 | -16.64 | -17.75 |
T = 1500 K | 5000 | -4.86 | -7.71 | -6.70 |
(gaslike regime) | 1000 | -5.43 | -5.85 | -5.56 |
| 10 | -5.05 | -6.77 | -6.38 |
H = 10 Å | 0.5 | -6.11 | -5.26 | -6.98 | -3.27
| 0.05 | -4.67 | -6.35 | -6.13 |
| 0.005 | -7.62 | -11.27 | -11.50 |

($\epsilon_{f-f}$) interaction strengths
$\frac{\epsilon_{w-f}}{\epsilon_{f-f}}$ = $1$. A value $\geq$ $1$ for this
ratio represents hydrophilic boundary conditions Morciano  _et al._ (2017);
Voronov, Papavassiliou, and Lee (2006). Studies were undertaken explicitly
varying the ratios of fluid-wall and fluid-fluid interaction strengths,
mimicking both hydrophobic ($\frac{\epsilon_{w-f}}{\epsilon_{f-f}}$ $<$ $1$)
and hydrophilic ($\frac{\epsilon_{w-f}}{\epsilon_{f-f}}$ $\geq$ $1$) boundary
conditions. It is observed that the key findings regarding the structure and
dynamics remain unaltered, thus making the results presented here quite
generic and fairly independent of the affinity between the fluid and the wall
particles.

In summary, soon after the ballistic regime, the highly dense SCF at T = $300$
K exhibits persisting in-plane sub-diffusive features under rigid-wall
confinement enabling amorphous-like structural features to emerge. Over long-
time scales, the structure stabilizes as evidenced by the appearance of
multiple distorted peaks in the radial distribution function and the
associated large negative excess entropy denoting ordering. The structure
causes a nearly $90$ $\%$ reduction in the diffusion coefficient from the
bulk. Perhaps not surprisingly, soft-wall confinement leads to a more liquid-
like structure to emerge with only a $30$ $\%$ reduction in the diffusion
coefficient from the bulk.

At T = $1500$ K, the gas-like SCF, although still dense, exhibits mild in-
plane sub-diffusive features which are influenced significantly by the soft-
wall confinement in sharp contrast with the behaviour at T = $300$ K. The
diffusion coefficient under soft-wall confinement is reduced by nearly $70$
$\%$ from the bulk, while that for rigid-wall confinement is reduced only by
$20$ $\%$ from the bulk.

The dynamics and structure of SCF under rigid and soft walls at two phase
space points lying on either side of the FL, emerges as a unique system,
showcasing features of almost all types of fluids: (a) liquid-like packing
with prominent layering across the confined width before crossing FL, (b) gas-
like diffusion dominated motion after crossing FL, (c) glass-like slowdown in
dynamics at intermediate timescale and (d) amorphous-like structural features
for rigid wall confinement in the liquid-like regime.

##  I Supplementary Material

In the supplementary material, detailed descriptions of velocity
autocorrelation function (VACF) have been discussed for supercritical LJ fluid
under both rigid and soft walls at two ends of the FL (300K and 1500 K).

We acknowledge the help of HPCE, IIT Madras for high performance computing. KG
expresses his gratitude to Department of Science and Technology, Government of
India for providing INSPIRE fellowship.

## References

  * Kumaran and Muralikrishnan (2000) V. Kumaran and R. Muralikrishnan, Phys. Rev. Lett. 84, 3310 (2000).
  * Dinesh and Pushpavanam (2018) B. Dinesh and S. Pushpavanam, Phys. Fluids 30, 104103 (2018).
  * Srinivas and Kumaran (2017) S. S. Srinivas and V. Kumaran, J. Fluid Mech. 812, 1076–1118 (2017).
  * Krott and Bordin (2013) L. B. Krott and J. R. Bordin, J. Chem. Phys. 139, 154502 (2013).
  * Bordin, Krott, and Barbosa (2014) J. R. Bordin, L. B. Krott, and M. C. Barbosa, J. Phys. Chem. C 118, 9497 (2014).
  * Mittal, Errington, and Truskett (2006) J. Mittal, J. R. Errington, and T. M. Truskett, Phys. Rev. Lett. 96, 177804 (2006).
  * Mittal  _et al._ (2008) J. Mittal, T. M. Truskett, J. R. Errington, and G. Hummer, Phys. Rev. Lett. 100, 145901 (2008).
  * Vanderlick and Davis (1987) T. K. Vanderlick and H. T. Davis, J. Chem. Phys. 87, 1791 (1987).
  * Schoen, Diestler, and Cushman (1987) M. Schoen, D. J. Diestler, and J. H. Cushman, J. Chem. Phys. 87, 5464 (1987).
  * Schoen  _et al._ (1988) M. Schoen, J. H. Cushman, D. J. Diestler, and C. L. Rhykerd, J. Chem. Phys. 88, 1394 (1988).
  * Pagonabarraga  _et al._ (1999) I. Pagonabarraga, M. H. J. Hagen, C. P. Lowe, and D. Frenkel, Phys. Rev. E 59, 4458 (1999).
  * Kumar  _et al._ (2005) P. Kumar, S. V. Buldyrev, F. W. Starr, N. Giovambattista, and H. E. Stanley, Phys. Rev. E 72, 051503 (2005).
  * Varadarajan, Dasgupta, and Ayappa (2018) V. Varadarajan, C. Dasgupta, and K. G. Ayappa, J. Chem. Phys. 149, 064503 (2018).
  * Liu, Harder, and Berne (2004) P. Liu, E. Harder, and B. J. Berne, J. Phys. Chem. B 108, 6595 (2004).
  * J.J.Magda, M.Tirrell, and H.T.Davis (1985) J.J.Magda, M.Tirrell, and H.T.Davis, J. Chem. Phys. 83, 1888 (1985).
  * Sievers  _et al._ (2016) C. Sievers, Y. Noda, L. Qi, E. M. Albuquerque, R. M. Rioux, and S. L. Scott, ACS Cat. 6, 8286–8307 (2016).
  * Minton (2001) A. P. Minton, J. Biol. Chem. 276, 10577–10580 (2001).
  * Patra  _et al._ (2014) A. Patra, T. Q. Luong, R. K. Mitra, and M. Havenith, Phys.Chem.Chem.Phys. 16, 12875 (2014).
  * Ghosh and Krishnamurthy (2018a) K. Ghosh and C. V. Krishnamurthy, Phys. Rev. E 98, 052115 (2018a).
  * Ghosh and Krishnamurthy (2018b) K. Ghosh and C. V. Krishnamurthy, Phys. Rev. E 97, 012131 (2018b).
  * Brazhkin  _et al._ (2012) V. V. Brazhkin, Y. D. Fomin, A. G. Lyapin, V. N. Ryzhov, and K. Trachenko, Phys. Rev. E 85, 031203 (2012).
  * Brazhkin  _et al._ (2013) V. V. Brazhkin, Y. D. Fomin, A. G. Lyapin, V. N. Ryzhov, E. N. Tsiok, and K. Trachenko, Phys. Rev. Lett. 111, 145901 (2013).
  * Prescher  _et al._ (2017) C. Prescher, Y. D. Fomin, V. B. Prakapenka, J. Stefanski, K. Trachenko, and V. V. Brazhkin, Phys. Rev. B 95, 134114 (2017).
  * (24) “Ovito-open visualization tool,” http://ovito.org/.
  * Stukowski (2010) A. Stukowski, Modelling Simul. Mater. Sci. Eng. 18, 015012 (2010).
  * Plimpton (1995) S. Plimpton, J. Comput. Phys 117, 1 (1995).
  * (27) “Lammps software package,” http://lammps.sandia.gov.
  * Nettleton and Green (1958) R. E. Nettleton and M. S. Green, J. Chem. Phys. 29, 1365 (1958).
  * Truskett, Torquato, and Debenedetti (2000) T. M. Truskett, S. Torquato, and P. G. Debenedetti, Phys. Rev. E 62, 993 (2000).
  * Morciano  _et al._ (2017) M. Morciano, M. Fasano, A. Nold, C. Braga, P. Yatsyshin, D. N. Sibley, B. D. Goddard, E. Chiavazzo, P.Asinari, and S. Kalliadasis, J. Chem. Phys. 146, 244507 (2017).
  * Voronov, Papavassiliou, and Lee (2006) R. S. Voronov, D. V. Papavassiliou, and L. L. Lee, J. Chem. Phys. 124, 204701 (2006).