# Can the double-slit experiment distinguish between quantum interpretations?

Ali Ayatollah Rafsanjani<EMAIL_ADDRESS>Department of Physics,
Sharif University of Technology, Tehran, Iran  School of Physics, Institute
for Research in Fundamental Sciences (IPM), Tehran, Iran  MohammadJavad Kazemi
<EMAIL_ADDRESS>Research center for quantum engineering and photonics
technology, Sharif University of Technology, Tehran, Iran  Alireza Bahrampour
Department of Physics, Sharif University of Technology, Tehran, Iran  Research
center for quantum engineering and photonics technology, Sharif University of
Technology, Tehran, Iran  Mehdi Golshani  School of Physics, Institute for
Research in Fundamental Sciences (IPM), Tehran, Iran

###### Abstract

Despite the astonishing successes of quantum mechanics, due to some
fundamental problems such as the measurement problem and quantum arrival time
problem, the predictions of the theory are in some cases not quite clear and
unique. Especially, there are various predictions for the joint spatiotemporal
distribution of particle detection events on a screen, which are derived from
different formulations and interpretations of the quantum theory. Although the
differences are typically small, our studies show that these predictions can
be experimentally distinguished by an unconventional double-slit
configuration, which is realizable using present-day single-atom
interferometry.

##  I Introduction

In textbook quantum theory, time is a parameter in the Schrödinger equation,
not a self-adjoint operator, hence there is no unique and unambiguous way to
compute the temporal probability distribution of events from the first
principles (i.e. the Born rule) 111As shown by Pauli [105], if the Hamiltonian
spectrum is discrete or has a lower bound, then there is no self-adjoint time
operator canonically conjugate to the Hamiltonian.. Nonetheless, since clocks
exist and time measurements are routinely performed in quantum experiments [2,
3], a complete quantum theory must be able to predict the temporal statistics
of detection events. For example, in the famous double slit experiment, each
particle is detected at a random time as same as at a random position on the
detection screen [4, 5, 6, 7, 8]. Therefore, one can ask: What is the
position-time joint probability density $\mathbb{P}(\mathbf{x},t)$ on the
screen? Although this question is very old [9, 10, 11, 12], it is still open
[13, 14, 15, 16, 17, 18]. In fact, the ambiguity in the arrival time
distribution even prevents a clear prediction of cumulative arrival position
distribution, $\int\mathbb{P}(\mathbf{x},t)dt$, which is typically measured in
a non-time-resolved double-slit experiment 222Note that, the Heisenberg
position operator describes position measurement at a specific time, not
position measurements at random times [106, 107]. In the other words,
$|\psi_{t}(\mathbf{x})|^{2}$ is just the conditional position probability
density $\mathbb{P}(\mathbf{x}|t)$ [14, 15, 104], not the position-time joint
probability density $\mathbb{P}(\mathbf{x},t)$ [78, 79, 98].

Nonetheless, usual experiments are performed in the far-field (or scattering)
regime, where a semiclassical analysis is often sufficient [13, 19]. In this
analysis, it is assumed that particles move along classical trajectories, and
the arrival time distribution is computed using the quantum momentum
distribution [8, 20, 21]. However, because of the quantum backflow effect
[22], even in free space, the quantum mechanical time evolution of position
probability density is not consistent with the underlying uniform motion
assumption, especially in near-field interference phenomena [23]. In fact, due
to recent progress in the ultra-fast detectors technology (e.g. see [24, 25,
26, 27]), it will be soon possible to investigate the near-field regime, where
the semiclassical approximation breaks down and a deeper analysis would be
demanded [13, 28, 29].

To remedy this problem, based on various interpretations and formulations of
quantum theory, several attempts have been made to introduce a suitable
arrival time distribution. On the one hand, according to the (generalized)
standard canonical interpretation, the arrival distribution is considered as a
generalized observable, which is described by a positive-operator-valued
measure (POVM), satisfying some required symmetries [10, 11, 30, 31]. On the
other hand, in the realistic-trajectory-based formulations of quantum theory,
such as the Bohmian mechanics [32], Nelson stochastic mechanics [33], and many
interacting worlds interpretation [34], the arrival time distribution could be
obtained from particles trajectories [35, 7, 36, 18]. Moreover, in other
approaches, the arrival time distribution is computed via phenomenological
modeling of the detection process, such as the (generalized) path integral
formalism in the presence of an absorbing boundary [37, 38, 12, 39],
Schrödinger equation with complex potential or absorbing boundary [40, 41, 42,
43, 44], and so on [45, 46, 47].

In principle, the results of these approaches are different. However, in most
of the experimental situations, the differences are typically slight, and so
far as we know, in the situation where differences are significant, none of
the proposals have been backed up by experiments in a strict manner [36, 8].
An experiment that can probe these differences would undoubtedly enrich our
understanding of the foundations of quantum mechanics. The purpose of the
present paper is to make it evident, via numerical simulations, that the
famous two-slit experiment could be utilized to investigate these differences
if we simply use a horizontal screen instead of a vertical one: see Fig. 1.
Using current laser cooling and magneto-optical trapping technologies, this
type of experiment can be realized by Bose-Einstein condensates, as a
controllable source of coherent matter waves [48, 49, 50]. Moreover, our
numerical study shows that the required space-time resolution in particle
detection is achievable using fast single-atom detectors, such as the recent
delay-line detectors described in [51, 52] or the detector used in [53, 6].

The structure of this paper is as follows: In Section II, we study the main
proposed intrinsic arrival distributions. Then, in section III we compare them
in the double-slit setup with vertical and horizontal screens and in different
detection schemes. In Section IV, we study the screen back-effect, and we
summarize in section V.

##  II “Intrinsic” arrival distributions

In this section, we first review the semi-classical approximation and then
scrutinize two main proposed intrinsic arrival time distributions [36, 16] and
their associated screen observables. In these approaches, the effect of the
detector’s presence on the wave function evolution, before particle detection,
is not considered. We discuss this effect in section IV.

###  II.1 Semiclassical approximation

As mentioned, in the experiments in which the detectors are placed far away
from the support of the initial wave function (i.e. far-field regime), the
semiclassical arrival time distribution is routinely used to the description
of the particle time-of-flight [54, 21, 55, 56, 57]. In this approximation, it
is assumed that particles move classically between the preparation and
measurement. In this approach, the arrival time randomness is understood as a
result of the uncertainty of momentum, and so the arrival time distribution is
obtained from momentum distribution [13, 58, 17, 36]. In the one-dimensional
case, the classical arrival time is given by

$t=m(L-x_{0})/p_{0},$ (1)

which is applicable for a freely moving particle of mass $m$ that at the
initial $t\\!=\\!0$ had position $x_{0}$ and momentum $p_{0}$ arriving at a
distant point $L$ on a line. Hence, for a particle with the momentum wave
fuction ${\tilde{\psi\mkern 3.0mu}\mkern-3.0mu}{}_{0}(p)$, assuming $\Delta
x_{0}\\!\ll\\!|L-\langle x\rangle_{0}|$, the semiclassical arrival time
distribution reads [58]

$\Pi_{\text{SC}}(t|x\\!=\\!L)=\frac{mL}{t^{2}}|{\tilde{\psi\mkern
3.0mu}\mkern-3.0mu}{}_{0}(mL/t)|^{2}.$ (2)

This analysis could be generalized in three-dimensional space. Then, the
distribution of arrival time at a screen surface $\mathbb{S}$ is given by [36]

$\Pi_{\text{SC}}(t|\mathbf{x}\\!\in\\!\mathbb{S})=\frac{m^{3}}{t^{4}}\int_{\mathbb{S}}|{\tilde{\psi\mkern
3.0mu}\mkern-3.0mu}{}_{0}(m\mathbf{x}/t)|^{2}\,\mathbf{x}\cdot d\mathbf{S},$
(3)

where the $d\mathbf{S}$ is the surface element directed outward. The other
main distribution that should be demanded is the joint position-time
probability distribution on the screen, which is also called ”screen
observable” [11]. Using the conditional probability definition, the joint
probability of finding the particle in $dS$ and in a time interval
$[t,t\\!+\\!dt]$ could be written as
$\mathbb{P}(\mathbf{x},t|\mathbf{x}\\!\in\\!\mathbb{S})dSdt=\left[\Pi(t|\mathbf{x}\\!\in\\!\mathbb{S})dt\right]\times\left[\mathbb{P}(\mathbf{x}|\mathbf{x}\\!\in\\!\mathbb{S},t)dS\right].$
In this regard, one can use the fact that $\psi_{t}(\mathbf{x})$ is the state
of the system, conditioned on the time being $t$ in the Schrödinger picture.
This implies that $|\psi_{t}(\mathbf{x})|^{2}$ refers to the position
probability density conditioned at a specific time $t$ [14, 15, 59].
Therefore, in the semiclassical approximation, the joint spatiotemporal
probability density reads as

$\mathbb{P}_{\text{SC}}(\mathbf{x},t|\mathbf{x}\\!\in\\!\mathbb{S})=N_{\text{SC}}\Pi_{SC}(t|\mathbf{x}\\!\in\\!\mathbb{S})\left|\psi_{t}(\mathbf{x})\right|^{2}$
(4)

in which
$N_{\text{SC}}\\!\equiv\\!1/\int_{\mathbb{S}}dS\left|\psi_{t}(\mathbf{x})\right|^{2}$
is the normalization constant, and $dS\\!=\\!\mathbf{n}\cdot d\mathbf{S}$,
where ${\mathbf{n}}$ is the outward unit normal vector at
$x\\!\in\\!\mathbb{S}$.

###  II.2 “Standard” approach

The first attempts to investigate the arrival time problem, based on the
standard rules of quantum theory, were made at the beginning of the 1960s by
Aharonov and Bohm [60], and also Paul [61]. This approach starts with a
symmetric quantization of classical arrival time expression (1), as follows
[62]:

$\hat{t}_{AB}=mL\,\hat{p}^{\,-1}-\frac{m}{2}\,(\hat{p}^{\,-1}\,\hat{x}\,+\,\hat{x}\,\hat{p}^{\,-1}),$
(5)

where $\hat{x}$ and $\hat{p}\\!=\\!-i\,\partial/\partial x$ are the usual
position and momentum operators, respectively, and $\hat{t}_{AB}$ is called
the Aharonov-Bohm time operator. This operator satisfies the canonical
commutation relation with the free Hamiltonian operator,
$[\hat{t}_{AB},\hat{p}^{2}/2m]\\!=\\!i\hbar$, which has been used to establish
the energy-time uncertainty relation [63, 64]. However, although
$\hat{t}_{AB}$ is Hermitian (or symmetric in mathematics literature), it is
not a self-adjoint operator [65]—a fact that is in agreement with Pauli’s
theorem [1]. The origin of this non-self-adjointness can be understood as a
result of the singularity at $p\\!=\\!0$ in the momentum representation,
$\hat{t}_{AB}\to(i\hbar m/2)(p^{-2}-2p^{-1}\partial_{p})$ [65]. Nevertheless,
although the (generalized) eigenfunctions of $\hat{t}_{AB}$ are not
orthogonal, they constitute an over-complete set and provide a POVM, which are
used to define the arrival-time distribution as follows [63, 65]:

$\\!\\!\Pi_{\text{STD}}(t|x\\!=\\!L)\\!=\\!\frac{1}{2\kern
1.00006pt\pi\hbar}\\!\sum_{\alpha\kern 0.70004pt=\kern
0.70004pt\pm}\\!\left|\int_{-\infty}^{\infty}\\!\\!\\!\\!\\!\\!dp\,\theta(\alpha
p)\,\sqrt{\frac{|p|}{m}}{\tilde{\psi\mkern
3.0mu}\mkern-3.0mu}{}_{t}(p)e^{\frac{i}{\hbar}Lp}\right|^{2}\\!\\!\\!\\!,$ (6)

where $\theta(\cdot)$ is Heaviside’s step function and ${\tilde{\psi\mkern
3.0mu}\mkern-3.0mu}{}_{t}(p)$ is the wave function in the momentum
representation which could be obtained from the initial wave function
${\tilde{\psi\mkern 3.0mu}\mkern-3.0mu}{}_{0}(p)$, as ${\tilde{\psi\mkern
3.0mu}\mkern-3.0mu}{}_{t}(p)={\tilde{\psi\mkern
3.0mu}\mkern-3.0mu}{}_{0}(p)\exp[-\,itp^{2}/2m\hbar].$ The distribution
$\Pi_{\text{STD}}$ and its generalization in the presence of interaction
potential have been referred to as the ”standard arrival-time distribution” by
some authors [66, 67, 68, 69, 16]. In fact, Grot, Rovelli, and Tate treated
the singularity of (5) by symmetric regularization and obtained equation (6)
via the standard Born rule [64]. The generalizations of equations (5) and (6)
in the presence of interaction potential have been investigated in various
works [70, 71, 72, 73, 74, 75, 31, 16]. Using these developments, it has been
shown that the non-self-adjointness of the free arrival time operator can also
be lifted by spatial confinement [76, 71], and the above arrival time
distribution could be derived from the limit of the arrival time distribution
in a confining box as the length of the box increases to infinity [72].
Furthermore, recently, the distribution (6) is derived from a space-time-
symmetric extension of non-relativistic quantum mechanics [77].

The three-dimensional generalization of (6) is derived by Kijowski’s [10] via
an axiomatic approach. The assumed axioms are implied by the principle of the
probability theory, the mathematical structure of standard quantum mechanics,
and the Galilei invariance [78]. Based on these axioms, Kijowski constructed
the following arrival time distribution for a free particle that passes
through a two-dimensional plane $\mathbb{S}$ as

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\Pi_{\text{STD}}(t|\mathbf{x}\in\mathbb{S})\\!\\!\\!\\!\\!\\!$
$\displaystyle=\frac{1}{2\kern 1.00006pt\pi\hbar\kern
1.00006pt}\\!\sum_{\alpha\kern 0.70004pt=\kern
0.70004pt\pm}\int_{\mathbb{R}^{2}}\\!\\!d^{2}\mathbf{p}_{\|}$ (7)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\times\,\left|\int_{-\infty}^{\infty}\\!\\!\\!dp_{\bot}\,\theta(\alpha\mathbf{p}.\mathbf{n})\sqrt{\frac{|\mathbf{p}_{\bot}|}{m}}~{}{\tilde{\psi\mkern
3.0mu}\mkern-3.0mu}{}_{t}(\mathbf{p})e^{\frac{i}{\hbar}\mathbf{x}.\mathbf{p}_{\bot}}\kern
1.00006pt\right|^{2}\\!\\!,$

where $\mathbf{p}_{\bot}\\!\equiv\\!\\!(\mathbf{p}\,.\,\mathbf{n})\mathbf{n}$
and $\mathbf{p}_{\|}\\!\equiv\\!\mathbf{p}-\mathbf{p}_{\bot}$ are
perpendicular and parallel components of $\mathbf{p}$ relative to $\mathbb{S}$
respectively, and $\mathbf{n}$ is the outward normal of plane $\mathbb{S}$. In
fact, he first proves the above expression for the wave functions whose
supports lie in the positive (or negative) amounts of ${p}_{\bot}$. Then he
uniquely derives the following self-adjoint variant of the (three-dimensional
version of) Aharonov-Bohm arrival time operator, by demanding that the time
operator be self-adjoint and leads to (7) for these special cases via the Born
rule [10, 78]:

$\hat{t}_{L}=\text{sgn}(\hat{p}_{\bot})\left[mL\hat{p}_{\bot}^{-1}-\frac{m}{2}(\hat{p}_{\bot}^{-1}\hat{x}_{\bot}+\hat{x}_{\bot}\hat{p}_{\bot}^{-1})\right],$
(8)

where $\hat{x}_{\bot}\equiv\hat{\mathbf{x}}.\mathbf{n}$ and $L$
($\equiv\mathbf{x}.\mathbf{n}$) represent the distance between the detection
surface and the origin 333The presence of $\text{sgn}(\hat{p}_{\bot})$
operator ensures the self-adjointness of this time operator, however, leads to
a modified commutation relation, i.e. $[\hat{t}_{K},\hat{H}]=i\hbar\
\text{sgn}(\hat{p}_{\bot})$. . Finally, for an arbitrary wave function, the
equation (7) could be derived from this self-adjoint operator. Moreover,
considering this time operator, besides the components of the position
operator in the detection plane,
$\hat{\mathbf{x}}_{\|}\equiv\hat{\mathbf{x}}-(\hat{\mathbf{x}}.\mathbf{n})\mathbf{n}$,
Kijowski obtains the following expression as the joint position-time
distribution on the detection screen via the Born rule [78]:

$\mathbb{P}_{\text{STD}}(\mathbf{x},t|\mathbf{x}\\!\in\\!\mathbb{S})=\sum_{\alpha=\pm}|\psi^{\alpha}_{\mathbb{S}}(\mathbf{x},t)|^{2},$
(9)

in which $\psi^{\pm}_{\mathbb{S}}(\mathbf{x},t)$ is the wave function on the
basics of joint eigenstates of the operators $\hat{t}_{L}$ and
$\hat{\mathbf{x}}_{\|}$. Explicitly

$\\!\\!\\!\\!\\!\psi^{\pm}_{\mathbb{S}}(\mathbf{x},t)=\frac{1}{(2\pi\hbar)^{3/2}}\int
d^{3}\mathbf{p}\,\theta(\pm\mathbf{p}.\mathbf{n})\,\sqrt{\frac{|\mathbf{p}_{\bot}|}{m}}{\tilde{\psi\mkern
3.0mu}\mkern-3.0mu}{}_{t}(\mathbf{p})e^{\frac{i}{\hbar}\mathbf{x}.\mathbf{p}}.$
(10)

Note that, the arrival time distribution (7) could be reproduced by taking the
integral of (9) over the whole of the screen plane. The joint space-time
probability distribution (9), and its generalization for the particles with
arbitary spin, have been also derived by Werner in another axiomatic manner
[11]. Moreover, it is easy to see that the results (7) and (9) can be obtained
from a regularized version of the (three-dimensional generalization of)
Aharonov-Bohm time operator, which is the same as the procedure used by Grot,
Rovelli and Tate in one-dimensional cases [64].

###  II.3 Quantum flux and Bohmian approach

Inspiring by classical intuition, another proper candidate for screen
observables is the perpendicular component of the quantum probability current
to the screen surface, $\mathbf{J}(\mathbf{x},t).\mathbf{n}$, where

$\mathbf{J}(\mathbf{x},t)=-\frac{\hbar}{m}\imaginary\left[\psi_{t}^{*}(\mathbf{x})\gradient\psi_{t}(\mathbf{x})\right],$
(11)

and $\mathbf{n}$ is the outward normal to the screen $\mathbb{S}$. This
proposal is applicable for a particle in a generic external potential and a
generic screen surface, not necessarily an infinite plane. There are several
attempts to derive this proposal in various approaches, such as Bohmian
mechanics for the scattering case in [79], decoherent histories approach in
[80] as an approximation, or in [81] as an exact formula using the concept of
extended probabilities, and so on [45, 46, 82]. Howover, even if the wave
function contains only momentum in the same direction as $\mathbf{n}$, the
$\mathbf{J}(\mathbf{x},t)\cdot\mathbf{n}$ could be negative due to the
backflow effect [22]. This property is incompatible with the standard notion
of probability. Nevertheless, this problem could be treated from the Bohmian
point of view: Using Bohmian trajectories, it can be shown that the positive
and negative values of $\mathbf{J}(\mathbf{x},t)\cdot\mathbf{n}$ correspond to
the particles that reach the point $\mathbf{x}$ at $\mathbb{S}$ in the same
direction of $\mathbf{n}$ or the opposite direction of it, respectively [83,
84]. In this regard, through the Bohmian mechanics in one-dimension, Leavens
demonstrates that the time distribution of arrival to $x\\!=\\!L$ from both
sides could be obtained from the absolute form of probability flux as [85, 35]

$\Pi_{\text{QF}}(t|x\\!=\\!L)=\frac{|J(L,t)|}{\int dt\,|J(L,t)|},$ (12)

which is free from the aforementioned problem.

The three-dimensional justification of
$\mathbf{J}(\mathbf{x},t)\cdot\mathbf{n}$ as an operational formulation of the
arrival time model has been made in [82]. Also, the generalization of (12) for
arrival to the surface $\mathbb{S}$ is given by [13, 16, 7, 86]

$\Pi_{\text{QF}}(t|\mathbf{x}\\!\in\\!\mathbb{S})=\frac{\int_{\mathbb{S}}dS|\mathbf{J}(\mathbf{x},t)\\!\cdot\\!\mathbf{n}|}{\int
dt\int_{\mathbb{S}}dS|\mathbf{J}(\mathbf{x},t).\mathbf{n}|},$ (13)

with $dS\\!=\\!\mathbf{n}\\!\cdot\\!d\mathbf{S}$ the magnitude of the surface
element $d\mathbf{S}$ which is directed outward at
$\mathbf{x}\\!\in\\!\mathbb{S}$. To illustrate (13) and to generalize it to
the case of joint arrival distribution, we can use the Bohmian point of view.
In this theory, each particle has a specific trajectory, depending on the
initial position, and so the rate of passing particles through an area element
$d\mathbf{S}$ centered at $\mathbf{x}\\!\in\\!\mathbb{S}$, in the time
interval between $t$ and $t+dt$, is proportional to
$\rho_{t}(\mathbf{x})|\mathbf{v}(\mathbf{x},t)\\!\cdot\\!d\mathbf{S}|dt$,
where
$\mathbf{v}(\mathbf{x},t)\\!=\\!\mathbf{J}(\mathbf{x},t)/|\psi_{t}(\mathbf{x})|^{2}$
is the Bohmian velocity of the particle. Hence, using quantum equilibrium
condition [87, 88], $\rho_{t}(\mathbf{x})\\!=\\!|\psi_{t}(\mathbf{x})|^{2}$,
and accomplishing normalization, the joint arrival distribution could be
represented by the absolute value of the current density as

$\mathbb{P}_{\text{QF}}(\mathbf{x},t|\mathbf{x}\\!\in\\!\mathbb{S})=\frac{|\mathbf{J}(\mathbf{x},t)\\!\cdot\\!\mathbf{n}|}{\int
dt\int_{\mathbb{S}}dS|\mathbf{J}(\mathbf{x},t)\\!\cdot\\!\mathbf{n}|}.$ (14)

Now, by integrating (14) over all $\mathbf{x}\\!\in\\!\mathbb{S}$, we arrive
at the three-dimensional arrival time distribution (13) for the screen surface
$\mathbb{S}$. It should be noted that Eq. (14) is not necessarily followed for
an ensemble of classical particles because a positive or negative current at a
space-time point, $(\mathbf{x},t)$, can in general have contributions from all
the particles arriving to $\mathbf{x}$ at $t$ from any direction. Nonetheless,
since the Bohmian velocity field is single-valued, the particle trajectories
cannot intersect each other at any point of space-time and so only a single
trajectory contributes to the current density $\mathbf{J}(\mathbf{x},t)$ at
the particular space-time point $(\mathbf{x},t)$. Moreover, this fact implies
that when $\mathbf{v}(\mathbf{x},t)\cdot\mathbf{n}\\!>\\!0$ we can say that
the trajectory and consequently the particle has passed through the screen
from the inside and vice versa for
$\mathbf{v}(\mathbf{x},t)\cdot\mathbf{n}\\!<\\!0$. Hence, one can define the
joint probability distribution for the time of arrival to each side of
$\mathbb{S}$ as

$\mathbb{P}^{\pm}_{\text{QF}}(\mathbf{x},t|\mathbf{x}\\!\in\\!\mathbb{S})=\frac{\mathbf{J}^{\pm}(\mathbf{x},t)\\!\cdot\\!\mathbf{n}}{\int
dt\int_{\mathbb{S}}dS\,\,\mathbf{J}^{\pm}(\mathbf{x},t)\\!\cdot\\!\mathbf{n}},$
(15)

where
$\mathbf{J}^{\pm}({\mathbf{x}},t)=\pm\,\theta(\pm\mathbf{J}\\!\cdot\\!\mathbf{n})\mathbf{\,}\mathbf{J}(\mathbf{x},t)$.
In addition, note that there may be some trajectories which cross $\mathbb{S}$
more than once—and we have multi-crossing trajectories (see the typical
Bohmian trajectory in Fig. 1). The course of the above inference to Eq. (14)
was in such a manner that multi-crossing trajectories could contribute several
times (see Fig. 2 (a)). However, one could assume the detection surface as a
barrier that does not allow the crossed particle to return inside (see Fig. 2
(c)). In this case, it is suggested to use the truncated current defined as

${\tilde{\mathbf{J}}}({\mathbf{x}},t)\coloneqq\left\\{\begin{aligned}
&{\mathbf{J}}({\mathbf{x}},t)&&\mbox{if $({\mathbf{x}},t)$ is a first exit
through $\mathbb{S}$}\\\ &0&&\mbox{otherwise}\end{aligned}\right.$ (16)

where $({\mathbf{x}},t)$ is a first exit event through the boundary surface
$\mathbb{S}$, if the trajectory passing through $\mathbf{x}$ at time $t$
leaves inside $\mathbb{S}$ at this time, for the first time since $t\\!=\\!0$
[79, 89, 13]. The limiting condition in (16), imposes that the joint
probability distribution based on it should be computed numerically using
trajectories:

$\tilde{\mathbb{P}}_{\text{QF}}(\mathbf{x},t|\mathbf{x}\\!\in\\!\mathbb{S})=\frac{\tilde{\mathbf{J}}(\mathbf{x},t)\\!\cdot\\!\mathbf{n}}{\int
dt\int_{\mathbb{S}}dS\,\,\tilde{\mathbf{J}}(\mathbf{x},t)\\!\cdot\\!\mathbf{n}}.$
(17)

Of course, the detection screen is not always a barrier-like surface (see Fig.
2 (b)), and one could assume that there is a point-like detector that lets the
multi-crossing trajectories to contribute to the distribution and we can use
(14) in such cases.

Figure 1: Schematic double-slit experiment setup. The center of two slits is
considered as the coordinate origin, and the vertical and horizontal screens
are placed at $x\\!=\\!L_{x}$ and $y\\!=\\!L_{y}$, respectively. The dashed
black line shows a typical Bohmian trajectory that arrives at the horizontal
screen. A suitable single-particle detector, in addition to particle arrival
position, can record the arrival time using a proper clock. Figure 2:
Different schemes of particle detection on the screen surface $\mathbb{S}$. In
the Bohmian point of view, particles could have a recursive motion on surface
$\mathbb{S}$ and cross it more than once (e.g. see the trajectory that plotted
in Fig. 1). Assuming different detector types, one can prob variant possible
observables on the screen. In panel (a) a conceivable particle trajectory is
depicted, which crosses $\mathbb{S}$ three times. In this panel, a movable
point-like detector is placed on $\mathbb{S}$, which can survey the whole
screen and detect particles that arrive only from one side, while in panel (b)
a two-sided point detector is placed on $\mathbb{S}$, which can move along it
and detect particles that arrive from up and down. In addition, one can assume
there is (c) an array of side-by-side detectors covering the entire screen
surface $\mathbb{S}$. The last configuration blocks the trajectory and does
not allow the crossed particle to return. In this scheme, we only detect
first-arrivals from one side.

##  III “Intrinsic” screen observable in two-slit experiment

In this section, we study the discussed proposals in the previous section for
the double-slit experiment. We compare the results of these proposals in the
cases of vertical and horizontal screens (see Fig. 1), and also in different
detection schemes. The main motivation for the study of the horizontal screen
is the non-classical particles’ motions along the $y$-direction, in the
Bohmian perspective; see a typical Bohmian trajectory in Fig. 1. This behavior
is due to changing the sign of the probability current’s component in the
$y$-direction. This behavior does not occur for $x$-component of $\mathbf{J}$
and consequently for the Bohmian motion of a particle along the $x$-direction.

As shown in Fig. 1, the setup contains two identical slits at $y\\!=\\!\pm s$,
and screens are placed at $x\\!=\\!L_{x}$ and $y\\!=\\!L_{y}$ correspond to
the vertical and horizontal screens, respectively. To avoid the mathematical
complexity of Fresnel diffraction at the sharp-edge slits, it is supposed that
the slits have soft edges that generate waves having identical Gaussian
profiles in the $y$-direction. So, for each slit, we can take the wave
function as an uncorrelated two-dimensional Gaussian wave packet, which in
each dimension has the form

$\displaystyle\psi^{(i)}_{G}$ $\displaystyle(x,t)=(2\pi
s_{t}^{2})^{\textrm{-}\frac{1}{4}}\,\text{exp}\left\\{\frac{(x-x^{(i)}_{0}-u_{x}t)^{2}}{4\sigma_{0}s_{t}}\right\\}$
$\displaystyle\times\text{exp}\left\\{\frac{i}{\hbar}mu_{x}(x-x_{0}^{(i)}-\frac{u_{x}t}{2})\right\\}\
\ (i=1,2),$ (18)

with $m$ the particle’s mass, $\sigma_{0}$ the initial dispersion, $u_{x}$ the
wave packet’s velocity, $x_{0}^{(i)}$ the initial position of wave packet or
in other words the location of $i$-th slit, and $s_{t}=\sigma_{0}(1+i\hbar
t/(2m\sigma_{0}^{2}))$. Therefore, when the particle passes through the slits,
we have the total wave function as

$\psi(x,y,t)=\frac{1}{\sqrt{2}}[\psi_{G}^{(1)}(x,t)\psi_{G}^{(1)}(y,t)+\psi_{G}^{(2)}(x,t)\psi_{G}^{(2)}(y,t)],$
(19)

where superscripts (1) and (2) correspond to upper and lower slits,
respectively. This form of Gaussian superposition state is commonly used in
the literature [90, 91, 7, 92, 93] and is feasible to implement by quantum
technologies because such a state could be produced and controlled readily
[94, 95], even without using slits [49]. In this paper, we have chosen the
metastable helium atom, with mass $m=6.64\times 10^{-27}$ kg, as the
interfering particle, and the parameters as $s=10\,\mu$m,
$\sigma_{x}=0.04\,\mu$m, $\sigma_{y}=0.5\,\mu$m, $u_{x}=3$ m/s, and $u_{y}=0$
m/s. These values are feasible according to the performed experiments [96].
Moreover, the meta-stable helium atom could be detected with high efficiency
because of its large internal energy [52, 97].

###  III.1 Vertical screen

The arrival time distribution for the vertical screen placed at different
distances from the two-slit is shown in Fig. 3. As one can see this
distribution is the same for all methods, and their average arrival time is
close to the corresponding quantity in classical uniform motion. To calculate
the mean time of arrival to the screen, we use the arrival time distribution
of each method presented in sec II, i.e., Eq. (3), (7) and (13), and we have

$\bar{t}_{\mathbb{S}}=\int_{0}^{\infty}\\!\\!\\!dt\,\,\Pi(t|\mathbf{x}\\!\in\\!\mathbb{S})\,t,$
(20)

as the mean arrival time at the surface $\mathbb{S}$. Furthermore, we can
compute the average arrival time to each point on the screen using the joint
probability distribution as

$\bar{t}_{\mathbf{x}}=\frac{\int_{0}^{\infty}dt\,\,\mathbb{P}(\mathbf{x},t|\mathbf{x}\\!\in\\!\mathbb{S})\,t}{\int_{0}^{\infty}dt\,\,\mathbb{P}(\mathbf{x},t|\mathbf{x}\\!\in\\!\mathbb{S})}.$
(21)

This observable is depicted in Fig. 4-b for a vertical screen placed at
$L_{x}\\!=\\!300$ mm. Apparently, the results of the standard and quantum flux
methods are the same and similar to one that resulted in [7] by Nelson’s
mechanics. Nevertheless, they are different from the semiclassical
approximation. However, when the interference pattern is calculated by either
method, we see that their predicted cumulative position distributions do not
differ much from the others (Fig. 4-a). This observable can be calculated by
using the joint distribution as

Figure 3: Arrival time distributions of particles that arrive at the vertical
screen of the double-slit experiment at different screen distances. Figure 4:
(a) The cumulative arrival position distribution, Eq. (22), for the vertical
screen at $L_{x}\\!=\\!300$ mm, and (b) the average arrival time at each point
of the screen, Eq. (21).

$P(\mathbf{x}|\mathbf{x}\\!\in\\!\mathbb{S})\\!=\\!\frac{\int_{0}^{\infty}dt\,\,\mathbb{P}(\mathbf{x},t|\mathbf{x}\\!\in\\!\mathbb{S})}{\int_{0}^{\infty}dt\int_{\mathbb{S}}dS\,\,\mathbb{P}(\mathbf{x},t|\mathbf{x}\\!\in\\!\mathbb{S})}.$
(22)

As mentioned, it should be noted that, $|\psi_{t}(\mathbf{x})|^{2}$ is just
the conditional position probability density at the specific time $t$, not the
position-time joint probability density and so the accumulated interference
pattern, $P(\mathbf{x}|\mathbf{x}\\!\in\\!\mathbb{S})$, is not given by $\int
dt|\psi_{t}(\mathbf{x})|^{2}$ [98].

###  III.2 Horizontal screen

In this section, we are going to compare the mentioned proposals in the
double-slit setup with a horizontal detection screen (see Fig. 1). In this
regard, in Fig. 5, the arrival time distributions at the screen are plotted
for some horizontal screens which are located at $L_{y}\\!=\\!15$, $\,20$,
$\,25$, and $\,30\,\mu$m. In this figure, solid-black, dashed-green, and dash-
dotted-blue curves represent the distributions $\Pi_{STD}$, $\Pi_{QF}$ and
$\Pi_{SC}$ respectively. Also, the vertical lines show the average time of
arrival to the screen, $\bar{t}_{\mathbb{S}}$, associated with these arrival
time distributions. From this figure, one can see that, although the averages
almost coincide, the distributions are distinct. Moreover, as expected, when
the screen’s distance from the center of the two slits $L_{y}$ decreases, the
difference between distributions increases. Most of these differences occur in
the early times, which are associated with the particles that arrive at the
$\mathbb{S}$ in the near field. Furthermore, we observe that the $\Pi_{SC}$
behaves quite differently from $\Pi_{QF}$ and $\Pi_{STD}$. The distributions
$\Pi_{QF}$ and $\Pi_{STD}$ are more or less in agreement, however, for the
screen that is located at $L_{y}\\!=\\!15\,\mu$m, a significant difference
between the standard and quantum flux distributions occurs around
$t\\!\approx\\!0.2$ ms.

Figure 5: Arrival time distributions of particles that arrive on the
horizontal screen at four different distances from the center of two slits.
The vertical lines show the average arrival time.

To have a more comprehensive insight, we can look at the joint spatiotemporal
arrival distributions in Fig. 6. In this figure, joint distributions,
$\mathbb{P}_{\text{SC}}$, $\mathbb{P}_{\text{STD}}$ and
$\mathbb{P}_{\text{QF}}$ are plotted in three panels, for the horizontal
screen surface located at $L_{y}=15\,\mu$m. These density plots clearly
visualize differences between the mentioned arrival distribution proposals. In
these plots, we can see separated fringes with different shapes, which this
fact imply that the particles arrive at the screen in some detached space-time
regions. In the insets, one can see that the shapes of these regions are
different for each proposal. In the joint density of the semiclassical
approximation (Fig.6-a), fringes are well-separated, while the standard
distribution (Fig. 6-b) exhibits more continuity in its fringes. In addition,
in the pattern of the quantum flux proposal (Fig. 6-c) there are grooves
between every two fringes which is due to changing the sign of
$\mathbf{J}(\mathbf{x},t)\cdot\mathbf{n}$ in (14). In all panels of Fig.6, the
duration of “temporal no-arrival windows” between every two typical fringes
variate in the range between $0.01$ and $0.2$ ms which has a spatial extension
of about $0.3$ to $2$ mm. These space-time scales are utterly amenable
empirically by current technologies [96, 53], which could be used to test
these results.

Figure 6: Density plots of joint arrival position-time distributions for
particles that arrive at the horizontal screen of the double-slit experiment.
Panels (a), (b), and (c) represent $\mathbb{P}_{\text{SC}}$,
$\mathbb{P}_{\text{STD}}$ and $\mathbb{P}_{\text{QF}}$, respectively. Insets:
Magnified contour plots of the joint distributions. Figure 7: The space-time
arrival statistics for the double-slit experiment with a horizontal screen
placed at $L_{y}\\!=\\!15\,\mu$m. Panel (a) represents the average time of
arrival at each point of the screen, $\bar{t}_{\mathbf{x}}$. Panel (b)
represents the cumulative position probability density. The panels (c)-(f)
show the local arrival time probability densities,
$\Pi_{\mathbf{x}}(t|\mathbf{x}\\!\in\\!\mathbb{S})$, at the at the points
$x\\!=\\!16.2,\,17.4,\,18.4,\,19.2$ mm on the screen, which are chosen from
the gray region in panel (b). The vertical lines in these panels represent the
average arrival times. Panel (g) is Histograms associated with probability
densities of panel (f), which are generated by $10^{4}$ numerical random
sampling.

The average time of arrival to each point of the screen and cumulative
position interference pattern could be calculated as in the vertical screen
case by Eqs. (21) and (22). In Fig. 7(a)-(b), these two quantities are shown
for the horizontal screen which is placed at $y\\!=\\!15\,\mu$m. In contrast
to the vertical screen, the cumulative position distribution of the
semiclassical approximation is entirely separate from the two other proposals.
The cumulative position distribution resulting from standard and quantum flux
approaches have obvious differences from each other, as well. As one can see
in Fig. 7(b), the average arrival times are the same for all three methods at
first and begin to deviate from each other at $x\\!\approx\\!5$ mm; then
again, these curves converge to each other at $x\\!\approx\\!25$ mm,
approximately. The maximum deviation between the standard and quantum flux
average arrival time occurs at $x\\!\approx\\!19$ mm, which is quite in the
far-field regime—the width of the initial wave function is $\sim O(10^{-3})$mm
which is smaller than $19$ mm. Therefore one can suggest the average arrival
time in the gray region of Fig. 7(b) as a practical target for comparing these
approaches experimentally. To this end, we study arrival time distributions at
some points of this region as local arrival distributions. The arrival time
distribution conditioned at a specific point $\mathbf{x}$ on the screen can be
obtained as follow

$\Pi_{\mathbf{x}}(t|\mathbf{x}\\!\in\\!\mathbb{S})=\frac{\mathbb{P}(\mathbf{x},t|\mathbf{x}\\!\in\\!\mathbb{S})}{\int_{0}^{\infty}dt\,\,\mathbb{P}(\mathbf{x},t|\mathbf{x}\\!\in\\!\mathbb{S})}.$
(23)

Using the associated joint distribution of each proposal, we have plotted Fig.
7(c)-(f) that show $\Pi_{\mathbf{x}}(t|\mathbf{x}\\!\in\\!\mathbb{S})$ at the
positions $x\\!=\\!16.2,\,17.4,\,18.4,\,19.2$ mm, on the screen placed at
$L_{y}\\!=\\!15\,\mu$m. The broken black curves in Fig. 7 (c)-(f), resulting
from the quantum flux proposal, against the smooth curves of the other two
methods could be understood as the result of the changing the signature of the
$y$-component of the probability current: Note that, quantum flux distribution
is given by the absolute value of the probability current. The origin of
distinctions between the local average arrival times is more perceptible from
these local arrival distributions. In principle, these distributions could be
probed using fast and high-resolution single-atom detectors [53, 97]. In
particular, the delay-line detector that is recently developed by Keller et
al. [51] seems suitable for our purpose: It has the capability to resolve
single-atom detection events temporally with $220$ ps and spatially with
$177\mu$m at rates of several $10^{6}$ events per second.

We estimate by a numerical investigation that these local arrival
distributions could be well reconstructed from about $10^{4}$ number of
detection events. As an example, in Fig. 7, the histograms associated with the
probability densities of the panel (f) are plotted in panel (g), using
$10^{4}$ numerical random sampling. It is easy to estimate that the recording
of $10^{4}$ particle detection events can determine the local average arrival
time with a statistical error of about $10^{-2}$ms, while the differences
between local average arrival times of various proposals are almost bigger
than $10^{-1}$ms. Using cumulative position distribution, Fig. 7(b), one can
estimate that, if the total number of particles that arrived at the screen is
about $10^{8}$, we have about $10^{4}$ particles around $x=19.2$ mm, in the
spacial interval $(19.1,19.3)$. Using recent progress in laser cooling and
magneto-optical trapping [97], the preparation of a coherent ensemble of
metastable helium atoms with this number of particles is quite achievable
[51].

One might be inclined to think that the difference between the quantum flux
and standard average arrival times is just due to changing the signature of
$\mathbf{J}(\mathbf{x},t)\cdot\mathbf{n}$, but in the following, we show that
even without the contribution of the negative part of
$\mathbf{J}(\mathbf{x},t)\cdot\mathbf{n}$, these proposals are significantly
distinguishable: see Fig. 8.

Figure 8: The space-time Bohmian arrival statistics for the double-slit
experiment with a horizontal screen placed at $L_{y}\\!=\\!15\,\mu$m. The
interior curves in the central figure are the averaged times of arrival
obtained by different detection schemes: see Fig. 2. The Left and top plots
are marginal arrival time distributions and marginal arrival position
distributions, respectively. The scatter plot is generated using $2\times
10^{6}$ Bohmian trajectories, and the black, blue, and green points of the
scatter plot represent the first, second, and third arrivals of Bohmian
particles to the screen, respectively. The inset is a zoom-in of the dashed
rectangle.

###  III.3 Detection schemes

As we mentioned in section II.3, according to the Bohmian deterministic point
of view, there are several possible schemes to detect arrived particles,
especially for the horizontal screen surface which we have recursive motions
on it (see Fig. 1 and 2). One can assume that the horizontal screen is swept
with a point-like detector that surveys all arrived particles at the surface
$\mathbb{S}$, which we call spot-detection scheme. In this scheme, one option
is to use a unilateral detector to detect arrived particles at the top or
bottom of $\mathbb{S}$. In this case, the positive and negative parts of the
quantum probability current have respectively corresponded to particles that
arrive at the top or bottom of $\mathbb{S}$ (as shown in Fig. 2 (a)), and we
must use Eq. (15) to calculate the screen observables. Additionally, we can
choose a bilateral detector (or two unilateral detectors) that prob all
particles that arrive from both sides of $\mathbb{S}$, along the time with
several repeats of the experiment (as shown in Fig. 2 (b)). In these
circumstances (i.e. spot-detection scheme), there is no barrier in front of
the particles before they reach the point of detection and we can use Eq. (14)
to obtain the screen observables as in the two previous subsections.

As we have already shown in section II.3, whether the particles arrive from
the top or bottom of $\mathbb{S}$, the absolute value of the quantum
probability current yield the trajectories’ density and consequently give the
joint distribution of the total arrival at each point of $\mathbb{S}$. This
fact is the case for the standard method, as well, however, there is a subtle
difference between the two proposals in the spot-detection scheme. When we
talk about the spot-detection in the Bohmian approach, it would be considered
the possibility of multi-crossing and the distribution includes all-arrivals
at $\mathbb{S}$. Although, in the standard method there is an interpretation
for $\psi^{+}_{\mathbb{S}}(\mathbf{x},t)$ and
$\psi^{-}_{\mathbb{S}}(\mathbf{x},t)$ in Eq. (10), which relates them to the
particles arrive at $\mathbb{S}$ in a direction which is the same or opposite
with the direction of outward normal of the screen $\mathbf{n}$, respectively
[10, 64], nevertheless, since there are no defined paths in this approach, it
is obscure whether it counts only the first-arrivals to each side of the
screen or includes recursive movements of particles.

Alternatively, along with the spot-detection scheme, it could be assumed that
there is a continuous flat barrier in front of the particle’s paths as the
detection surface or screen surface that does not allow particles to cross
this surface. Depending on the screen’s length and position, there are several
possibilities for the detection process. In each case, a specific number of
particle paths contribute to the distribution of arrival time. In the simplest
case, the screen blocks all the trajectories that reach the horizontal surface
$\mathbb{S}$, and we only detect the first-arrivals. In such a setup, we can
no longer use the quantum flux method to represent Bohmian trajectories’ first
encounter with the surface; hence, the screen observables must be obtained by
numerical analysis, due to the definition of truncated current as in Eq. (16)
and its corresponding joint distribution,
$\tilde{\mathbb{P}}_{\text{QF}}(\mathbf{x},t|\mathbf{x}\\!\in\\!\mathbb{S})$,
defined in Eq. (17). By computing the Bohmian trajectories, we can find
positions and times of the first-arrivals to the screen, and consequently
calculate the arrival time distribution which mathematically could be defined
as

$\tilde{\Pi}_{\text{QF}}(t|\mathbf{x}\\!\in\\!\mathbb{S})=\int_{\mathbb{S}}\tilde{\mathbb{P}}_{\text{QF}}(\mathbf{x},t|\mathbf{x}\\!\in\\!\mathbb{S})dS.$
(24)

Also, other observable quantities such as the cumulative spatial distribution
and averaged arrival time over the detection surface could be defined and
calculated numerically in a similar way—by substituting
$\tilde{\mathbb{P}}_{\text{QF}}(\mathbf{x},t|\mathbf{x}\\!\in\\!\mathbb{S})$
in Eqs. (21) and (22). Furthermore, we can complete the computations to find
the second and third encounters to the surface (regardless of the barrier).

In Fig. 8, we show our numerical results of Bohmian trajectories simulation.
The background scatter plot is the position and time of arrivals of $2\times
10^{6}$ trajectories. In this plot, the second and third arrivals are shown in
blue and green, respectively. Here, it is more clear why we interpret the
grooves of the quantum flux density plot (Fig. 6 (c)) as a result of the
multi-crossing of Bohmian trajectories. The three middle graphs are the
average time of the first and all-arrivals, which are simulation results of
$10^{8}$ trajectories, and are compared by the quantum flux method. As
expected, the average time of all-arrivals fits on the quantum flux curve.
However, the average time of first-arrivals deviates from all-arrivals in the
area discussed in the previous section (between $x=16.2$ mm and $x=19.2$ mm).

To scrutinize the deviation zone of Fig. 8 (the gray region), Fig. 9 is drawn
to show the arrival time distributions of screen positions
$x\\!=\\!16.2,\,17.4,\,18.4,\,19.2$ mm. As one can see, at the first recursive
points of quantum flux distribution, the first-arrival distributions raise
down to zero. This implies that in the presence of a barrier-like screen,
there would be a big temporal gap between arrived particles. These gaps could
be investigated as a result of the non-intersection property of Bohmian
trajectories that cause a unilateral motion of particles along the direction
of the probability current field.

Figure 9: Arrival time distribution at the horizontal screen positions
$x\\!=\\!16.2,\,17.4,\,18.4,\,19.2$ mm, and $L_{y}\\!=\\!15\,\mu$m, which are
in the gray region of Fig (8). The width of sampling in each point is about
$\delta x=0.25$ mm, and $10^{8}$ Bohmian trajectories are simulated to obtain
these distributions.

##  IV Screen back-effect

In principle, the presence of the detector could modify the wave function
evolution, before the particle detection, which is called detector back-
effect. To have a more thorough investigation of detection statistics, we
should consider this effect. Howsoever, due to the measurement problem and the
quantum Zeno effect [9], a complete investigation of the detector effects is
problematic at the fundamental level, and it is less obvious how to model an
ideal detector. Nonetheless, some phenomenological non-equivalent models are
proposed, such as the generalized Feynman path integral approach in the
presence of absorbing boundary [37, 38, 12, 39], Schrödinger equation with a
complex potential [44], Schrödinger equation with absorbing (or complex Robin)
boundary condition [40, 41, 42, 43, 44], and so on. The results of these
approaches are not the same, and a detailed study of the differences is an
interesting topic. In this section, we provide a brief review of the absorbing
boundary rule (ABR) and path-Integral with absorbing boundary (PAB) models,
then we compare them in the double-slit setup with the horizontal screen.

###  IV.1 Absorbing Boundary Rule

Among the above-mentioned phenomenological models, the absorbing boundary
condition approach has the most compatibility with Bohmian mechanics [42]. The
application of absorbing boundary condition in arrival time problem was first
proposed by Werner [40], and recently it is re-derived and generalized by
Tumulka and others using various methods [41, 42, 43, 44]. Especially, it is
recently shown that in a suitable (non-obvious) limit, the imaginary potential
approach yields the distribution of detection time and position in agreement
with the absorbing boundary rule [44]. According to this rule, the particle
wave function $\psi$ evolves according to the free Schrödinger equation, while
the presence of a detection screen is modeled by imposing the following
boundary conditions on the Detection screen, $\mathbf{x}\in\mathbb{S}$,

$\mathbf{n}\cdot\nabla\psi=i\kappa\psi,$ (25)

where $\kappa\\!>\\!0$ is a constant characterizing the type of detector, in
which $\hbar\kappa/m$ represents the momentum that the detector is most
sensitive to. This boundary condition ensures that waves with wave number
$\kappa$ are completely absorbed while waves with other wave numbers are
partly absorbed and partly reflected [41, 99]. In the absorbing boundary rule,
the joint spatiotemporal distribution of the detection event is given by
quantum flux. Considering (25), this distribution reads

$\mathbb{P}_{\text{ABR}}(t,\mathbf{x}|\mathbf{x}\\!\in\\!\mathbb{S})=\frac{|\psi_{\text{ABC}}|^{2}}{\int
dt\int_{\mathbb{S}}dS|\psi_{\text{ABC}}|^{2}},$ (26)

where $\psi_{\text{ABC}}$ represent the solution of the free Schrödinger
equation satisfying the aforementioned absorbing boundary condition. This
distribution can be understood in terms of Bohmian trajectories. The Bohmian
particle equation of motion,
$\dot{\mathbf{X}}=(\hbar/m)\text{Im}\left[\nabla\psi_{\text{ABC}}/\psi_{\text{ABC}}\right]$,
together with the boundary condition (25), imply that trajectories can cross
the boundary $\mathbb{S}$ only outwards and so there are no multi-crossing
trajectories. If it is assumed that the detector clicks when and where the
Bohmian particle reaches $\mathbb{S}$, the probability distribution of
detection events is given by (26), because the initial distribution of the
Bohmian particle is $|\psi_{\text{ABC}}(x,0)|^{2}$ [41].

###  IV.2 Path-Integral with Absorbing Boundary

In several papers [37, 38, 12, 39], Marchuwka and Schuss develop an
interesting method to calculate the detection effect of absorbing surface
using the Feynman path integral method. They postulate a separation principle
for the wave function in which we could consider the (bounded wave function)
as a sum of two parts,
$\psi(\mathbf{x},t)=\psi_{1}(\mathbf{x},t)+\psi_{2}(\mathbf{x},t)$, such that
$\psi_{1}(\mathbf{x},t)$ corresponds to the survival part of the wave which is
orthogonal to $\psi_{2}(\mathbf{x},t)$ at a time $t$ and evolve independently
[38]. So, we can obtain the probability of survival of the particle, denoted
$S(t)$, which is the probability of the particle not being absorbed by the
time $t$, as $\int_{\mathbb{D}}d^{3}\mathbf{x}|\psi_{1}(\mathbf{x},t)|^{2}$,
where the integral is over the domain $\mathbb{D}$, outside the absorbing
region. By discretizing the path integral in a time interval $[0,t]$ and
eliminating the trajectories that, in each time interval
$[t^{\prime},t^{\prime}+\Delta t^{\prime}]$ for all $t^{\prime}<t$, are
reached to the absorbing surface $\mathbb{S}$, the survival and consequently
absorbing probability would be obtained. Based on this analysis, we could
define a unidirectional probability current into the surface as
$\frac{d}{dt}[1-S(t)]$, which yields a normal component of the
multidimensional probability current density at any point on $\mathbb{S}$ as

$\displaystyle\mathbf{J}(\mathbf{x},t)\\!\cdot\\!\mathbf{n}$
$\displaystyle\\!\\!=\\!\\!$
$\displaystyle\frac{\lambda\hbar}{m\pi}\,|\mathbf{n}\\!\cdot\\!\gradient\psi(\mathbf{x},t)|^{2}$
$\displaystyle\\!\times\,\,\text{exp}\left\\{-\frac{\lambda\hbar}{m\pi}\int_{0}^{t}\\!\\!dt^{\prime}\\!\oint_{\mathbb{S}}dS|\mathbf{n}\\!\cdot\\!\gradient\psi(\mathbf{x}^{\prime},t^{\prime})|^{2}\right\\},$

where $dS\\!=\\!\mathbf{n}\cdot d\mathbf{S}$ is the magnitude of the surface
element $d\mathbf{S}$, $\mathbf{n}$ is the unit outer normal to the absorbing
surface $\mathbb{S}$, and $\lambda$ is a proportionality factor with the
dimension of length [37, 62]. Also, $\psi(\mathbf{x},t)$ is the solution of
Schrödinger equation bounded and normalized in the domain $\mathbb{D}$.
Moreover, the normal component $\mathbf{J}(\mathbf{x},t)\\!\cdot\\!\mathbf{n}$
is supposed to be the probability density for observing the particle at the
point $\mathbf{x}$ on the screen at time $t$ [12, 39].

###  IV.3 Screen back-effect in two-slit experiment

In order to complete the investigations carried out in section III, we are
going to study the screen back-effect in the double-slit experiment with a
horizontal screen. In this regard, we compare the arrival distributions which
are resulted from the absorbing boundary rule (ABR), path-Integral with
absorbing boundary (PAB), and Bohmian truncated current (BTC).

Figure 10: Density plots of joint probability distributions of position and
time (screen observable) for the horizontal screen placed at
$y\\!=\\!15\,\mu$m in the double-slit experiment. These densities are
calculated by the three methods which take the screen effects into account.

We continue with the same initial conditions as in section III, and choose
$\kappa\\!=\\!1\,\mu\text{m}^{-1}$ for ABR. This value of $\kappa$ leads to
the maximum absorption probability—which is almost $0.4$—for the chosen
initial wave function. In addition, for a more meaningful comparison, we
consider $\lambda\\!=\\!1\,\mu$m in the PAB method, which leads to the same
absorption probability as ABR. The resulting joint arrival time-position
distributions of the three methods are depicted in Fig. 10. As one can see,
the distributions of the ABR and PAB methods—i.e., panels (a) and (b) in Fig.
10—have more compatibility with each other than the result of the BTC method.
However, there are differences between them which are more obvious in the
zoomed areas. The joint density of the ABR is more

Figure 11: Averaged time of arrival at each point of the screen (central
figure), cumulative interference pattern (upper figure), and distribution of
time of arrival to the horizontal screen of the double-slit experiment placed
at $y\\!=\\!15\,\mu$m (right-hand figure). Figure 12: Arrival time
distribution at the horizontal screen positions
$x\\!=\\!16.2,\,17.4,\,18.4,\,19.2$ mm, and $L_{y}\\!=\\!15\,\mu$m, which are
calculated for the three methods which take the screen effects into account.

uniformly distributed than of the PAB method. The empty areas between the
fringes of the panel (c) of Fig. 10 are due to the elimination of the
recursive trajectories—or in other words, are due to the elimination of second
and third arrivals in Fig. 8.

For a more detailed comparison, in Fig. 11 the spatial and temporal marginal
distributions are shown. In addition, the associated local average arrival
times are compared in the central panel of this figure. The PAB method leads
to significant discrepancies in marginal distributions; The maximum difference
is about $40\%$ that occurs around $x\\!\approx\\!0.8$ mm, which seems
testable clearly. In contrast to the previous results on intrinsic
distributions, in which the difference between average arrival times was
significant, there is a good agreement in this observable for the ABR and PAB
methods. However, there is a significant difference between the average
arrival time in these two methods and BTC around $x\\!=\\!6$ mm. In Fig. 12,
the local arrival time distributions at some points on the screen are plotted,
which show similar behavior.

##  V Summary and Discussion

When and where does the wave function collapse? How one can model a detector
in quantum theory? These are the questions that we investigated in this work.
We tried to show that there is no agreed answer for these questions, even for
the double-slit experiment that has in it the heart of quantum mechanics
[100]. This is a practical encounter with the measurement problem [73]. In
this regard, we numerically investigated and compared the main proposed
answers to these questions for a double-slit setup with a horizontal detection
screen. It is shown that these proposals lead to experimentally
distinguishable predictions, thanks to the current single-atom detection
technology.

In this work, we suggest the meta-stable helium atom as a proper coherent
source of the matter wave, however, other sources may lead to some practical
improvements. For example, using heavier condensate atoms can lead to more
clear discrepancies. Moreover, it is worth noting that although the experiment
with photons may have some practical advantages, there are more complications
in its theoretical analysis. This is partially because of the relativistic
localization-causality problem [101, 102, 103, 104]. The theoretical
investigation of a proposed experiment for photons would be an interesting
extension of the present work, which has been left for future studies.

###### Acknowledgements.

We sincerely thank Mohammad Hossein Barati for carefully reviewing the
manuscript, and Sheldon Goldstein for his helpful comments.

## References

  * Note [1] As shown by Pauli [105], if the Hamiltonian spectrum is discrete or has a lower bound, then there is no self-adjoint time operator canonically conjugate to the Hamiltonian.
  * Zimmermann  _et al._ [2016] T. c. v. Zimmermann, S. Mishra, B. R. Doran, D. F. Gordon, and A. S. Landsman, Phys. Rev. Lett. 116, 233603 (2016).
  * Kataoka  _et al._ [2016] M. Kataoka, N. Johnson, C. Emary, P. See, J. P. Griffiths, G. A. C. Jones, I. Farrer, D. A. Ritchie, M. Pepper, and T. J. B. M. Janssen, Phys. Rev. Lett. 116, 126803 (2016).
  * Kolenderski  _et al._ [2014] P. Kolenderski, C. Scarcella, K. D. Johnsen, D. R. Hamel, C. Holloway, L. K. Shalm, S. Tisa, A. Tosi, K. J. Resch, and T. Jennewein, Scientific reports 4, 1 (2014).
  * Frabboni  _et al._ [2012] S. Frabboni, A. Gabrielli, G. Carlo Gazzadi, F. Giorgi, G. Matteucci, G. Pozzi, N. S. Cesari, M. Villa, and A. Zoccoli, Ultramicroscopy 116, 73 (2012).
  * Kurtsiefer  _et al._ [1997] C. Kurtsiefer, T. Pfau, and J. Mlynek, Nature 386, 150 (1997).
  * Nitta and Kudo [2008] H. Nitta and T. Kudo, Phys. Rev. A 77, 014102 (2008).
  * Das  _et al._ [2019] S. Das, M. Nöth, and D. Dürr, Phys. Rev. A 99, 052124 (2019).
  * Allcock [1969] G. Allcock, Annals of Physics 53, 253 (1969).
  * Kijowski [1974] J. Kijowski, Rep. Math. Phys. 6, 361 (1974).
  * Werner [1986] R. Werner, J. Math. Phys. 27, 793 (1986).
  * Marchewka and Schuss [2001] A. Marchewka and Z. Schuss, Phys. Rev. A 63, 032108 (2001).
  * Vona  _et al._ [2013] N. Vona, G. Hinrichs, and D. Dürr, Phys. Rev. Lett. 111, 220404 (2013).
  * Maccone and Sacha [2020] L. Maccone and K. Sacha, Phys. Rev. Lett. 124, 110402 (2020).
  * Dias and Parisio [2017a] E. O. Dias and F. Parisio, Phys. Rev. A 95, 032133 (2017a).
  * Das and Nöth [2021] S. Das and M. Nöth, Proceedings of the Royal Society A 477, 20210101 (2021).
  * Das and Struyve [2021] S. Das and W. Struyve, Phys. Rev. A 104, 042214 (2021).
  * Kazemi and Hosseinzadeh [2022] M. Kazemi and V. Hosseinzadeh, arXiv preprint arXiv:2208.01325 (2022).
  * Note [2] Note that, the Heisenberg position operator describes position measurement at a specific time, not position measurements at random times [106, 107]. In the other words, $|\psi_{t}(\mathbf{x})|^{2}$ is just the conditional position probability density $\mathbb{P}(\mathbf{x}|t)$ [14, 15, 104], not the position-time joint probability density $\mathbb{P}(\mathbf{x},t)$ [78, 79, 98].
  * Shucker [1980] D. S. Shucker, Journal of Functional Analysis 38, 146 (1980).
  * Wolf and Helm [2000] S. Wolf and H. Helm, Phys. Rev. A 62, 043408 (2000).
  * Bracken and Melloy [1994] A. J. Bracken and G. F. Melloy, Journal of Physics A: Mathematical and General 27, 2197 (1994).
  * Hofmann [2017] H. F. Hofmann, Phys. Rev. A 96, 020101 (2017).
  * Korzh  _et al._ [2020] B. Korzh, Q.-Y. Zhao, J. P. Allmaras, S. Frasca, T. M. Autry, E. A. Bersin, A. D. Beyer, R. M. Briggs, B. Bumble, M. Colangelo, _et al._ , Nature Photonics 14, 250 (2020).
  * Steinhauer  _et al._ [2021] S. Steinhauer, S. Gyger, and V. Zwiller, Applied Physics Letters 118, 100501 (2021).
  * Azzouz  _et al._ [2012] H. Azzouz, S. N. Dorenbos, D. De Vries, E. B. Ureña, and V. Zwiller, AIP Advances 2, 032124 (2012).
  * Rosticher  _et al._ [2010] M. Rosticher, F. Ladan, J. Maneval, S. Dorenbos, T. Zijlstra, T. Klapwijk, V. Zwiller, A. Lupaşcu, and G. Nogues, Applied Physics Letters 97, 183106 (2010).
  * Delgado  _et al._ [2006] F. Delgado, J. G. Muga, and G. García-Calderón, Phys. Rev. A 74, 062102 (2006).
  * Note [3] The presence of $\text{sgn}(\hat{p}_{\bot})$ operator ensures the self-adjointness of this time operator, however, leads to a modified commutation relation, i.e. $[\hat{t}_{K},\hat{H}]=i\hbar\ \text{sgn}(\hat{p}_{\bot})$.
  * Hegerfeldt and Muga [2010] G. C. Hegerfeldt and J. G. Muga, Journal of Physics A: Mathematical and Theoretical 43, 505303 (2010).
  * Hegerfeldt  _et al._ [2010] G. C. Hegerfeldt, J. G. Muga, and J. Muñoz, Phys. Rev. A 82, 012113 (2010).
  * Bohm [1952] D. Bohm, Phys. Rev. 85, 166 (1952).
  * Nelson [1966] E. Nelson, Phys. Rev. 150, 1079 (1966).
  * Hall  _et al._ [2014a] M. J. W. Hall, D.-A. Deckert, and H. M. Wiseman, Phys. Rev. X 4, 041013 (2014a).
  * Leavens [1998a] C. R. Leavens, Phys. Rev. A 58, 840 (1998a).
  * Das and Dürr [2019] S. Das and D. Dürr, Sci. Rep. 9, 2242 (2019).
  * Marchewka and Schuss [1998] A. Marchewka and Z. Schuss, Physics Letters A 240, 177 (1998).
  * Marchewka and Schuss [2000] A. Marchewka and Z. Schuss, Phys. Rev. A 61, 052107 (2000).
  * Marchewka and Schuss [2002] A. Marchewka and Z. Schuss, Phys. Rev. A 65, 042112 (2002).
  * Werner [1987] R. Werner, in  _Annales de l’IHP Physique théorique_, Vol. 47 (1987) pp. 429–449.
  * Tumulka [2022a] R. Tumulka, Annals of Physics 442, 168910 (2022a).
  * Tumulka [2022b] R. Tumulka, Phys. Rev. A 106, 042220 (2022b).
  * Dubey  _et al._ [2021] V. Dubey, C. Bernardin, and A. Dhar, Phys. Rev. A 103, 032221 (2021).
  * Tumulka [2022c] R. Tumulka, Communications in Theoretical Physics (2022c).
  * Damborenea  _et al._ [2002] J. A. Damborenea, I. L. Egusquiza, G. C. Hegerfeldt, and J. G. Muga, Phys. Rev. A 66, 052104 (2002).
  * Muga  _et al._ [1995] J. Muga, S. Brouard, and D. Macias, Annals of Physics 240, 351 (1995).
  * Halliwell [2008] J. J. Halliwell, Phys. Rev. A 77, 062103 (2008).
  * Andrews  _et al._ [1997] M. Andrews, C. Townsend, H.-J. Miesner, D. Durfee, D. Kurn, and W. Ketterle, Science 275, 637 (1997).
  * Shin  _et al._ [2004] Y. Shin, M. Saba, T. A. Pasquini, W. Ketterle, D. E. Pritchard, and A. E. Leanhardt, Phys. Rev. Lett. 92, 050405 (2004).
  * Cronin  _et al._ [2009] A. D. Cronin, J. Schmiedmayer, and D. E. Pritchard, Rev. Mod. Phys. 81, 1051 (2009).
  * Keller  _et al._ [2014] M. Keller, M. Kotyrba, F. Leupold, M. Singh, M. Ebner, and A. Zeilinger, Phys. Rev. A 90, 063607 (2014).
  * Khakimov  _et al._ [2016] R. I. Khakimov, B. Henson, D. Shin, S. Hodgman, R. Dall, K. Baldwin, and A. Truscott, Nature 540, 100 (2016).
  * Kurtsiefer and Mlynek [1996] C. Kurtsiefer and J. Mlynek, Applied Physics B 64, 85 (1996).
  * Gliserin  _et al._ [2016] A. Gliserin, M. Walbran, and P. Baum, Review of Scientific Instruments 87, 033302 (2016).
  * Kurtsiefer  _et al._ [1995] C. Kurtsiefer, T. Pfau, C. R. Ekstrom, and J. Mlynek, Applied Physics B 60, 229 (1995).
  * Copley and Udovic [1993] J. R. Copley and T. J. Udovic, Journal of research of the National Institute of Standards and Technology 98, 71 (1993).
  * Kothe  _et al._ [2013] A. Kothe, J. Metje, M. Wilke, A. Moguilevski, N. Engel, R. Al-Obaidi, C. Richter, R. Golnak, I. Y. Kiyan, and E. F. Aziz, Review of Scientific Instruments 84, 023106 (2013).
  * Vona and Dürr [2015] N. Vona and D. Dürr (Springer, 2015) pp. 95–112.
  * Arce [2012] J. C. Arce, Phys. Rev. A 85, 042108 (2012).
  * Aharonov and Bohm [1961] Y. Aharonov and D. Bohm, Physical Review 122, 1649 (1961).
  * Paul [1962] H. Paul, Annalen der Physik 464, 252 (1962).
  * Muga and Leavens [2000] J. Muga and C. Leavens, Physics Reports 338, 353 (2000).
  * Giannitrapani [1997] R. Giannitrapani, International Journal of Theoretical Physics 36, 1575 (1997).
  * Grot  _et al._ [1996] N. Grot, C. Rovelli, and R. S. Tate, Phys. Rev. A 54, 4676 (1996).
  * Egusquiza and Muga [1999] I. L. Egusquiza and J. G. Muga, Phys. Rev. A 61, 012104 (1999).
  * Muga  _et al._ [2007] G. Muga, R. S. Mayato, and I. Egusquiza,  _Time in quantum mechanics_, Vol. 734 (Springer Science & Business Media, 2007).
  * Egusquiza  _et al._ [2003] I. Egusquiza, J. Muga, B. Navarro, and A. Ruschhaupt, Physics Letters A 313, 498 (2003).
  * Leavens [2005] C. Leavens, Physics Letters A 338, 19 (2005).
  * Leavens [2007] C. Leavens, Physics Letters A 362, 256 (2007).
  * Galapon [2004] E. A. Galapon, Journal of mathematical physics 45, 3180 (2004).
  * Galapon  _et al._ [2004] E. A. Galapon, R. F. Caballar, and R. T. B. Jr, Phys. Rev. Lett. 93, 180406 (2004).
  * Galapon  _et al._ [2005] E. A. Galapon, F. Delgado, J. G. Muga, and I. n. Egusquiza, Phys. Rev. A 72, 042107 (2005).
  * Galapon [2009] E. A. Galapon, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 465, 71 (2009).
  * Flores and Galapon [2019] P. C. M. Flores and E. A. Galapon, Phys. Rev. A 99, 042113 (2019).
  * Galapon [2012] E. A. Galapon, Phys. Rev. Lett. 108, 170402 (2012).
  * Galapon [2002] E. A. Galapon, Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences 458, 2671 (2002).
  * Dias and Parisio [2017b] E. O. Dias and F. Parisio, Physical Review A 95, 032133 (2017b).
  * Kijowski [1999] J. Kijowski, Phys. Rev. A 59, 897 (1999).
  * Daumer  _et al._ [1997] M. Daumer, D. Dürr, S. Goldstein, and N. Zanghì, Journal of Statistical Physics 88, 967 (1997).
  * Halliwell and Yearsley [2009] J. Halliwell and J. Yearsley, Physics Letters A 374, 154 (2009).
  * Boonchui and Hutem [2013] S. Boonchui and A. Hutem, Journal of Physics A: Mathematical and Theoretical 46, 105305 (2013).
  * Hannstein  _et al._ [2005] V. Hannstein, G. C. Hegerfeldt, and J. G. Muga, Journal of Physics B: Atomic, Molecular and Optical Physics 38, 409 (2005).
  * Leavens [1993] C. Leavens, Physics Letters A 178, 27 (1993).
  * McKinnon and Leavens [1995] W. R. McKinnon and C. R. Leavens, Phys. Rev. A 51, 2748 (1995).
  * Leavens [1998b] C. Leavens, Superlattices and Microstructures 23, 795 (1998b).
  * Ali  _et al._ [2003] M. M. Ali, A. S. Majumdar, D. Home, and S. Sengupta, Phys. Rev. A 68, 042105 (2003).
  * Dürr  _et al._ [1992] D. Dürr, S. Goldstein, and N. Zanghi, Journal of Statistical Physics 67, 843 (1992).
  * Valentini and Westman [2005] A. Valentini and H. Westman, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 461, 253 (2005).
  * Grübl and Rheinberger [2002] G. Grübl and K. Rheinberger, Journal of Physics A: Mathematical and General 35, 2907 (2002).
  * Hall  _et al._ [2014b] M. J. W. Hall, D.-A. Deckert, and H. M. Wiseman, Phys. Rev. X 4, 041013 (2014b).
  * Viale  _et al._ [2003] A. Viale, M. Vicari, and N. Zanghì, Phys. Rev. A 68, 063610 (2003).
  * Paul and Qureshi [2017] T. Paul and T. Qureshi, Phys. Rev. A 95, 042110 (2017).
  * Mishra  _et al._ [2019] S. Mishra, A. Venugopalan, and T. Qureshi, Phys. Rev. A 100, 042122 (2019).
  * Fang  _et al._ [2010] A.-p. Fang, Y.-l. Chen, F.-l. Li, H.-r. Li, and P. Zhang, Phys. Rev. A 81, 012323 (2010).
  * Laurat  _et al._ [2005] J. Laurat, G. Keller, J. A. Oliveira-Huguenin, C. Fabre, T. Coudreau, A. Serafini, G. Adesso, and F. Illuminati, Journal of Optics B: Quantum and Semiclassical Optics 7, S577 (2005).
  * Barnea  _et al._ [2018] A. R. Barnea, O. Cheshnovsky, and U. Even, Phys. Rev. A 97, 023601 (2018).
  * Vassen  _et al._ [2012] W. Vassen, C. Cohen-Tannoudji, M. Leduc, D. Boiron, C. I. Westbrook, A. Truscott, K. Baldwin, G. Birkl, P. Cancio, and M. Trippenbach, Rev. Mod. Phys. 84, 175 (2012).
  * Das  _et al._ [2022] S. Das, D.-A. Deckert, L. Kellers, and W. Struyve, arXiv preprint arXiv:2211.13362 (2022).
  * Fevens and Jiang [1999] T. Fevens and H. Jiang, SIAM Journal on Scientific Computing 21, 255 (1999).
  * Aharonov  _et al._ [2017] Y. Aharonov, E. Cohen, F. Colombo, T. Landsberger, I. Sabadini, D. C. Struppa, and J. Tollaksen, Proceedings of the National Academy of Sciences 114, 6480 (2017).
  * Hegerfeldt [1985] G. C. Hegerfeldt, Phys. Rev. Lett. 54, 2395 (1985).
  * Sebens [2019] C. T. Sebens, Foundations of Physics 49, 365 (2019).
  * Kazemi  _et al._ [2018] M. J. Kazemi, H. Hashamipour, and M. H. Barati, Phys. Rev. A 98, 012125 (2018).
  * Terno [2014] D. R. Terno, Phys. Rev. A 89, 042111 (2014).
  * Pauli [1958] W. Pauli, in  _in Encyclopedia of Physics_ , Vol. 5/1, edited by S. Flugge (Springer,Berlin, 1958) p. 60.
  * Dürr  _et al._ [2004] D. Dürr, S. Goldstein, and N. Zanghì, Journal of Statistical Physics 116, 959 (2004).
  * Dürr and Teufel [2004] D. Dürr and S. Teufel, in  _Multiscale Methods in Quantum Mechanics_ (Springer, 2004) pp. 41–58.


BS

    base station
CMS

    canonical minimum scattering
DMA

    dynamic metasurface antenna
LIS

    large intelligent surface
LoS

    line-of-sight
mMIMO

    massive multiple-input multiple-output
RF

    radio-frequency

# Tensor-based Channel Tracking for RIS-Empowered Multi-User MIMO Wireless
Systems

Jide Yuan1, George C. Alexandropoulos2, Eleftherios Kofidis3,
Tobias Lindstrøm Jensen1, and Elisabeth De Carvalho1

1Department of Electronic Systems, Aalborg University, Aalborg, Denmark
2Department of Informatics and Telecommunications, National and Kapodistrian
University of Athens, Greece  3Department of Statistics and Insurance Science,
University of Piraeus, Greece

###### Abstract

The accurate estimation of Channel State Information (CSI) is of crucial
importance for the successful operation of Multiple-Input Multiple-Output
(MIMO) communication systems, especially in a Multi-User (MU) time-varying
environment and when employing the emerging technology of Reconfigurable
Intelligent Surfaces (RISs). Their predominantly passive nature renders the
estimation of the channels involved in the user-RIS-base station link a quite
challenging problem. Moreover, the time-varying nature of most of the
realistic wireless channels drives up the cost of real-time channel tracking
significantly, especially when RISs of massive size are deployed. In this
paper, we develop a channel tracking scheme for the uplink of RIS-enabled MU
MIMO systems in the presence of channel fading. The starting point is a tensor
representation of the received signal and we rely on its PARAllel FACtor
(PARAFAC) analysis to both get the initial estimate and track the channel time
variation. Simulation results for various system settings are reported, which
validate the feasibility and effectiveness of the proposed channel tracking
approach.

###### Index Terms:

Channel estimation, tracking, multi-user MIMO, PARAFAC, reconfigurable
intelligent surface.

##  I Introduction

The technology of Reconfigurable Intelligent Surfaces (RISs) is currently
considered as one of the most promising technologies for future wireless
communications, due to its inherent capability to offer programmable signal
propagation conditions [1]. RISs are artificial surfaces containing massive
numbers of almost zero power consumption elements that have tunable impedance
properties, whose joint configuration can dynamically adjust the phase of the
impinging signals and thus shape the channel incurred [2]. The tunable signal
reflections can realize more directive beams towards the desired receiver (and
cancel the interference to undesired ones), resulting in improved Signal-to-
Interference-plus-Noise Ratio (SINR) [1], or enable spatio-temporal signal
focusing in rich scattering environments [2].

For the RIS-based configuration to deploy its full potential, Channel State
Information (CSI) for the end-to-end wireless link is required. Its estimation
is, in general, a challenging task due to the absence of baseband processing
capability of the (baseline) RISs [3]; most available RIS designs only contain
basic circuitry for low-resolution phase control, that is, they lack the
ability to also transmit training signals [2]. Moreover, the size of the
Channel Estimation (CE) problem increases with the number of the RIS elements,
and of course the number of antennas and users in a Multi-User (MU) Multiple-
Input Multiple-Output (MIMO) system, rendering it a very challenging task in
practically large RIS-empowered networks. This challenge has already motivated
intensive research and the related literature is quite extensive [4]. An
additional challenge comes from the fact that the wireless propagation
environment is, in general, time-varying in practice, resulting in channel
matrices that are correlated in the temporal dimension. The latter aspect has
only recently started to be considered in the RIS literature, also for channel
estimation _and tracking_ purposes. Related works address the problem via
Doppler effect estimation and mitigation [5, 6, 7, 8] and/or state-space
modeling and associated Kalman filtering [9, 10].

Exploiting the tensor representation and a respective decomposition for the
received signal constitutes a recently developed research line in the area of
CE for RIS-empowered systems, inspired from the high representation and
recovery potential of the tensor decomposition models [11], which are highly
relevant here in view of the intrinsically high-dimensional nature of the
signals involved in such systems. Indeed, the received signal — in an uplink
setting — can be viewed as having three dimensions, namely the spatial one,
corresponding to the Base Station (BS) antennas, another dimension induced by
the RIS phase configuration patterns, and of course the temporal one. Methods
based on the PARAllel FACtor (PARAFAC) decomposition have thus recently been
proposed for estimating all involved channel matrices in RIS-empowered MU MIMO
networks; see, e.g., [12, 13, 14, 15, 16] and references therein. However, the
assumption of stationarity in time was made for the channels of the system in
all these works.

In this paper, we consider the channel tracking problem in the uplink of a MU
MIMO system empowered by a passive RIS, assuming time-varying wireless
channels and relying on an adaptive PARAFAC decomposition scheme for the
signal at the BS to both get the initial estimates of the User Equipment
(UE)-RIS channels and track their variations in time. Moreover, the sparse
nature of the incurred channels is exploited with the aid of a Generalized
Approximate Message Passing (GAMP) algorithm [17]. Simulation results for
various system settings are reported, which validate the feasibility and
effectiveness of the proposed channel tracking approach. To the best of our
knowledge, this is the first tensor-based approach to the problem. We should
stress at this point that this work is only a preliminary step in this
direction and is to be improved and extended through related on-going work.

Notation: Vectors and matrices are denoted by bold lower- and upper-case
symbols, respectively. Calligraphic letters are used to denote higher-order
tensors. $(\cdot)^{\mathrm{T}}$, $(\cdot)^{\mathrm{H}}$, and
$(\cdot)^{{\dagger}}$ stand for the transposition, Hermitian transposition,
and pseudo-inversion, respectively. $\\{\mathbf{A}\\}_{j,k}$ is the $(j,k)$-th
element of the matrix $\mathbf{A}$ and similarly for higher-order tensors. The
notation ${\text{diag}}(\mathbf{a})$ is used to denote the diagonal matrix
with the elements of the vector $\mathbf{a}$ on its main diagonal. The
identity matrix of order $n$ is denoted by $\mathbf{I}_{n}$. The symbol
$\diamond$ denotes the Khatri-Rao (columnwise Kronecker) product.
$j=\sqrt{-1}$ is the imaginary unit. The expectation operator is denoted by
$\mathbb{E}\\{\cdot\\}$. $\|\cdot\|_{\mathrm{F}}$ stands for the Frobenius
norm. $\mathbb{C}$ is the field of complex numbers.

##  II System and Channel Models

###  II-A System model

We consider an RIS-empowered MU communication system in the uplink direction,
as illustrated in Fig. 1, where $M$ single-antenna UEs communicate with one BS
equipped with $N_{\mathrm{r}}$ antennas, with the assistance of a passive RIS.

Figure 1: The considered RIS-empowered uplink multi-user MIMO system.

All nodes are considered synchronized and the RIS configuration control is
handled by the BS [1]. The RIS contains $K$ unit elements, each with a tunable
impedance capable of configuring its reflection coefficient [2]. Moreover,
channels are assumed to undergo flat fading. We neglect any direct link
between the BS and any of the UEs and focus on the problem of estimating and
tracking the RIS-UE channels. In cases where a direct link also exists,
traditional methods can be deployed while keeping the RIS deactivated [13, 4].

We divide the time in _slots_ and consider the following transmission
protocol:

  1. i)

Each UE periodically sends its own pilot sequence, of length $S$, $L$ times
per slot, with $L$ being the number of different phase profiles implementable
in the RIS.

  2. ii)

The RIS maintains a phase profile for the duration of a pilot sequence and
changes to another when a whole pilot sequence has been received at the BS.

We refer to the time needed to transmit one complete pilot sequence as the
duration of a _block_. Let
$\mathbf{X}=\left[{{\mathbf{x}}_{1}\,\cdots\,{\mathbf{x}}_{M}}\right]^{\mathrm{T}}\in{\mathbb{C}^{M\times
S}}$ be the matrix of the pilot sequences sent from all UEs at each block,
where ${\mathbf{x}}_{m}\in{\mathbb{C}^{S\times 1}}$ represents the pilot
sequence of the $m$-th UE. The received pilot signal
${\mathbf{Y}}[{i,l}]\in\mathbb{C}^{N_{\text{r}}\times S}$ at the $l$-th block
of the $i$-th slot can be expressed as

$\displaystyle{\mathbf{Y}}[{i,l}]=\mathbf{G}[i]{\text{diag}}(\bm{\phi}[l]){\mathbf{H}}[i]{\mathbf{X}}+{\mathbf{W}}[{i,l}],$
(1)

where
${\mathbf{H}}[i]=\left[{{{\mathbf{h}}_{1}}[i],\ldots,{{\mathbf{h}}_{M}}[i]}\right]\in{\mathbb{C}^{K\times
M}}$ contains the channels of the RIS-UEs links and
${\mathbf{G}}[i]\in{\mathbb{C}^{{N_{\text{r}}}\times K}}$ is the BS-RIS
channel at the $i$-th slot.
$\bm{\phi}[l]\triangleq{[{{e^{j{\theta_{1}}[l]}}\,\cdots\,{e^{j{\theta_{K}}[l]}}}]^{\mathrm{T}}}\in{\mathbb{C}^{K\times
1}}$ is the chosen RIS phase profile at the $l$th block of the $i$th slot,
with ${\theta_{k}}[l]\in[0,2\pi)$, $k=1,2,\ldots,K$.
${\mathbf{W}}[i,l]\in{\mathbb{C}^{{N_{\text{r}}}\times S}}$ represents the
zero-mean spatially and temporally uncorrelated noise at the BS receiver.

###  II-B Channel model

Considering that both the BS and the RS can be described as uniform linear
arrays (normally a planar array for the RIS), we can write the following for
the RIS-BS and UE-RIS channels [18]

$\displaystyle{\mathbf{G}}[i]$
$\displaystyle=\sum\nolimits_{p=1}^{P}{{\alpha_{p,i}}{{\mathbf{a}}_{{N_{\text{r}}}}}\left({{\psi_{p,i}}}\right){\mathbf{a}}_{K}^{\mathrm{H}}\left({{\omega_{p,i}}}\right)},$
(2) $\displaystyle{{\mathbf{h}}_{m}}[i]$
$\displaystyle=\sum\nolimits_{j=1}^{{J_{m}}}{{\beta_{m,j,i}}{{\mathbf{a}}_{K}}\left({{\varphi_{m,j,i}}}\right)},$
(3)

where $P$ and $J_{m}$ denote the number of the propagation paths of the RIS-BS
link and the link from the $m$-th UE to the RIS, respectively, and
${\alpha_{p,i}}$ and ${\beta_{m,j,i}}$ are the corresponding complex path
gains. The generic vector ${{\mathbf{a}}_{U}}\left(u\right)$ denotes the array
steering vector ${{\mathbf{a}}_{U}}\left(u\right)=\left[{1,{e^{-j2\pi
u}},\ldots,{e^{-j2\pi\left({U-1}\right)u}}}\right]^{\mathrm{T}}$, where
$U\in\\{N_{\text{r}},K\\}$ and
$u\in\\{{{\psi_{p,i}}},{{\omega_{p,i}}},{{\varphi_{m,j,i}}}\\}$, with
${{\psi_{p,i}}},{{\omega_{p,i}}}$, and ${{\varphi_{m,j,i}}}$ being the
corresponding directional cosines. We make the assumption that
${\mathbf{G}}[i]$ changes much more slowly than ${\mathbf{H}}[i]$ due to the
fact that the RIS and the BS locations are fixed, rendering their wireless
channel relatively stationary. We therefore assume that ${\beta_{m,j,i}}$ and
${{\varphi_{m,j,i}}}$ vary independently across slots, while ${\alpha_{p,i}}$,
${{\psi_{p,i}}}$ and ${{\omega_{p,i}}}$ may vary only every $I$ time slots.

##  III Tensor-based Channel Tracking

###  III-A PARAFAC signal model

Assuming the pilot sequences are long enough ($S\geq M$) to be orthogonal
among the UEs, $\mathbf{X}$ can be assumed to have orthonormal rows, which
then yields the following equivalent of (1):

$\displaystyle{\mathbf{\tilde{Y}}}[{i,l}]=\mathbf{G}[i]{\text{diag}}(\bm{\phi}[l]){\mathbf{H}}[i]+{\mathbf{\tilde{W}}}[{i,l}],$
(4)

with
$\mathbf{\tilde{Y}}[{i,l}]\triangleq\mathbf{Y}[{i,l}]\mathbf{X}^{\mathrm{H}}\in\mathbb{C}^{N_{\mathrm{r}}\times
M}$ and similarly for $\mathbf{\tilde{W}}[i,l]$. It is readily verified [11]
that this is the approximate rank-$K$ PARAFAC decomposition of the
$N_{\mathrm{r}}\times M\times L$ tensor with frontal slices as in (4). The
mode-1, 2, and 3 PARAFAC factors are $\mathbf{G}[i]$,
$\mathbf{H}[i]^{\mathrm{T}}$, and
$\bm{\Phi}\triangleq\left[\bm{\phi}[1],\bm{\phi}[2],\ldots,\bm{\phi}[L]\right]^{\mathrm{T}}\in\mathbb{C}^{L\times
K}$, respectively. This formulation can be made use of to estimate the
channels at the first slot (as in, e.g., [13]) and subsequently feed those
estimates to a channel tracking scheme that pursues the slot-by-slot variation
of the UE-RIS channels.

However, the above requires the pilot sequence length, $S$, to be large enough
to allow for the (right) inversion of $\mathbf{X}$. In a MU context, with a
large number of users, $M$, this amounts to an unrealistically strict
requirement for the training overhead. We therefore stick instead to (1), that
we re-write in a more compact form as

$\displaystyle{\mathbf{Y}}[i,l]=\mathbf{G}[i]{\text{diag}}(\bm{\phi}[l])\mathbf{Z}[i]+{\mathbf{W}}[i,l],\;\;l=1,2,\ldots,L,$
(5)

where

$\displaystyle\mathbf{Z}[i]\triangleq\mathbf{H}[i]\mathbf{X}$ (6)

is the $K\times S$ matrix of the pilot signals received at the RIS elements.
Defining the $N_{\mathrm{r}}\times S\times L$ tensor $\mathcal{Y}$ with
frontal slices as in (5), we then arrive at a PARAFAC model with factors
$\mathbf{G}[i],\mathbf{Z}[i]^{\mathrm{T}}$, and $\bm{\Phi}$, that we can use
to first get initial estimates of the unknown factors,
$\mathbf{G}[1],\mathbf{Z}[1]$, before we track the time-varying
$\mathbf{Z}[i]$, $i=2,3,\ldots,I-1$, with the aid of a recursive scheme to be
developed later in this section. Then, capitalizing on the sparse nature of
$\mathbf{H}[i]$, we recover it from (6) with the aid of a GAMP approach.

###  III-B Uniqueness analysis

Clearly, if the known matrix $\bm{\Phi}$ is of full column rank, which
requires that $L\geq K$, the above PARAFAC decomposition is unique [19].
Indeed, one can then estimate the Khatri-Rao product of the two unknown
factors, whereby these can be found via least-squares (LS) Khatri-Rao
factorization (KRF) [19, Proposition 3.1].

Of course, the requirement that $L\geq K$ can be unrealistic for large-scale
RISs. An alternative uniqueness condition, valid for the case that $\bm{\Phi}$
is not necessarily of full column rank, can be deduced from [19, Proposition
3.2] and, in its generic form, and for a realistic setting where $L<K$ and
$S<N_{\mathrm{r}},K$, states that $L+S-2\geq K$.

It should be stressed that these are only sufficient, not necessary
conditions.

###  III-C The proposed channel tracking scheme

####  III-C1 Initial channel estimation

Estimates for $\mathbf{G}$ and $\mathbf{Z}$ at the fist slot can be computed
in a way similar to [13], namely with the aid of the Alternating Least Squares
(ALS) algorithm. In fact, since one of the PARAFAC factors is known, this
should be called Bilinear ALS (BALS). The idea is to alternatingly solve the
following two LS problems, fixing the one of the factors when computing the
other:

$\displaystyle\hat{\mathbf{G}}[1]$
$\displaystyle\triangleq\mathop{\text{argmin}}\limits_{\mathbf{G}}\left\|\mathbf{Y}_{(1)}[1]-\mathbf{G}\left(\bm{\Phi}\diamond\mathbf{Z}^{{T}}[1]\right)^{{T}}\right\|_{{F}}^{2},$
(7) $\displaystyle{\mathbf{\hat{Z}}}[1]$
$\displaystyle\triangleq\mathop{\text{argmin}}\limits_{\mathbf{Z}}\left\|\mathbf{Y}_{(2)}[1]-\mathbf{Z}^{{T}}\left(\bm{\Phi}\diamond\mathbf{G}[1]\right)^{{T}}\right\|_{{F}}^{2},$
(8)

where
$\mathbf{Y}_{(1)}[i]\triangleq\left[\mathbf{Y}[i,1]\,\cdots\,\mathbf{Y}[i,L]\right]$
and
$\mathbf{Y}_{(2)}[i]\triangleq\left[\mathbf{Y}^{{T}}[i,1]\,\cdots\,\mathbf{Y}^{{T}}[i,L]\right]^{{T}}$
are the mode-1 and mode-2 unfoldings of $\mathcal{Y}$, respectively. This
results in the iterative procedure of alternating between the following until
convergence has been achieved:

$\displaystyle{\mathbf{\hat{G}}}[1]$
$\displaystyle=\mathbf{Y}_{(1)}[1]\left[\left(\bm{\Phi}\diamond\mathbf{\hat{Z}}^{\mathrm{T}}[1]\right)^{\mathrm{T}}\right]^{\dagger},$
(9) $\displaystyle{\mathbf{\hat{Z}}}[1]$
$\displaystyle={\left({{\bm{\Phi}}\diamond{\mathbf{\hat{G}}}[1]}\right)^{\dagger}}{\mathbf{Y}}_{(2)}^{\mathrm{T}}[1].$
(10)

Alternatively, LS KRF may be employed.

####  III-C2 Recursive channel tracking

Recall here that $\mathbf{G}[i+1]\approx\mathbf{G}[i]$, over a period of $I$
slots. However, $\mathbf{H}[i]$ may be changing and so will $\mathbf{Z}[i]$.
Once $\mathbf{G}[i]\approx\mathbf{G}[1]$ has been estimated, an estimate for
$\mathbf{Z}[i+1]$ can be simply obtained from the analog of (10) for slot
$i+1$, namely

$\displaystyle{\mathbf{\hat{Z}}}[i+1]$
$\displaystyle=\mathbf{F}^{\dagger}\mathbf{Y}_{(2)}^{\mathrm{T}}[i+1],$ (11)

where $\mathbf{F}$ is the $N_{\mathrm{r}}L\times K$ matrix
$\bm{\Phi}\diamond\mathbf{G}[1]$. However, this is a formidable task, in view
of the generally large size of this inversion problem.

It is thus preferable to follow a recursive approach to computing
$\mathbf{Z}[i+1]$ based on its previous estimates. Specifically, we consider
the minimization with respect to $\mathbf{Z}[i+1]$ of the following
exponentially weighted squared error function at slot $i+1$:

$\displaystyle
J[i+1]\triangleq\sum_{\tau=1}^{i+1}\lambda^{i+1-\tau}\left\|\mathbf{Y}_{(2)}^{\mathrm{T}}[\tau]-\mathbf{F}[i+1]\mathbf{\hat{Z}}[\tau]\right\|_{\mathrm{F}}^{2},$
(12)

where $\lambda\in(0,1]$ is the forgetting factor. This yields a Recursive LS
(RLS) method, which resembles the classical RLS algorithm and hence its
details are omitted here. It can be initialized with $\mathbf{Z}[1]$, computed
in the BALS step. This two-step procedure will henceforth be referred to as
BALS-RLS.

####  III-C3 Resolution of ambiguities

Note that the permutation ambiguity in our PARAFAC model is trivially resolved
with the aid of the knowledge of $\bm{\Phi}$. Regarding the scaling
ambiguities that affect the rows of $\hat{\mathbf{Z}}[i]$ and the columns of
$\hat{\mathbf{G}}[1]$, they can be addressed via appropriate normalization of
the channel estimates, which is an essential step for the GAMP-based channel
recovery presented next.

###  III-D GAMP-based channel recovery

In order to exploit the sparsity of $\mathbf{H}$ when estimating it from (6),
we resort to a GAMP algorithm [17]. This way, the channel can be recovered
with fewer pilots than needed in sparsity-unaware methods. This algorithm
belongs to a family of methods that are based on Gaussian and quadratic
approximations of loopy belief propagation, a unified methodology of
estimating sparse/dense matrices. The problem can be stated as follows:

$\displaystyle{\mathbf{\hat{H}}}[i]=\mathop{\text{argmin}}\limits_{\mathbf{H}}\left\|{{\hat{\mathbf{Z}}}[i]-{\mathbf{H}}[i]{\mathbf{X}}}\right\|_{\mathrm{F}}^{2}.$
(13)

To fit the GAMP formalism, this optimization problem can be equivalently
transformed into [20]:

$\displaystyle\hat{\mathbf{H}}_{\text{a}}[i]=\mathop{\text{argmin}}\limits_{\mathbf{H}_{\text{a}}}\left\|{{\hat{\mathbf{Z}}}^{T}[i]\mathbf{U}-{\mathbf{X}}^{T}}{\mathbf{H}_{\text{a}}}[i]\right\|_{\mathrm{F}}^{2},$
(14)

where
${\mathbf{H}_{\text{a}}}\left[i\right]\triangleq{\mathbf{H}^{T}}\left[i\right]\mathbf{U}$
is a sparse matrix representing the channel in the angular domain, with
$\mathbf{U}$ being the $K\times K$ normalized Discrete Fourier Transform (DFT)
matrix. This manipulation explores the row-sparsity of the channel, making it
more efficient to apply GAMP for computing $\hat{\mathbf{H}}_{\text{a}}[i]$.
The sought for estimate $\hat{\mathbf{H}}[i]$ can then be computed as
$\hat{\mathbf{H}}[i]=\mathbf{U}{\hat{\mathbf{H}}_{a}^{\mathrm{T}}}[i]$.

GAMP performs better than other sparsity-aware methods in this task. The
reason is that ${\mathbf{H}}_{a}[i]$ is not a strictly sparse matrix due to
the angular spread and the resolution limit of the DFT matrix. GAMP, however,
can cope with both sparse and dense matrices, and this gives it an advantage
over other algorithms [21]. Its worst-case complexity per iteration is
$\mathcal{O}(MS)$, yielding an overall complexity $\mathcal{O}(TMS)$, with $T$
being the (maximum) number of iterations.

##  IV Simulation Results

We adopt the Normalized Mean Squared Error (NMSE) at each slot $i$,

$\text{NMSE}_{\mathbf{A}[i]}=\mathbb{E}\left\\{\frac{\left\|\mathbf{\hat{A}}[i]-\mathbf{A}[i]\right\|_{\mathrm{F}}^{2}}{\left\|\mathbf{A}[i]\right\|_{\mathrm{F}}^{2}}\right\\},$

with ${\mathbf{A}}[i]\in\\{{\mathbf{G}[i]\mathbf{Z}[i]},{\mathbf{H}[i]}\\}$,
as the performance metric in our comparative study, where the expectation is
computed with averaging over 100 Monte Carlo runs.

All angles are chosen uniformly at random. The ${{\varphi_{m,j,i}}}$s are left
to vary from one slot to the other, while the ${{\psi_{p,i}}}$s and
${{\omega_{p,i}}}$s are kept constant over $I$ slots. We have adopted
truncated DFT matrices for both the pilot matrix $\mathbf{X}$ and the RIS
phase profile matrix $\bm{\Phi}$. There are $M=20$ users, with $J_{m}=4$ for
all of them, while the BS is equipped with $N_{\text{r}}=16$ antennas. The
forgetting factor was set to $\lambda=0.5$.

To evaluate the tracking ability of the BALS-RLS algorithm, we first check the
NMSE for ${\mathbf{A}}[i]={\mathbf{G}[1]\mathbf{Z}[i]}$. Fig. 2 depicts the
NMSE performance as a function of the Signal-to-Noise Ratio (SNR) for three
algorithms: BALS-RLS, recursive tracking with random initialization (RLS), and
BALS for all time slots.

Figure 2: The NMSE performance of the BALS, RLS, and BALS-RLS algorithms
versus SNR, at the $3^{\text{rd}}$ and 50-th slots. $K=64$ and $I=100$.

The NMSE at the 3rd and the 50-th time slots is plotted, assuming a relatively
stable RIS-BS channel, with $I=100$. We can see that the BALS and BALS-RLS
algorithms provide nearly the same performance in both time slots. On the
contrary, the RLS algorithm performs poorly at the 3rd slot and is still not
comparable to the other two algorithms at $i=50$. The NMSE performance of the
two recursive schemes, RLS and BALS-RLS, is plotted in Fig. 3 as a function of
the slot index, for $I=200$ time slots in total, and for various RIS sizes.

Figure 3: The NMSE performance of the RLS and BALS-RLS algorithms for $I=100$
slots. SNR=10 dB.

In this experiment, the RIS-BS channel changes every 100 slots. We observe
that RLS and BALS-RLS achieve the same level of accuracy when $K=16$ and 64.
However, RLST requires observations from more time slots to converge with
larger $K$, and it does not manage to converge when $K=256$, which seriously
limits its application. In contrast, the BALS-RLS algorithm converges within a
few time slots in all cases.

Fig. 4 plots the computational cost of the three algorithms with respect to
$K$.

Figure 4: Comparison of the run times of the BALS, RLS, and BALS-RLS
algorithms versus the number of the RIS elements, $K$.

The curves represent the average execution time per time slot. We first
observe that the cost of BALS is nearly two orders of magnitude higher than
those of the other algorithms. The complexities of the two recursive schemes
are significantly lower and very similar, especially at high SNRs.

We finally illustrate the recovery performance of the GAMP-based algorithm. In
Fig. 5, we depict the NMSE performance of $\mathbf{H}[i]$ as a function of the
pilot sequence length, $S$, for different numbers of propagation paths per UE.
The orthogonal pilot scheme is included as a benchmark.

Figure 5: Comparison of GAMP and the orthogonal pilot scheme for different
pilot sequence lengths. $K=64$ and SNR=0 dB.

It can be observed that a relatively stable estimation accuracy is obtained by
GAMP, when the number of pilots is approximately equal to $\frac{M}{2}$. Due
to the angular spread effect, the channel is not a strictly sparse matrix in
the angular domain. This imposes a limit on the recovery ability of GAMP and
hence it performs somewhat worse than the orthogonal pilot scheme when $S=M$.

##  V Conclusion

In this paper, we proposed a low complexity channel tracking algorithm for
RIS-enabled multi-user SIMO systems operating in time-varying wireless
channels. Capitalizing on the PARAFAC decomposition of the received signal
model, we first proposed a BALS-RT algorithm incorporating the RLST approach
for tracking the channel variations, which is initialized by BALS CE to
accelerate convergence. We then adopted GAMP to recover the RIS-UE channels
from the BALS-RT estimations with reduced the pilot overhead. Our numerical
simulations showcased that the presented tracking algorithm can achieve the
same level of estimation accuracy as the BALS algorithm, while exhibiting
similar complexity to RLST with halving of overheads.

## References

  * [1] E. C. Strinati _et al._ , “Reconfigurable, intelligent, and sustainable wireless environments for 6G smart connectivity,” _IEEE Commun. Mag._ , vol. 59, no. 10, pp. 99–105, Oct. 2021.
  * [2] G. C. Alexandropoulos _et al._ , “Reconfigurable intelligent surfaces for rich scattering wireless communications: Recent experiments, challenges, and opportunities,” _IEEE Commun. Mag._ , vol. 59, no. 6, pp. 28–34, Jun. 2021.
  * [3] S. Lin _et al._ , “Reconfigurable intelligent surfaces with reflection pattern modulation: Beamforming design, channel estimation, and achievable rate analysis,” _IEEE Trans. Wireless Commun._ , vol. 20, no. 2, pp. 741–754, Feb. 2021.
  * [4] F. Danufane, P. Mursia, and J. Liu, “Channel estimation in RIS-aided networks,” in _Enabling 6G Mobile Networks_ , J. Rodriguez, C. Verikoukis, J. S. Vardakas, and N. Passas, Eds. Springer, 2022.
  * [5] T. Zhou _et al._ , “Aerial intelligent reflecting surface-enhanced cell-free massive MIMO for high-mobility communication: joint Doppler compensation and power optimization,” _EURASIP J. Adv. Signal Process._ , 2021, https://doi.org/10.1186/s13634-021-00781-z.
  * [6] Z. Huang, B. Zheng, and R. Zhang, “Transforming fading channel from fast to slow: IRS-assisted high-mobility communication,” arXiv:2011.03147v2 [cs.IT], Mar. 2021.
  * [7] S. Sun and H. Yan, “Channel estimation for reconfigurable intelligent surface-assisted wireless communications considering Doppler effect,” _IEEE Wireless Commun. Lett._ , vol. 10, no. 4, pp. 790–794, Apr. 2021\.
  * [8] K. Wang, C.-T. Lam, and B. K. Ng, “IRS-aided predictable high-mobility vehicular communication with Doppler effect mitigation,” in _Proc. VTC (Spring)_ , Helsinki, Finland, Apr. 2021.
  * [9] Z. Mao, M. Peng, and X. Liu, “Channel estimation for reconfigurable intelligent surface assisted wireless communication systems in mobility scenarios,” _China Commun._ , pp. 29–38, Mar. 2021.
  * [10] P. Cai _et al._ , “Downlink channel tracking for intelligent reflecting surface-aided FDD MIMO systems,” _IEEE Trans. Veh. Technol._ , vol. 70, no. 4, pp. 3341–3353, Apr. 2021.
  * [11] N. D. Sidiropoulos _et al._ , “Tensor decomposition for signal processing and machine learning,” _IEEE Trans. Signal Process._ , vol. 65, no. 13, pp. 3551–3582, Jul. 2017.
  * [12] G. T. de Araújo _et al._ , “Channel estimation for intelligent reflecting surface assisted MIMO systems: A tensor modeling approach,” _IEEE J. Sel. Topics Signal Process._ , vol. 15, no. 3, pp. 789–802, Apr. 2021.
  * [13] L. Wei _et al._ , “Channel estimation for RIS-empowered multi-user MISO wireless communications,” _IEEE Trans. Commun._ , vol. 69, no. 6, pp. 4144–4157, Jun. 2021.
  * [14] Y. Lin _et al._ , “Tensor-based algebraic channel estimation for hybrid IRS-assisted MIMO-OFDM,” _IEEE Trans. Wireless Commun._ , vol. 20, no. 6, pp. 3770–3784, Jun. 2021.
  * [15] K. Ardah _et al._ , “Double-RIS versus single-RIS aided systems: Tensor-based MIMO channel estimation and design perspectives,” arXiv:2109.09099v1 [eess.SP], Sep. 2021.
  * [16] S. Gherekhloo _et al._ , “Tensor-based channel estimation and reflection design for RIS-aided millimeter-wave MIMO communication systems,” arXiv:2107.13851v2 [eess.SP], Dec. 2021.
  * [17] S. Rangan, “Generalized approximate message passing for estimation with random linear mixing,” in _Proc. ISIT_ , Saint-Petersburg, Russia, Jul. 2011\.
  * [18] A. Sayeed and V. Raghavan, “Maximizing MIMO capacity in sparse multipath with reconfigurable antenna arrays,” _IEEE J. Sel. Topics Signal Process._ , vol. 1, no. 1, pp. 156––166, Jun. 2007.
  * [19] M. Sørensen and L. De Lathauwer, “New uniqueness conditions for the canonical polyadic decomposition of third-order tensors,” _SIAM J. Matrix Anal. Appl._ , vol. 36, no. 4, pp. 1381––1403, 2015.
  * [20] C. Hu _et al._ , “Super-resolution channel estimation for mmwave massive MIMO with hybrid precoding,” _IEEE Trans. Veh. Technol._ , vol. 67, no. 9, pp. 8954–8958, 2018.
  * [21] D. L. Donoho _et al._ , “Message-passing algorithms for compressed sensing,” _Proc. Nat. Acad. Sci._ , vol. 106, no. 45, 2009.