\section{Introduction}

One of the cornerstones of the celebrated general theory of relativity is undoubtedly represented by the equivalence principle (EP). According to Einstein, the (locally) indistinguishable equivalence between a gravitational field and an acceleration of the reference frame has been the main guiding light for the development of a relativistic treatment of gravitation. Ever since this archetypical notion was introduced, the equivalence principle has been thoroughly investigated from several standpoints, thus allowing for a plethora of distinct, brand-new interpretations. Indeed, when the equivalence principle is addressed nowadays there is a certain degree of ambiguity, as it may come in different formulations, each of which has its own requirements and range of validity. For a complete and detailed overview on this topic, the interested reader can consult Refs.~\cite{libe,new}.

For the purpose of the present paper, we restrict our attention to the weak equivalence principle, which establishes the equality between the inertial and the gravitational mass, henceforth labeled as $m_i$ and $m_g$, respectively. However, the aforementioned equality can be heavily influenced by the physical system and environment under consideration, as under specific circumstances the two masses may no longer have the same magnitude. For instance, when studying an electrically charged particle that interacts with a thermal bath of photons, radiative corrections affect $m_i$ and $m_g$ separately and unevenly, thereby giving rise to a violation of EP in its weak formulation~\cite{new,don,our}. The same occurrence can be found in the analysis of non-relativistic neutrinos, where the flavor mixing is responsible for a redefinition of the inertial mass, which thus becomes distinguishable from $m_g$~\cite{neutrino}. On a broader scale, the interplay between neutrino physics and EP violation is much deeper than the above feature, as the vast literature on the subject confirms~\cite{wlw,ms,wwg,sck,nmf,sah,ash,bha,ours,gn,hl}.

Here, we investigate a different scenario where EP violation is manifest. Specifically, we consider a relativistic quantum mechanical framework whose underlying Heisenberg uncertainty principle (HUP) is modified so as to accommodate the existence of a minimal momentum scale~\cite{kempf,kempf2,cbl,mz,bambi,oy,mignemi,symm}; a similar generalization goes by the name of extended uncertainty principle (EUP). The justification for the above deformation sinks its roots in cosmological arguments, since the curvature of spacetime at large distances is so extreme that there is a fundamental limit to the resolution with which momentum can be described~\cite{kempf,kempf2,mz}. As an additional piece of evidence for this observation, it is no coincidence that, in accordance with several findings~\cite{bambi,oy,mignemi}, a well-known formulation of the EUP is based upon the employment of the cosmological constant. Nevertheless, it has also been shown~\cite{schur,schur2,fab,fab2,fab3} that another version of EUP is deduced from a straightforward application of quantum mechanical laws in a general relativistic setting, thus letting the deviation from HUP depend exclusively on the intrinsic aspects of the spacetime under examination. In this work, we will analyze the most common setting, according to which the characteristic parameter associated to EUP must be fixed by requiring consistency with experimental data. In this regard, it is worth mentioning that the most stringent bounds on the EUP parameter are discussed in Refs.~\cite{mz,symm}.

In order to comply with the aforementioned purposes, the paper is structured as follows: in Sec.~II, the main features of the extended uncertainty principle are summarized, focusing in particular on the correction it entails at the level of the momentum operator. Section~III is devoted to the non-relativistic limit of the Dirac equation with the modified momentum, whilst the same study in conjunction with the Klein-Gordon equation is the topic of Sec.~IV. Finally, Sec.~V contains concluding remarks and future perspectives.

Throughout the whole manuscript, we use the natural units $\hbar=c=1$ but keep the speed of light explicit only in Sec.~IV, since the non-relativistic regime for a scalar field can be readily derived from a perturbative treatment of the Klein-Gordon equation in powers of $c$ (as seen in Ref.~\cite{kiefer}).

\section{Extended uncertainty principle}

In simple terms, the Heisenberg uncertainty principle states that arbitrarily small spatial distances can be resolved with sufficiently high-energy probes. To mathematically express this concept, we know that
\begin{equation}\label{hup}
\Delta x\Delta p\geq\frac{1}{2}\,.
\end{equation}
However, if a minimal momentum scale is accounted for, this is no longer the case; indeed, the previous equation must be suitably modified. Drawing inspiration from Refs.~\cite{kempf,kempf2,cbl,mz,bambi,oy,mignemi}, we argue that the influence of spacetime curvature at large distances induces a position-dependent correction into Eq.~\eqref{hup}, which then becomes
\begin{equation}\label{eup}
\Delta X\Delta P\geq\frac{1}{2}\left}\def\lrar{\leftrightarrow(1+\alpha\Delta X^2\right}\def\ti{\tilde)\,,
\end{equation}
where $\alpha$ is the so-called deformation parameter and we assume to work in the regime in which the dimensionless quantity $\alpha\Delta X^2$ can be regarded as a perturbation. Since the deformation parameter is the inverse of a squared length, the close relation between $\alpha$ and the cosmological constant $\Lambda$ can be easily inferred. As a matter of fact, it was shown in Ref.~\cite{park} that $\alpha=\pm\Lambda/3$, where the minus (plus) sign is associated to a (anti-)de Sitter spacetime. Therefore, in compliance with the experimental observations, we deduce that $\alpha<0$. Nonetheless, as already anticipated, we analyze the more general scenario where $\alpha$ may not be necessarily related to $\Lambda$, in line with Refs.~\cite{kempf,kempf2,mz,symm}.

At this point, from Eq.~\eqref{eup} it is straightforward to compute the value of the minimal uncertainty in momentum, which turns out to be $\Delta P_{\mathrm{min}}=\sqrt{|\alpha|}$. Furthermore, from the same equation, it is also possible to trace back to a general expression for the deformed commutator between the position and momentum operators. Indeed, starting from
\begin{equation}\label{comm}
\left}\def\lrar{\leftrightarrow[X,P\right}\def\ti{\tilde]=i\left}\def\lrar{\leftrightarrow(1+\alpha X^2\right}\def\ti{\tilde)\,,
\end{equation}
one can make use of the Robertson-Schr\"odinger prescription to deduce Eq.~\eqref{eup}.

The deformed canonical commutator plays a relevant role, as it provides an easy representation of $X$ and $P$ in terms of the auxiliary operators $x$ and $p$ for which the usual relation $\left}\def\lrar{\leftrightarrow[x,p\right}\def\ti{\tilde]=i$ holds. To see this, it is sufficient to replace the following identities
\begin{equation}\label{op}
X=x\,, \qquad P=\left}\def\lrar{\leftrightarrow(1+\alpha x^2\right}\def\ti{\tilde)p-i\alpha x\,,
\end{equation}
into Eq.~\eqref{comm} and verify that it is always satisfied. The second term in the definition of the physical momentum has been added to make the operator self-adjoint.

The identification~\eqref{op} is valid in one spatial dimension and as long as all the orders higher than and up to $\mathcal{O}(\alpha^2)$ are neglected. Should we perform the full expansion with more than one spatial dimension, we would have to deal with an emerging non-commutativity for different momentum components, because $\left}\def\lrar{\leftrightarrow[P_i,P_j\right}\def\ti{\tilde]\propto\mathcal{O}(\alpha^2)$. However, for the purposes of the present paper, we can safely restrict the attention to the lowest-order terms in the deformation parameter without harming the generality of our results. For further details on higher-order EUP, see Refs.~\cite{kempf,kempf2,mignemi}.

Before concluding this Section, we point out that calculations will be carried out in the position representation; hence, by generalizing Eq.~\eqref{op} to higher dimensions, we obtain
\begin{equation}\label{op2}
X_i=x_i\,, \qquad P_i=-i(1+\alpha x^2)\partial}\def\ran{\rangle_i-i\alpha x_i\,,
\end{equation}
with $x^2=\sum_ix_i^2$.

\section{Equivalence principle violation for spinor fields}

In order to approach the modified Dirac equation, one has to modify
\begin{equation}\label{deq}
\left}\def\lrar{\leftrightarrow(i\gamma^\mu\partial}\def\ran{\rangle_\mu-m\right}\def\ti{\tilde)\psi=0\,,
\end{equation}
by properly generalizing the auxiliary momentum operator to the physical one in the context of EUP. This is achieved by exploiting the relation~\eqref{op2} to cast the extended version of the Dirac equation as
\begin{equation}\label{mdeq}
i\gamma^0\partial}\def\ran{\rangle_t\psi=-i\left}\def\lrar{\leftrightarrow(1+\alpha x^2\right}\def\ti{\tilde)\vec{\gamma}\cdot\vec{\nabla}\psi-i\alpha\vec{\gamma}\cdot\vec{x}\,\psi+m\psi\,.
\end{equation}
A similar field equation has already been addressed in Ref.~\cite{gine}, where the authors discuss a completely unrelated issue.

The gamma matrices appearing in~\eqref{mdeq} can be written in the Dirac representation, namely
\begin{equation}\label{gamma}
\gamma^0=\begin{pmatrix} \textbf{1} & 0 \\ 0 & -\textbf{1} \end{pmatrix}\,, \qquad \vec{\gamma}=\begin{pmatrix} 0 & \vec{\si} \\ -\vec{\si} & 0 \end{pmatrix}\,,
\end{equation}
where $\textbf{1}$ denotes the $2\times2$ identity matrix and $\vec{\si}=\left}\def\lrar{\leftrightarrow(\si_x,\si_y,\si_z\right}\def\ti{\tilde)$ the Pauli matrices.

To provide a comprehensive analysis of the subject, one shall incorporate the effects of an external gravitational field in such a way to have access to both the inertial and the gravitational mass from the non-relativistic limit. Indeed, along the line of Refs.~\cite{don}, it is known that the two different masses can only be distinguished within the non-relativistic Schr\"odinger equation for the particle sector. This identification is achieved by looking at the kinetic operator to deduce $m_i$ and the interaction with the gravitational potential to derive $m_g$.

In order to study the dynamics of a spinor field in the presence of gravity, Eq.~\eqref{deq} must be modified by introducing the spin connection $\Gamma_\mu$ and vierbein fields $e^\mu_a$ as follows:
\begin{equation}\label{gdeq}
\Big[i\gamma^ae^\mu_a\left}\def\lrar{\leftrightarrow(\partial}\def\ran{\rangle_\mu+\Gamma_\mu\right}\def\ti{\tilde)-m\Big]\psi=0\,.
\end{equation}
The nature of $\Gamma_\mu$ and $e^\mu_a$ is closely intertwined with the shape of the metric tensor; for a detailed review, see Ref.~\cite{gravitation}. Here, since we seek the non-relativistic limit and an easy-to-interpret interaction with the gravitational field, we require that the metric tensor describes a weak-field approximation of the Schwarzschild solution in isotropic coordinates, i.e.,
\begin{equation}\label{metric}
ds^2=\left}\def\lrar{\leftrightarrow(1+2\phi\right}\def\ti{\tilde)dt^2-\left}\def\lrar{\leftrightarrow(1-2\phi\right}\def\ti{\tilde)d\textbf{r}^2\,,
\end{equation}
where the potential $\phi=-GM/r$ is assumed to be a small quantity. From Eq.~\eqref{metric}, one can promptly compute the spin connection and the vierbeins; a simple calculation (which can be found in Ref.~\cite{neutrino}) shows that
\begin{equation}\label{spinc}
\Gamma_\mu=\frac{1}{8}\left}\def\lrar{\leftrightarrow[\gamma^a,\gamma^b\right}\def\ti{\tilde]e^\lambda_a\left}\def\lrar{\leftrightarrow(\eta_{\mu\lambda}\partial}\def\ran{\rangle_\rho\phi-\eta_{\mu\rho}\partial}\def\ran{\rangle_\lambda\phi\right}\def\ti{\tilde)e^\rho_b\,,
\end{equation}
\begin{equation}\label{tetrad}
e^t_0=1-\phi\,, \qquad e^x_1=e^y_2=e^z_3=1+\phi\,.
\end{equation}
At this point, one can further require that the spatial gradients of the gravitational potential are negligible, so that $\partial}\def\ran{\rangle_i\phi\approx0$, $\forall i$. Although such a choice does not affect the validity of our results, it turns out to be useful, as it significantly streamlines the intermediate steps. As a matter of fact, by virtue of this assumption, the spin connection~\eqref{spinc} is identically vanishing, which means that the influence of gravity only enters the Dirac equation through the vierbein fields. In light of this, the extension of Eq.~\eqref{mdeq} which is inclusive of the gravitational interaction reads
\begin{equation}\label{mdeq2}
i\gamma^0\left}\def\lrar{\leftrightarrow(1-\phi\right}\def\ti{\tilde)\partial}\def\ran{\rangle_t\psi=-i\left}\def\lrar{\leftrightarrow(1+\alpha x^2+\phi\right}\def\ti{\tilde)\vec{\gamma}\cdot\vec{\nabla}\psi-i\alpha\vec{\gamma}\cdot\vec{x}\,\psi+m\psi\,.
\end{equation}
We can now proceed by separating the Dirac spinor into the upper and lower components
\begin{equation}\label{spinor}
\psi=\begin{pmatrix} \varphi \\ \chi \end{pmatrix}\,,
\end{equation}
so as to explicitly deal with the equations for $\varphi$ and $\chi$ independently. In so doing, one has
\begin{equation}\label{varphi}
i\left}\def\lrar{\leftrightarrow(1-\phi\right}\def\ti{\tilde)\partial}\def\ran{\rangle_t\varphi=-i\left}\def\lrar{\leftrightarrow(1+\alpha x^2+\phi\right}\def\ti{\tilde)\vec{\si}\cdot\vec{\nabla}\chi-i\alpha\vec{\si}\cdot\vec{x}\,\chi+m\varphi\,,
\end{equation}
\begin{equation}\label{chi}
-i\left}\def\lrar{\leftrightarrow(1-\phi\right}\def\ti{\tilde)\partial}\def\ran{\rangle_t\chi=i\left}\def\lrar{\leftrightarrow(1+\alpha x^2+\phi\right}\def\ti{\tilde)\vec{\si}\cdot\vec{\nabla}\varphi+i\alpha\vec{\si}\cdot\vec{x}\,\varphi+m\chi\,.
\end{equation}
In order to reach the non-relativistic regime, we first must note that, in this limit, the most important contribution to the energy comes from the rest mass. Therefore, if we single out a fast oscillating factor $\exp(-imt)$ from the Dirac spinor, we have that~\cite{neutrino,drell}
\begin{equation}\label{nrl}
\psi=e^{-imt}\tilde{\psi}\,,
\end{equation}
but since $\tilde{\psi}$ still depends on the time coordinate, it is immediate to realize $|i\partial}\def\ran{\rangle_t\tilde{\psi}|\ll|m\tilde{\psi}|$, which provides a non-relativistic approximation~\cite{drell}. In light of this result, it is possible to cast $\varphi$ and $\chi$ in terms of $\tilde{\varphi}$ and $\tilde{\chi}$ and enforce the above condition; the outcome of such a step yields
\begin{equation}\label{varphi2}
i\partial}\def\ran{\rangle_t\tilde{\varphi}=m\phi\,\tilde{\varphi}-i\left}\def\lrar{\leftrightarrow(1+\alpha x^2+\phi\right}\def\ti{\tilde)\vec{\si}\cdot\vec{\nabla}\tilde{\chi}-i\alpha\vec{\si}\cdot\vec{x}\,\tilde{\chi}\,,
\end{equation}
\begin{equation}\label{chi2}
\left}\def\lrar{\leftrightarrow(\phi-2\right}\def\ti{\tilde)m\tilde{\chi}=i\left}\def\lrar{\leftrightarrow(1+\alpha x^2+\phi\right}\def\ti{\tilde)\vec{\si}\cdot\vec{\nabla}\tilde{\varphi}+i\alpha\vec{\si}\cdot\vec{x}\,\tilde{\varphi}\,.
\end{equation}
The latter equation can be employed to find an explicit dependence of $\tilde{\chi}$ on $\tilde{\varphi}$. By recalling the weak-field approximation, we obtain
\begin{equation}\label{chi3}
\tilde{\chi}=-\frac{i}{2m}\left}\def\lrar{\leftrightarrow(1+\alpha x^2+\frac{3}{2}\phi\right}\def\ti{\tilde)\vec{\si}\cdot\vec{\nabla}\tilde{\varphi}-\frac{i\alpha}{2m}\vec{\si}\cdot\vec{x}\,\tilde{\varphi}\,.
\end{equation}
As it also occurs for the standard Dirac equation, in the non-relativistic regime the lower component $\chi$ is sub-dominant with respect to the upper component $\varphi$~\cite{drell}. Hence, we can focus on Eq.~\eqref{varphi2} only and derive the Schr\"odinger-like equation for the particle sector. To this aim, the expression~\eqref{chi3} must be substituted in \eqref{varphi2}, so as to get
\begin{equation}\label{find1}
i\partial}\def\ran{\rangle_t\tilde{\varphi}=-\left}\def\lrar{\leftrightarrow(1+2\alpha x^2\right}\def\ti{\tilde)\frac{\nabla^2\tilde{\varphi}}{2m}+m\phi\,\tilde{\varphi}-\frac{5}{4m}\phi\,\nabla^2\tilde{\varphi}-\frac{2\alpha}{m}\left}\def\lrar{\leftrightarrow(\vec{\si}\cdot\vec{x}\right}\def\ti{\tilde)\vec{\si}\cdot\vec{\nabla}\tilde{\varphi}-\frac{3\alpha}{2m}\tilde{\varphi}\,.
\end{equation}
A further manipulation of Eq.~\eqref{find1} allows us to reach the final expression, where the non-relativistic Schr\"odinger equation associated with $\tilde{\varphi}$ is manifest
\begin{equation}\label{find}
i\partial}\def\ran{\rangle_t\tilde{\varphi}=-\left}\def\lrar{\leftrightarrow(1+2\alpha x^2\right}\def\ti{\tilde)\frac{\nabla^2\tilde{\varphi}}{2m}+m\phi\,\tilde{\varphi}-\frac{5}{4m}\phi\,\nabla^2\tilde{\varphi}-\frac{2i\alpha}{m}\vec{\si}\cdot\left}\def\lrar{\leftrightarrow(\vec{x}\wedge\vec{\nabla}\tilde{\varphi}\right}\def\ti{\tilde)-\frac{2\alpha}{m}\vec{x}\cdot\vec{\nabla}\tilde{\varphi}-\frac{3\alpha}{2m}\tilde{\varphi}\,.
\end{equation}
The relevant factors to be considered are the first two appearing in the r.h.s., which represent the kinetic operator and the interaction with the external gravitational potential, respectively. The third term can be interpreted as a post-Newtonian correction, whilst the fourth one may be viewed as an EUP-induced spin-orbit interaction. The last two quantities are distinctive of the EUP quantum mechanical treatment, but the study of their implications lies beyond the purpose of the present paper.

Now, we need to extrapolate the shape of the inertial and gravitational mass connected with the spinor field; concerning $m_i$, its definition can be directly deduced from the kinetic operator, and turns out to be
\begin{equation}\label{mi}
m_i=m\left}\def\lrar{\leftrightarrow(1-2\alpha x^2\right}\def\ti{\tilde)\,.
\end{equation}
On the other hand, the gravitational mass is easily recognizable from the gravitational interaction, which tells us that $m_g=m$. Therefore, $m_i$ and $m_g$ do not undergo the same redefinition, which in turn signals a violation of the weak equivalence principle. This violation can be quantified by the ratio
\begin{equation}\label{ratio}
\frac{m_g}{m_i}=1+2\alpha x^2\,,
\end{equation}
which smoothly recovers the fulfillment of EP as long as $\alpha\to0$.

In the next Section, we exhibit that a similar feature emerges not only in the context of spinor fields, but also for scalar particles. Remarkably, the degree of EP violation will be quantified by the same function as in Eq.~\eqref{ratio}, thus suggesting that the present achievement does not depend on the nature of the system under investigation.

\section{Equivalence principle violation for scalar fields}

The starting point for the analysis of scalar fields is the Klein-Gordon equation
\begin{equation}\label{kg}
\left}\def\lrar{\leftrightarrow(\partial}\def\ran{\rangle_\mu\partial}\def\ran{\rangle^\mu+m^2c^2\right}\def\ti{\tilde)\psi=0\,,
\end{equation}
where the speed of light is kept explicit for later convenience. In the presence of gravity, the above equation is generalized as~\cite{gravitation}
\begin{equation}\label{gkg}
\frac{1}{\sqrt{-g}}\partial}\def\ran{\rangle_\mu\left}\def\lrar{\leftrightarrow(\sqrt{-g}g^{\mu\nu}\partial}\def\ran{\rangle_\nu\psi\right}\def\ti{\tilde)+m^2c^2\psi=0\,,
\end{equation}
where $g=\mathrm{det}(g_{\mu\nu})$. In order to establish a faithful comparison with the spinor case, $g_{\mu\nu}$ must be the one described by Eq.~\eqref{metric}. However, by restoring $c$ appropriately, the line element~\eqref{metric} now reads
\begin{equation}\label{metric2}
ds^2=\left}\def\lrar{\leftrightarrow(1+2\frac{\phi}{c^2}\right}\def\ti{\tilde)c^2dt^2-\left}\def\lrar{\leftrightarrow(1-2\frac{\phi}{c^2}\right}\def\ti{\tilde)d\textbf{r}^2\,,
\end{equation}
with $\phi$ left untouched.

By embedding the above physical setting in the framework of EUP and neglecting contributions that go like $\mathcal{O}(\alpha^2)$, Eq.~\eqref{gkg} becomes
\begin{equation}\label{gkg2}
\frac{1}{c^2}\left}\def\lrar{\leftrightarrow(1-2\frac{\phi}{c^2}\right}\def\ti{\tilde)\partial}\def\ran{\rangle^2_t\psi-\left}\def\lrar{\leftrightarrow(1+2\alpha x^2+2\frac{\phi}{c^2}\right}\def\ti{\tilde)\nabla^2\psi-4\alpha\vec{x}\cdot\vec{\nabla}\psi-3\alpha\psi+m^2c^2\psi=0\,.
\end{equation}
Equation~\eqref{gkg2} has been thoroughly investigated in Ref.~\cite{chung}, where the scalar particle couples to an external magnetic field rather than the gravitational one.

In order to consider the non-relativistic limit, we can perform an expansion in powers of $c$. In accordance with Ref.~\cite{kiefer}, a possible solution for the Klein-Gordon equation can be written as
\begin{equation}\label{sol}
\psi\left}\def\lrar{\leftrightarrow(\vec{x},t\right}\def\ti{\tilde)=e^{i\varphi\left}\def\lrar{\leftrightarrow(\vec{x},t\right}\def\ti{\tilde)}\,,
\end{equation}
where $\varphi$ is assumed to depend on $c$ as
\begin{equation}\label{exp}
\varphi=c^2\varphi_0+\varphi_1+c^{-2}\varphi_2+\dots\,.
\end{equation}
The series is infinite, but for our purposes we can focus the attention on the three contributions appearing in Eq.~\eqref{exp}. Clearly, the above working hypothesis entails that any derivative acting on $\psi$ results in the emergence of further powers of $c$ in Eq.~\eqref{gkg2}. Specifically, by denoting with $\partial}\def\ran{\rangle$ a generic derivative (which can be either spatial or temporal), one has
\begin{equation}\label{fd}
\partial}\def\ran{\rangle\psi=i\left}\def\lrar{\leftrightarrow(c^2\partial}\def\ran{\rangle\varphi_0+\partial}\def\ran{\rangle\varphi_1+c^{-2}\partial}\def\ran{\rangle\varphi_2\right}\def\ti{\tilde)e^{i\varphi}\,,
\end{equation}
whilst second derivatives imply
\begin{equation}\label{sd}
\partial}\def\ran{\rangle^2\psi=\left}\def\lrar{\leftrightarrow\{i\Big[c^2\partial}\def\ran{\rangle^2\varphi_0+\partial}\def\ran{\rangle^2\varphi_1+c^{-2}\partial}\def\ran{\rangle^2\varphi_2\Big]-\left}\def\lrar{\leftrightarrow[c^4\left}\def\lrar{\leftrightarrow(\partial}\def\ran{\rangle\varphi_0\right}\def\ti{\tilde)^2+2c^2\partial}\def\ran{\rangle\varphi_0\partial}\def\ran{\rangle\varphi_1+2\partial}\def\ran{\rangle\varphi_0\partial}\def\ran{\rangle\varphi_2+\left}\def\lrar{\leftrightarrow(\partial}\def\ran{\rangle\varphi_1\right}\def\ti{\tilde)^2+\frac{2}{c^2}\partial}\def\ran{\rangle\varphi_1\partial}\def\ran{\rangle\varphi_2+\frac{\left}\def\lrar{\leftrightarrow(\partial}\def\ran{\rangle\varphi_2\right}\def\ti{\tilde)^2}{c^4}\right}\def\ti{\tilde]\right}\def\ti{\tilde\}e^{i\varphi}\,.
\end{equation}
We can now apply the aforementioned scheme to Eq.~\eqref{gkg2} in order to write a Klein-Gordon equation which explicitly depends on $\varphi_0$, $\varphi_1$ and $\varphi_2$. To achieve the sought non-relativistic limit, the analysis must be truncated at $c^0$ starting from the highest order in the expansion parameter, which turns out to be $c^4$. Despite the seeming complexity, computations are instead straightforward, because in order to satisfy Eq.~\eqref{gkg2} it must necessarily be imposed that all the functions in front of any power of $c$ must vanish identically, order by order~\cite{kiefer}.

The first step is to look at the order $c^4$, whose contribution yields
\begin{equation}\label{c4}
\left}\def\lrar{\leftrightarrow(1+2\alpha x^2\right}\def\ti{\tilde)(\vec{\nabla}\varphi_0)^2=0\,,
\end{equation}
which tells that $\varphi_0$ is a function of time only. The explicit expression for $\varphi_0$ can be deduced from the order $c^2$, which reads
\begin{equation}\label{c2}
-\left}\def\lrar{\leftrightarrow(\partial}\def\ran{\rangle_t\varphi_0\right}\def\ti{\tilde)^2+m^2=0\,.
\end{equation}
Hence, up to $c^2$ the solutions to the Klein-Gordon equation are simply $\psi=\exp(\pm i m c^2t)$, which represent the particle/antiparticle wave functions with $E$ being the rest energy $m c^2$. In the non-relativistic regime, no particle creation occurs, so the two scenarios can be treated separately; however, since they would give the exactly same outcome at the next order, we can arbitrarily choose to work with only one of them. Thus, by selecting the minus sign (i.e., particle wave function) and moving on to the order $c^0$, we obtain
\begin{equation}\label{c0}
2m\partial}\def\ran{\rangle_t\varphi_1-i\left}\def\lrar{\leftrightarrow(1+2\alpha x^2\right}\def\ti{\tilde)\nabla^2\varphi_1-4i\alpha\vec{x}\cdot\vec{\nabla}\varphi_1-3\alpha+2m^2\phi+\left}\def\lrar{\leftrightarrow(1+2\alpha x^2\right}\def\ti{\tilde)\left}\def\lrar{\leftrightarrow(\vec{\nabla}\varphi_1\right}\def\ti{\tilde)^2=0\,.
\end{equation}
By introducing $\psi_1=\exp(i\varphi_1)$, the previous equation can be cast in a Schr\"odinger-like form, that is
\begin{equation}\label{finkg}
i\partial}\def\ran{\rangle_t\psi_1=-\left}\def\lrar{\leftrightarrow(1+2\alpha x^2\right}\def\ti{\tilde)\frac{\nabla^2\psi_1}{2m}+m\phi\psi_1-\frac{2\alpha}{m}\vec{x}\cdot\vec{\nabla}\psi_1-\frac{3\alpha}{2m}\psi_1\,.
\end{equation}
Higher-order terms in $c$ would only have added relativistic corrections to Eq.~\eqref{finkg}, which thereby represents the non-relativistic limit of Eq.~\eqref{gkg2}.

As for Eq.~\eqref{find}, we recover the two characteristic terms which are reminiscent of the minimal momentum scale (namely, the last two quantities of the r.h.s.). Furthermore, from the first two factors of Eq.~\eqref{finkg} the definitions of $m_i$ and $m_g$ can be immediately deduced; it turns out that the two masses are identical to the ones already found for the spinor field. In light of this, we conclude that the equivalence principle does not hold for scalar fields in the presence of the extended uncertainty principle, and the degree of its violation is precisely equal to the one evaluated with the Dirac equation, i.e., Eq.~\eqref{ratio}.

\section{Concluding remarks}

In this paper, we have shown how the existence of a minimal momentum scale induces a violation of the weak equivalence principle. To demonstrate the generality of our statement, we have separately considered spinor and scalar fields to verify that, in both cases, the degree of EP violation is still the same, thus suggesting that a similar occurrence may not depend on the nature of the examined quantum system. To unambiguously recognize the magnitude of the inertial and gravitational mass, we have performed a non-relativistic limit onto the Dirac and the Klein-Gordon equations, which have been written so as to simultaneously account for the presence of EUP and a weak and static gravitational field. In so doing, we have then identified $m_i$ with the mass term appearing in the kinetic part of the Schr\"odinger Hamiltonian and $m_g$ with the factor that multiplies the gravitational potential $\phi$. The same strategy has already been used in different frameworks with success~\cite{don,neutrino}, and it still represents the most reliable option, given that in the fully relativistic regime such an identification is a challenging task.

On a final note, we want to remark that, because of the physical origins of the extended uncertainty principle, the implications of the achieved result become relevant at cosmological scales. This concept can be quantified by comparing the EP-violating factor $2\alpha x^2$ with the available experimental constraint for the term $|m_g-m_i|/m_i$. As a matter of fact, by recalling that the above ratio is bounded by a factor $10^{-14}$~\cite{neutrino}, we immediately obtain
\begin{equation}\label{phys}
x\simeq\frac{10^{-7}}{\sqrt{|\alpha|}}\,.
\end{equation}
Since the magnitude of $\alpha$ is expected to be extremely small, the characteristic length scale $x$ at which the EP violation becomes significant is huge. By resorting to the formulation of EUP where $|\alpha|\simeq\Lambda\simeq10^{-52}$ m$^{-2}$, a straightforward calculation would yield $x_\Lambda\simeq10^{19}$ m, which is roughly $10^3$ light years. On the other hand, a more conservative approach that relies on experimental bounds for the EUP deformation parameter provides $\sqrt{|\alpha|}\lesssim10^{-12}\,\,\mathrm{m}^{-1}$~\cite{symm}, which instead would give a lower bound, that is $x_{\mathrm{min}}\gtrsim10^5\,\,\mathrm{m}$. Thence, the most promising phenomenological perspectives for the EUP-induced equivalence principle violation lies in the interval $[x_{\mathrm{min}},x_\Lambda]$. To unravel further pieces of information, more work is inevitably required along this direction.

\section*{Acknowledgement}

The author acknowledges networking support by the COST Action CA18108.


\section{Introduction}
\IEEEPARstart{T}{he} increased popularity of streaming, entertainment and social media content has garnered significant interest towards the development of methods that can improve perceptual video quality and provide better consumer experiences. Some of the approaches that have been explored in the past include higher spatial resolutions (4K/8K) \cite{ge2017toward}, high dynamic range (HDR) \cite{mai2010optimizing}, immersive content (virtual reality/augmented reality) \cite{smolic2007coding} and increasing temporal resolutions by employing high frame rate (HFR) videos. Among these, improving quality via HFR has been relatively less explored, as the majority of the currently available content is limited to 60 frames-per-second (fps). This has begun to change however, given the availability of high-speed cameras and displays, and the desire to stream high motion content, such as live sports.

A variety of factors have contributed to the limited adoption of HFR videos. In the past, capturing HFR content required expensive cameras which were not designed for television formats, thereby reducing the availability of HFR sequences. Further, display devices capable of supporting HFR content were less commonly available. Further, the perceptual benefits associated with higher frame rate were not clearly understood. However, recently advancements in both capture and display technologies have made HFR videos much more viable, thus increasing their potential use.

As the video industry transitions towards HFR regimes, it is important to systematically investigate and analyze the advantages associated with HFR videos. Further, perceptually optimizing modern bit rate ladders to accommodate frame rate changes, much as spatial downscaling is used now, could significantly improve streaming efficiency. Of course, adopting HFR implies significant potential increases in bandwidth consumption, hence is necessary to understand the merits of HFR. A common notion is that HFR videos tend to have better quality as the number of frames is increased, since it can often replicate motion more smoothly. However, very few studies have systematically evaluated these notions. The databases which are directed towards understanding and modeling the perceptual implications of frame rate include Waterloo HFR \cite{nasiri2015perceptual}, BVI-HFR \cite{mackin2018study} and LIVE-YT-HFR \cite{madhusudana2021liveythfr}, all of which contain Variable Frame Rate (VFR) and/or HFR videos.

Objective Video Quality Assessment (VQA) models are designed to predict human perception of video quality. VQA models can be broadly categorized into three types : Full Reference (FR) models which have access to pristine high quality reference signals \cite{wang2004image,wang2003multiscale,VMAF2016}, Reduced Reference (RR) models which use partial reference information \cite{soundararajan2012rred,soundararajan2012video}, and No Reference (NR) models which only have the distorted videos to make quality predictions. VFR-VQA problem is concerned with comparing reference and distorted videos at possibly difference frame rates, thus can be categorized as either FR or RR VQA models.

Existing FR VQA models typically require the pristine and distorted sequences being compared to share same frame rate, thus limiting their direct applicability to VFR-VQA problem. Although these methods can be extended to VFR videos by employing suitable preprocessing techniques, such as temporal upsampling of the distorted sequence, it has been observed that these operations often introduced additional artifacts, leading to inferior correlations against human judgments \cite{madhusudana2021liveythfr}. Moreover, these FR models are mainly intended to account for commonly observed artifacts such as compression, blur, scaling etc. but not those arising due to changes in frame rate.

Frame rate distortions are spatio-temporal, but the temporal aspects are more evident since these artifacts arise due to temporal subsampling. Thus, it is imperative to formulate effective, perceptually relevant temporal representations to achieve better quality predictions on VFR videos. Existing methods generally place less emphasis on temporal aspects of visual quality, usually employing frame differences \cite{soundararajan2012video,bampis2017speed,VMAF2016,bampis2018spatiotemporal}, which are generally insufficient to capture the complex variations, often spanning multiple frame instants, that can arise due to frame rate changes. Towards advancing progress on these problems, here we describe and implement an effective way to incorporate temporal features into existing VQA models, making them more sensitive to frame rate induced artifacts.

Our proposed model uses a fusion framework, whereby we combine spatial features from any existing FR-VQA model with temporal features from the recently developed GREED \cite{madhusudana2021greed} model. The GREED model has been shown to achieve superior correlations against subjective judgments of VFR video quality, when evaluated on videos of varied frame rates. In \cite{madhusudana2021pcs}, a preliminary version of this current work, we experimented with fusing GREED and Video Multi-method Fusion (VMAF) \cite{VMAF2016} features towards capturing frame rate dependent artifacts. The results we achieved were provocative enough to motivate the current greatly expanded study. Here, we extend this idea and show that this type of feature fusion can be generalized to any existing FR VQA models, including FR Image Quality Assessment (IQA) models resulting in significant performance gains in capturing the perceptual effects of frame rate distortions. The design framework is quite simple, it effectively promotes the advantages of both the fused models, resulting in robust performance on both HFR/VFR and non-VFR VQA standard databases.

\begin{algorithm}[t] \caption{GST Feature Fusion}
\begin{algorithmic}[1]
\Require reference video $R$, distorted video $D$, VQA model F
\Ensure GST-F features
\State Calculate spatial index from model F, $\text{S} = \text{F}(R,D)$
\State Temporally subsample $R$ to get $PR$
\State Scale $S = \{4,5\}$, band-pass filter bank $b_k: k \in \{1,\ldots 7\}$
\ForEach {$s \in S $}
\ForEach {$b_k$}
\State Calculate TGREED$_k$ from (\ref{eqn:GTI})
\EndFor
\EndFor
\State Concatenate S and TGREED to obtain GST-F features.
\end{algorithmic}
\label{alg1}
\end{algorithm}

\section{Background and Method}
\label{sec:GREED}
Our proposed method involves employing the spatial components from existing FR VQA models, and the temporal component of GREED. In our experiments we use four leading VQA models : SSIM \cite{wang2004image}, MS-SSIM \cite{wang2003multiscale}, ST-RRED \cite{soundararajan2012video} and SpEED \cite{bampis2017speed}. The framework is fairly generic in nature and can easily be extended to other FR VQA models. In this section we first provide a brief overview of temporal GREED, followed by a discussion of our proposed feature fusion methodology.

\begin{table}[t]
\caption{Performance comparison of FR algorithms on the LIVE-YT-HFR Database. The reported numbers are median values from 200 iterations of randomly chosen train-test sets. In each column the first and second best models are boldfaced. We obtained numbers from \cite{madhusudana2021greed,madhusudana2021pcs} for the results of GREED and GREED-VMAF.}
    \label{Table:MOS_comparison}
    \centering
    \begin{tabular}{|c||c|c|c|c|}
        \hline
        ~    & SROCC $\uparrow$ & KROCC $\uparrow$ & PLCC $\uparrow$ & RMSE $\downarrow$ \\ \hline \hline

        SSIM \cite{wang2004image} & 0.5566 & 0.4042 & 0.5418 & 9.99 \\
        GST SSIM & 0.7576 & 0.5648 & 0.7700 & 7.40 \\ \hline
        MS-SSIM \cite{wang2003multiscale} & 0.5742 & 0.4135 & 0.5512 & 10.01 \\
        GST MS-SSIM & 0.8128 & 0.6204 & 0.8179 & 6.75 \\ \hline
        ST-RRED \cite{soundararajan2012video} & 0.6394 & 0.4516 & 0.6073 & 9.58 \\
        GST ST-RRED & 0.8029 & 0.5937 & 0.8144 & 6.91 \\ \hline
        SpEED \cite{bampis2017speed} & 0.6051 & 0.4437 & 0.5206 & 10.28 \\
        GST SpEED & 0.8276 & 0.6343 & 0.8346 & 6.51 \\ \hline
        VMAF \cite{VMAF2016}& 0.7782 & 0.5918 & 0.7419 & 8.10 \\
        GREED & \textbf{0.8822} & \textbf{0.7046} & \textbf{0.8869} & \textbf{5.48} \\
        GREED-VMAF & \textbf{0.8658} & \textbf{0.6840} & \textbf{0.8723} & \textbf{5.89} \\
        \hline
    \end{tabular}
\end{table}

\vspace{-7pt}
\subsection{Temporal GREED (TGREED)}
GREED model is motivated by deviations that occur in the statistical distributions of temporal band-pass video coefficients that arise as frame rates are varied. In particular TGREED, captures the temporal aspects of frame rate distortions by employing a bank of temporal band-pass filters. The video $V(\mathbf{x},t)$ ($\mathbf{x} = (x,y)$ is subjected to temporal band-pass decomposition by a bank of $K$ 1D filters denoted by $b_k$ for $k \in \{1,\ldots K\}$
\begin{align}
    B_k(\mathbf{x},t) = V(\mathbf{x},t)*b_k(t) \text{\hspace{10pt}} \forall k \in \{1,\ldots K\},
\end{align}
where $B_k$ denotes the band-pass response of a video to the $k^{th}$ filter. The distributions of $B_k$ reliably exhibit heavy tailed shapes and can be modeled as following a Generalized Gaussian Distribution (GGD). Deviations from these model distributions caused by distortion are quantified by comparing the entropies of the reference and distorted sequences. The entropy of the GGD is given by
\begin{align}
    h(X) = \frac{1}{\beta} - \log \left(\frac{\beta}{2\alpha\Gamma(1/\beta) }\right).
    \label{eqn:ggd_entropy}
\end{align}
where $\beta$ and $\alpha$ are parameters of the GGD, and $\Gamma$ denotes the Gamma function:
\begin{align*}
    \Gamma(a) = \int_0 ^{\infty}x^{a-1}e^{-x} dx .
\end{align*}
The band-pass coefficients $B_k$ are partitioned into non-overlapping patches of size $\sqrt{M} \times \sqrt{M}$, which are indexed by $p \in \{ 1,2, \ldots P \}$. In our experiments we employed spatial patches of size $5 \times 5$ ({i}.{e}.{ }$\sqrt{M} = 5$). The GGD parameters $\beta$ and $\alpha$ are computed on each patch using the kurtosis matching procedure detailed in \cite{soury2015new,pan2012exposing}.

\begin{table*}[t]
	\caption{Performance comparison of FR methods for specific frame rates on the LIVE-YT-HFR Database. The reported numbers are median values over 200 iterations of randomly chosen train-test sets. In each column the first and second best values are marked in boldface. We obtained numbers from \cite{madhusudana2021greed,madhusudana2021pcs} for the results of GREED and GREED-VMAF.}
	\label{Table:FPS_comparison}
	\centering
	\scriptsize
	\scalebox{1.02}{
		\begin{tabular}{|c||c|c|c|c|c|c|c|c|c|c|c|c|}
			\hline
			& \multicolumn{2}{|c|}{24 fps} & \multicolumn{2}{|c|}{30 fps} & \multicolumn{2}{|c|}{60 fps} & \multicolumn{2}{|c|}{82 fps} & \multicolumn{2}{|c|}{98 fps} & \multicolumn{2}{|c|}{120 fps} \\
			\cline{2-13}
			~ & SROCC$\uparrow$ & PLCC$\uparrow$ & SROCC$\uparrow$ & PLCC$\uparrow$ & SROCC$\uparrow$ & PLCC$\uparrow$ & SROCC$\uparrow$ & PLCC$\uparrow$ & SROCC$\uparrow$ & PLCC$\uparrow$ & SROCC$\uparrow$ & PLCC$\uparrow$ \\ \hline \hline
			SSIM \cite{wang2004image} & 0.266 & 0.222 & 0.283 & 0.189 & 0.382 & 0.302 & 0.371 & 0.362 & 0.537 & 0.497 & 0.867 & 0.833 \\
			GST SSIM & 0.386 & 0.635 & 0.461 & 0.715 & 0.516 & 0.722 & 0.698 & 0.835 & 0.743 & 0.833 & 0.797 & 0.817 \\ \hline
			MS-SSIM \cite{wang2003multiscale} & 0.305 & 0.260 & 0.296 & 0.238 & 0.416 & 0.338 & 0.439 & 0.393 & 0.578 & 0.561 & 0.706 & 0.696 \\
			GST MS-SSIM & 0.505 & 0.682 & 0.495 & 0.769 & 0.579 & 0.796 & 0.704 & 0.844 & 0.775 & 0.841 & 0.832 & 0.838 \\ \hline
			ST-RRED \cite{soundararajan2012video} & 0.305 & 0.275 & 0.296 & 0.206 & 0.612 & 0.613 & 0.584 & 0.513 & 0.650 & 0.604 & 0.755 & 0.696 \\
			GST ST-RRED & 0.518 & 0.698 & 0.421 & 0.664 & 0.580 & 0.814 & 0.684 & 0.833 & 0.752 & 0.816 & 0.888 & 0.885 \\ \hline
			SpEED \cite{bampis2017speed} & 0.432 & 0.273 & 0.410 & 0.233 & 0.439 & 0.292 & 0.546 & 0.390 & 0.578 & 0.471 & 0.758 & 0.739 \\
			GST SpEED & 0.645 & 0.744 & 0.616 & 0.759 & 0.577 & 0.745 & 0.723 & 0.789 & 0.787 & 0.827 & 0.860 & 0.867 \\ \hline
			VMAF \cite{VMAF2016} & 0.250 & 0.368 & 0.362 & 0.471 & 0.630 & 0.680 & 0.734 & 0.793 & \textbf{0.860} & 0.868 & 0.818 & 0.816 \\
			GREED & \textbf{0.727} & \textbf{0.822} & \textbf{0.702} & \textbf{0.843} & \textbf{0.732} & \textbf{0.840} & \textbf{0.818} & \textbf{0.896} & \textbf{0.864} & \textbf{0.891} & \textbf{0.888} & \textbf{0.895} \\
			GREED-VMAF & \textbf{0.748} & \textbf{0.805} & \textbf{0.743} & \textbf{0.833} & \textbf{0.773} & \textbf{0.836} & \textbf{0.786} & \textbf{0.880} & \textbf{0.860} & \textbf{0.899} & \textbf{0.881} & \textbf{0.903} \\
			\hline
		\end{tabular}
	}
\end{table*}

To lend a more local nature to the entropy estimates, and to also provide numerical stability in regions having low variances, while preferring more interesting regions of greater activity, the entropy estimates $h(\Tilde{B}_{kpt})$ are scaled by their local variances as
\begin{align*}
    \epsilon_{kpt} = [\log(1+\sigma^2(\Tilde{B}_{kpt}))]h(\Tilde{B}_{kpt}),
\end{align*}
where $B_{kpt}$ the denotes band-pass coefficients corresponding to frame $t$, patch $p$ and band-pass filter $k$, and $\sigma^2$ denotes variance. The changes in entropy values can arise due to two lossy operations : temporal subsampling and compression. To achieve effective quality prediction, a quality model must account for a combination of these operations. However, in \cite{madhusudana2021greed} it was observed that entropy changes due to variations in frame rate often dominated those arising from compression, {i}.{e}., are biased towards VFR distorions. Thus, simple comparisons between the entropies of reference and distorted sequences is often insufficient to capture the overall quality. To overcome this shortcoming, a "debiasing" stage is implemented, which requires defining an additional sequence "Pseudo Reference" video $PR$ in addition to the reference video $R$ and distorted video $D$ . $PR$ is obtained by temporally subsampling $R$ to have the same frame rate as $D$. For cases where $R$ and $D$ have the same frame rates, then $PR$ is the same as $R$. Using these signals TGREED is defined as

\begin{align}
    \text{TGREED}_{kt} = \frac{1}{P}\sum_{p=1} ^P \Bigg|\Big(1 + |\epsilon_{kpt} ^D - \epsilon_{kpt} ^{PR}|\Big) \frac{\epsilon_{kpt} ^R + 1}{\epsilon_{kpt} ^{PR} + 1} - 1 \Bigg|.
    \label{eqn:GTI}
\end{align}
This expression can be understood by decomposing it into two components - absolute difference term, and the ratio term. The ratio term is dependent only on the frame rate of distorted video $D$, and is invariant to compression effects, since $D$ is not required for its calculation. On the other hand, the difference term mainly accounts for compression artifacts, since $PR$ can be considered as a lossless version of $D$. The entire TGREED index can be considered as a weighted entropic difference.
\vspace{-8pt}
\subsection{Feature Integration}
\label{sec:feat_integration}

The main goal of this work is to show that using existing FR VQA models that perform poorly when VFR distortions are present, can be transformed into more powerful models able to address frame rate distortions. We accomplish this using a fusion framework, combining spatial features from existing VQA methods, and TGREED features which account for temporal artifacts. Current state-of-the-art (SOTA) VQA models, like those studied here exhibit excellent performance when evaluated on legacy databases like LIVE-VQA \cite{seshadrinathan2010study}, CSIQ-VQA \cite{vu2014vis3} etc. These datasets contain reference and distorted videos having the same frame rates and, include common distortions such as compression, scaling, packet loss etc. However, when evaluated on VFR sequences, these models under-perform due to the presence of frame rate artifacts, as evidenced in Table \ref{Table:MOS_comparison}.

Our feature fusion method is motivated by the successes of earlier fusion-based VQA models \cite{VMAF2016,bampis2018spatiotemporal,madhusudana2021pcs}. Combining features exploits the advantages offered by both models into a single framework. We retain the spatial components present in the current models, and replace the temporal components (if any) with TGREED features. We refer to these modified VQA versions using the Generalized Space Time (GST) prefix. TGREED features are computed at two scales $s=4,5$, where the frames are spatially subsampled by $2^s$ times before calculating the entropy values. The above subsampling was observed to have the optimal trade-off in terms of computational complexity and performance. A biorthogonal-2.2 wavelet filter with three levels of wavelet packet decomposition \cite{coifman1992entropy} was used to obtain temporal band-pass coefficients. Ignoring the low pass subband, 7 TGREED features from the respective temporal subbands were calculated at each scale, yielding a total of 14 features from both scales. For each tested VQA model, a single value corresponding to the spatial component was computed, then concatenated with the TGREED features, resulting in a total of 15 features. These features are then fed to a Support Vector Regressor (SVR), and a mapping from feature space to quality scores is learned using ground-truth subjective opinion scores. The entire feature fusion methodology is summarized in Algorithm \ref{alg1}. The code implementations used in this work have been made available online\footnote{\url{https://github.com/pavancm/GST_VQA}}.

\begin{table*}[t]
\caption{Performance comparison of various FR methods on multiple VQA databases. These databases do \textbf{not} contain variable frame rates. The reported numbers are median values over multiple iterations of randomly chosen train-test sets. In each column the first and second best values are boldfaced. We obtained numbers from \cite{madhusudana2021greed,madhusudana2021pcs} for the results of GREED and GREED-VMAF.}
\label{table:MOS_VQA}
\centering
    \begin{tabular}{|c||cc|cc|cc|cc|cc|}
        \hline
        \multirow{2}{*}{} & \multicolumn{2}{|c|}{LIVE-VQA \cite{seshadrinathan2010study}}  & \multicolumn{2}{|c|}{LIVE-mobile \cite{moorthy2012video}} & \multicolumn{2}{|c|}{CSIQ-VQA \cite{vu2014vis3}} & \multicolumn{2}{|c|}{ETRI-LIVE STSVQ \cite{lee2021subjective}} & \multicolumn{2}{|c|}{AVT-VQDB-UHD-1 \cite{rao2019avt}}\\ \cline{2-11}
        ~ & SROCC$\uparrow$ & PLCC$\uparrow$ & SROCC$\uparrow$ & PLCC$\uparrow$ & SROCC$\uparrow$ & PLCC$\uparrow$ & SROCC$\uparrow$ & PLCC$\uparrow$ & SROCC$\uparrow$ & PLCC$\uparrow$ \\ \hline \hline

        SSIM \cite{wang2004image} & 0.802 & 0.791 & 0.798 & 0.870 & 0.705 & 0.724 & 0.512 & 0.549 & 0.777 & 0.693 \\
        GST SSIM & 0.747 & 0.754 & 0.875 & 0.896 & 0.836 & 0.821 & 0.820 & 0.860 & 0.804 & 0.820 \\ \hline
        MS-SSIM \cite{wang2003multiscale} & \textbf{0.830} & \textbf{0.804} & 0.800 & 0.874 & 0.757 & 0.755 & 0.498 & 0.473 & 0.784 & 0.715 \\
        GST MS-SSIM & 0.745 & 0.730 & 0.879 & 0.888 & 0.873 & 0.847 & 0.831 & 0.845 & 0.807 & 0.820\\ \hline
        ST-RRED \cite{soundararajan2012video} & \textbf{0.826} & \textbf{0.828} & 0.882 & 0.935 & 0.813 & 0.793 & 0.514 & 0.562 & 0.837 & 0.778 \\
        GST ST-RRED & 0.793 & 0.794 & 0.880 & 0.925 & \textbf{0.890} & \textbf{0.927} & \textbf{0.873} & \textbf{0.893} & 0.820 & 0.822 \\ \hline
        SpEED \cite{bampis2017speed} & 0.801 & 0.815 & 0.886 & \textbf{0.929} & 0.743 & 0.743 & 0.466 & 0.552 & 0.834 & 0.785 \\
        GST SpEED & 0.782 & 0.774 & 0.890 & \textbf{0.927} & 0.793 & 0.776 & 0.842 & 0.860 & 0.821 & 0.834 \\ \hline
        VMAF \cite{VMAF2016}& 0.794 & 0.794 & \textbf{0.897} & 0.902 & 0.618 & 0.623 & 0.660 & 0.651 & \textbf{0.872} & \textbf{0.873} \\
        GREED & 0.750 & 0.755 & 0.863 & 0.908 & 0.780 & 0.766 & 0.730 & 0.760 & 0.823 & 0.823 \\
        GREED-VMAF & 0.784 & 0.786 & \textbf{0.904} & 0.923 & \textbf{0.907} & \textbf{0.899} & \textbf{0.860} & \textbf{0.866} & \textbf{0.842} & \textbf{0.853} \\
        \hline
    \end{tabular}
\end{table*}

\section{Experiments and Results}
\label{sec:experiments}

\subsection{Experimental Settings}
\textbf{Compared Methods} : For performance comparison we choose 5 popular VQA models : SSIM \cite{wang2004image}, MS-SSIM \cite{wang2003multiscale}, ST-RRED \cite{soundararajan2012video}, SpEED \cite{bampis2017speed} and VMAF \cite{VMAF2016}. Of these, the first two are IQA (frame based) models, which are computed on every frame, then temporally averaged to obtain final video quality scores. The above models require the reference and distorted sequences to have same frame rate, thus, when comparing videos with differing frame rates, the distorted sequences were temporally upsampled using frame duplication to have the same frame rate as the reference. In addition to the above models, we also included the recently proposed GREED-VMAF model \cite{madhusudana2021pcs}, which combines GREED and VMAF features to obtain better performance. We employ commonly used evaluation metrics for VQA models, such as Spearman's rank order correlation coefficient (SROCC), Kendall's rank order correlation coefficient (KROCC), Pearson's linear correlation coefficient (PLCC), and the root mean squared error (RMSE).

\vspace{-7pt}
\subsection{Performance Comparison on the LIVE-YT-HFR Dataset}
The first experiment we performed was to analyze the impact of adding TGREED features to VQA models in regards to capturing frame rate artifacts. In the analysis, we used the LIVE-YT-HFR dataset \cite{madhusudana2021liveythfr}, which contains a total of 480 videos spanning 6 different frame rates. For training purposes, the LIVE-YT-HFR dataset was randomly partitioned into 70\%/15\%/15\% sets for training, validation and testing, respectively. We also ensured that there was no overlap of contents between the subsets. To avoid bias towards any specific train-test combination, the experiment was repeated over 200 different random train-test splits, and the median performance across these combinations is reported. The performances of the different models is compared in Table \ref{Table:MOS_comparison}. For ease of comparison, the original VQA models and corresponding GST versions are clustered together. From the Table, it can be observed that inclusion of the TGREED features provides significant boosts in performance of all the compared models.

In order to analyze model correlations on each individual frame rate, the LIVE-YT-HFR dataset was divided into subsets containing videos having a single frame rate with the resulting performances on each subset reported in Table \ref{Table:FPS_comparison}. From the Table it is evident that there was significant enhancement in performance due to the inclusion of TGREED features on each frame rate, especially low frame rates.

\vspace{-7pt}
\subsection{Performance Comparison on other VQA Databases}
To analyze the generalizability of the proposed design, we conducted performance comparison of the compared models on traditional VQA databases containing a variety of spatial and temporal artifacts, but single, shared frame rates. For evaluation purposes, we used 5 such legacy VQA datasets, LIVE-VQA \cite{seshadrinathan2010study}, LIVE-mobile \cite{moorthy2012video}, CSIQ-VQA \cite{vu2014vis3}, ETRI-LIVE STSVQ \cite{lee2021subjective} and AVT-VQDB-UHD-1 \cite{rao2019avt}. These datasets span a variety of spatial resolutions, contain various types of distortions such as scaling, compression, white noise etc., and employ multiple compression schemes such as MPEG, H264, H265, VP9 etc. Thus an exhaustive evaluation on them should validate the generalizability of a GREED-enhanced objective VQA model. Note that all the videos in the above datasets have reference and corresponding distorted videos of identical frame rates. Although some datasets contain distorted videos at different frame rates, the authors reported that appropriate temporal processing was applied so that the videos viewed by the human subjects had the same frame rate as the reference.

Each dataset was divided into 80\% and 20\% non-overlapping subsets based on content for training and testing, respectively. This procedure was repeated over all possible train-test splits, and the median performance values reported in Table \ref{table:MOS_VQA}. ETRI-LIVE STSVQ and AVT-VQDB-UHD-1 datasets contain many more source videos, thus for these two datasets the train-test splits were repeated 200 times, similar to LIVE-YT-HFR. From the Table it may be observed that including TGREED features leads to better correlations on most of the compared datasets, and a few of the modified versions also achieve top performance on certain datasets. These results highlight the generalization capability of the TGREED index as well as the efficacy of including effective temporal features for objective video quality prediction.

\section{Conclusion}
\label{sec:conclusion}
We proposed a feature fusion method to improve existing VQA models to be able to capture frame rate distortions. The design approach exploits the superior prediction power of the temporal components present in GREED. This flexible framework is quite generic and can be applied to any existing FR VQA models, including IQA models. We performed a comprehensive evaluation on the LIVE-YT-HFR dataset, and showed that fusing TGREED with existing SOTA models is highly beneficial for improving quality predictions on videos afflicted with frame rate distortions. We also evaluated the generalizability of the proposed method on legacy fixed frame rate databases and observed that the TGREED-enhanced models often achieved better performance than original models. Our results indicate that using powerful temporal representations can lead to better VQA models for both VFR and non-VFR scenarios. These VQA models can be beneficial in determining frame rate dynamically in a constrained bandwidth environment.

\ifCLASSOPTIONcaptionsoff
  \newpage
\fi

\bibliographystyle{IEEEtran}


\section{Introduction}
The future sixth generation of communication standards is envisioned to support nanomachine communication \cite{tripathi2021}. For intra-body communication, molecular communication (MC) is a promising solution for reliable communication among bio-nanomachines (nanomachines with biological functional parts)  \cite{nakano2013}. MC uses molecules as the carrier of information, termed information molecules (IMs). The out of the body communication can be effectively performed with the help of electromagnetic (EM) waves like tera-hertz waves\cite{tripathi2021}. A hybrid system utilizing both MC and EM-based communication can help realize efficient healthcare systems. Here, MC systems can gather information from inside the human body, and  EM-based systems convey the gathered information to a healthcare facility to process the information and take necessary actions.

In intra-body communication systems, multiple receiver bio-nanomachines may coexist in the same communication channel. For example, multiple malicious receivers may coexist with MC systems inside a human body. The malicious receivers can act as eavesdroppers or hinder the triggering of MC detection systems by absorbing the IMs \cite{eved}. Such actions have profound security implications in e-health systems. Multiple receivers can also be utilized to create single-input multiple-output (SIMO) and multiple-input multiple-output (MIMO) communication links with improved reliability and performance.

An important metric in characterizing a molecular channel is hitting probability which denotes the probability that an IM hits the intended receiver (or equivalently the fraction of IMs hitting the receiver) within time $t$.
For a diffusion-based MC system with a point transmitter and a fully-absorbing receiver (FAR) in three-dimension (3-D), the hitting probability was discussed in \cite{yilmaz2014}.
When there are other FARs present in the same medium, they can capture the IMs transmitted to the intended receiver.
Therefore, the hitting probability decreases with the number of additional receivers, thus degrading the communication reliability. It is important to derive the hitting probability for general $N$ FARs to study the extent of the impact of additional receivers.
However, there is no analytical channel model for an MC system with multiple FARs. Many of the works \cite{koo2016, damrath2018,lee2017} relied on fitting models, machine learning-based approaches, and particle-based simulations to analyze systems with multiple FARs. For a system with two FARs, an approximate analytical expression for the hitting probability on \textit{each of} the FAR was derived in our past work \cite{sabu2020a}. Later this work was extended to derive the approximate expression for hitting probability of an IM on FARs with different sizes \cite{sabu2021ab}. An approximate analytical equation for the hitting probability of an IM on \textit{any one} of the FAR in a system with multiple FARs was derived in \cite{sabu2020}. However, analysis for an arbitrary value of $N$ is still missing in the current literature.
 Analyzing MC systems with multiple FARs using the methods employed in  \cite{koo2016, damrath2018,lee2017} are time-consuming and inconvenient. Even though obtaining exact analytical expressions for the hitting probability of MC systems with multiple FARs can be difficult, approximate analytical expressions with good accuracy are of a greater value, which is the focus of the work.

This work considers a diffusion-based MC system in 3-D with a point transmitter and $N$ FARs. We first consider $N=3$ and derive  the hitting probability of an IM on one FAR in presence of the rest of the FARs, which is validated using simulations. We then provide the approach for the general $N$ case.
The contributions of this work are listed below.
\begin{itemize}[leftmargin=10pt]
	\item We consider an MC system with three FARs located in arbitrary positions in 3-D. We derive an approximate yet accurate analytical expression for the hitting probability of an IM on each of the FARs while considering the impact of the rest of the FARs. The derived equation is then validated using particle-based simulations.
	\item Using the derived equation for the 3 FAR system, a very simple expression is obtained when the FARs are symmetrically positioned (\textit {i.e., } mutually equidistant and equidistant from the transmitter).

	\item We also study the impact of malicious receivers on the intended receiver's performance by  quantifying the mutual influence of FARs and discuss when the impact is minimal.
	\item We also study the gain in detection performance when FARs can cooperate with each other.
	\item We also provide a framework to derive the  hitting probability expression for an MC system with  $ N $  FARs for an arbitrary value of $N$.

\end{itemize}
\textbf{Notation}:
 $ \mathcal{L}\left[.\right] $ denotes the Laplace transform and
 $ \mathcal{L}^{-1}\left[.\right] $ denotes the inverse Laplace transform.

\section{System Model}
	\begin{figure}[ht]
	\centering
	\includegraphics[width=0.65\linewidth]{Figures/3rxSM}
	\caption{An illustration showing a 3-D MC  system with a point transmitter at the origin and three FARs located at positions $\pmb{A}_{1},\ \pmb{A}_{2},\ \text{and},\ \pmb{A}_{3}$ in $ \mathbb{R}^3 $. }
	\label{Fig:sm}
\end{figure}
We consider a diffusion-limited MC system with a point transmitter located at origin $\mathbf{o}$ and $N=3$ FARs located at $\pmb{A}_i$'s ($ i\in\{1,2,3\} $) in a 3-D medium as shown in Fig. \ref{Fig:sm}. The radius of all the FARs is assumed to be the same and is denoted by $ a $.
The IM emitted by the point transmitter propagates in the medium via Brownian motion\cite{einstein2011}. The propagation medium is assumed to be homogeneous, and the diffusion coefficient (denoted by $ D $) of the IM is assumed to be constant over space and time. An IM is absorbed and is detected with a probability $ 1 $ when it hits the surface of FAR.

The considered setup can be applied to various scenarios, including cooperative and non-cooperative settings. For example, the multiple FARs can be cooperating where they work together as a SIMO/MIMO system to decode the same message transmitted from $\mathbf{o}$. In a different scenario, these FARs can act as malicious or non-cooperative receivers interfering with the other receivers' communication. Without loss of generality, we assume that the first FAR FAR$_1$ is the intended FAR and others are cooperating or non-cooperating FARs. The probability that an IM hits the FAR$_1$ within time $t$  is known as the hitting probability.

Let the distance between the transmitter and the FAR at $ \pmb{A}_i $ be denoted by $ r_i $, where $ r_i=\|\pmb{A}_i\|,\ i\in\mathbb{N}$. The angle between $ i $th FAR located at $ \pmb{A}_i $ (denoted by FAR$ _i $) and FAR$_j $ is represented by $ \phi_{ij},\ i,j\in\mathbb{N},\ i\neq j $. Let us denote the closest point on the surface of FAR$ _i $ from the transmitter by $ \pmb{B}_{i} $. Therefore the distance between the closest point of FAR$ _i $ from the transmitter  and the center of FAR$ _j $ is
$ 	R_{ij}=\sqrt{(r_i-a)^2+r_{j}^2-2(r_i-a)r_{j}\cos\phi_{ij}}. $

Note that the hitting probability expression is known  for $N=1$ and $N=2$, which are given respectively as  \cite{yilmaz2014,sabu2020a}
	\begin{align}
		\overline{\mathsf{p}}(t,r_1)&=\frac{a}{r_1}\mathrm{erfc}\left(\frac{r_1-a}{\sqrt{4Dt}}\right),  \text{ and}\label{eq1rx}\\

		\hp{}(t,\pmb{A}_1{\mid} \pmb{A}_2)&=\sum_{n=0}^{\infty} \frac{a^{2 n}}{R_{12}^{n} R_{21}^{n}}\left[\frac{a}{r_1} \mathrm{erfc}\left(\frac{r_1{-}a+n\left(R_{21}{-}a\right)+n\left(R_{12}{-}a\right)}{\sqrt{4 D t}}\right)\right. \nonumber\\
		&\left.\quad -\frac{a^{2}}{r_2 R_{21}} \mathrm{erfc}\left(\frac{r_2{-}a+(n+1)\left(R_{21}{-}a\right)+n\left(R_{12}{-}a\right)}{\sqrt{4 D t}}\right)\right].
		\label{eq2rx}
	\end{align}
However,  the expression for $N\ge3$ is not available. We will first consider  three FARs (\textit {i.e., } $N=3$) located at $ \pmb{A}_1,\pmb{A}_2 $ and $ \pmb{A}_3 $ and later extend the analysis for the general case.

\section{Channel Model for an MC System with 3-FARs}
In this section, we will analyze an MC system with 3 FARs in terms of hitting probability.
\subsection{Hitting probability}
We now derive the hitting probability of an IM on the FAR in a system with three FARs, located at $\pmb{A}_1$, $\pmb{A}_2$ and $\pmb{A}_3$. It can be seen from \eqref{eq2rx} that $ \hp{}(t,\pmb{A}_1{\mid} \pmb{A}_2)\leq \overline{\mathsf{p}}(t,r_1)$, which shows that the hitting probability decreases when an additional receiver is present.
It is expected that for an MC system with three FARs, the hitting probability would be less than   $ \hp{1}(t,\pmb{A}_1{\mid} \pmb{A}_2)$.
Let us denote the hitting probability until time $ t $ of an IM emitted by a point transmitter located at the origin on FAR$ _1 $ (located at $ \pmb{A}_1 $) in the presence of the other two FARs by $ \hp{}\left(t,\pmb{A}_1{\mid} \pmb{A}_2,\pmb{A}_3\right)$,
or just $ \hp{1}(t) $ (when there is no ambiguity)
which is given in the theorem below.

\begin{theorem}\label{t3rx}
	In an MC system with a point transmitter at the origin and three FARs located at $\pmb{A}_1,\ \pmb{A}_2$ and $\pmb{A}_3$, the probability that an IM emitted from the point transmitter hits  $\mathrm{FAR} _1 $ within time $ t $ in the presence of other two FARs is approximately given by
	\begin{align}
	\hp{1}(t)&=\mathcal{L}^{-1}\left[\lhp{1}(s)\right],\label{eq3rx}
\end{align}
{ where } $\lhp{1}(s)$ is the Laplace transform of $\hp{1}(t)$ and is given by
\begin{align}
	\lhp{1}(s)= \frac{\overline{\mathcal{P}}(s,r_1) -s\alpha_1(s)+s^2\beta_1(s)}
	{1-s^2\gamma(s)+s^3\delta(s)}\label{lp3rx}.
\end{align}
 Here,
 \begin{align*}
 \alpha_1(s)&=\overline{\mathcal{P}}(s,r_2)\overline{\mathcal{P}}(s,R_{21}){+}\overline{\mathcal{P}}(s,r_{13})\overline{\mathcal{P}}(s,R_{31})\nonumber\\
 	\beta_1(s)&=-\overline{\mathcal{P}}(s,r_1)\overline{\mathcal{P}}(s,R_{23})\overline{\mathcal{P}}(s,R_{32})+\overline{\mathcal{P}}(s,r_2)\overline{\mathcal{P}}(s,R_{23})\overline{\mathcal{P}}(s,R_{31})\nonumber\\
 	 &\qquad+\overline{\mathcal{P}}(s,r_{3})\overline{\mathcal{P}}(s,R_{32})\overline{\mathcal{P}}(s,R_{21})\nonumber\\
 	\gamma(s)&=\overline{\mathcal{P}}(s,R_{12})\overline{\mathcal{P}}(s,R_{21}) + \overline{\mathcal{P}}(s,R_{32})\overline{\mathcal{P}}(s,R_{23})
	+\overline{\mathcal{P}}(s,R_{13}) \overline{\mathcal{P}}(s,R_{31})\nonumber\\
 	\text{and
 	\delta(s)&=\overline{\mathcal{P}}(s,R_{12}) \overline{\mathcal{P}}(s,R_{23})\overline{\mathcal{P}}(s,R_{31})
	+ \overline{\mathcal{P}}(s,R_{13})\overline{\mathcal{P}}(s,R_{32})\overline{\mathcal{P}}(s,R_{21})\nonumber.
 \end{align*}
Here, $\overline{\mathcal{P}}(s,x
$ is the Laplace transforms of  $\overline{\mathsf{p}}(t,x)$
given as
\begin{align*}
	\overline{\mathcal{P}}(s,x)&=\laplace{}{\overline{\mathsf{p}}(t,x)}=\frac{a}{sx}\exp\left(-\left(x-a\right)\sqrt{\frac{s}{D}}\right)\nonumber\hspace{3.5in}\ \quad
\end{align*}
with $\ x\in\{r_j,R_{ij}\}$ for all $i,j\in\{1,2,3\}, i\neq j
$.
\begin{proof}
	See Appendix \ref{a3rx}.
\end{proof}
\end{theorem}

The source of approximation in the above theorem is the approximation of the hitting point of an IM on the surface of the FAR with the closest point on FAR as given in Appendix \ref{a3rx}. This approximation is of good accuracy and is verified in Section \ref{acc}.\par

\subsection{Evaluation of the hitting probability for $ 3-$FAR case}

Fig. \ref{fig:3rxf1} shows the variation of the hitting probability within time $ t $ with time $ t $ for an MC system with three FARs. The validity of the Theorem \ref{t3rx} is verified using particle-based simulation in MATLAB. In the particle-based simulation for the considered system, an IM is generated at the transmitter at the origin at time $ t=0 $. The IM undergoes Brownian motion with diffusion coefficient $ D $. The position of the IM is tracked for every time step $ \Delta t $. If the distance between the IM and the center of the FAR is less than the radius of the FAR, the IM is considered to be hit at the surface of the FAR. The hitting process is accurate when $ \Delta t $ is too small. This procedure is continued for many iterations, and the mean of the hitting events is taken to get the hitting probability.
 Fig. \ref{fig:3rxf1} verifies that the derived equations are in good match with the simulation results. We can also observe that the hitting probability of FARs increases when it is near the transmitter.

\begin{figure}[ht]
	\centering
	\includegraphics[width=0.65\linewidth]{Figures/threerx}
	\caption{Variation of $\hp{i}(t)$ with time $t$ for a MC system with three FARs. Parameters: $\pmb{A}_1= [25,0,0],\ \pmb{A}_2=[-25,5,0],\ \pmb{A}_3=[20,-15,10],\ a=5\ \mu m,\ D=100\mu m^2/s,\ \Delta t=10^{-5}\ s $.}
	\label{fig:3rxf1}
	\vspace{-0.5cm}
\end{figure}

\subsection{Accuracy of the hitting probability equations}\label{acc}

Fig.  \ref{fig:aberr} shows the variation of the absolute error in the derived hitting probability expression of FAR$ _1 $
defined as
\begin{align*}
\mathrm{ Absolute ~error} = \left| \mathrm{Analytical ~result} -\mathrm{ Simulation ~result}\right| ,
\end{align*}
when it is moved in the $y,z$-plane, while FAR$ _2 $ and FAR$ _3 $ are positioned at fixed locations. The shaded regions in gray represents the excluded locations as  the FARs overlap with each other in this region. Fig. \ref{fig:aberr} confirms that the derived equation has good accuracy in most of the regions. However, similar to \cite{sabu2020a}, the accuracy of the derived equation can be bad if the FARs in the medium are too close to each other and the FAR of interest hinders the other FARs from the line of sight of the point transmitter.\par

\begin{figure}[ht]
	\centering
	\includegraphics[width=.65\linewidth]{Figures/AbErrnew}
	\caption{Variation of the absolute error of the hitting probability of FAR$ _3 $ in a $ 3- $FAR system. Parameters: $ a=5\mu m,\ D=100\mu m^2/s,\ t=1s,\ \Delta t=10^{-4}s,\ \pmb{A}_2=[10,\	14.14,\	14.14],\  \pmb{A}_3=[10,\	14.14,\	-14.14],\  \pmb{A}_1=[10,\ y,\	z], \text{where } y,z$ is varied.
	}
	\label{fig:aberr}
	\vspace{-0.4cm}
\end{figure}

The analytical expression given in Theorem \ref{t3rx} simplifies when the FARs are equidistant from the transmitter and equidistant from each other, which we consider next.

\subsection{Special case: Symmetric system with equidistant FARs}\label{secuca}
\begin{figure}[ht]
	\centering
	\includegraphics[width=.65\linewidth]{Figures/UCA1}
	\caption{An MC system with three FARs equidistant from each other and the transmitter.}
	\label{fig:uca1}
\end{figure}

Now, we consider a symmetric system where all three FARs are equidistant  the transmitter and their mutual distance is also the same as shown in Fig. \ref{fig:uca1}. The point transmitter is at the origin and the centers of the FARs are located in a circle with center $ [w,0,0] $ and radius $ d $. Therefore the location of the FAR$ _i $ is $\pmb{A}_i= [w,d\cos\left(2\pi (i-1)/3\right),d\sin\left(2\pi (i-1)/3\right)], \ i\in\{1,2,3\} $. This type of arrangement of FARs is called uniform circular array (UCA) \cite{gursoy2019,tang2021}.
The hitting probability for this special case
is given in the following theorem.

\begin{theorem}\label{c3rxeq}

	In an MC system with three FARs that are equidistant from the transmitter at the origin and equidistant from each other (\textit {i.e., } $ r_i=r \text { and } R_{ij}=R,\ i,j\in{1,2,3},\ i\neq j$), the probability that an IM emitted from the point source at the origin  hits  $\mathrm{FAR} _1 $ within time $ t $ in the presence of other two FARs is

	\begin{align}
		\hp{1}(t)=\frac{a}{r}\sum_{n=0}^{\infty}\frac{(-2a)^n}{R^n}\mathrm{erfc}\left(\frac{r-a+n(R-a)}{\sqrt{4Dt}}\right).\label{eq3rxeq}
	\end{align}
\end{theorem}
\begin{proof}
	See Appendix \ref{a3rxeq}.
\end{proof}
\begin{coro}\label{corolim}
	The fractions of IMs that hits each of the FARs eventually (\textit {i.e., } $ t\rightarrow\infty $) in an MC system with three FARs that are equidistant from the transmitter at the origin and equidistant from each other is given by
	\begin{align}
		\hp{\infty}=\frac{a/r}{1+2a/R}.
	\end{align}
\begin{proof}
	Applying the final value theorem for the Laplace transform (\textit {i.e., } $ \lim_{t\rightarrow\infty} f(t)=\lim_{s\rightarrow0}sF(s)$) \cite{spiegel1992} in \eqref{eq3rxeq} gives Corollary \ref{corolim}.
\end{proof}
\end{coro}

Fig. \ref{fig:htd} shows the variation of hitting probability (\textit {i.e., } eq. \eqref{eq3rxeq}) to the diffusion coefficient for different values of FAR radius. The equation \eqref{eq3rxeq} is validated using particle-based simulations. As $ D $ increases, the hitting probability increases due to the faster motion of the IMs. The increase in the hitting probability with the increase in $ a $ is due to the increase in the surface area of the FARs and thereby the increase in the probability of capturing the IM.

  \begin{figure}[ht]
	\centering
	\includegraphics[width=0.65\linewidth]{Figures/HPvsD}
	\caption{Variation of the hitting probability $\hp{1}(t)$ with diffusion coefficient $ D $ for different values of $ a $ for the symmetric case. Here, $d=20 \mu m, w=10\mu m,\ r=22.36\mu m,\

		t=10^{-4}\ s,\ t=1s $.}
	\label{fig:htd}
	\vspace{-0.3cm}
\end{figure}

We now discuss two specific scenarios to demonstrate the applicability of the desired result in the next two subsections.

\subsection{Impact of malicious receivers}
Let us consider a scenario where FAR$_1$ is the intended receiver for the transmitter, and the rest of the two FARs are acting as malicious receivers.
The relative change in the IMs that were supposed to hit the intended FAR$ _1 $ within time $ t $, but are hitting the other two FARs  first is given by
	\begin{align}
		\hpdiff{1}(t)=\frac{\overline{\mathsf{p}}(t,r_1)-\hp{1}(t)}{\overline{\mathsf{p}}(t,r_1)}.\label{ediff}
	\end{align}
Note that, $ \hpdiff{1}(t) $ quantifies the relative influence of the malicious receivers on the desired receiver in capturing the IMs. By useful IMs, we mean the IMs that would reach FAR$ _1 $ and are responsible for detection. The more the $ \hpdiff{1}(t) $ is, the more the IMs captured by the malicious FARs would be. This will result in the degradation of hitting probability of the intended receiver, degrading its communication performance. To avoid the malicious FARs to degrade the desired communication, $ \hpdiff{1}(t) $ has to be minimal.

\begin{coro}
The relative change in the IMs that are supposed to hit the intended FAR eventually, but are hitting the other two FARs in the symmetric system (\textit {i.e., } all the FARs that are equidistant from the transmitter at the origin and equidistant from each other) is given
	by\begin{align}
		\hpdiff{\infty}&=

		\frac{1}{1+R/2a}.
	\end{align}
\end{coro}

Therefore, the malicious receivers' influence is minimal when $a$ is small and $ r $ and $ R $ are large. For a secure MC system, the eavesdropping by the malicious FARs is minimal when the FARs are of small radii and located far from the transmitters and each other. This can also be verified from Fig. \ref{fig:qfig} which shows the variation of this influence with respect to the angular separation between the receivers. With sufficient angular separation, effect of malicious receivers can be reduced signicantly.

  \begin{figure}[ht]
	\centering
	\includegraphics[width=0.65\linewidth]{Figures/q1}
	\caption{Variation of $  \hpdiff{1}(t) $ with angle between the intended and malicious receivers (\textit {i.e., } $ \phi_{12}=\phi_{13}=\theta $). Here, the intended receiver is located at $ \pmb{A}_1=[r,0,0] $ and the malicious receivers are located $ \pmb{A}_2=[r\cos(\theta),r\sin(\theta),0] $ and $ \pmb{A}_3=[r\cos(-\theta),r\sin(-\theta),0] $ with  $D=100 \mu m^2/s$.}
	\label{fig:qfig}.
	\vspace{-0.5cm}
\end{figure}

  \begin{figure}
	\centering
	\includegraphics[width=0.65\linewidth]{Figures/N_HP}
	\caption{Variation of the hitting probability $\hp{1}(t)$ with $ t $ in the presence of $m$ malicious receivers. Here, receivers are located according to Fig. \ref{fig:uca1} with $d=20 \mu m, w=10\mu m,\ r=22.36\mu m,\ a=4\mu m,\

		t=10^{-4}\ s,\ t=1s $.}
	\label{fig:nuca}
	\vspace{-0.5cm}
\end{figure}
Fig. \ref{fig:nuca} shows the variation of the hitting probability within time $ t $ with $ t $ for $ m$ number of malicious receivers with values $m=0,1,2$. Here, $m=0$ refers to the case when no malicious receiver is present. When the number of malicious receivers is increased, the hitting probability decreases due to the increase in their influence. The widening gap between the curves shows that the degrading effect of these receivers grows with time. Therefore, it is essential to consider the effect of this influence while designing MC systems with multiple malicious FARs.

\subsection{Impact of cooperating receivers}
Let us consider a scenario where FAR$_1$ is the intended receiver, and the rest of the two FARs cooperate in the receiving. In this case, the receiver can detect the transmission if any of the FARs have received IMs. The fractions of IMs  hitting any of the FARs within time $t$ to the fraction of IMs hitting a single FAR in the absence of other FARs  is given by
	\begin{align}
		\mathsf{s}_1(t)&=
		\frac{\sum_i\hp{i}(t)}
		{\overline{\mathsf{p}}(\infty,r_1)}.
	\end{align}
$ \mathsf{s}_1(t) $ can be considered as the array gain when multiple FARs are employed instead of a single FAR.

\begin{coro}
	The fractions of IMs eventually hitting any of the FARs  to the fraction of IMs hitting a single FAR in the absence of other FARs  in a symmetric system (\textit {i.e., } all the FARs that are equidistant from the transmitter at the origin and equidistant from each other) is given by
	\begin{align}
		\mathsf{s}_\infty&=\frac{3\hp{\infty}}{\overline{\mathsf{p}}(\infty,r)}
		=\frac{3}{1+2a/R}.
	\end{align}
\end{coro}

Note that the asymptotic array gain $\mathsf{s}_\infty$ is less than 3. Further, it decreases when FARs are placed closer to each other (\textit {i.e., } $ R $ decreases). When FARs are far apart (\textit {i.e., } $ R\rightarrow\infty $), $ \mathsf{s}\rightarrow 3 $.

Having derived the hitting probability for $N=3$ FARs, we  will now show how we can extend the proposed analysis for the general value of $N$ FARs.

\section{Channel Model for a  MC system with $N$ FARs
}
We now consider that there are $ N $ FARs coexisting in the same communication channel.
Let us first focus on FAR$_i$.
We assume that the point transmitter at the origin emits an IM at time $ t=0 $. The probability that this IM  first hits the surface of FAR$_{j}$
($j\neq i $) in the time interval $[\tau,\tau+\mathrm{d}\tau] $ is  $\displaystyle \frac{\partial \hp{j}(\tau)}{\partial \tau} \mathrm{d}\tau$. If the IM were not absorbed by FAR$ _j $, it would have eventually hit the FAR$ _i $ in the remaining $ t-\tau $ time with probability $\overline{\mathsf{p}}(t-\tau,R_{ji})$, where $R_{ji}$ is the distance between the nearest point ($ \pmb{B}_j $) on the surface of FAR$_{j} $  from the transmitter and the center of FAR$_{i} $ ($ \pmb{A}_i $). Note that, the hitting point of IM at FAR$_j$ is approximated  by $\pmb{B}_j$.

Similar to Appendix \ref{a3rx}, the hitting probability of an IM that was supposed to hit FAR$ _{i} $ within time $ t $ but is hitting the other $ N-1 $ FARs is given by
\begin{align}
	\overline{\mathsf{p}}(t,r_i) - \hp{i}(t)   =&  \sum_{j=1,\ j\neq i}^{N}
	\int_{0}^{t} \frac{\partial \hp{j}(\tau)}{\partial \tau}\overline{\mathsf{p}}(t-\tau,R_{ji}) 	\mathrm{d}\tau.\label{eqtnxr}
\end{align}
Now taking the Laplace transform of \eqref{eqtnxr} gives
\begin{align}
	\overline{\mathcal{P}}(s,r_i) -\lhp{i}(s)  =  \sum_{j=1,\ j\neq i}^{N}
	s\lhp{j}(s) \overline{\mathcal{P}}(s,R_{ji}) , \ \ \forall\ i.\label{eqsnxr}
\end{align}
Note that, \eqref{eqsnxr} forms a system of $ N $ linear equations with $ N $ unknowns for $ i\in\{1,2,\cdots,N\} $. Solving the system of linear equations gives the value of hitting probability at the $ i $th FAR, \textit {i.e., } $ \hp{i}(t) $. The solution to the equations and further analysis is left for the future work.
\vspace{-0.2cm}
\section{Conclusions}
Multiple absorbing receivers may coexist in the same communication channel in a molecular communication system. It is essential to obtain the analytical channel model for analyzing and understanding such systems. However, exact channel models are not available in the literature due to the mathematical intractability in solving the diffusion equations. In this work, we bridged this gap by deriving an approximate analytical expression for the hitting probability of an IM emitted by a point transmitter of each FARs. The derived hitting probability expression was simplified by considering that the FARs are equidistant from each other and the transmitter. Using the derived expressions, we discussed several design insights like quantifying the mutual influence of FARs, the fraction of IMs eventually hitting the FARs, and the fractions of IMs eventually hitting all of the FARs to the fraction of IMs hitting a single FAR.  We also provided an analytical framework to derive the approximate hitting probability expression for an MC system with $ N $ FARs, and further analysis will be done as future work. The derived analytical expressions can be used to analyze and design SIMO and MIMO systems to improve the reliability and data rate. It also helps design communication systems that secure the desired receiver bio-nanomachines from eavesdroppers.
\vspace{-0.2cm}
\appendices
	\section{Derivation of Fraction of IMs Absorbed at Each FAR in a $3-$FARs System}\label{a3rx}

Consider that the point transmitter at the origin emits an IM at time $ t=0 $. The probability that this IM  first hits the surface of FAR$_{2} $ in the time interval $[\tau,\tau+\mathrm{d}\tau] $ is  $\frac{\partial \hp{2}(\tau)}{\partial \tau} \mathrm{d}\tau$. If not absorbed by FAR$ _2 $, the IM would have eventually hit the FAR$ _1 $ in the remaining $ t-\tau $ time with probability $\overline{\mathsf{p}}(t-\tau,R_{21})$, where $R_{21}$ is the distance between the nearest point ($ \pmb{B}_2 $) on the surface of FAR$_{2} $  from the transmitter and the center of FAR$_{1} $ ($ \pmb{A}_1 $). Here, the hitting point of IM at FAR$_2$ is approximated  by $\pmb{B}_2$. Similarly, the IM emitted from the point source at the origin first hits the surface of FAR$ _3 $ in the time interval $[\tau,\tau+\mathrm{d}\tau] $ with probability  $\frac{\partial \hp{3}(\tau)}{\partial \tau} \mathrm{d}\tau$. If not absorbed by FAR$_3 $, the IM would have eventually hit the FAR$ _1 $ in the remaining $ t-\tau $ time  with probability $\overline{\mathsf{p}}(t-\tau,R_{31})$, where $R_{31}$ is the distance between the nearest point ($ \pmb{B}_3 $) on the surface of FAR$_{3} $  from the transmitter and the center of FAR$_{1} $ ($ \pmb{A}_1 $). Therefore,
the hitting probability of an IM that was supposed to hit FAR$ _{1} $ within time $ t $ but is hitting either FAR$ _{2} $ or FAR$ _{3} $ is given by
\begin{align}
	\overline{\mathsf{p}}(t,r_{1}) - \hp{1}(t)   =&
	\int_{0}^{t} \frac{\partial \hp{2}(\tau)}{\partial \tau}\overline{\mathsf{p}}(t-\tau,R_{21}) 	\mathrm{d}\tau
	+
	\int_{0}^{t} \frac{\partial \hp{3}(\tau)}{\partial \tau}\overline{\mathsf{p}}(t-\tau,R_{31}) \mathrm{d}\tau , \label{eqn1}
\end{align}
Similarly, the hitting probability of an IM that was supposed to hit FAR$ _{2} $ within time $ t $ but is hitting either FAR$ _{1} $ or FAR$ _{3} $, and the hitting probability of an IM that is supposed to hit FAR$ _{3} $ within time $ t $ but is hitting either FAR$ _{1} $ or FAR$ _{2} $  first is given by
\begin{align}
	\overline{\mathsf{p}}(t,r_{2}) - \hp{2}(t)   = &
	\int_{0}^{t} \frac{\partial \hp{3}(\tau)}{\partial \tau}\overline{\mathsf{p}}(t-\tau,R_{32}) \mathrm{d}\tau
	+ \int_{0}^{t} \frac{\partial \hp{1}(\tau)}{\partial \tau}\overline{\mathsf{p}}(t-\tau,R_{12}) \mathrm{d}\tau ,
	\label{eqn2}\\
\text{and }
	\overline{\mathsf{p}}(t,r_{3}) -  \hp{3}(t)   = &
	\int_{0}^{t} \frac{\partial \hp{1}(\tau)}{\partial \tau}\overline{\mathsf{p}}(t-\tau,R_{13}) \mathrm{d}\tau
	+ \int_{0}^{t} \frac{\partial \hp{2}(\tau)}{\partial \tau}\overline{\mathsf{p}}(t-\tau,R_{23}) \mathrm{d}\tau
	\label{eqn3}  ,
\end{align}
respectively. Taking Laplace transform on both sides of \eqref{eqn1}, \eqref{eqn2} and \eqref{eqn3} gives the following set of equations:
\begin{align}
	\overline{\mathcal{P}}(s,r_1) -\lhp{1}(s)  =
	s\lhp{2}(s) \overline{\mathcal{P}}(s,R_{21}) +
	s\lhp{3}(s)\overline{\mathcal{P}}(s,R_{31}),\label{eqn4}\\
	\overline{\mathcal{P}}(s,r_2) -\lhp{2}(s)   =
	s\lhp{1}(s) \overline{\mathcal{P}}(s,R_{12})
	+ s \lhp{3}(s) \overline{\mathcal{P}}(s,R_{32}),\label{eqn5}\\

	\overline{\mathcal{P}}(s,r_3) - \lhp{3}(s)   =
	s\lhp{1}(s)\overline{\mathcal{P}}(s,R_{13})
	+ s\lhp{2}(s) \overline{\mathcal{P}}(s,R_{23}).	\label{eqn6}
\end{align}
Here $\lhp{i}(s)$ and $\overline{\mathcal{P}}(s,x),\ x\in\{r_i,R_{ij}\}\ i,j\in\{1,2,3\}, i\neq j
$ are the Laplace Transforms of $\hp{i}(t)$ and $\overline{\mathsf{p}}(t,x)$ respectively. Note that,
$ \overline{\mathcal{P}}(s,x)=\frac{a}{sx}\exp\left(-\left(x-a\right)\sqrt{\frac{s}{D}}\right).\label{lapeq} $

Solving \eqref{eqn4}, \eqref{eqn5} and \eqref{eqn6} for the value of $\lhp{1}(s)$ and taking the inverse Laplace transform gives \eqref{eq3rx}.
\section{Derivation of Theorem \ref{c3rxeq}}\label{a3rxeq}
 When three FARs that are equidistant from the transmitter at the origin and equidistant from each other, $ r_i=r,\ \forall i  $ and $ R_{ij}=R,\ \forall i,j,\  i\neq j $. Therefore \eqref{lp3rx} simplifies to
\begin{align}
	\lhp{1}(s)   =
	\frac{\overline{\mathcal{P}}(s,r) -s\mathsf{a}(s)+s^2\mathsf{b}(s)}
	{1-s^2\mathsf{c}(s)+s^3\mathsf{d}(s)}, \label{abeq19}
\end{align}\\
where $\mathsf{a}(s){=}2\overline{\mathcal{P}}(s,r)\overline{\mathcal{P}}(s,R) ,\ \mathsf{b}(s){=}\overline{\mathcal{P}}(s,r)\overline{\mathcal{P}}^2(s,R),\ \mathsf{c}(s)=3\overline{\mathcal{P}}^2(s,R),\  \mathsf{d}(s){=}2\overline{\mathcal{P}}^3(s,R) $.
Simplifying \eqref{abeq19} further gives
\begin{align}
	\lhp{1}(s)  &
	=\overline{\mathcal{P}}(s,r)\frac{1}{1+2s\overline{\mathcal{P}}(s,R)}\label{l3rx}\\
	&\stackrel{(a)}{=}\overline{\mathcal{P}}(s,r)\times\sum_{n=0}^{\infty}(-2)^n\left(s\overline{\mathcal{P}}(s,R)\right)^n,
	\label{abeq20}
\end{align}
where  $ (a) $ is due to the identity $ (1+2x)^{-1}=\sum_{n=0}^{\infty}(-2)^nx^n $.
Substituting \eqref{lapeq} in \eqref{abeq20} gives
{\small\begin{align}
	\lhp{1}(s)=\frac{a}{r}\sum_{n=0}^{\infty}\frac{(-2a)^n}{sR^n}\exp\left(-(r-a+n(R-a)\sqrt{\frac{s}{D}})\right).\label{eqcor}
\end{align}}
Applying the inverse Laplace transform on both sides of \eqref{eqcor} gives Theorem \ref{c3rxeq}.

	\bibliographystyle{IEEEtran}