# Quantum Equilibrium in Stochastic de Broglie-Bohm-Bell Quantum Mechanics

Jeroen C. Vink  Shell Global Solutions International B.V., Gasweg 31, 1031 HW
Amsterdam, The Netherlands. Email<EMAIL_ADDRESS>
###### Abstract

This paper investigates dynamical relaxation to quantum equilibrium in the
stochastic de Broglie-Bohm-Bell formulation of quantum mechanics. The time-
dependent probability distributions are computed as in a Markov process with
slowly varying transition matrices. Numerical simulations, supported by exact
results for the large-time behavior of sequences of (slowly varying)
transition matrices, confirm previous findings that indicate that de Broglie-
Bohm-Bell dynamics allows an arbitrary initial probability distribution to
relax to quantum equilibrium; i.e., there is no need to make the ad-hoc
assumption that the initial distribution of particle locations has to be
identical to the initial probability distribution prescribed by the system’s
initial wave function. The results presented in this paper moreover suggest
that the intrinsically stochastic nature of Bell’s formulation, which is
arguable most naturally formulated on an underlying discrete space-time, is
sufficient to ensure dynamical relaxation to quantum equilibrium for a large
class of quantum systems without the need to introduce coarse-graining or any
other modification in the formulation.

##  1 Introduction

The formulation of quantum mechanics due to de Broglie, Bohm and Bell (dBBB)
has many advantages over the conventional Copenhagen interpretation (cf. for
example refs. [1, 2, 3]). These advantages and benefits hold for the original,
causal formulation originally proposed by de Broglie [4] and rediscovered and
perfected by Bohm [5] as well as for the stochastic version introduced by Bell
[6] and further developed in refs. [7, 8, 9]. The benefits and attractive
features notwithstanding, these Bohm-type formulations continue to be
criticized. One of the more fundamental challenges concerns a seemingly ad-hoc
assumption that is required to reconcile the computed probabilities in dBBB
formulations with the established probability interpretation and experimental
results of quantum mechanics.

In dBBB formulations, particles always have a well-defined position and move
(either causally or stochastically) along trajectories that are guided by the
system’s wave function. If the probability distribution of the initial
particle locations is the same as the probability distribution defined by the
system’s initial wave function, then the particle dynamics is such that for
all future times, the probability distribution of the particle locations is
identical to the distribution defined by the time-dependent wave function.
Recognizing this ad-hoc assumption on the initial particle distribution as a
weakness in the formulation, already Bohm and Hiley proposed [10] that it
could be avoided if an arbitrary initial probability distributions would
dynamically relax to the quantum “equilibrium” distribution. This was further
explored and substantiated in work by Valentini and others [11, 12, 13, 14]
(see also ref. [15] for an alternative approach to understand the role of
probability in dBBB formulations).

Such a dynamical relaxation to quantum equilibrium is not at all obvious. In
the causal Bohm formulation, the same dynamics that ensures that the
probability distribution of the ensemble of particles continues to reproduce
the probabilities computed with the system’s wave function, also implies that
any deviation from this distribution in the initial state will be preserved
for all later times. Hence, the exact dynamics cannot accommodate a
spontaneous relaxation to quantum equilibrium and one has to invoke a coarse-
graining procedure where particle locations are evaluated as averaged values,
for example on a spatial lattice. Owing to underlying chaotic particle
dynamics, these coarse grained locations could be, and in 2D examples studied
in refs. [12, 13, 14] in fact are, found to relax to quantum equilibrium -
with a relaxation time that Valentini argues [11, 12] scales as

$t^{eq}\propto m\hbar^{2}/\varepsilon(\Delta P)^{3}.$ (1)

Here, $m$ is the particle mass, $\varepsilon$ the coarsening scale and $\Delta
P$ the momentum spread in the system’s wave function.

Unlike the causal Bohm formulation, Bell’s formulation is stochastic, and in
the version of refs. [7, 9] space is discrete with particle locations limited
to the sites of a spatial lattice. Here, the dynamics does not imply that an
initial non-equilibrium must persist - in fact, the similarity of this
stochastic formulation to a Markov process might suggest that equilibration
could happen without the need for further modifications of the dynamics. This
is explored in the present paper using numerical simulations of a 2D quantum
system describing a free particle in a box. In these simulations there is no
need to compute large ensembles of trajectories from which the probability
distribution of the particle location must be determined. Provided the systems
are reasonably sized, the time evolution of the probability distributions can
be computed directly, applying the time-dependent transition matrices of the
system to the evolving probability distribution. The definition of these
transition matrices and their features, along with the similarities and
difference between a Markov and a Bell process will be further discussed in
section 2 below.

The results of these numerical simulations show rapid and robust relaxation to
quantum equilibrium, which by the nature of the (weak) ergodicity of the
process automatically applies to all (non-equilibrium) initial probability
distributions. These simulations can be performed for various values of the
particle mass, lattice distance and momentum spread in the wave functions,
with results for the equilibrium time that suggest a slightly different
scaling behavior than proposed in (1):

$t^{eq}\propto mL\hbar/a(\Delta P)^{2}.$ (2)

Here, $a$ is the lattice distance and $L$ the box dimension. As in Eq. (1)
this time scale is linear in the particle mass and linear in the inverse
lattice distance - which acts as a natural substitute for the coarsening scale
required in the causal dynamics. The dependence on momentum spread is less
severe, with one (inverse) factor replaced by a factor $L$.

These results confirm that dBBB formulations do not need the additional
assumption that the initial distribution of particle locations has to be
identical to the initial probability distribution prescribed by the system’s
initial wave function. They moreover suggest that the intrinsically stochastic
nature of Bell’s formulation, which arguable is most naturally formulated on
an underling discrete space-time, is sufficient to ensure dynamical relaxation
to quantum equilibrium for a large class of quantum systems - without the need
to introduce coarse-graining or any other modification in the formulation.

The remainder of the paper is organized as follows. The next section starts
with a review of the discrete space-time version of the dBBB formulation [9],
followed by a brief summary of Markov processes and their convergence, after
which the discussion is extended to explore properties of the Bell-process in
which the Markovian transition matrices are time-dependent. Section 3 begins
with the description of the 2D quantum systems, followed by an overview of the
simulation results for the lattice distance, particle mass and momentum spread
dependent equilibrium times and, finally, section 4 contains a summary with
concluding remarks.

##  2 Markov Processes with Time-dependent Transition Probabilities

###  2.1 Stochastic Bell Dynamics

Bell’s stochastic version [6] of Bohm mechanics in the discrete time version
of ref. [9] can be formulated as follows. The system’s quantum state evolves
according to Schrödinger dynamics as

$|{\psi(t+\epsilon)}\rangle=\hat{U}|{\psi(t)}\rangle,$ (3)

where $t$ is the discrete time111As was briefly discussed in ref. [9] the time
step size $\epsilon$ could be time (index) dependent, with a magnitude that is
self-consistently determined by the system’s dynamics. Here, such a time
dependence will be ignored, as it does not materially impact the results
presented below. that progresses in steps $\epsilon$ and

$\hat{U}=e^{-i\epsilon\hat{H}}$ (4)

is the evolution operator with its usual dependence on the system’s (time-
independent) Hamiltonian $\hat{H}$. Also space is discrete and finite such
that particle configurations can be labeled with indices $n=1,\dots,N$; in
such a location representation, Eq. (3) takes the form

$\psi_{n}(t+\epsilon)=\sum_{m}U_{nm}\psi_{m}(t),$ (5)

with $\psi_{n}(t)=\langle{n}|\psi(t)\rangle$ and
$U_{nm}=\langle{n}|\hat{U}|{m}\rangle$. As was shown in ref. [9], the time-
dependent configuration probabilities obey a discrete-time continuity
equation,

$P_{n}(t+\epsilon)=P_{n}(t)+\sum_{m}J_{nm}(t),$ (6)

with

$P_{n}(t)=|\psi_{n}(t)|^{2}$ (7)

and

$J_{nm}(t)={\rm
Re}(\psi^{*}_{n}(t+\epsilon)U_{nm}\psi_{m}(t))-(n\leftrightarrow m).$ (8)

Instead of using the recursion defined in Eq. (6), the probability
distribution $P_{n}(t)$, which will also be referred to as “quantum
equilibrium distribution”, at time $t=k\epsilon$ given an initial state
$\psi(0)$ can equally well be directly computed from the Schrödinger equation
(3) as

$P_{n}(t)=|\sum_{m}(U^{k})_{nm}\psi_{m}(0)|^{2}.$ (9)

These time-dependent probabilities can also be generated from an ensemble of
stochastically evolving configuration trajectories $\\{n_{i}(t)\\}_{i=1}^{M}$.
Here, $n(t)$ is the index for the particle configuration at time $t$, and $i$
labels the configuration trajectories in the ensemble of size $M$. The
trajectories can be generated using transition probabilities $T_{nm}(t)$
defined as [6, 7]

$T_{nm}(t)=\theta(J_{nm}(t))J_{nm}(t)/P_{m}(t),$ (10) $T_{mm}(t)=1-\sum_{n\neq
m}T_{nm}(t),$ (11)

where $\theta(x)$ is the Heaviside step function and $T_{nm}(t)$ is the
probability that configuration $m=n_{i}(t)$ at time step $t$ changes to
configuration $n=n_{i}(t+\epsilon)$ at time $t+\epsilon$. The configuration
probabilities can then be computed from the ensemble as

$P_{m}(t)=\lim_{M\rightarrow\infty}\sum_{i=1}^{M}\delta_{m,n_{i}(t)}/M.$ (12)

Equivalently, the time dependence of the configuration probabilities (12) can
be computed using the following time dependent Master Equation,

$P_{n}(t+\epsilon)=P_{n}(t)+\sum_{m}\left(T_{nm}(t)P_{m}(t)-T_{mn}(t)P_{n}(t)\right).$
(13)

Since the transition matrices $T(t)$ are stochastic matrices,

$T_{nm}(t)\geq 0,\forall n,m,t;\;\;\;\;\sum_{n}T_{nm}(t)=1,\forall m,t,$ (14)

it follows that $\sum_{m}T_{mn}(t)P_{n}(t)=P_{n}(t)$ and hence Eq. (13) can be
rewritten in a more explicitly Markovian form as

$P(t+\epsilon)=T(t)P(t).$ (15)

This in turn implies that the probability at time $t$ can be computed directly
from the initial probabilities as

$P(t)=\tilde{T}(t)P(0),$ (16)

where the cumulative transition matrix $\tilde{T}(t)$ is defined as the
backward product of the preceding transition matrices $T(t)$,

$\tilde{T}(t)\equiv T(t-\epsilon)T(t-2\epsilon)\dots T(0).$ (17)

As was shown in ref. [9], $P(t)$ computed from Eq. (16) is identical to the
quantum distribution defined in Eq. (9) computed using the Schrödinger
dynamics (3) if the distributions at time $0$ are the same,
$P_{m}(0)=|\psi^{*}_{m}(0)\psi_{m}(0)|^{2}$, and provided the time step size
$\epsilon$ in Eq. (4) is sufficiently small to ensure that $T_{mm}(t)$ defined
in Eq. (11) is always non-negative.

###  2.2 Markov Processes

To set the stage for exploring if $P(t)$ defined in Eq. (16) also converges to
the quantum equilibrium distribution (9) for arbitrary $P(0)$, it is
worthwhile to review some properties of a (fixed) stochastic matrix $A$ and
the associated Markov process,

$\pi^{(k+1)}=A\pi^{(k)}.$ (18)

Since the sum of the matrix elements in every column of $A$ is one, it follows
that $u^{T}A=u^{T}$ for a unit-elements vector $u=(1,1,\dots,1)^{T}$ and hence
$u$ is a left eigenvector of $A$ with eigenvalue $1$. The Gershgorin circle
theorem (see e.g. [16]) applied to the rows of the transposed matrix $A^{T}$
then implies that all other eigenvalues have norm less or equal to one. If $A$
is also primitive, i.e., for a sufficiently large power $k$ all matrix
elements of $A^{k}$ are positive, then the Perron-Frobenius theorem (see,
e.g., [17, 18]) states that the eigenvalue $\lambda^{0}=1$ is unique and
strictly larger than all other eigenvalues $\lambda^{s}$ (i.e.,
$|\lambda^{s}|<1-\delta$ with $\delta>0$ and $s>0$). If this is the case, only
the dominant eigenmode remains in the spectral decomposition of $A^{k}$ for
sufficiently large $k$:

$\lim_{k\rightarrow\infty}A^{k}=\pi^{eq}u^{T},$ (19)

where $\pi^{eq}$ is the probability-normed right eigenvector with eigenvalue
$1$ associated with $u$,

$A\pi^{eq}=\pi^{eq},\;\;\;\;u^{T}\pi^{eq}=1.$ (20)

This implies that the Markov process with such a transition matrix converges
to the unique equilibrium distribution $\pi^{eq}$,

$\pi^{eq}=\lim_{k\rightarrow\infty}A^{k}\pi^{(0)},$ (21)

for any initial probability distribution $\pi^{(0)}$.

###  2.3 Markov Process with Time-Dependent Transition Matrix

The Markov process (18) resembles the Bell process (15), but besides the
formal similarity, there is the important difference that the transition
matrices in (15) are time dependent. It is therefore not clear how much one
can learn from the spectrum of individual matrices. However, since every
$T(t)$ is a stochastic matrix, it remains the cases that for all $t$, $T(t)$
has (at least) one eigenvalue $1$ with left eigenvector $u=(1,\dots,1)^{T}$
and all other eigenvalues have an absolute value less than or equal to $1$. It
also remains the case that the properties of the spectral decomposition of the
cumulative transition matrix,

$\tilde{T}(t)=\sum_{s=0}^{N-1}\lambda^{s}(t)v^{s}(t)\overline{v}^{sT}(t),$
(22)

determine if the Bell process (15) will equilibrate to a late-time behavior
that is independent of the initial state. Here, $v^{s}(t)$ and
$\overline{v}^{s}(t)$ are the right and left eigenvectors of $\tilde{T}(t)$,
with eigenvalue $\lambda^{s}(t)$. This equilibration will be established if
the process is “weakly ergodic”222See ref. [18] for the formal definition of
weak (and strong) ergodicity., i.e., only a single eigenmode with eigenvalue
$1$ dominates in the cumulative transition matrix for times (much) larger than
an equilibrium time $t^{eq}$:

$\tilde{T}(t)\rightarrow\lambda^{0}(t)v^{0}(t)\overline{v}^{0T}(t)\equiv\pi(t)u^{T}\;\;\;\;{\rm
for}\;\;\;\;t\gg t^{eq},$ (23)

with $\overline{v}^{0}(t)=u$ the constant left eigenvector with eigenvalue
$\lambda^{0}=1$ and $v^{0}(t)=\pi(t)$ the corresponding time-dependent right
eigenvector. If this is the case, this right eigenvector $\pi(t)$ necessarily
will be equal to the quantum distribution $P(t)$, because
$\tilde{T}(t)\pi(0)=\pi(t)$ for any initial distribution $\pi(0)$ and by
construction $\tilde{T}$ produces the quantum equilibrium distribution $P(t)$
if the initial distribution $\pi(0)=P(0)$.

This shows that it is sufficient that $\tilde{T}(t)$ satisfies the large time
behavior (23) to obtain a Bell process (17) that automatically converges to
reproduce quantum dynamics irrespective of the initial probability
distribution. The late-time behavior (23) is also necessary, since the
contribution from other modes, i.e., terms
$\lambda^{s}v^{s}(t)(\overline{v}^{sT}(t)\pi(0))$ with $s>0$, can only be
independent of $\pi(0)$ if $(\overline{v}^{sT}(t)\pi(0))$ is the same for
arbitrary $\pi(0)$, which can only be the case333More precisely (dropping the
argument $(t)$): Suppose $\tilde{T}\pi$ is invariant under an infinitesimal
shift of $\pi$ of the form, $\pi\rightarrow\pi+\delta^{j}-\delta^{j^{\prime}}$
with $\delta^{j}$ a vector which is zero everywhere, except on location $j$
where it has value $\epsilon$. Then the infinitesimal change of $\tilde{T}$
must be zero for all pairs $j,j^{\prime}$, i.e.,
$\delta(\sum_{m}\tilde{T}_{nm}\pi_{m})=\sum_{m}\tilde{T}_{nm}(\delta^{j}_{m}-\delta^{j^{\prime}}_{m})=\epsilon\sum_{s}\lambda^{s}v^{s}_{n}(\overline{v}^{s}_{j}-\overline{v}^{s}_{j^{\prime}})=\epsilon(\tilde{T}_{nj}-\tilde{T}_{nj^{\prime}})=0$,
which shows that $\tilde{T}$ must have identical columns, as in Eq. (23).  if
$\overline{v}^{s}(t)\propto u$. However, there cannot be additional left
eigenvectors proportional to $u$, since they would have to be orthogonal to
the right eigenvector $v^{0}(t)$.

To establish the result (23), it is therefore sufficient to show that the
Markov process generating the cumulative transition matrix $\tilde{T}(t)$ is
weakly ergodic. In practical terms, it is then sufficient to show that for
large $t$ the next-largest eigenvalue $\lambda^{1}(t)$ (and hence every
subsequent eigenvalue) approaches zero, $|\lambda^{1}(t)|\rightarrow 0$ for
$t\rightarrow\infty$. Slightly more specifically and in analogy with the $k$
dependence of the eigenvalue spectrum of $A^{k}$ in a Markov process (19), it
is sufficient that the time dependence of the eigenvalues can be expressed as,

$|\lambda^{s}(t)|\approx e^{-c(t)w_{s}},$ (24)

with $c(t)$ a positive (real) value that increases with $t$ and $w_{s}$ a
positive weight factor that increase with increasing eigenvalue index $s$ and
$w_{0}=0$. Somewhat remarkably, the results of the numerical experiments
discussed below can be described with weights $w_{s}\approx s$ and
coefficients $c(t)$ that grow (approximately) linearly with time, $c(t)\approx
c_{0}t$. The coefficient $c_{0}$ then can be taken as a measure for the
inverse equilibrium time, $c_{0}\approx 1/t^{eq}$.

Before proceeding to the next section, which discusses numerical evaluations
of the eigenvalue spectra of $\tilde{T}(t)$, it is worthwhile to briefly
contemplate how much one could learn from the properties of the individual
matrices $T(t)$.

Since for a normal Markov process the matrices $A$ must be primitive, it may
be reasonable to assume that also the $T(t)$ matrices (at least predominantly)
must be primitive, or equivalently, have a single eigenvalue $1$ and sub-
leading eigenvalues that are well-separated from $1$. Since matrices
$T(t+\epsilon)$ and $T(t)$ are “almost the same” (because the time step size
$\epsilon$ is very small), it is then not unreasonable to assume that
subdominant eigenmodes will be suppressed in each iteration by a factor
roughly equal to the modulus of their eigenvalue $|\lambda^{s}(t)|$.
Interestingly, sequences of slowly varying non-negative primitive matrices
have been studied in ref. [19]. There, Artzrouni proves that for a sequence of
slowly varying non-negative primitive matrices444Technically, the matrices
must furthermore be bounded and have non-zero elements well-separated from
zero. $T(t)$ with $\|T(t+\epsilon)-T(t)\|<\varepsilon_{0}$, there is an
$\varepsilon_{0}>0$ such that the backward product $\tilde{T}(t)$ is weakly
ergodic. Weak egodicity in this more generic situation implies that for all
$t$ the product matrix $\tilde{T}(t)$ can be written as,

$\tilde{T}_{nm}(t)=(w_{m}+\varepsilon_{nm}(t))L_{n}(t),$ (25)

with constant, positive weight vector $w$ and time-dependent, positive vector
$L(t)$ and $\varepsilon_{nm}(t)\rightarrow 0$ for $t\rightarrow\infty$. For
matrices that are furthermore stochastic, the columns of $\tilde{T}(t)$ always
sum up to one and hence the weights must be equal, $w_{m}=w$, such that
$wL(t)$ becomes the dominant right eigenvector of $\tilde{T}(t)$,
$wL(t)=v^{0}(t)$.

This suggests that under suitable, not necesserily very stringent conditions,
the Bell process will obtain weak ergodicity (i.e., the cumulative transition
matrix $\tilde{T}(t)$ will assume the form Eq. (23)) and hence it will support
the desired relaxation to quantum equilibrium for arbitrary initial
probability distributions. Unfortunately, it is not straightforward to upfront
determine if these conditions (primitivity, sufficiently slowly varying) apply
to matrices $T(t)$ associated with a specific quantum system. Therefore
numerical simulations will still be needed to further explore and quantify the
large-time behavior of $\tilde{T}(t)$.

Figure 1: Contour maps of an initial probability distribution with low $\Delta
P$, $N_{k}=4$ (left) and high $\Delta P$, $N_{k}=7$ (right), displayed in the
2D plane. The box dimension $L$ is set to 1 and is discretized with 45 lattice
sites in each direction.

##  3 Numerical Experiments

###  3.1 Model Definition

As in refs. [13, 14] the quantum systems investigated in the present paper
describe a free particle in a 2D box. Unlike in this previous work, the space
within the box is a square lattice of size $L\times L$ with $N$ lattice sites
in each direction. The wave functions $\psi_{x}$ on this lattice have
Dirichlet boundary conditions $\psi=0$, as in a box with impenetrable walls.
In order to assess the impact of boundary conditions, also periodic boundary
conditions will be briefly considered.

The Hamiltonian for a free particle on this 2D lattice takes the form,555
Units are such that $\hbar=c=1$ and the scale is set using the linear box
dimension $L$.

$H_{x,y}=\frac{1}{2ma^{2}}\sum_{\mu=1,2}(2\delta_{x,y}-\delta_{x,y+a\hat{\mu}}-\delta_{x,y-a\hat{\mu}}),$
(26)

where $a=L/N$ is the lattice distance, $m$ the particle mass; $x$ and $y$ are
2D vectors $(x_{1},x_{2})^{T}$ and $(y_{1},y_{2})^{T}$ with $x_{1,2}$ and
$y_{1,2}\in\\{0,a,\dots,Na-a\\}$; the $\hat{\mu}$ are unit vectors,
$\hat{1}=(1,0)^{T}$ and $\hat{2}=(0,1)^{T}$ and
$\delta_{x,y}\equiv\delta_{x_{1},y_{1}}\delta_{x_{2},y_{2}}$ etc., with
Kronecker $\delta_{x_{\mu},y_{\mu}}$.

For a box with impenetrable walls, the wave function at $x_{1,2}=-a$ and $aN$
must be zero. These Dirichlet boundary conditions imply that Eq. (26) must be
modified at the boundary, $\delta_{x,y+a\hat{\mu}}\rightarrow\delta_{x,y}$ for
$y_{\mu}=aN-a$ and $\delta_{x,y-a\hat{\mu}}\rightarrow\delta_{x,y}$ for
$y_{\mu}=0$; With periodic boundary conditions, the Kronecker deltas at the
boundaries are modified by identifying $y_{\mu}+a\hat{\mu}$ with $0$ if
$y_{\mu}=aN-a$ and $y_{\mu}-a\hat{\mu}$ with $aN-1$ if $y_{\mu}=0$.

With Dirichlet boundary conditions the eigenvectors of $H$ are,

$\psi^{(k)}_{x}=\sin((x_{1}+a)k_{1}\pi/(L+a))\sin((x_{2}+a)k_{2}\pi/(L+a)),$
(27)

with eigenvalues

$E^{(k)}=[2-\cos(k_{1}\pi/(N+1))-\cos(k_{2}\pi/(N+1))]/ma,$ (28)

where $k=(k_{1},k_{2})^{T}$ and $k_{1},k_{2}\in\\{1,\dots,N\\}$. As in refs.
[13, 14] the initial wave function consists of a superposition of the
$N_{k}^{2}$ lowest energy eigenvectors multiplied with a random phase factor,
such that the time dependent wave function is

$\psi_{x}(t)=\sum_{k_{1},k_{2}=1}^{N_{k}}e^{i\phi_{k}-itE^{(k)}}\psi_{x}^{(k)},$
(29)

with $\phi_{k}$ a random phase between $0$ and $2\pi$. As an example, Figure 1
shows contour maps of the probability distribution $|\psi(0)|^{2}$ in the 2D
box for a low $\Delta P$ and high $\Delta P$ initial state ($N_{k}=4$ and
$N_{k}=7$ respectively).

In most experiments below, the wave function is a superposition of the 16
lowest-energy states ($N_{k}=4$); when exploring the impact of the momentum
spread in the wave function, different ranges of the wave number values will
also be used, $N_{k}=3,\dots,8$. Given limitations on compute capacity, it was
not possible to explore systems in which both the lattice distance in mass
units, $ma$, would be very small and simultaneously the box dimension in mass
units, $mL$, very large. Particle masses will range from $mL=5$ to $40$, with
most computations done at $mL=20$ and $mL=5$; lattices with different
resolution will be employed with $N$ ranging from $15$ to $60$.

Figure 2: Eigenvalue spectra of $\tilde{T}(t)$ evaluated on a system with
Dirichlet boundary conditions, low $L\Delta P=8.55$ and $N=30$. The dots in
the plot show $\log(|\lambda^{s}|)$ for $t/L=40,120,240$ and $360$, in order
of increasing (negative) slope. The straight lines are least-square fits to
the first 25 data points.

Figure 3: Density plots of the probability distribution $P(t)$ (left plots),
the dominant right-eigenvalue $v^{0}(t)$ (middle plots) and the difference of
$P(t)$ and $v^{0}(t)$ (right plots). The top, middle and bottom rows are at
$t/L=0,120$ and $360$ respectively.

For a given quantum system it is straightforward to compute $\tilde{T}(t)$:
Given the solution of the Schrödinger equation shown in Eq.(29), the time-
dependent transition matrices $T(t)$ are computed from Eqs. (10) and (11),
after which $\tilde{T}(t)$ follows from Eq. (17). Both $\psi(t)$ and
$\tilde{T}(t)$ are evaluated at discrete times $t=k\epsilon$, with
$\epsilon/L=0.02/N$, such that the time step size $\epsilon/L$ and lattice
distance $a/L$ both scale as $1/N$.666To test the impact of the time step
size, a few computations were repeated with $\epsilon/L=0.01/N$, which
produced essentially identical results. The number of time steps varies per
simulation, ranging from $N_{t}=1000$ (for systems with fast equilibration) to
$50,000$ (for systems with slow equilibration).

In the numerical evaluations of the evolving transition matrix it may happen
that the constraint $\sum_{n\neq m}T_{nm}\leq 1$ is violated for one or more
columns $m$. When this happens (very rarely), the diagonal element for the
offending column is put to zero and the off-diagonals are normalized to sum up
to 1. Since this could lead to inaccuracies, it is always checked that the
final probabilities computed using the Schrödinger equation (9) and using the
Markov process (16) on the same initial probability are the same. In all cases
the two final probabilities are very close:
$\sum_{n}|\psi^{*}_{n}(t)\psi_{n}(t)-(\tilde{T}(t)P(0))_{n}|<10^{-8}$ (with
$P(0)=|\psi(0)|^{2}$).

###  3.2 Eigenvalue Spectrum of the Cumulative Transition Matrix

The objective of this subsection is to explore if the large-time eigenvalue
spectrum of the cumulative transition matrix $\tilde{T}(t)$ defined in Eq.
(17) exhibits the desired dominance of the eigenvector with eigenvalue $1$
that is shown in Eq. (24).

Figure 2 shows an example of the time-dependence of the eigenvalues of
$\tilde{T}(t)$, computed on a lattice with $N=30$. The dots represent
$\log(|\lambda^{s}|)$ for $s=0,1,\dots,40$ and the straight lines are least-
square fits of the form $b_{0}+b_{1}s$, using the first $25$ eigenvalues. From
top to bottom, the spectra are for $t/L=40,120,240$ and $360$. These results
confirm first of all that the spectra only have a single eigenvalue $1$.
Furthermore, the absolute values of the eigenvalues decrease roughly
exponentially with the eigenvalue index (as in Eq. (24) with $w_{s}\approx
s$). This suppression of eigenvalues grows with time and the steepness of the
lines in Fig. 2 should be a measure for the level of equilibration: the
steeper the slope, the closer the system is to quantum equilibrium.

The gradual approach to quantum equilibrium can also be seen from the
evolution of the dominant right eigenvector. As shown in section 2.2, this
eigenvector $v^{0}(t)$ must become equal to the quantum probability
distribution $P(t)$ when $\tilde{T}(t)$ obtains the form shown in Eq. (23).
This gradual evolution can be seen in Fig. 3, which shows a density plot of
the probability distribution $P(t)$ (left plots), the dominant right-
eigenvector $v^{0}(t)$ (middle plots) and their difference (right plots). The
top row is at the initial state, $t/L=0$, the middle and bottom row are at
$t/L=120$ and $360$ respectively.

###  3.3 Quantum Equilibrium Times

Having established that the eigenvalues $\lambda^{s}$ with $s>0$ are
(approximately) exponentially suppressed,

$\lambda^{s}(t)\approx e^{-c(t)s},$ (30)

this subsection will explore the time-dependence of this suppression. This
will lead to estimates of the equilibration time scale as a function of the
lattice distance. The dependence of the time scale on particle mass and
momentum spread in the wave function will be discussed in the next subsection.

Figure 4 shows the time dependence of the slope-coefficients $c(t)$ for the
closed box system with $mL=5$ and $L\Delta P=8.55$ evaluated on lattices with
$N=15,17,20,25,30$ and $45$. As can be seen from the quality of the fit (which
always have $R^{2}$ values larger than $0.96$), the suppression of the non-
dominant eigenvalue modes progresses to a good approximation linearly in time.
Hence the slope-coefficient can be parameterized as

$c(t)\approx c_{0}+t/t^{eq},$ (31)

where the coefficient of $t$ is identified with the inverse of the equilibrium
time scale $t^{eq}$. The intercept $c_{0}$ is small, typically less than $0.1$

As was mentioned in the introduction, the coarsening scale that is needed to
demonstrate equilibration in the causal Bohm formulation appears to play the
same role as the lattice distance in the stochastic Bell formulation.
Consistency between these two alternative formulations implies that the
equilibrium time should diverge towards the continuum limit $a\rightarrow 0$.
Adopting the expression for the equilibrium time, Eq. (1) proposed by
Valentini [11, 12, 13] and equating the coarsening scale $\varepsilon$ with
the lattice distance $a$, $t^{eq}$ is expected to grow linearly with $N=L/a$.
The results in Figure 5 indeed show that $1/t^{eq}$ is roughly proportional to
$1/N$.

The data points in this figure were obtained from linear fits of $c(t)$ vs.
$t$ in three different versions of the 2D quantum system: a closed box with
Dirichlet boundary conditions and particle mass $mL=20$ (dots with solid
line), a closed box with Dirichlet boundary conditions and particle mass
$mL=5$ (diamonds with dashed line), and an open box with periodic boundary
condition with $mL=20$ (squares with dotted line). The wave function in the
open box consists of a superposition of 9 plane wave energy eigenvectors,
similarly randomized as the eigenvectors in the closed box, with roughly
similar momentum spread.

Figure 4: Time dependence of the slope-coefficients $c(t)$ in Eq. (30), for
the closed box model with $mL=5$ and $L\Delta P=8.55$, evaluated on lattices
with different lattice distances $a=L/N$. The straight lines are least-square
fits to the data. From top to bottom the lines are for $N=15,17,20,25,30$ and
$45$.
Figure 5: Discretization scale dependence of the equilibrium times. The plots
show $L/t^{eq}$ vs. $a/L$ for three quantum systems: closed box with $mL=20$
(dots with solid line), closed box with $mL=5$ (diamonds with dashed line) and
periodic box with $mL=20$ (squares with dotted line). The results of the fits
are: $L/t^{eq}=0.006(5)+3.0(1)(a/L)$, $L/t^{eq}=0.004(5)+6.6(1)(a/L)$ and
$L/t^{eq}=0.01(2)+5.4(4)(a/L)$ respectively. For the periodic box, the
superposition of plane-wave energy states has approximately the same $L\Delta
P\approx 9$ as the wave functions in the closed box, which have $L\Delta
P=8.55$. The lines are least-square fits through the data, including the point
at the origin.

The straight lines in Figure 5 are least square fits to the data points, where
the error bars indicate the standard error777As can be seen in Figure 2, the
exponential suppression of the eigenvalues of $\tilde{T}(t)$ is not strictly
linear in the eigenvalue index $s$: the slope is not quite constant and there
are clusters of near degenerate eigenvalues. This introduces additional
uncertainty in the fitted results for $c(t)$ at different values of $1/N$,
which is not included in the size of the error bars. on the fitted coefficient
of $t$ (in Eq. (31)). The slopes of the lines are slightly biased, since also
the point at the origin has been included in the fits. However, the trend in
the data supports the expected behavior that $L/t^{eq}\propto a/L$, and hence
$t^{eq}\rightarrow\infty$ for $a\rightarrow 0$. This is also the case for the
model in which the box has periodic boundary conditions. Here, the equilibrium
time scales are actually systematically shorter than in the box with hard wall
boundaries (the dotted line in Figure 4 lies above the solid line). This is
somewhat surprising, since the (rapid) establishment of quantum equilibrium in
the causal Bohm approach has been linked to the level of chaotic behavior of
the trajectories computed. One might intuitively think that reflecting
boundaries would be important to promote chaotic behavior.

###  3.4 Mass and Momentum Spread Dependence of Equilibrium Times

As is expressed in the scaling relations (1) and (2), the equilibrium time
scale is expected to grow with increasing particle mass and decrease with
increasing momentum spread in the system’s wave function. Figure 5 indeed
shows that $L/t^{eq}$ values for $mL=20$ are smaller than those for $mL=5$. To
further explore the mass dependence of the equilibrium time, $t^{eq}$ is
evaluated for a range of masses on the closed-box system with $L\Delta P=8.55$
and $N=30$. The results in Figure 6 convincingly show that $t^{eq}/L$ grows
linearly with $mL$, as expected.

Figure 6: Mass dependence of the equilbrium times $t^{eq}/L$ for the closed
box model with $N=30$ and $L\Delta P=8.55$. The straight line,
$t^{eq}=-0.6(3)+0.89(2)mL$, is a least-squares fit of the data, including the
point at the origin.
Figure 7: Momentum spread dependence of the equilbrium times $t^{eq}/L$ for
the closed box model with $N=30$ and $mL=20$. The wave functions are composed
of the $N_{k}^{2}$ lowest energy eigenvectors, with $N_{k}=4,5,\dots,8$
($L\Delta P=8.55,10.5,12.4,14.3$ and $16.2$). The straight line ,
$L/t^{eq}=0.004(4)+0.00076(3)(L\Delta P)^{2}$, is a least-squares fit of the
data including the point at the origin.

In analogy with thermal equilibration, the quantum equilibrium time is
expected to also depend on the momentum or energy spread in the system’s wave
function [11, 12]. This is confirmed by the results in Figure 7, which shows
the $\Delta P$ dependence of $L/t^{eq}$. The momentum spread $\Delta P$ is
computed from the standard deviation of the energy: $L\Delta P=L(2m\Delta
E)^{1/2}$ with $\Delta E=(\langle H^{2}\rangle-\langle H\rangle^{2})^{1/2}$. A
least-squares fit $L/t^{eq}$ vs. $(L\Delta P)^{3}$ (inspired by Eq. (1)) gives
a relatively poor result ($R^{2}=0.94$); a linear fit of $\log(L/t^{eq})$ vs.
$\log(\Delta P)$ indicates a lower power of $L\Delta P$, $1.91(9)$ (and
$2.0(1)$ if the largest $\Delta P$ value is left out). This suggest to use a
power 2 instead of 3, which leads to the result shown in Figure 7, where
$R^{2}=0.99$.

Combining this result with that of the previous subsection, the scaling
relation for the quantum equilibrium time using dimensionless parameters takes
the form,

$t^{eq}/L\propto mL/(a/L)(L\Delta P)^{2},$ (32)

which in dimensionful parameters is the result anticipated in Eq. (2) of the
introduction,

$t^{eq}\propto mL\hbar/a(\Delta P)^{2}.$ (33)

The factor $m$ is easy to understand, because for a free particle the
evolution matrix only depends on the combination $\epsilon/m$, i.e., time
scales are naturally proportional to the particle mass. Or, in the particle
picture, the probability to move to a neighboring location (i.e., its typical
velocity) is proportional to $\epsilon/am\propto 1/Lm$.

The appearance of a factor $L/a=N$ is not unreasonable as one would expect
that relaxation requires the particles to travel some system-specific
distance, a mean free path lengths in a multi-particle system or the distance
between two low-probability regions in this single-particle example. Since a
particle will hop at most to a neighboring location, $N$ sets the scale for
the number of particle moves that are needed to cover such a distance. Or,
given the sparse nature of the individual transition matrices, which only have
non-zero off-diagonal elements for nearest neighbors (in at most one
direction), it is also clear that at least $N$ multiplications are needed to
obtain a cumulative matrix in which all elements are larger than zero.

Even though it seems reasonable to expect that an initial state with a large
momentum or energy spread equilibrates faster than one with a smaller spread,
it is not true that initial states with zero momentum spread cannot
equilibrate. Obviously, an energy eigenstate that is real will never
equilibrate, since the transition matrix for such a system has no off-diagonal
elements. As is the case in the causal Bohm formulation, particles do not move
when the system is in a real-valued energy eigenstate. However, on a space
with periodic boundary conditions, where energy eigenstates can be take as
complex-valued plane waves, particles do move and perhaps surprisingly, even a
system with a single plane wave as initial state is found to equilibrate. This
can happen owing to the discrete time formulation in which the probability
current (8) is split over two time steps. When expanding the transition
probabilities in orders of the time step size $\epsilon$, one will find
subleading terms that allow transitions between neighboring locations (in the
direction orthogonal to the direction the particle moves in), even though the
wave function values at these locations evaluated at the same time are equal.
In the example of a single plane wave, the equilibrium time scale is found to
increase with decreasing time step size - as expected.

As the above comments already suggest, it is not easy to find initial states
that cannot support quantum equilibration, i.e., states for which the
transition matrix is non-primitive. Wave functions composed of random
superpositions of energy states, but such that they possess a rotation or
mirror symmetry that is preserved during time evolution, still lead to
primitive transition matrices and these systems readily equilibrate. Besides
real-valued energy eigenstates, it turns out that initial (Gaussian) states
that are sufficiently localized start out with a transition matrix that has
multiple eigenvalues $1$. However, once the wave packet spreads out to cover
enough of the box volume, such that the transition matrix develops off-
diagonal elements connecting neighboring locations everywhere in the system,
also here transition matrices become primitive, only have a single isolated
eigenvalue $1$, and also this system readily equilibrates.

##  4 Discusion

The stochastic version of the de Broglie-Bohm formulation of quantum mechanics
introduced by Bell, can be rigorously defined on discrete space-time [9]. As
was shown in section 2, the dBBB dynamics for the probability distribution of
particle locations resembles a Markov process with transition matrices $T(t)$
that are slowly varying in time. This analogy with a Markov process, supported
by proven properties of the backward product of such matrices [19] suggests
that the dBBB dynamics (for sufficiently small time step sizes) will be weakly
ergodic if the transition matrices are primitive (i.e., if they are such that
all matrix elements of $T^{k}(t)$ are positive for a sufficiently large power
$k$). Weak ergodicity then implies that for large enough times, the cumulative
transition matrix $\tilde{T}(t)$ is dominated by a single eigenmode, as shown
in Eqs. (23) and (25), which in turn implies that the system will relax to
quantum equilibrium irrespective of the initial probabilities attributed to
the particle locations in the ensemble.

Direct calculations of the spectral decomposition of $\tilde{T}(t)$ for a free
particle in a 2D box confirm that this dynamical relaxation to quantum
equilibrium indeed happens for many variations of the system’s wave function -
in fact, it proves to be difficult to find (non-static) wave functions for
which equilibrium relaxation does not happen. The example of a spreading
Gaussian wave package, which starts to equilibrate once the package has
sufficiently spread out over the box volume, illustrates that a necessary
condition for equilibration (i.e., for the system’s transition matrix to be
primitive) is that the (non-static) wave function is non-zero everywhere, such
that there is a non-zero transition probability for the particle to move, at
every location in the system – it may well be that this condition is also
sufficient.

The numerical results furthermore confirm that the equilibrium time scale
diverges when the lattice distance approaches zero. This has to be the case,
since in this limit the stochastic Bell dynamics must reproduce [7] the causal
Bohm dynamics for which quantum equilibration cannot happen without
introducing coarse-graining [11, 12, 13]. Since the quantum world in the
discrete dBBB formulation is naturally discrete (and finite), such an
additional coarse-grainig step, or interpretation, is not needed.

The computed values of $t^{eq}$ are sufficiently accurate to establish the
scaling behavior of equilibrium times shown in Eq. (2), which has the expected
dependence on particle mass and lattice distance. Unlike in Valentini’s
proposal shown in Eq. (1), there is an additional dependence on the system
size, with a correspondingly weaker dependence on the momentum spread.

Even though the simulations in this paper were performed on a very simple
quantum system, the proven properties of slowly evolving sequences of
transition matrices, corroborated by these simulation results are sufficiently
encouraging to expect that, for typical quantum systems of interest, the
stochastic dBBB dynamics will, after a system-dependent relaxation time,
reproduce the probability distributions computed with quantum mechanics
irrespective of the details of the initial particle distribution. This
dynamical relaxation to quantum equilibrium elegantly counters one of the more
persistent objections against Bohm-type interpretations of quantum mechanics.

### Acknowledgements

I would like to thank Marc Artzrouni for helpful comments and discussion.

## References

  * [1] D. Bohm and B.J. Hiley, “The Undivided Universe”, Routledge Publ. (1993).
  * [2] S. Goldstein, “Bohmian Mechanics”, The Stanford Encyclopedia of Philosophy (Fall 2021 Edition), Edward N. Zalta (ed.),
URL = $<$https://plato.stanford.edu/archives/fall2021/entries/qm-bohm/$>$.

  * [3] J. Barrett, “The Conceptual Foundations of Quantum Mechanics”, Oxford Univ. Press (2019).
  * [4] L. de Broglie, “Tentative d’Interpretation Causale et Non-linéaire de la Mécanique Ondulatoire”, Gauthier-Villars, Paris (1956).
  * [5] D. Bohm, Phys. Rev. 85 (1952) 166 and 180.
  * [6] J.S. Bell, “Speakable and unspeakable in quantum mechanics”, Cambridge University Press (1987), ch. 19.
  * [7] J.C. Vink, Phys. Rev. A48 (1993) 1808.
  * [8] J.C. Vink, Found. Phys. 48 (2018) 209.
  * [9] J.C. Vink, “Spin and Contextuality in Extended de Broglie-Bohm-Bell Quantum Mechanics”, arXiv:2203.14034 [quant-ph].
  * [10] D. Bohm and B.J. Hiley, “The Undivided Universe”, Routledge Publ. (1993), ch. 9.
  * [11] A. Valentini, On the Pilot-Wave Theory of Classical, Quantum and Subquantum Physics, PhD. Thesis, International School for Advanced Studies, Trieste (1992).
  * [12] A. Valentini, in: Chance in Physics: Foundations and Perspectives, eds. J. Bricmont et al. (Springer, Berlin, 2001) [arXiv:quant-ph/0104067].
  * [13] A. Valentini and H. Westman, Proc. R. Soc. A 461 (2005) 253 [arXiv:quantph/ 0403034].
  * [14] S. Colin and W. Struyve, New Journal of Physics 12 (2010) 043008.
  * [15] D. Dürr, S. Goldstein and N. Zanghì, J. Stat. Phys. 67 (1992) 843 [arXiv:quant- ph/0308039].
  * [16] G.H. Golub, and F. Van Loan, “Matrix Computations”, The Johh Hopkins University Press, third edition, 1996.
  * [17] S. Sternberg, “Dynamical Systems”, Dover Publications; Illustrated edition (June 10, 2014), Ch. 9.
  * [18] E. Seneta, “Non-negative Matrices and Markov Chains”, 2nd Edition, Springer-Verlag New York Inc. (1981).
  * [19] M. Artzrouni, Linear Algebra Appl. 145(5) (1991) 33.


#  How Generative Models Improve LOS Estimation in 6G Non-Terrestrial Networks

Saira Bano  Achilles Machumilane  Pietro Cassarà  Alberto Gotta A. Machumilane
and S.Bano are with the Department of Information Engineering (DII) of the
University of Pisa and the Institute of Information Science and Technologies
(ISTI), CNR, Pisa - e-mails<EMAIL_ADDRESS>A. Gotta and P.
Cassarà are with the Institute of Information Science and Technologies (ISTI),
CNR, Pisa - e-mails<EMAIL_ADDRESS>
###### Abstract

With the advent of 5G and the anticipated arrival of 6G, there has been a
growing research interest in combining mobile networks with Non-Terrestrial
Network platforms such as low earth orbit satellites and Geosynchronous
Equatorial Orbit satellites to provide broader coverage for a wide range of
applications. However, integrating these platforms is challenging because
Line-Of-Sight (LOS) estimation is required for both inter satellite and
satellite-to-terrestrial segment links. Machine Learning (ML) techniques have
shown promise in channel modeling and LOS estimation, but they require large
datasets for model training, which can be difficult to obtain. In addition,
network operators may be reluctant to disclose their network data due to
privacy concerns. Therefore, alternative data collection techniques are
needed. In this paper, a framework is proposed that uses generative models to
generate synthetic data for LOS estimation in non-terrestrial 6G networks.
Specifically, the authors show that generative models can be trained with a
small available dataset to generate large datasets that can be used to train
ML models for LOS estimation. Furthermore, since the generated synthetic data
does not contain identifying information of the original dataset, it can be
made publicly available without violating privacy.

###### Index Terms:

NTNs, Satellites, Channel Modeling, Generative Models.

##  1 Introduction

In recent years, the Third Generation Paternship Project (3GPP) has envisioned
the integration of 5G and 6G mobile networks with Non-Terrestrial Network
(NTN) technologies such as Low Earth Orbit (LEO) satellites, Unmanned Aerial
Systems, and High Altitude Platforms as a promising solution for providing
ubiquitous coverage in inaccessible areas [1]. This integration will
significantly improve network connectivity, accessibility and data rates and
will also support a wide range of applications and services, including rescue
missions, remote monitoring, and goods delivery [1]. The main communications
challenge with this integration is modeling the Line-of-Sight (LOS)
availability of the link between the satellite and terrestrial segments since
satellite communications require a clear LOS that can be blocked by obstacles
such as buildings and vegetation, resulting in signal blockage, diffraction,
and reflection. The elevation angle of the satellite also affects the LOS,
with lower angles less likely to result in a LOS that yields to weak or no
signal. Existing 3GPP and International Telecommunication Union (ITU) models
define channel parameters based on elevation angle, frequency, and deployment
scenarios [2]. However, certain critical parameters, such as LOS probability,
have no temporal correlation and therefore do not account for satellite
motion. Consequently, changes in LOS/ Non-Line-of-Sight (NLOS) states may be
inaccurately represented, which could complicate the impact of radio mobility
on 5G-based satellite networks. Therefore, it is important to model these
changes more accurately, considering the impact of satellite and user
mobility. Statistical modeling is one potential solution for LOS prediction
and transmission of data only when conditions are favorable. However, this
approach can be tedious and time-consuming. Recently, researchers have
explored the use of Machine Learning (ML) techniques for LOS estimation.
Although these methods have shown promising results, they typically require
large amounts of data to train the ML models, which can be difficult and
expensive to obtain. In addition, privacy concerns of network operators for
sensitive data have made it difficult to obtain large data sets for ML model
training. To address this problem, generative models such as Generative Neural
Networks (GNNs) have been proposed to generate synthetic datasets using small
available real datasets. These networks have numerous applications in various
fields, including communications, manufacturing, and healthcare. They are
mainly used to create new images, music, text, and videos. The two most
commonly used GNNs are Generative Adversarial Network (GAN) and Variational
Autoencoder (VAE), which are used to generate synthetic data that mimics the
quality and statistical distribution of real-world data. This study aims to
demonstrate the effectiveness of generating models in generating new satellite
channel data for LOS estimation. To accomplish this, the study utilizes the
channel models provided by ITU and creates datasets of LOS and NLOS traces for
various satellite elevation angles that can account for satellite mobility.
These datasets are then used to train the generative models, namely GAN and
VAE. By using these models, this work shows that synthetic data can be
generated that closely resembles the original data and retains the statistical
distribution, thereby providing a solution to the limited availability of
training data in ML models for LOS estimation. Furthermore the generated data
would be free of legal, privacy, and security issues. People can use them for
academic and research purposes and mitigate the difficulties of obtaining
real-world data. We comprehensively evaluate, compare, and analyze the
performance levels of the proposed generative models using various statistical
parameters. The experimental results show that the proposed models are robust
for LOS estimation using the generated dataset.

##  2 Related Work

This section provides various techniques proposed in the literature for LOS
estimation and synthetic data generation. For LOS estimation, the authors of
[3] provide a theoretical model that calculates the likelihood of having LOS
when there are no clouds. They showed that the LOS availability depends on the
height of the ground station and the satellite elevation angle. A ML-based
method for NLOS identification is proposed in [4] , which uses two ML
algorithms Support Vector Machine (SVM) and Logistic Regression (LR) to
predict the NLOS using data from the global navigation system. Authors in [5]
use time-varying angular information of a channel to train ML algorithms for
LOS identification in Vehicle to Vehicle (V2V) communication. The problem of
these ML-based models is that they require a huge amount of measurement data
that may not be available or too small to train a model. Moreover, the
measurement campaigns to obtain real channel data may be costly and time-
consuming. Our work proposes a solution by showing that it is possible to use
synthetic data instead of real channel data for LOS estimation. The synthetic
data can be obtained by using generative models using a small amount of
available real channel data or from publicly available synthetic data of the
channel characteristics.

Numerous data-driven methods for generating synthetic time series data are
also studied in the literature. One example is the approach proposed by [6],
which uses GAN to generate energy consumption data. GAN uses a discriminator
to indirectly train the generator network, enabling it to generate synthetic
data. In [7], a general modeling approach is presented based on training a
generative neural network with data. The proposed generative model includes
two stages: Prediction of the link state (LOS, NLOS, or no available link) and
subsequent use of a conditional VAE to generate path losses, delays, and
arrival and departure angles for all propagation paths given the previously
predicted link state. In [8], the authors used GAN for path loss prediction
for satellite images using the dataset from raytracing simulations.

##  3 System Model

###  3.1 Channel Model

The 3GPP has defined two main architectures for integrating cellular networks
and NTN [1]: the transparent and regenerative architecture. In the transparent
architecture, the satellite acts like a Radio Frequency (RF) repeater,
transparently forwarding traffic between the UE and the New Generation (NG)-
Radio Access Network (RAN) of the mobile network. In the regenerative
architecture, on the other hand, the satellite has gNB capabilities with
onboard gNB- Distributed Unit (DU) and the gNB- Centralized Unit (CU) deployed
in the ground segment of the mobile network, as shown in Figure 2. In this
work, we use regenerative architecture. As explained earlier, the link between
the satellite and the UE can be in LOS or NLOS.

Figure 1: Reference Scenario.

In this work, we investigate the generation of synthetic data that can be used
to estimate the LOS probability for the link between the satellite and the UE.
We use the channel model provided by the ITU [9] and modify it to a simplified
Lutz model [10] , [3]. According to these models, a channel between a
satellite and any land mobile terminal can be in either a good (G) or bad (B)
state, as the signal power varies as a result of shadowing and multipath
caused by signal obstructions and reflections from obstacles such as
buildings, vegetation, and the ground. We assume that a channel is in a good
state if there is a LOS, and in a bad state otherwise. The ITU [9]
recommendations provide several statistical parameters that can be used to
calculate the average duration of each of the two states in different
environments, including the type of terrain (urban, suburban, and rural),
elevation angles, and frequencies. These parameters were collected in a city
in France and are shown in Table I. They include each state’s mean, standard
deviation, and minimum state lengths in meters. For this study, we utilized
the parameters for a dense urban environment at 2.2 GHz and computed the
transition probabilities, which denote the likelihood of transitioning from
one state to another given an initial state. Table II lists the transition
probabilities we derived, and we modeled the state transitions using a two-
state Markov process. To train our generative models, we generated LOS-NLOS
traces from the statistical model. The generative models learned the latent
space and distributions of the training data, allowing them to produce
synthetic traces with similar statistical characteristics to the original
dataset. Alternatively, without access to a channel model, the training data
could be obtained from traces acquired through real-time transmissions using
feedback mechanisms, as demonstrated in [11]

TABLE I: Satellite link parameters and transition probabilities for the dense
urban environment in France at 2.2 GHz.

Elevation | $\mu_{G,B}$ | $\sigma_{G,B}$ | $dur_{minG,B}$ |  $P(B\to\,G)$ (g) |  $P(G\to\,B)$ (b)
---|---|---|---|---|---
$20^{\circ}$ | 2.0042, 3.6890 | 1.2049, 0.9796 | 3.9889, 10.3114 | 0.00014310 | 0.00047466
$30^{\circ}$ | 2.7332, 2.7582 | 1.1030, 1.2210 | 7.3174, 5.7276 | 0.00024460 | 0.00027570
$45^{\circ}$ | 3.0639, 2.9108 | 1.6980, 1.2602 | 10.0, 6.0 | 0.00020318 | 0.00007556
$60^{\circ}$ | 2.8135, 2.0211 | 1.9595, 0.6568 | 10.0, 1.9126 | 0.00105161 | 0.00010797
$70^{\circ}$ | 4.2919, 2.1012 | 2.4703, 1.0341 | 118.3312, 4.8569 | 0.00052923 | $2.76683\cdot 10^{-}6$

###  3.2 Generative Neural Networks

####  3.2.1 Generative Adversarial Networks (GANs)

In [12] Ian Goodfellow introduces GAN, a type of ML algorithm that has gained
popularity because of its ability to generate high-quality synthetic images,
videos, and even sounds that closely resemble real-world data. GANs use an
unsupervised learning approach to detect patterns in the input data and
generate new samples with the same distribution as the original data. GANs
train two neural networks, the generator network and the discriminator. The
generator network generates the fake data by taking samples from a random
distribution and converting them into data that resembles real data. In
contrast, the discriminator network tries to distinguish between the real data
and the fake data generated by the generator. The generator does not have
direct access to real samples and learns only through interaction with the
discriminator, which has access to synthetic and real samples. Training GANs
is about finding the parameters of a discriminator that maximize its
classification accuracy and finding the parameters of a generator that
maximally confuse the discriminator. As the generator network improves, more
realistic data is generated, making it increasingly difficult for the
discriminator network to distinguish between real and fake data. This work
uses GANs to generate LOS and NLOS traces for 6G-NTN channels.

####  3.2.2 VAE - Variational Auto-Encoders

Variational autoencoders are a powerful type of Deep Learning (DL) algorithm
that can be used for unsupervised learning and generative modeling. They are
particularly useful for learning a compressed dataset representation that
captures the data’s underlying structure in a low-dimensional space. VAEs have
many practical applications, such as image and speech recognition and natural
language processing. VAEs are capable of generating new data that is similar
to training data. VAEs use a two-part architecture consisting of an encoder
and a decoder. The encoder maps the input data to a distribution in latent
space, while the decoder maps points in latent space back to the original data
space. By varying the sampled points in the latent space, new data that is
similar to, but not identical to, the training data can be generated. The VAE
architecture includes an encoder that transforms input data into a Gaussian
distribution over the latent space using convolutional or dense layers. The
decoder takes samples from the latent space and maps them back to the original
data space using similar layers. VAE aims to minimize the difference between
input data and the decoder’s output while maximizing the Gaussian
distribution’s parameter likelihood, which defines the latent space.

####  3.2.3 Conditional Tabular GAN (CTGAN) and Tabular VAE (TVAE)

This work uses Tabular Conditional GAN and Tabular VAE, both presented by Lei
Xu et al. in [13], which are part of the Synthetic Data Vault (SDV) package in
TensorFlow. These models were chosen because of their ability to process
tabular data and allow training of a single model that can generate synthetic
data for any available channel between UE and satellite. CTGAN and TVAE are
capable of capturing the distribution of each column in the tabular data. In
our case, each column contains LOS-NLOS traces for each channel or elevation
angle. We consider three elevation angles:
$70^{\circ}$,$60^{\circ}$,$45^{\circ}$. Knowing the distribution of each
column, the models used in this work can generate synthetic data for each
column simultaneously, based on the distribution of each column, which saves
time compared to training a model for each channel. However, since the LOS
probability changes with the elevation angle of the satellite, different
channel models are needed for each angle. For LEO satellites, this means that
a different model is required for each elevation angle. Since satellites are
visible from certain locations at certain angles, the terminal must switch to
different satellites as they become visible. However, the proposed approach
uses a single model that can generate synthetic data for all channels or
satellite elevation angles that the UE connects to the satellite.

##  4 Performance Evaluation

In this part, we evaluate and compare the accuracy of the generative models in
generating synthetic data that closely resemble real data. First, we train the
models to generate synthetic data for the LOS and NLOS traces. Then we
evaluate the trained models and compare the similarity between the generated
data and the real data. We use Wasserstein Distance, Kolmogorov-Smirnov (KS)
test, and Kullback-Leibler (KL) divergence as measures for comparison and
evaluation. These measures allow us to determine how well the probability
distribution of the synthetic data matches that of the real data.

###  4.1 Model Training

As explained previously, we trained our models using the training dataset
consisting of the LOS-NLOS -traces obtained using the state transition
probabilities in I. For this work, we used the probabilities at three
elevation angles: $45^{\circ}$, $60^{\circ}$, $70^{\circ}$. We assume the
satellite is visible from our reference UE at these three elevation angles.
The trained models should produce synthetic data that can estimate or predict
the presence of the LOS at these three angles. The training data set is a
table with 100,000 rows of LOS and NLOS traces in three columns, with the
traces of each angle in each column. We trained each model for 100 epochs with
a batch size of 50 and a learning rate of $2e-4$.

###  4.2 Performance Metrics

We used the following performance metrics to evaluate the effectiveness of our
geneative models for LOS estimation:

  * •

Wasserstein distance is a metric for measuring the distance between two
probability distributions, that is, for measuring the similarity of the
probability distribution of the synthetic data to that of the real data. The
smaller the Wasserstein distance between the synthetic and real data, the more
similar the probability distributions, which means that the generative model
has produced synthetic data of high quality.

  * •

The KS test is a statistical measure that gets the distance between two
empirical Cumulative Distributions Function (CDF), a popular non-parametric
measure used in statistics. We use the KS test to measure how far the CDF of
the synthetic data is from the real data. It is usually presented as a
complementary measure, i.e., a 1-difference in CDF. Thus, the higher the
value, the more similar the synthetic data is to the real data.

  * •

The KL divergence measures how much two probability distributions differ from
each other. The KL divergence between two probability distributions, P and Q,
is calculated as the sum of the log difference between the probabilities of
each value in P and Q multiplied by the probability of that value in P. The KL
divergence is always non-negative and equal to zero if and only if the two
distributions are identical. A low KL divergence indicates that the generative
model has produced synthetic data of high quality that are similar to the real
data.

###  4.3 Results

TABLE II: Distance between real and synthetic data for CTGAN.

Metric | $70^{\circ}$ | $60^{\circ}$ | $45^{\circ}$
---|---|---|---
| mean | variance | mean | variance | mean | variance
KS-test | 0.9697 | 7.84e-08 | 0.9646 | 8.41e-08 | 0.9632 | 6.5025e-08
Wasserstein Distance | 0.0591 | 4.096e-07 | 0.0700 | 1.21e-08 | 0.0743 | 7.48e-08
KL-Divergence | 0.0199 | 7.85e-08 | 0.0101 | 1.42e-09 | 0.0034 | 6.87e-10

TABLE III: Distance between real and synthetic data for TVAE.

| $70^{\circ}$ | $60^{\circ}$ | $45^{\circ}$
---|---|---|---
| mean | variance | mean | variance | mean | variance
KS-test | 0.9651 | 9.6e-08 | 0.9763 | 8.1e-09 | 0.9570 | 1.156e-07
Wasserstein Distance | 0.0591 | 4.096e-07 | 0.0701 | 1.21e-08 | 0.0743 | 7.84e-08
KL-Divergence | 0.0247 | 8.969e-09 | 0.0039 | 4.242e-08 | 0.0052 | 7.904e-08

Figure 2: Distribution of real vs synthetic data for CTGAN and TVAE at
different satellite elevation angles.

####  4.3.1 Mean distance between real and synthetic data

In Tables II and III, we show the mean and variance distance between synthetic
and real data at different elevation angles using the selected evaluation
metrics. To determine the mean distance, we first created a test dataset with
100,000 samples from LOS-NLOS, the same as the one created for the training
dataset. Then we use the trained CTGAN and TVAE models to create synthetic
datasets of the same size. We repeat this procedure 50 times and calculate the
KS test, Wasserstein distance, and KL divergence distances between each
synthetic data set and the real data, obtaining a set of fifty distances for
each metric for each model. We then calculated the mean and variance of these
distances, and the results are shown in Table II for CTGAN and Table III for
TVAE. As can be seen from the tables, both the Wasserstein and KL divergences
are very low for both models at all three elevation angles. Similarly, the
values for the KS test are very high for both models, ranging from 0.9763 to
0.9570, indicating that our models produce high-quality synthetic data with
distributions that are very close to the real data.

####  4.3.2 Variance of the distance between real and synthetic data

The variance of the distance between real and synthetic data shows whether the
models can produce similar synthetic data given different specifications. The
variance shows how stable and robust the models are in generating the
synthetic data. The lower the variance, the higher the stability and
robustness of the models. The results presented in II and III show low
variances, as low as 1.42 x $10^{-}9$, which means that the distances between
the different synthetic data sets generated in different instances vary
little. This shows that our models are stable, robust, and reliable in
generating LOS/NLOS estimates at different instances and different elevation
angles.

####  4.3.3 Distributions of real data and synthetic data

In Figure 2 we have shown the distribution of real and synthetic data for
CTGAN and TVAE at different satellite elevation angles. In generating our
dataset, LOS and NLOS were coded as 1 and -1, respectively. As expected, the
LOS probability decreases with decreasing elevation angle and is highest at
$70^{\circ}$ and lowest at $45^{\circ}$. The results show that this variation
is similar for real and synthetic data, implying that our models can correctly
estimate the LOS and NLOS probabilities at different elevation angles.

####  4.3.4 Comparative Analysis

The following briefly compares the generative models used in this paper. As
mentioned earlier, both GANs and VAEs are DL models used to generate synthetic
data. However, they differ in the way they generate synthetic data. GANs use a
generator that learns to generate synthetic data that is indistinguishable
from real data. In contrast, VAEs use a probabilistic encoder and decoder
network to generate synthetic data by mapping the real data to a low-
dimensional latent space and then mapping the latent space back to the
original space. In general, GANs are known to be more difficult to train and
may suffer from ”mode collapse” compared to VAE, where the generator produces
only a limited amount of synthetic data. The results in Figure 3 show that
VAEs exhibit faster stability with fewer runs over the training dataset, while
GANs require longer training periods to achieve stability. The Figure shows
the KL divergence and Wasserstein distance between the test dataset and the
generated dataset for the elevation angle of $70^{\circ}$, indicating that the
data distribution generated by the TVAEs is very similar to that of the real
data set, even after only a few runs over the training dataset. Thus, the
TVAEs perform better than the GAN for the LOS estimation for the given
dataset. We also consider the training times of the CTGAN and TVAE. CTGAN
requires more training time than TVAE. For example, in our simulations, CTGAN
required 1.18 hours in the given training environment, while TVAE required 47
minutes. This shows that for the given problem of LOS estimation and with the
given dataset, TVAE is the best option in terms of performance and training
time.

Figure 3: Comparison of CTGAN and TVAE in terms of KL and Wasserstein distance
against number of epochs for the channel at $70^{\circ}$

##  5 Conclusion

In this study, a DL technique was used to generate synthetic data for LEO
satellites operating in non-terrestrial 6G networks, considering both LOS and
NLOS scenarios. Two generative network variants, CTGAN and TVAE, were used
because they are well suited for tabular data. The simulation results showed
that the generative models mimicked the real dataset very well and
successfully estimated the LOS probabilities over multiple satellite channels.
The statistical metrics used to measure the performance of the models showed
that both models were comparable. However, TVAE outperformed CTGAN in terms of
training time, while CTGAN exhibited some oscillations in training due to the
min-max-adversarial effect.

## References

  * [1] 3GPP, “Technical specification group radio access network; solutions for nr to support non-terrestrial networks (ntn): Tr 38.821 v16.1.0 (2021-05), (release 16).”
  * [2] E. Juan, M. Lauridsen, J. Wigard, and P. E. Mogensen, “A time-correlated channel state model for 5g new radio mobility studies in leo satellite networks,” in _2021 IEEE 93rd Vehicular Technology Conference (VTC2021-Spring)_ , 2021, pp. 1–5.
  * [3] H. Bischel, M. Werner, and E. Lutz, “Elevation-dependent channel model and satellite diversity for ngso s-pcns,” in _Proceedings of Vehicular Technology Conference-VTC_ , vol. 2. IEEE, 1996, pp. 1038–1042.
  * [4] Y. Sun and L. Fu, “Stacking ensemble learning for non-line-of-sight detection of global navigation satellite system,” _IEEE Transactions on Instrumentation and Measurement_ , vol. 71, pp. 1–10, 2022.
  * [5] C. Huang, A. F. Molisch, R. He, R. Wang, P. Tang, B. Ai, and Z. Zhong, “Machine learning-enabled los/nlos identification for mimo systems in dynamic environments,” _IEEE Transactions on Wireless Communications_ , vol. 19, no. 6, pp. 3643–3657, 2020.
  * [6] C. Zhang, S. R. Kuppannagari, R. Kannan, and V. K. Prasanna, “Generative adversarial network for synthetic time series data generation in smart grids,” in _2018 IEEE international conference on communications, control, and computing technologies for smart grids (SmartGridComm)_. IEEE, 2018, pp. 1–6.
  * [7] W. Xia, S. Rangan, M. Mezzavilla, A. Lozano, G. Geraci, V. Semkin, and G. Loianno, “Generative neural network channel modeling for millimeter-wave uav communication,” _IEEE Transactions on Wireless Communications_ , vol. 21, no. 11, pp. 9417–9431, 2022.
  * [8] A. Marey, M. Bal, H. F. Ates, and B. K. Gunturk, “Pl-gan: Path loss prediction using generative adversarial networks,” _IEEE Access_ , vol. 10, pp. 90 474–90 480, 2022.
  * [9] I. Recommendation, “Propagation data required for the design of earth-space land mobile telecommunication systems,” _International Telecommunication Union_ , pp. 681–10, 2017.
  * [10] E. Lutz, D. Cygan, M. Dippold, F. Dolainsky, and W. Papke, “The land mobile satellite communication channel-recording, statistics, and channel model,” _IEEE transactions on vehicular technology_ , vol. 40, no. 2, pp. 375–386, 1991.
  * [11] A. Machumilane, A. Gotta, P. Cassarà, and M. Bacco, “A path-aware scheduler for air-to-ground multipath multimedia delivery in real time,” _IEEE Communications Magazine_ , vol. 60, no. 9, pp. 54–58, 2022.
  * [12] I. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio, “Generative adversarial networks,” _Communications of the ACM_ , vol. 63, no. 11, pp. 139–144, 2020.
  * [13] L. Xu, M. Skoularidou, A. Cuesta-Infante, and K. Veeramachaneni, “Modeling tabular data using conditional gan,” _Advances in Neural Information Processing Systems_ , vol. 32, 2019.

  *[3GPP]: Third Generation Paternship Project
  *[NTN]: Non-Terrestrial Network
  *[LEO]: Low Earth Orbit
  *[LOS]: Line-of-Sight
  *[ITU]: International Telecommunication Union
  *[NLOS]: Non-Line-of-Sight
  *[ML]: Machine Learning
  *[GNNs]: Generative Neural Networks
  *[GAN]: Generative Adversarial Network
  *[VAE]: Variational Autoencoder
  *[SVM]: Support Vector Machine
  *[LR]: Logistic Regression
  *[V2V]: Vehicle to Vehicle
  *[RF]: Radio Frequency
  *[NG]: New Generation
  *[RAN]: Radio Access Network
  *[DU]: Distributed Unit
  *[CU]: Centralized Unit
  *[DL]: Deep Learning
  *[CDF]: Cumulative Distributions Function