# One- and two-dimensional solitons in spin-orbit-coupled Bose-Einstein
condensates with fractional kinetic energy

Hidetsugu Sakaguchi1 and Boris A. Malomed2,3 1Department of Applied Science
for Electronics and Materials, Interdisciplinary Graduate School of
Engineering Sciences, Kyushu University, Kasuga, Fukuoka 816-8580, Japan
2Department of Physical Electronics, School of Electrical Engineering, Faculty
of Engineering, and Center for Light-Matter Interaction, Tel Aviv University,
P.O. Box 39040 Tel Aviv, Israel
3Instituto de Alta Investigación, Universidad de Tarapacá, Casilla 7D, Arica,
Chile

###### Abstract

We address effects of spin-orbit coupling (SOC), phenomenologically added to a
two-component Bose-Einstein condensate composed of particles moving by Lévy
flights, in one- and two-dimensional (1D) and (2D) settings. The corresponding
system of coupled Gross-Pitaevskii equations includes fractional kinetic-
energy operators, characterized by the Lévy index, $\alpha<2$ (the normal
kinetic energy corresponds to $\alpha=2$). The SOC terms, with strength
$\lambda$, produce strong effects in the 2D case: they create families of
_stable solitons_ of the semi-vortex (SV) and mixed-mode (MM) types in the
interval of $1<\alpha<2$, where the supercritical collapse does not admit the
existence of stable solitons in the absence of the SOC. At $\lambda\rightarrow
0$, amplitudes of these solitons vanish $\sim\lambda^{1/(\alpha-1)}$.

##  I Introduction

The concept of derivatives of fractional orders was proposed in 1823 by Niels
Henrik Abel Abel . In physics, the implementation of the fractionality was
proposed by N. Laskin Laskin , who had introduced a linear fractional
differential operator (actually, it is defined by an integral expression, see
Eq. (2) below) which replaces the usual quantum-mechanical kinetic energy. The
result is the fractional Schrödinger equation (FSE) for wave function
$\psi(x,t)$, with the Lévy index $\alpha$, time $t$, spatial coordinate $x$,
and potential $U(x)$. The normalized form of the one-dimensional (1D) FSE is

$i\frac{\partial\psi}{\partial
t}=\frac{1}{2}\left(-\frac{\partial^{2}}{\partial
x^{2}}\right)^{\alpha/2}\psi+U(x)\psi.$ (1)

It was derived, by means of the Feynman’s path-integral method (tantamount to
stochastic quantization), for a particle whose stochastic motion is realized
not by the usual Brownian regime, but rather by Lévy flights. This means that
the average distance of the corresponding randomly walking classical particle
from the initial position growth with time as $t^{1/\alpha}$ Levy flights . In
the case of $\alpha=2$, this is the usual random-walk law for a Brownian
particle, and, accordingly, Eq. (1) amounts to the usual Schrödinger equation.
The Lévy-flight regime, corresponding to $\alpha<2$, implies diffusive walk
faster than Brownian, which is performed, at the classical level, by random
leaps. The quantum counterpart of this regime, represented by Eq. (1) in the
1D case, is the basis of fractional quantum mechanics Laskin book .

While there are different formal definitions of fractional derivatives, the
one which is relevant in the context of quantum mechanics is the Riesz
derivative, which is based on the juxtaposition of the direct and inverse
Fourier transforms Riesz-original ; Riesz :

$\left(-\frac{\partial^{2}}{\partial
x^{2}}\right)^{\alpha/2}\psi(x)=\frac{1}{2\pi}\int_{-\infty}^{+\infty}|p|^{\alpha}dp\int_{-\infty}^{+\infty}\psi(\xi)\exp\left(ip(x-\xi)\right)d\xi.$
(2)

Similarly, in the two-dimensional (2D) version of the FSE the fractional
kinetic-energy operator is defined as a fractional power of the Laplacian,
i.e.,

$\displaystyle\left(-\nabla^{2}\right)^{\alpha/2}\psi(x,y)\equiv\left(-\frac{\partial^{2}}{\partial
x^{2}}-\frac{\partial^{2}}{\partial y^{2}}\right)^{\alpha/2}\psi(x,y)$
$\displaystyle=\frac{1}{(2\pi)^{2}}\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}\left(p^{2}+q^{2}\right)^{\alpha/2}dpdq\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}\psi(\xi,\eta)\exp\left[ip(x-\xi)+iq(y-\eta)\right]d\xi
d\eta,$ (3)

The similarity of the quantum-mechanical Schrödinger equation to the wave-
propagation equation under the action of the paraxial diffraction in optics
suggests a possibility to simulate FSE in optical setups, with time $t$
replaced by the propagation distance, $z$. This option was proposed by S.
Longhi Longhi , who elaborated a scheme based on a Fabry-Perot resonator, in
which the 4f setup (two lenses with four focal lengths 4f ) is employed to
approximate the fractional diffraction by transforming the spatial structure
of the light beam into its Fourier counterpart, applying the fractional
diffraction in its straightforward form in the Fourier layer (as proposed in
early works early1 ; early2 ), and then transforming the result back into the
spatial domain. This physical picture also corresponds to the definition of
the Riesz derivative in Eqs. (2) and (3). In other physical contexts, it was
proposed to create a Lévy crystal, which provides discretization of the FSE in
a dynamical lattice Stickler , and to realize an effective FSE in a condensate
of polaritons in a semiconductor cavity Pinsker .

The emulation of the fractional diffraction in optics suggests to take into
regard the natural Kerr nonlinearity of optical media, and thus to address a
possibility of the existence of fractional solitons. Various aspects of this
topic were a subject of many recent theoretical works Zhong -Zeng 2 , see also
a review in Ref. review . In this case, the FSE in free space (without the
external potential) is replaced by one containing the self-focusing cubic
term:

$i\frac{\partial\psi}{\partial
z}=\frac{1}{2}\left(-\nabla^{2}\right)^{\alpha/2}\psi-|\psi|^{2}\psi,$ (4)

where the fractional-diffraction operator
$\left(-\nabla^{2}\right)^{\alpha/2}$ may be one- or two-dimensional, as
defined in Eqs. (2) and (3), respectively. In 1D, Eq. (4) gives rise to stable
solitons at $\alpha>1$ and to the critical collapse at $\alpha=1$, if the norm
of a localized input exceeds a critical value collapse ; review . In 2D, the
critical collapse takes place at $\alpha=2$ (it is the usual critical collapse
in the 2D nonlinear Schrödinger equation with normal diffraction Sulem ;
Fibich ), and the supercritical collapse (for which the critical value is
zero) at $\alpha<2$.

Another natural possibility is to consider a Bose-Einstein condensate (BEC) in
an ultracold gas of quantum particles moving by Lévy flights. In this case,
Eq. (4), with $z$ replaced by time $t$, represents the scaled form of the
corresponding free-space Gross-Pitaevskii equation (GPE), assuming attractive
interaction between particles. As well as in the realization in optics, a
crucially important issue is a possibility of the stabilization of solutions
for matter-wave solitons, created by Eq. (4), against the collapse. In this
connection, it is relevant to mention that the spin-orbit coupling (SOC)
efficiently stabilizes 2D matter-wave solitons of two types, viz., semi-
vortices (SVs) and mixed modes (MMs), against the critical collapse, in the
framework of the two-component system of GPEs with the normal kinetic-energy
operators and cubic self- and cross-attraction terms we ; Sherman , see also a
review in Ref. EPL . SVs are composite states with zero vorticity in one
component, and vorticity $1$ in the other, while MMs mix terms with zero and
nonzero vorticities in both components. Thus, a relevant issue is to consider
the SOC system with the fractional kinetic-energy operators, and construct
soliton states in them, with emphasis on possibilities to stabilize them
against the collapse. This is the subject of the present work.

The character of the motion of the particles in the condensate (Lévy flights)
does not affect the form of the nonlinearity in the GPE. As concerns the SOC
term in the system’s Hamiltonian, its systematic derivation, starting from the
consideration of transitions between two intrinsic quantum states in the
particle, driven by the Raman laser fields, and applying the unitary
transformation to the spinor wave function Zhai ; Busch , produces a
cumbersome result, in the case of the fractional kinetic energy. In this work,
we adopt a phenomenological model, in which the SOC is represented by the same
term of the Rashba type in the system’s Hamiltonian as in the usual situation
usual . Actually, the results reported below demonstrate that the SOC produces
a minor effect in the 1D setup (therefore, the 1D case is considered below in
a brief form in Section 2), while in the most interesting 2D case, considered
in Section 3, the effect is conspicuous: it shows partial stabilization of the
solitons of the SV and MM types at $1<\alpha\leq 2$ (the range of $\alpha<1$
is considered too, but it produces no soliton states).

It is plausible that the results will be quite similar beyond the framework of
the phenomenological model. Indeed, the stabilization of the 2D solitons by
SOC is accounted for by its spatial scale. It breaks the scaling invariance
which underlies the onset of the collapse Sulem ; Fibich . This fact is
equally true for the full and phenomenological forms of the system which
includes the fractional kinetic energy and SOC interaction. On the other hand,
it is expected that the cumbersome form of the full system will not allow one
to clearly distinguish the SV and MM species of 2D solitons, adding extra
vortical terms to both components of each soliton.

In the analysis following below, stationary shapes of both 1D and 2D solitons
are produced in a numerical form, and, in parallel, analytically by means of
the variational approximation (VA). Stability of the solitons is tested with
the help of direct simulations of their perturbed evolution.

##  II The 1D system

The system of 1D GPEs for the spinor wave function,
$\boldsymbol{\phi}=(\phi_{+},\phi_{-})$, of the binary BEC with attractive
contact interactions and the phenomenologically introduced SOC terms of the
Rashba type usual , with strength $\lambda$, under the action of the
fractional kinetic-energy operator (2) is written, in the scaled form, as

$\displaystyle i\frac{\partial\phi_{+}}{\partial t}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(-\frac{\partial^{2}}{\partial
x^{2}}\right)^{\alpha/2}\phi_{+}-(|\phi_{+}|^{2}+\gamma|\phi_{-}|^{2})\phi_{+}+\lambda\frac{\partial\phi_{-}}{\partial
x},$ $\displaystyle i\frac{\partial\phi_{-}}{\partial t}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(-\frac{\partial^{2}}{\partial
x^{2}}\right)^{\alpha/2}\phi_{-}-(|\phi_{-}|^{2}+\gamma|\phi_{+}|^{2})\phi_{-}-\lambda\frac{\partial\phi_{+}}{\partial
x},$ (5)

where $\gamma\geq 0$ is the relative strength of the inter-component
attraction, while the strength of the self-attraction is scaled to be $1$.
Stationary soliton solutions of Eq. (5) with chemical potential $\mu<0$ are
looked for as

$\phi_{\pm}\left(x,t\right)=e^{-i\mu t}u_{\pm}(x),$ (6)

with real functions $u_{\pm}(x)$ satisfying the following equations:

$\displaystyle\mu u_{+}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(-\frac{\partial^{2}}{\partial
x^{2}}\right)^{\alpha/2}u_{+}-(u_{+}^{2}+\gamma
u_{-}^{2})u_{+}+\lambda\frac{du_{-}}{dx},$ $\displaystyle\mu u_{-}$
$\displaystyle=$ $\displaystyle\frac{1}{2}\left(-\frac{\partial^{2}}{\partial
x^{2}}\right)^{\alpha/2}u_{-}-(u_{-}^{2}+\gamma
u_{+}^{2})u_{-}-\lambda\frac{du_{+}}{dx}.$ (7)

In agreement Eqs. (7), it is possible to set $u_{+}(x)$ and $u_{-}(x)$ to be,
respectively, even and odd functions of $x$.

Note that the system of 1D stationary equations (7) can be derived from the
Lagrangian,

$\displaystyle
L_{\mathrm{1D}}=\frac{\mu}{2}\int_{-\infty}^{+\infty}dx\left[u_{+}^{2}(x)+u_{-}^{2}(x)\right]$
$\displaystyle-\frac{1}{4\pi}\int_{0}^{+\infty}p^{\alpha}dp\int_{-\infty}^{+\infty}d\xi\int_{-\infty}^{+\infty}dx\cos\left(p(x-\xi\right))\cdot\left[u_{+}(x)u_{+}(\xi)+u_{-}(x)u_{-}(\xi)\right]$
$\displaystyle+\int_{-\infty}^{+\infty}\left[\lambda
u_{-}\frac{du_{+}}{dx}+\frac{1}{4}\left(u_{+}^{4}+u_{-}^{4}\right)+\frac{\gamma}{2}u_{+}^{2}u_{-}^{2}\right]dx.$
(8)

Taking into regard the parities of functions $u_{\pm}(x)$, the Lagrangian can
be additionally simplified:

$\displaystyle
L_{\mathrm{1D}}=\frac{\mu}{2}\int_{-\infty}^{+\infty}dx\left[u_{+}^{2}(x)+u_{-}^{2}(x)\right]$
$\displaystyle-\frac{1}{4\pi}\int_{0}^{+\infty}p^{\alpha}dp\int_{-\infty}^{+\infty}d\xi\int_{-\infty}^{+\infty}dx\left\\{\cos(px)\cos(p\xi)\left[u_{+}(x)u_{+}(\xi)\right]+\sin(px)\sin(p\xi)\left[u_{-}(x)u_{-}(\xi)\right]\right\\}$
$\displaystyle+\int_{-\infty}^{+\infty}dx\left[\lambda
u_{-}\frac{du_{+}}{dx}+\frac{1}{4}\left(u_{+}^{4}+u_{-}^{4}\right)+\frac{\gamma}{2}u_{+}^{2}u_{-}^{2}\right].$
(9)

The Lagrangian is used to look for solitons in the framework of the VA, which
may be based on the Gaussian-shaped ansatz,

$u_{+}(x)=\sqrt{\frac{N}{\sqrt{\pi}W}}\left(\cos\theta\right)\exp\left(-\frac{x^{2}}{2W^{2}}\right),u_{-}(x)=\sqrt{\frac{2N}{\sqrt{\pi}W}}\left(\sin\theta\right)\frac{x}{W}\exp\left(-\frac{x^{2}}{2W^{2}}\right).$
(10)

Free parameters $\theta$ and $W$ in Eq. (10) determine the distribution of the
norm between the components and the width of the soliton, for fixed total
norm,

$N=\int_{-\infty}^{+\infty}dx\left(u_{+}^{2}+u_{-}^{2}\right).$ (11)

The substitution of ansatz (10) in Lagrangian (9) yields an expressions for
the effective Lagrangian,

$\displaystyle\left(L_{\mathrm{eff}}\right)_{\mathrm{1D}}$ $\displaystyle=$
$\displaystyle\frac{\mu}{2}N-\Gamma\left(\frac{\alpha+1}{2}\right)\frac{N}{4\sqrt{\pi}W^{\alpha}}\left[\cos^{2}\theta+\left(\alpha+1\right)\sin^{2}\theta\right]-\frac{\lambda
N}{2\sqrt{2}W}\sin\left(2\theta\right)$ (12)
$\displaystyle+\frac{N^{2}}{4\sqrt{2\pi}W}\left(\cos^{4}\theta+\frac{3}{4}\sin^{4}\theta\right)+\frac{\gamma
N^{2}}{16\sqrt{2\pi}W}\sin^{2}(2\theta),$

where $\Gamma$ is the Euler’s Gamma-function. Then, for given $N$ the VA
predicts values $\theta$ and $W$ as solutions of the Euler-Lagrange equations,

$\frac{\partial\left(L_{\mathrm{eff}}\right)_{\mathrm{1D}}}{\partial
W}=\frac{\partial\left(L_{\mathrm{eff}}\right)_{\mathrm{1D}}}{\partial\theta}=0.$
(13)

In the case of $\lambda=0$ in Eq. (5) and, accordingly, $\theta=0$ in Eq.
(10), equation $\partial\left(L_{\mathrm{eff}}\right)_{\mathrm{1D}}/\partial
W=0$ is identical to that derived in Ref. Guangzhou . In addition to Eq. (13),
a relation between $\mu$ and $N$ is determined, in the framework of the VA, by
equation

$\frac{\partial\left(L_{\mathrm{eff}}\right)_{\mathrm{1D}}}{\partial N}=0.$
(14)

Note that the equations which correspond to the critical collapse, i.e., the
1D version of Eq. (4) with $\alpha=1$, or its 2D version with $\alpha=2$,
produce degenerate families of (unstable) soliton solutions, with the single
(critical) value of the norm, $N_{c}$, that does not depend on $\mu$. In
particular, in the 2D case with $\alpha=2$, this is the well-known family of
Townes solitons Townes ; Sulem ; Fibich , whose single (numerically found)
value of the norm is

$\left(N_{c}\right)_{\mathrm{2D}}\equiv\left(N_{\mathrm{Townes}}^{\mathrm{(num)}}\right)_{\mathrm{2D}}\approx
5.85.$ (15)

The simple Gaussian-based VA for the Townes solitons predicts an approximate
value of the norm Anderson ,

$\left(N_{\mathrm{Townes}}^{\mathrm{(VA)}}\right)_{\mathrm{2D}}=2\pi,$ (16)

the relative error, in comparison to Eq. (15) being $\simeq 7\%$. In the case
of the 1D version of Eq. (4) with $\alpha=1$, similar results for the
degenerate soliton family are Zeng 2

$~{}\left(N_{c}\right)_{\mathrm{1D}}\equiv
N_{\mathrm{1D,}\alpha=1}^{\mathrm{(num)}}\approx
1.23,~{}N_{\mathrm{1D,}\alpha=1}^{\mathrm{(VA)}}=\sqrt{2},$ (17)

the relative error being $\simeq 13\%$, cf. Eqs. (15) and (16).

If the SOC terms are included ($\lambda>0$), width $W$ drops out from Eqs.
(13) and (12) at $\alpha=1$, hence this system too admits a solution at a
single (critical) value $N_{c}$ of $N$, i.e., the soliton family remains a
degenerate one in this case. This conclusion, which is a natural consequence
of the fact that operator (2) with $\alpha=1$ and SOC terms in Eq. (5) feature
the same scaling with respect to $x$, is confirmed by a numerically found
solution of Eq. (7), which was obtained by applying the imaginary-time
integration method to the full system of equations (5).

Dependences $N_{c}(\lambda)$ for the 1D system with $\alpha=1$, as produced by
the numerical solution and VA, are displayed in Fig. 1(a). In particular, a
straightforward consequence of Eqs. (13) and (12) is the value of the SOC
strength, $\lambda$, at which the VA-predicted critical norm vanishes:

$N_{c}^{\mathrm{(var)}}\left(\lambda=\left(\lambda_{\max}^{\mathrm{(var)}}\right)_{\alpha=1}\equiv
1/\sqrt{\pi}\approx 0.56\right)=0,$ (18)

with no solitons existing at $\lambda>\lambda_{\max}$. The exact value of
$\lambda_{\max}$ can be found too, without the use of the VA. Indeed, it
follows from the consideration of the scaling of soliton’s tails decaying at
$|x|\rightarrow\infty$, which are produced by Eq. (7) with $\alpha=1$, that
exponentially decaying tails are possible at

$\lambda<\left(\lambda_{\max}^{\mathrm{(exact)}}\right)_{\alpha=1}=1/2.$ (19)

At $\lambda>1/2$, there exist solutions for tails decaying algebraically,
rather than exponentially:

$u_{+,-}(x)\sim|x|^{-1/\left(4\lambda^{2}-1\right)}.$ (20)

The results demonstrate that slowly decaying tails (20) cannot be used for
constructing solitons (in addition, it is relevant to mention that the total
norm, determined by the asymptotic form (20) at $|x|\rightarrow\infty$,
diverges at $\lambda\geq\sqrt{3}/2$).

Figure 1: (a) The chain of rhombuses shows the numerically found critical norm
$N_{c}$ of the degenerate 1D soliton family as a function of the SOC strength,
$\lambda$, at $\alpha=1$. The numerical results are produced by imaginary-time
simulations of Eq. (5). The dashed line is the same value as predicted by the
VA, as per Eqs. (13) and (12). The numerical and variational values of $N_{c}$
vanish precisely at points (19) and (18), respectively. (b) Chains of
rhombuses and crosses show the amplitude of component $|\phi_{+}(x)|$ of the
1D solitons for two fixed values of the norm, $N=1.3$ and $0.8$, respectively,
as functions of $\alpha/2$ at fixed $\lambda=0.3$. Dashed lines represent the
respective VA-predicted results. (c) An example of a stable 1D soliton, found
in the numerical form at $\alpha=1.04$ and $\lambda=1/2$, with total norm
$N=1$, is displayed by means of the double-logarithmic plot of
$|\phi_{+}(x)|$. The dashed-line fit is $0.001/x^{2}$ (this fit is not valid
at still larger values of $|x|$, where the soliton’s tail decays
exponentially). In this figure and Fig. 2, the results are presented for
$\gamma=0$ in Eq. (5) (no attraction between the two components).

The norm-degenerate family of the 1D solitons existing at $\alpha=1$ and
$\lambda<1/2$ is completely unstable. Its relation to nondegenerate families
of stable solitons existing at $\alpha>1$ is illustrated by Fig. 1(b), which
shows the amplitude of the solitons as a function of $\alpha$ for a fixed SOC
strength, $\lambda=0.3$, which is smaller than the one given by Eq. (19) (at
$\lambda>1/2$, the solitons definitely do not exist in the limit of
$\alpha=1$, as shown above), and for two different fixed norms, $N=0.8$ and
$1.3$. Because both are different from the respective value
$N_{c}^{\mathrm{(num)}}(\lambda=0.3)\approx 0.9$, see Fig. 1(a), the 1D
solitons with these fixed values of $N$ do not exist in the limit of
$\alpha=1$. Accordingly, Fig. 1(b) demonstrates that the solitons with
$N<N_{c}^{\mathrm{(num)}}(\lambda=0.3)$ suffer delocalization at
$\alpha\rightarrow 1$, with the amplitude approaching zero in this limit,
while the solitons with $N>N_{c}^{\mathrm{(num)}}(\lambda=0.3)$ suffer
collapse (divergence of the amplitude).

Figure 1(c) shows that, at values of $\alpha$ slightly exceeding $1$ (here,
$\alpha=1.04$) and $\lambda=1/2$, the 1D soliton exists with tails close to
the algebraic shape (comparison to that expression is relevant as $\alpha$ is
close to $1$, although Eq. (20) cannot be directly applied because $\alpha$ is
not 1). Lastly, Fig. 2 illustrates the stability of a typical 1D soliton with
$N=1.3$ at $\alpha=1.4$ and $\lambda=0.3$.

Figure 2: The evolution of $|\phi_{+}(x)|$ and $|\phi_{-}(x)|$ (solid and
dashed lines, respectively), produced by simulations of Eq. (5) with
$\alpha=1.4$ and $\lambda=0.3$. The norm of this 1D soliton is $N=1.3$.

Thus, the effect of the SOC terms in the 1D system changes the critical norm
for the collapse but fails to stabilize the solitons in the critical case of
$\alpha=1$. At $\alpha>1$, the effect is not strong either. In particular,
additional calculations clearly demonstrate that variation of $\alpha$ does
not strongly affect stability of the 1D solitons (details are not displayed
here, as they do not convey something remarkable). It is shown below that the
SOC terms produce a much more conspicuous impact in the 2D system.

##  III The 2D system

Following Ref. we , the 2D system of GPEs with the SOC coupling of the Rashba
type and the fractional kinetic-energy operator (3) is written, in the scaled
form, as

$\displaystyle i\frac{\partial\phi_{+}}{\partial t}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(-\nabla^{2}\right)^{\alpha/2}\phi_{+}-(|\phi_{+}|^{2}+\gamma|\phi_{-}|^{2})\phi_{+}+\lambda\left(\frac{\partial\phi_{-}}{\partial
x}-i\frac{\partial\phi_{-}}{\partial y}\right),$ $\displaystyle
i\frac{\partial\phi_{-}}{\partial t}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(-\nabla^{2}\right)^{\alpha/2}\phi_{-}-(|\phi_{-}|^{2}+\gamma|\phi_{+}|^{2})\phi_{-}-\lambda\left(\frac{\partial\phi_{+}}{\partial
x}+i\frac{\partial\phi_{+}}{\partial y}\right),$ (21)

cf. Eq. (5). Note that, except for case of $\alpha=1$, rescaling

$\left(x,y\right)=\lambda^{-1/(\alpha-1)}\left(\tilde{x},\tilde{y}\right),t=\lambda^{-\alpha/(\alpha-1)}\tilde{t},\phi_{\pm}=\lambda^{\alpha/\left(2(\alpha-1)\right)}\tilde{\phi}_{\pm}$
(22)

makes it possible to set $\lambda\equiv 1$ in Eq. (21). Accordingly, the norm
of the 2D states scales with the variation of the SOC strength as

$N(\lambda)\equiv\int\int\left(\left|u_{+}\right|^{2}+\left|u_{-}\right|^{2}\right)dxdy=\lambda^{-(2-\alpha)/(\alpha-1)}N(\lambda=1),$
(23)

where $N(\lambda=1)$ is the norm at the reference point, $\lambda=1$.

Solutions of Eqs. (21) for 2D stationary states with chemical potential $\mu$
are looked for as

$\phi_{\pm}\left(x,y,t\right)=e^{-i\mu t}u_{\pm}(x,y),$ (24)

with complex functions $u_{\pm}(x,y)$, cf. the 1D counterpart given by Eq. (6)
with real functions $u_{\pm}(x)$. The substitution of ansatz (24) in Eqs. (21)
leads to the stationary equations for complex functions $u_{\pm}(x,y)$:

$\displaystyle\mu u_{+}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(-\nabla^{2}\right)^{\alpha/2}u_{+}-(|u_{+}|^{2}+\gamma|u_{-}|^{2})u_{+}+\lambda\left(\frac{\partial
u_{-}}{\partial x}-i\frac{\partial u_{-}}{\partial y}\right),$
$\displaystyle\mu u_{-}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(-\nabla^{2}\right)^{\alpha/2}u_{-}-(|u_{-}|^{2}+\gamma|u_{+}|^{2})u_{-}-\lambda\left(\frac{\partial
u_{+}}{\partial x}+i\frac{\partial u_{+}}{\partial y}\right),$ (25)

cf. the 1D stationary equations (7). Taking into regard definition (3), Eqs.
(25) can be derived from the respective Lagrangian:

$\displaystyle
L=\mu\int_{-\infty}^{+\infty}dx\int_{-\infty}^{+\infty}dy\left[\left|u_{+}(x,y)\right|^{2}+\left|u_{-}(x,y)\right|^{2}\right]$
$\displaystyle-\frac{1}{2\pi^{2}}\int_{0}^{+\infty}dp\int_{0}^{+\infty}dq(p^{2}+q^{2})^{\alpha/2}\int_{-\infty}^{+\infty}d\xi\int_{-\infty}^{+\infty}d\eta\int_{-\infty}^{+\infty}dx\int_{-\infty}^{+\infty}dy\cos\left[p(x-\xi)+q(y-\eta)\right]$
$\displaystyle\times\left[u_{+}^{\ast}(x,y)u_{+}(\xi,\eta)+u_{-}^{\ast}(x,y)u_{-}(\xi,\eta)\right]$
$\displaystyle-\lambda\int_{-\infty}^{+\infty}dx\int_{-\infty}^{+\infty}dy\left[u_{+}^{\ast}\frac{\partial
u_{-}}{\partial x}+u_{+}\frac{\partial u_{-}^{\ast}}{\partial
x}-i\left(u_{+}^{\ast}\frac{\partial u_{-}}{\partial y}-u_{+}\frac{\partial
u_{-}^{\ast}}{\partial y}\right)\right]$
$\displaystyle+\int_{-\infty}^{+\infty}dx\int_{-\infty}^{+\infty}dy\left[\frac{1}{2}\left(|u_{+}|^{4}+|u_{-}|^{4}\right)+\gamma|u_{+}|^{2}|u_{-}|^{2}\right],$
(26)

cf. Eq. (8). The integration with respect to $p$ and $q$ in Eq. (26) is
reduced from original domains $\left(-\infty,+\infty\right)$ to
$\left(0,\infty\right)$, using the symmetry of the domain.

As suggested by Ref. we  (which addressed the system with the normal kinetic-
energy operator, i.e., $\alpha=2$), we first aim to construct 2D soliton
solutions of the SV type. To this end, the following ansatz is adopted:

$u_{+}=A_{+}\exp(-\beta(x^{2}+y^{2})),u_{-}=A_{-}(x+iy)\exp(-\beta(x^{2}+y^{2})),$
(27)

which implies vorticities $0$ and $1$ in the components $u_{+}$ and $u_{-}$,
respectively, in accordance with the definition of the SV mode. Note that the
maximum of component $u_{-}$ of the ansatz,
$\left|u_{-}(r)\right|_{\max}=e^{-1/2}A_{-}/\sqrt{2\beta}$, is attained at
$r=1/\sqrt{2\beta}$.

The substitution of ansatz (27) in Lagrangian (26) produces the respective
effective Lagrangian, cf. Eq. (12):

$L_{\mathrm{eff}}=\mu
N-\frac{\pi\Gamma\left(1+\alpha/2\right)}{2\left(2\beta\right)^{1-\alpha/2}}\left[A_{+}^{2}+\left(1+\frac{\alpha}{2}\right)\frac{A_{-}^{2}}{2\beta}\right]-\frac{\pi\lambda}{\beta}A_{+}A_{-}+\frac{\pi}{8\beta}\left(A_{+}^{4}+\frac{A_{-}^{4}}{8\beta^{2}}+\frac{\gamma
A_{+}^{2}A_{-}^{2}}{2\beta}\right),$ (28)

where the total 2D norm of ansatz (27) is

$N=\frac{\pi}{2\beta}\left(A_{+}^{2}+\frac{A_{-}^{2}}{2\beta}\right).$ (29)

Then, for given $\mu$, SV’s parameters are predicted by the Euler-Lagrange
equations (cf. Eq. (13)),

$\frac{\partial L_{\mathrm{eff}}}{\partial A_{\pm}}=\frac{\partial
L_{\mathrm{eff}}}{\partial\beta}=0.$ (30)

Figure 3 presents the most essential results of the analysis for the 2D case
with $\gamma=0$ (no nonlinear interaction between components $\phi_{\pm}$ in
Eq. (21)). Namely, unlike the 2D version of Eq. (4) with $\alpha\leq 2$, which
does not include SOC and generates only completely unstable soliton solutions,
Eq. (21) gives rise to _stable SVs_ at $1<\alpha\leq 2$, with norms falling
below the respective critical value, viz.,
$N\,<N_{c}^{\mathrm{(SV)}}(\alpha)$. The dependence of $N_{c}^{\mathrm{(SV)}}$
on $\alpha$, as obtained from the numerical solution (using simulations of Eq.
(21) in imaginary time) and VA, is displayed in Fig. 3(a) for $\lambda=0.5$.
As it follows from the above argument which demonstrates the peculiarity of
the case of $\alpha=1$, $N_{c}^{\mathrm{(SV)}}$ vanishes at $\alpha=1$, which
is indeed demonstrated by the numerical results in Fig. 3(a) (while the VA
result is inaccurate at $\alpha\rightarrow 1$, as the ansatz (27) is not
relevant in this limit). Thus, SV states with a finite amplitude do not exist
at $\alpha=1$. On the other hand, for the same case of $\alpha=1$ the solution
with the vanishing amplitude can be obtained from the linearized version of
Eq. (21). This solution, generated by the imaginary-time simulations of the
linearized equations, with periodic boundary conditions fixed for domain
$|x|,|y|\leq 10$, is displayed by means of cross sections
$\left|\phi_{\pm}\left(x\right)\right|$, drawn through $y=0$, in Fig. 3(b).
The effective localization of the solution implies that it is, roughly
speaking, similar to periodic states represented by the Jacobi’s elliptic
functions with modulus close to $k=1$.

Figure 3: (a) The critical value of the norm of the SV solitons in the 2D
system, $N_{c}^{\mathrm{(SV)}}$, as a function of the half Lévy index,
$\alpha/2$, at $\lambda=0.5$ for $\gamma=0$ in Eq. (21). At $1<\alpha\leq 2$,
the system gives rise to stable SV (semi-vortex) solitons at
$N<N_{c}^{\mathrm{(SV)}}$, and to collapse at $N>N_{c}^{\mathrm{(SV)}}$. Here
and in panel (c), the chain of rhombuses and dashed line represent,
respectively, the numerical findings, obtained from the imaginary-time
simulations of Eq. (21), and VA results produced by Eqs. (28) and (30). (b)
The cross-section profiles of $|\phi_{+}|$ and $\left|\phi_{-}\right|$ (the
solid and dashed lines, respectively), produced at $\lambda=0.5$ and
$\alpha=1$ by the numerical solution of the _linearized system_ of size
$20\times 20$. Amplitudes of the linear solution are determined by an
(arbitrarily adopted) normalization condition $N=1$. (c) The amplitude of the
SV solitons (largest value of $|\phi_{+}(x)|$ vs. the SOC strength $\lambda$,
at $\alpha=1.5$, $\gamma=0$ and a fixed total norm, $N=0.5$. At small
$\lambda$, the dependence is explained by the scaling relation (31) with
$\alpha=1.5$, i.e., $A_{0}\sim\sqrt{N}\lambda^{2}$.

The build of the family of stable SV solitons, constructed by means of the
numerical solution and VA, is illustrated by Fig. 3(c) for a generic value of
the Lévy index, $\alpha=1.5$. In this figure, the SOC strength $\lambda$ is
not fixed but varied, to explicitly display the effect of the SOC, while the
scaling is used to keep a fixed value of the norm, $N=0.5$ (the mutual scaling
of $N$ and $\lambda$ is determined by Eq. (23). It is seen that the SV’s
amplitude, $A_{0}(\lambda)$, vanishes at $\lambda\rightarrow 0$ (in agreement
with the fact there are no stable 2D solitons in the absence of the SOC in the
system). Straightforward analysis of Eqs. (25) and its VA counterpart (30)
demonstrates the following scaling relation between the SV’s amplitude,
$A_{0}$, its width, $\beta^{-1/2}$ (see Eq. (27)), and $\lambda$, at
$\lambda\rightarrow 0$ (cf. Eq. (23):

$A_{0}\sim\sqrt{N}\lambda^{1/(\alpha-1)},\beta^{-1/2}\sim\lambda^{-1/(\alpha-1)}.$
(31)

For $\alpha=1.5$, Eq. (31) yields $A_{0}\sim\sqrt{N}\lambda^{2}$, which
obviously agrees with the curves plotted in Fig. 3(c).

Further, Fig. 3(c) shows that the amplitude increases with the increase of the
SOC strength, up to a maximum value, $\lambda_{\max}\approx 2.4$, above which
the SV soliton is destroyed by the collapse. While the VA produces, in the
same figure, the amplitude-vs.-$\lambda$ curve which is very close to the
numerical counterpart, it does not terminate at $\lambda$ close to
$\lambda_{\max}$. Nevertheless, the termination of the soliton family at large
$\lambda$ is a natural finding, as the existence of solitons in the system
with the SOC terms dominating over the kinetic energy requires the presence of
the Zeeman splitting between the components we2 , which is absent in Eqs. (5)
and (21).

It is plausible that the stable SV solitons, populating the area beneath the
right branch of the $N_{c}^{\mathrm{(SV)}}(\alpha)$ curve in Fig. 3(a), play
the role of the system’s ground state, which is missing (formally replaced by
the collapsing solution) in the absence of the SOC terms, cf. a similar
conclusion made in Ref. we  for the usual 2D system with $\alpha=2$. This
point calls for a deeper consideration, which is beyond the scope of the
present paper.

At values of the norm $N>N_{c}^{\mathrm{(SV)}}$ simulations of Eqs. (21)
demonstrate the collapse of the SV, no stationary solitons being possible. At
$\alpha<1$, the system does not give rise to solitons either.

The evolution and stability of typical SV solitons in real-time simulations of
Eq. (21) with $\gamma=0$ are illustrated in Fig. 4 for (a) $\lambda=1$,
$\alpha=1.5$, and $N=1$ and (b) $\lambda=0.4$, $\alpha=1.9$, and $N=5.15$. The
latter case represents the 2D system with the Lévy index close to $\alpha=2,$
which corresponds to the system with the usual (non-fractional) kinetic-energy
operator, that was considered in Ref. we .

Figure 4: Stable evolution of $|\phi_{+}(x,y)|$ and $|\phi_{-}(x,y)|$ in the
SV soliton, plotted in the cross section $y=0$ by solid and dashed lines,
respectively, (a) An example for a typical fractional system, with
$\lambda=1$, $\alpha=1.5$, and $\gamma=0$ in Eq. (21), and the soliton’s norm
$N=1$. (b) The stable evolution in a system close to the regular (non-
fractional) one, with $\lambda=0.4$, $\alpha=1.9$, $\gamma=0$, and $N=5.15$.

The above consideration was performed for $\gamma=0$ in Eq. (21), when (as
well as for $\gamma<1$, i.e., for the self-attraction of each component
stronger than the cross-attraction) the SV is expected to be the dominant
species of 2D solitons we . On the other hand, in the usual (non-fractional)
system, SVs are unstable at $\gamma>1$, while stable 2D solitons are MM
solution, which, as mentioned above, mix the zero-vorticity and vortex terms
in each component. The MM represents the ground state of the non-fractional
system with $\gamma>1$ we .

In the present case, the MMs may be approximated by the following variational
ansatz for the stationary components (see Eq. (24)),

$\displaystyle u_{+}$ $\displaystyle=$ $\displaystyle
A_{1}\exp(-\beta(x^{2}+y^{2}))-A_{2}(x-iy)\exp(-\beta(x^{2}+y^{2})),$
$\displaystyle u_{-}$ $\displaystyle=$ $\displaystyle
A_{1}\exp(-\beta(x^{2}+y^{2}))+A_{2}(x+iy)\exp(-\beta(x^{2}+y^{2})),$ (32)

which indeed mixes terms with zero and nonzero vorticities. The substitution
of ansatz (32) in Lagrangian (26) yields

$L_{\mathrm{eff}}=\mu
N-\frac{\pi\Gamma\left(2+\alpha/2\right)}{\left(2\beta\right)^{2-\alpha/2}}\left(2\beta
A_{1}^{2}+A_{2}^{2}\right)+(1+\gamma)\left(\frac{\pi
A_{1}^{4}}{4\beta}+\frac{\pi A_{2}^{4}}{32\beta^{3}}\right)+\frac{\pi
A_{1}^{2}A_{2}^{2}}{4\beta^{2}}-\frac{2\lambda\pi A_{1}A_{2}}{\beta},$ (33)

cf. Eq. (28), where the total norm is
$N=\pi[A_{1}^{2}/(\beta)+A_{2}^{2}/(2\beta^{2})]$, cf. Eq. (29). The
corresponding variational equations are then derived as per Eq. (30).

Figure 5: The critical value of the norm of the MM (mixed-mode) solitons in
the 2D system, $N_{c}^{\mathrm{(MM)}}$, as a function of the half Lévy index,
$\alpha/2$, for fixed $\lambda=0.5$ and $\gamma=2$ in Eq. (21). At
$1<\alpha\leq 2$, the system gives rise to stable MM solitons at
$N<N_{c}^{\mathrm{(MM)}}$, and to collapse at $N>N_{c}^{\mathrm{(MM)}}$. The
chain of rhombuses and dashed line represent, respectively, the numerical
findings, obtained from the imaginary-time simulations of Eq. (21), and VA
results produced by Eq. (33).

Similar to the SV family, stable 2D solitons of the MM type exist at
$\alpha>1$, below the respective critical value of the norm,
$N<N_{c}^{\mathrm{(MM)}}$. The numerically found and VA-predicted dependences
$N_{c}^{\mathrm{(MM)}}(\alpha)$ are shown in Fig. 5, with $\lambda=0.5$ and
$\gamma=2$ fixed in Eq. (21). Note that values of $N_{c}^{\mathrm{(MM)}}$ in
Fig. 5 are close to $2/3$ of the respective values of $N_{c}^{\mathrm{(SV)}}$,
which are presented, for the same values of $\lambda$ and $\alpha$, in Fig.
3(a). This relation is explained by the same argument which was presented, for
$\alpha=2$, in Ref. we : taking into regard that the dominant contribution to
the total norm of the SV soliton is produced by the zero-vorticity component,
and that the nonlinear interaction between the MM components may be, roughly
speaking, absorbed by means of rescaling,
$u_{\pm}\rightarrow\left(1+\gamma\right)^{-1/2}u_{\pm}$, an approximate
relation between the SV and MM norms is $N^{\mathrm{(MM)}}\approx
2\left(1+\gamma\right)^{-1}N^{\mathrm{(SV)}}$. For $\gamma=2$, it amounts to
the above-mentioned norm ratio, $N^{\mathrm{(MM)}}/N^{\mathrm{(SV)}}\approx
2/3$. A similar result, $N^{\mathrm{(MM)}}/N^{\mathrm{(SV)}}\approx
2/(1+\gamma)$, can be obtained for other values of $\gamma>1$, which maintain
the MMs as stable states.

If $N<N_{c}^{\mathrm{(MM)}}$ is fixed, while the SOC strength $\lambda$ is
varying, the dependence of the MM’s amplitude on $\lambda$ is similar to that
for the SVs plotted in Fig. 3(b), i.e., the amplitude vanishes in the limit of
$\lambda\rightarrow 0$ according to the same relation (31) which is produced
above for the SVs. Figure 6 displays the evolution of components
$|\phi_{\pm}(x,y)|$ at $\lambda=1$, $\alpha=1.5$, $\gamma=2$ and $N=0.8$, as
produced by real-time simulations of Eq. (21) (note that these parameters
correspond to a state located beneath the stability boundary in Fig. 5). The
plots clearly demonstrate the structure and stability of the generic MM.

Figure 6: The stable evolution of components $|\phi_{+}(x,y)|$ and
$|\phi_{-}(x,y)|$ of a typical MM soliton, shown by solid and dashed lines,
respectively, in the cross section of $y=0$. The parameters are $\lambda=1$,
$\alpha=1.5$, $\gamma=2$. The total norm of the soliton is $N=0.8$. The input
was produced by simulations of Eq. (21) in imaginary time.

At $N>N_{c}^{\mathrm{(MM)}}$, the MM-shaped input leads to the collapse. In
the range of $\alpha<1$, where SV solitons do not exist, MM solitons were not
found either.

##  IV Conclusion

This works aims to consider effects of SOC (spin-orbit coupling) in 1D and 2D
binary matter waves with the fractional kinetic-energy operator and the usual
cubic self- and cross-attractive nonlinearity. This is a model of BEC composed
of particles moving by Lévy flights, characterized by the value of the Lévy
index, $\alpha<2$ (the normal, non-fractional, kinetic energy corresponds to
$\alpha=2$). In the 1D setting, the effect of the SOC is not dramatic, leading
to decrease of the norm at which the collapse takes place in the system with
$\alpha=1$. Essential effects are predicted in the 2D setting. In that case,
the SOC creates regions of stable solitons of the SV (semi-vortex) and MM
(mixed-mode) solitons in the interval of $1<\alpha<2$, where the supercritical
collapse occurs and no stable modes exist in the absence of SOC. The stable
solitons exist at values of the norm below the respective critical values,
$N<N_{c}^{\mathrm{(SV,MM)}}$. Amplitudes of the stable solitons in these
regions vanish along with SOC strength, $\lambda$, as per Eq. (31).

As an extension of the present analysis, it may be relevant to consider moving
solitons and collisions between them. The 2D system may also be used to
predict other patterns, such as vortex lattices, cf. Refs. vort-latt -Fukuoka
.

## Acknowledgment

The work of B.A.M. was supported, in part, by the Israel Science Foundation
through grant No. 1286/17.

## References

  * (1) Podlubny I, Magin R L and Trymorush I 2017 Niels Henrik Abel and the birth of fractional calculus Fractional Calculus Appl. Anal. 20 1068-1075
  * (2) Laskin N 2000 Fractional quantum mechanics and Lévy path integrals Phys. Lett. A 268 298-305
  * (3) Dubkov A A, Spagnolo B and Uchaikin V V 2008 Lévy flight superdiffusion: an introduction Int. J. Bifurcat. Chaos 18 2649\.
  * (4) Laskin N 2018 Fractional Quantum Mechanics (World Scientific, Singapore)
  * (5) Riesz M 1949 L’intégrale de Riemann-Liouville et le problème de Cauchy Acta Mathematica 81 1–223, doi:10.1007/BF02395016
  * (6) Cai M Li C P 2019 On Riesz derivative Fractional Calculus Appl. Anal. 22 287–301
  * (7) Longhi S 2015 Fractional Schrödinger equation in optics Opt. Lett. 40 1117-1120
  * (8) Monmayrant A, Weber S and Chatel B 2010 A newcomer’s guide to ultrashort pulse shaping and characterization J. Phys. B: At. Mol. Opt. Phys. 43 103001\.
  * (9) Kasprzak H 1982 Differentiation of a noninteger order and its optical implementation Appl. Opt. 21 3287-3291
  * (10) Davis J A, Smith D A, McNamara D E, Cottrell D M and Campos J 2001 Fractional derivatives - analysis and experimental implementation Appl. Opt. 40 5943-5948
  * (11) Stickler B A 2013 Potential condensed-matter realization of space-fractional quantum mechanics: The one-dimensional Lévy crystal Phys. Rev. E 88 012120
  * (12) Pinsker F, Bao W, Zhang Y, Ohadi H, Dreismann A, and Baumberg J 2015 Fractional quantum mechanics in polariton condensates with velocity-dependent mass Phys. Rev. B 92 195310
  * (13) Zhong W P, Belić M R and Zhang Y 2016 Accessible solitons of fractional dimension Ann. Phys. 368 110–116
  * (14) Huang C and Dong L 2016 Gap solitons in the nonlinear fractional Schrödinger equation with an optical lattice Opt. Lett. 41 5636–5639
  * (15) Zhang L, He Z, Conti C, Wang Z, Hu Y, Lei D, Li Y and Fan A D 2017 Modulational instability in fractional nonlinear Schrödinger equation Commun. Nonlin. Sci. Numer. Simul. 48, 531–540
  * (16) Xiao J, Tian Z, Huang C and Dong L 2018 Surface gap solitons in a nonlinear fractional Schrödinger equation Opt. Exp. 26 2650–2658
  * (17) Zeng L and Zeng J 2019 One-dimensional gap solitons in quintic and cubic-quintic fractional nonlinear Schrödinger equations with a periodically modulated linear potential Nonlinear Dyn. 98, 985–995
  * (18) Molina M I 2020 The fractional discrete nonlinear Schrödinger equation Phys. Lett. A 384, 126180
  * (19) Qiu Y, Malomed B A, Mihalache D, Zhu X, Zhang L and He Y 2020 Soliton dynamics in a fractional complex Ginzburg-Landau model Chaos, Solitons and Fractals 131 109471
  * (20) Li P, Malomed B A and Mihalache D 2020 Vortex solitons in fractional nonlinear Schrödinger equation with the cubic-quintic nonlinearity Chaos Solitons & Fractals 137 109783
  * (21) Wang Q and Liang G 2020 Vortex and cluster solitons in nonlocal nonlinear fractional Schrödinger equation J. Optics 22 055501
  * (22) Li P, Malomed B A, and Mihalache D 2020 Metastable soliton necklaces supported by fractional diffraction and competing nonlinearities Opt. Exp. 28 34472-34488
  * (23) Qiu Y, Malomed B A, Mihalache D, Zhu X, Peng X and He Y 2020 Stabilization of single- and multi-peak solitons in the fractional nonlinear Schrödinger equation with a trapping potential Chaos, Solitons & Fractals 140 110222
  * (24) Li P, Malomed B A and Mihalache D 2020 Symmetry breaking of spatial Kerr solitons in fractional dimension Chaos Solitons & Fractals 132 109602
  * (25) Li P and Dai C 2020 Double loops and pitchfork symmetry breaking bifurcations of optical solitons in nonlinear fractional Schrödinger equation with competing cubic-quintic nonlinearities Ann. Phys. (Berlin) 532 2000048
  * (26) Li P, Li R and Dai C 2021 Existence, symmetry breaking bifurcation and stability of two-dimensional optical solitons supported by fractional diffraction Opt. Exp. 29, 3193–3210
  * (27) Zeng L, Mihalache D, Malomed B A, Lu X, Cai Y, Zhu Q, and Li J 2021 Families of fundamental and multipole solitons in a cubic-quintic nonlinear lattice in fractional dimension Chaos Sol. Fract. 144 110589
  * (28) Zeng L, Zhu Y, Malomed B A, Mihalache D, Wang Q, Long H, Cai Y, Lu X and Li J 2022 Quadratic fractional solitons Chaos, Solitons & Fractals 154 111586
  * (29) Malomed B A 2021 Optical solitons and vortices in fractional media: A mini-review of recent results Photonics 8 353
  * (30) Chen M N, Zeng S H, Lu D Q, Hu W and Guo Q 2018, Optical solitons, self-focusing, and wave collapse in a space-fractional Schrödinger equation with a Kerr-type nonlinearity Phys. Rev. E 98 022211
  * (31) Sulem C and Sulem P L 1999 The Nonlinear Schrödinger Equation: Self-focusing and Wave Collapse (Springer, Berlin).
  * (32) Fibich G 2015 The Nonlinear Schrödinger Equation: Singular Solutions and Optical Collapse (Springer, Heidelberg)
  * (33) Sakaguchi H, Li B, and Malomed B A 2014 Creation of two-dimensional composite solitons in spin-orbit-coupled self-attractive Bose-Einstein condensates in free space Phys. Rev. E 89 032920
  * (34) Sakaguchi H, Sherman E Ya, and Malomed B A 2016 Vortex solitons in two-dimensional spin-orbit coupled Bose-Einstein condensates: Effects of the Rashba-Dresselhaus coupling and the Zeeman splitting Phys. Rev. E 94 032202
  * (35) Malomed B A 2018 Creating solitons by means of spin-orbit coupling EPL 122 36001
  * (36) Zhai H 2015 Degenerate quantum gases with spin–orbit coupling: a review Rep. Prog. Phys. 78 026001
  * (37) Zhang Y, Mossman M E, Busch T, Engels P and Zhang C 2016 Properties of spin–orbit-coupled Bose–Einstein condensates Front. Phys. 11 118103
  * (38) Zhang Y, Mao L, and Zhang C 2012 Mean-field dynamics of spin-orbit coupled Bose-Einstein condensates Phys. Rev. Lett. 108 035302
  * (39) Chiao R Y, Garmire E and Townes C H 1964 Self-trapping of optical beams Phys. Rev. Lett. 13 479–482
  * (40) Desaix M, Anderson D and Lisak M 1991 Variational approach to collapse of optical pulses J. Opt. Soc. Am. B 8 2082-2086
  * (41) Sakaguchi H and Malomed B A 2016 One- and two-dimensional solitons in $\mathcal{PT}$-symmetric systems emulating spin–orbit coupling New J. Phys. 18, 105005
  * (42) Zhou X F, Zhou J, Wu, and Wu C J 2011 Vortex structures of rotating spin-orbit-coupled Bose-Einstein condensates Phys. Rev. A 84, 063624
  * (43) Xu Z F, Kawaguchi Y, You L, and Ueda M 2012 Symmetry classification of spin-orbit-coupled spinor Bose-Einstein condensates Phys. Rev. A 86, 033628
  * (44) Ruokokoski E, Huhtamaki J A M, and Mottonen M 2012 Stationary states of trapped spin-orbit-coupled Bose-Einstein condensates Phys. Rev. A 86, 051607
  * (45) H. Sakaguchi and B. Li 2013 Vortex lattice solutions to the Gross-Pitaevskii equation with spin-orbit coupling in optical lattices Phys. Rev. A 87, 015602


# Self-supervised Point Cloud Prediction
Using 3D Spatio-temporal Convolutional Networks

Benedikt Mersch Xieyuanli Chen Jens Behley Cyrill Stachniss
University of Bonn
{mersch, xieyuanli.chen, jens.behley<EMAIL_ADDRESS>

###### Abstract

Exploiting past 3D LiDAR scans to predict future point clouds is a promising
method for autonomous mobile systems to realize foresighted state estimation,
collision avoidance, and planning. In this paper, we address the problem of
predicting future 3D LiDAR point clouds given a sequence of past LiDAR scans.
Estimating the future scene on the sensor level does not require any preceding
steps as in localization or tracking systems and can be trained self-
supervised. We propose an end-to-end approach that exploits a 2D range image
representation of each 3D LiDAR scan and concatenates a sequence of range
images to obtain a 3D tensor. Based on such tensors, we develop an encoder-
decoder architecture using 3D convolutions to jointly aggregate spatial and
temporal information of the scene and to predict the future 3D point clouds.
We evaluate our method on multiple datasets and the experimental results
suggest that our method outperforms existing point cloud prediction
architectures and generalizes well to new, unseen environments without
additional fine-tuning. Our method operates online and is faster than the
common LiDAR frame rate of $10$ Hz.

###  1 Introduction

Sequential 3D point clouds obtained from Light Detection And Ranging (LiDAR)
sensors are widely used for numerous tasks in robotics and autonomous driving,
like localization [1], detection [2, 3], mapping [4], segmentation [5, 6], and
trajectory prediction [7]. The ability to forecast what the sensor is likely
to see in the future can enhance decision-making for an autonomous vehicle. A
promising application is to use the predicted point clouds for path planning
tasks like collision avoidance. In contrast to most approaches, which predict
for example future 3D bounding boxes of traffic agents, point cloud prediction
does not need any preceding inference steps such as localization, detection,
or tracking to predict a future scene. Running an off-the-shelf detection and
tracking system on the predicted point clouds yields future 3D object bounding
boxes as demonstrated by recent approaches for point cloud forecasting [8, 9].
From a machine learning perspective, point cloud prediction is an interesting
problem since the ground truth data is always given by the next incoming LiDAR
scans. Thus, we can train our point cloud prediction self-supervised without
the need for expensive labeling and evaluate its performance online in unknown
environments.

In this paper, we address the problem of predicting large and unordered future
point clouds from a given sequence of past scans as shown in Fig. 1. High
dimensional and sparse 3D point cloud data render point cloud prediction a
challenging problem that is not yet fully explored in the literature. A future
point cloud can be estimated by applying a predicted future scene flow to the
last received scan or by generating a new set of future points. In this paper,
we focus on the generation of new point clouds to predict the future scene. In
contrast to existing approaches [8, 9], which exploit recurrent neural
networks for modeling temporal correspondences, we use 3D convolutions to
jointly encode the spatial and temporal information. Our proposed approach
takes a new 3D representation based on concatenated range images as input and
jointly estimates a future range image and per-point scores for being a valid
or an invalid point for multiple future time steps. This allows us to predict
detailed future point clouds of varying sizes with a reduced number of
parameters to optimize resulting in faster training and inference times.
Furthermore, our approach is also fully self-supervised and does not require
any manual labeling of the data.

Pre-
---
diction
T-P+1
---
T
---
T+1
---
T+F
---
Figure 1: Given a sequence of $P$ past point clouds (left in red) at time $T$,
the goal is to predict the $F$ future scans (right in blue). Note that all
point clouds are in the sensor’s coordinate system. Best viewed in color.

The main contribution of this paper is a novel range image-based encoder-
decoder neural network to jointly process spatio-temporal information from
points clouds using 3D convolutions. Our method can obtain structural details
of the environment by using skip connections and horizontal consistency using
circular padding, and in the end, provides more accurate predictions than
other state-of-the-art approaches for point cloud prediction. We make three
main key claims: Our approach is able to (i) predict a sequence of detailed
future 3D point clouds from a given input sequence by a fast joint spatio-
temporal point cloud processing using temporal 3D convolutional networks
(CNNs), (ii) outperform state-of-the-art point cloud prediction approaches,
(iii) generalize well to unseen environments, and operate online faster than a
typical rotating 3D LiDAR sensor frame rate (10 Hz). These claims are backed
up by the paper and our experimental evaluation. The open-source code and the
trained models are available at https://github.com/PRBonn/point-cloud-
prediction.

###  2 Related Work

Extracting spatio-temporal information from point cloud sequences has been
exploited in the literature for different tasks related to point cloud
prediction. Liu et al. [10] extend PointNet++ [6], which extracts spatial
features from a single point cloud, and propose FlowNet3D for estimating the
3D motion field between two point clouds. FlowNet3D is used by Lu et al. [11]
to interpolate intermediate frames between two LiDAR scans. Other deep
learning-based scene flow estimation methods can be found in the literature
[12, 13]. In contrast to scene flow between two past frames, point cloud
prediction estimates future point clouds. In addition to that, scene flow
requires labels like point-wise correspondences or segmentation masks. This
makes learning costly. In contrast to that, we focus on the self-supervised
generation of new, future point clouds.

Predicting future trajectories of traffic participants based on their past
motion is one possible application of LiDAR-based future scene prediction.
Given LiDAR scans only, these methods rely on preceding detection and tracking
modules to infer past trajectories. For example, Luo et al. [7] develop a 3D
convolutional architecture based on a bird’s eye view (BEV) of voxelized LiDAR
scans transformed by the sensor’s ego-motion. Inspired by this, Casas et al.
[14] increase the prediction horizon with IntentNet and also estimate a high-
level driver behavior from HD maps. Liang et al. [15] integrate a tracking
module into the end-to-end pipeline to improve the temporal consistency of the
predictions while Meyer et al. [3] and Laddha et al. [16] explore the range
image representation to predict future 3D bounding boxes. A second method to
estimate a future scene representation without the need for prior object
detections is prediction performed at the sensor level. Hoermann et al. [17]
directly predict a future occupancy grid map from a sequence of past 3D LiDAR
scans and Song et al. [18] predict a motion flow for indoor 2D LiDAR maps. Wu
et al. [19] use a spatio-temporal pyramid network to estimate voxelized BEV
classification and future motion maps. Recently, Toyungyernsub et al. [20]
differentiate between dynamic and static cells for grid map prediction but
require tracking labels to segment moving cells. In contrast to the
aforementioned approaches, our method operates on the full-size point clouds
without the need for voxelization. As shown by Weng et al. [8], a reversed
trajectory prediction pipeline can take the point cloud predictions and detect
and track objects in future scans. This makes it possible to first train the
prediction of point clouds self-supervised and then use off-the-shelf
detection and tracking modules to get future trajectories.

There has been a large interest in vision-based prediction. Srivastava et al.
[21] propose an LSTM [22] encoder-decoder model that takes an input sequence
of image patches or feature vectors and outputs the sequence of future patches
or feature vectors, respectively. Shi et al. [23] propose ConvLSTMs that
employ convolutional structures in the temporal transitions imposed by the
LSTM. This makes it possible to jointly encode temporal and spatial
correlations in image sequences. The integration of 3D convolutions into LSTMs
has been researched by Wang et al. [24]. Aigner et al. [25] omit the use of
recurrent structures and develop a video prediction approach called FutureGAN
based on 3D convolutions combined with a progressively growing GAN. Using the
range image representation of LiDAR scans, we can transfer the techniques from
vision-based methods for solving point cloud prediction. However, there are
significant differences between video prediction and range image-based point
cloud prediction. The main difference is that range images contain 3D point
information rather than RGB values. In addition to that, range images have a
360 degrees field of view, which causes a strong correlation between the left
and right image borders because of the principle of a typical rotating LiDAR
sensor, e.g., a Velodyne or an Ouster scanner.

Though vision-based and voxel-based prediction methods have been well studied,
predicting full-scale future point clouds is still challenging and has not yet
been fully explored in existing works. To predict a future sequence of
unordered and unstructured point clouds, Fan et al. [26] propose a point
recurrent neural network (PointRNN). Their method models spatio-temporal local
correlations of point features and states. Lu et al. [9] introduce MoNet,
which integrates a point-based context and motion encoder network and
processes these features in a novel recurrent neural network. A 3D scene flow-
based point cloud prediction approach has been recently developed by Deng and
Zakhor [27]. All three methods are 3D point-based but only consider a limited
number of points. Weng et al. [8] propose an encoder-decoder network that
takes a highly compressed feature representation of a full-scale LiDAR scan
and models the temporal correlations with a recurrent auto-encoder. The
feature vectors are obtained either by a point-based architecture or 2D
convolutional neural networks processing 2D range image representations. Deep
learning-based generation of LiDAR scans from a 2D range image has been
investigated by Caccia et al. [28] as well. In contrast to existing
approaches, we propose a technique that concatenates the projected 2D range
images to a spatio-temporal 3D representation, which makes it possible to use
3D convolutions without the need for recurrent architectures for temporal
modeling. As discussed by Bai et al. [29] for general sequence modeling, using
convolutional networks results in an architecture that is easier to train and
faster at inference since the input range images are no longer processed
sequentially. Compared to MoNet [9], this allows us to predict more future
point clouds in a shorter amount of time.

###  3 Our Approach

The goal of our approach is to predict a future sequence of full-scale LiDAR
points given a sequence of past scans. As illustrated in the overview of our
method in Fig. 2, we first project each past 3D point cloud into a range image
representation and concatenate the 2D images to obtain a 3D spatio-temporal
tensor, see also Appendix A for details. Second, we pass this tensor to our
encoder-decoder architecture to extract the spatial and temporal correlations
with 3D convolutional kernels as described in Sec. 3.1. We use skip
connections and circular padding to maintain details and horizontal
consistency of the predicted range images, see Sec. 3.2. Details on the
training procedure are provided in Sec. 3.3.

Spatio-Temporal
---
3D CNN
Spherical
---
Projection
$\mathcal{S}_{T}$
---
$\mathcal{S}_{T{-}1}$
---
$\mathcal{S}_{T{-}\text{P}{+}1}$
---
Mask and
---
Re-project
$\hat{\mathcal{S}}_{T{+}2}$
---
$\hat{\mathcal{S}}_{T{+}\text{F}}$
---
Spherical
---
Projection
Range
---
Mask and
---
Re-project
Mask and
---
Re-project
Valid Point
---
Mask
Range
---
Spherical
---
Projection
$\hat{\mathcal{S}}_{T{+}1}$
---

Figure 2: Overview of our approach. At time $T$, the past point clouds are
first projected into a 2D range image representation and then concatenated.
After passing through our proposed spatio-temporal 3D CNN network, the
combined predicted mask and range tensors are re-projected to obtain the
future 3D point cloud predictions.

####  3.1 Spatio-Temporal Encoder-Decoder Architecture

Our approach first converts the 3D LiDAR points into spherical coordinates and
maps them to image coordinates resulting in a dense 2D range image from which
3D information can still be recovered. If no point is projected into a
specific pixel, its value is set to zero. If multiple points are mapped to the
same pixel, we keep the closest point to focus on the prediction of foreground
objects. For technical details on the range image projection see Appendix A.
The main task of the encoder-decoder architecture illustrated in Fig. 3 is to
jointly extract spatio-temporal features from the input sequence and to output
future range image predictions. Compared to a multilayer perceptron
architecture that disregards spatial and temporal dimensions, the use of
convolutions imposes an inductive bias of having local correlations along
these dimensions. Convolutional networks usually require less trainable
parameters and are less prone to overfitting. When using range images for
point cloud prediction with 3D convolutions, sufficient receptive fields along
the spatial and temporal dimensions between encoder and decoder are required
to capture the motion of points in the past frames and to propagate their
locations into the future range images.

The encoder takes the 3D input tensor of size $(P,H,W)$, containing $P$ range
images with height $H$ and width $W$, and first standardizes the range values
based on mean and standard deviation computed from the training data. We pass
the standardized tensor through a 3D convolutional input layer with $C$
kernels to get a $C$-channel feature representation of size $(C,P,H,W)$.
Inspired by FutureGAN [25], the encoder block receives a tensor of size
$(C_{l},P_{l},H_{l},W_{l})$ at each encoder stage $l$, applies a combination
of 3D convolution, 3D batch normalization (BN) [30], and leaky ReLU activation
function [31] while maintaining the size of the tensor. A strided 3D
convolution follows resulting in a downsampled tensor of size
$(C_{l+1},P_{l}{-}p_{l},\frac{H_{l}}{h_{l}},\frac{W_{l}}{w_{l}})$ where
$h_{l}$ and $w_{l}$ are predefined downsampling factors for the spatial
feature size and $p_{l}$ reduces the temporal feature dimension. We batch
normalize the resulting tensor and apply a leaky ReLU activation function. The
kernel size is $(p_{l},h_{l},w_{l})$ and the stride $(1,h_{l},w_{l})$ to
achieve the desired downsampling. The downsampling compresses the sequential
point cloud feature representation and forces the network to learn meaningful
spatio-temporal features.

To predict future range images for $F$ future time steps, the decoder
subsequently upsamples the feature tensor to the final output size
$(2,F,H,W)$. Note that the number of future range images is fixed in the
architecture, but longer prediction horizons can be achieved in an auto-
regressive manner. This is done by sequentially feeding back the predicted
range images as the input tensor. In this work, we will only focus on the
prediction of a fixed number of future point clouds. The decoder architecture
is a mirrored version of the encoder. First, a transposed 3D convolutional
layer with kernel size $(p_{l},h_{l},w_{l})$ and stride $(1,h_{l},w_{l})$
increase the feature map. We insert a 3D batch normalization and leaky ReLU
activation layer before a second 3D CNN, 3D BN, and leaky ReLU follows.
Finally, we pass the output tensor of size $(C,F,H,W)$ through a 3D CNN output
layer with two kernels of size $(1,1,1)$ and stride $(1,1,1)$ and apply the
final Sigmoid function to get normalized values between 0 and 1. The first
output channel maps to a predefined range interval resulting in the future
range predictions. As done by Weng et al. [8], the second channel contains a
probability for each range image point to be a valid point for re-projection.
The re-projection mask keeps all points with a probability above 0.5. This
makes it possible to mask out e.g. the sky since there are no ground truth
points available.

3D CNN K(3,3,3) S(1,1,1) P(1,1,1)
---
\+ 3D BN + LeakyReLU
3D CNN K({1,2},2,4) S(1,2,4)
---
\+ 3D BN + LeakyReLU
Concat. + 3D CNN K(3,3,3) S(1,1,1) P(1,1,1)
---
\+ 3D BN + LeakyReLU
ConvTranspose3D K({1,2},2,4) S(1,2,4)
---
\+ 3D BN + LeakyReLU
Standardization
---
\+ 3D CNN K(1,1,1) S(1,1,1)
3D CNN K(1,1,1) S(1,1,1)
---
\+ Sigmoid

(64,3,16,256)
---

(16,5,64,2048)
---

(16,5,64,2048)
---

(128,2,8,32)
---

(256,1,4,8)
---

(512,1,2,2)
---

(32,4,32,512)
---

(256,1,4,8)
---

(32,4,32,512)
---

(128,2,8,32)
---

(256,1,4,8)
---

(128,2,8,32)
---

(64,3,169,256)
---

(64,3,169,256)
---

(32,4,32,512)
---

Predicted Range
---

Predicted Mask Logits
---

Stacked Range Images
---

Figure 3: Our detailed spatio-temporal architecture using 3D convolutional
neural networks with feature map sizes $(C_{i},P_{i},H_{i},W_{i})$. Solid and
dashed arrows indicate information flow and skip connections, respectively.
Details of each block are in the colorized boxes with kernel sizes K, stride S
and padding P. The feature map is not padded in case P is not specified. To be
viewed in color.

####  3.2 Skip Connections and Horizontal Circular Padding

We use strided 3D convolutions and downsample the range images as described in
Sec. 3.1. The reduced size of the feature space however causes a loss of
details in the predicted range images. We address this problem by adding skip
connections [32] between the encoder and decoder to maintain details from the
input scene. As shown in Fig. 3, the feature maps bypass the remaining
encoding steps and the mirrored decoder stage concatenates them with the
previously upsampled feature volume along the channel dimension. Concatenation
enables the network to account for the temporal offset between encoder and
decoder feature maps. A combination of 3D convolutions, 3D batch
normalization, and leaky ReLU follows to merge the features back to the
original number of channels while maintaining the temporal and spatial
dimensions. We investigate the effect of skip connections in Sec. 4.1 and Sec.
4.4 and show that they maintain details in the predicted point clouds.

Another challenge when using range images for 3D point cloud prediction is to
maintain spatial consistency on the horizontal borders of the range images.
Range images obtained by rotating LiDAR sensors such as a Velodyne or an
Ouster sensor, are panoramic images and have strong horizontal correlations
between the image borders. If the ego vehicle is rotating along the vertical
$z$ axis, an object passing the left border of the range image will appear on
the right border. To take this property into account, we introduce circular
padding for the width dimension. We pad the left side of each feature map with
its right side and vice versa. The vertical dimension is padded with zeros. An
experiment on the effectiveness of this padding is provided in Sec. 4.1.

####  3.3 Training

When training the network, we project the ground truth point clouds into the
range images of size $(H,W)$, which allows for the computation of 2D image-
based losses. We slice the data into samples of sequences consisting of $P$
past and $F$ future frames, where subsequent samples are one frame apart. Our
architecture is trained with a combination of multiple losses. The average
range loss $\mathcal{L}_{\text{R},t}$ at time step $t$ between a predicted
range image $\left(\hat{r}_{i,j}\right)_{t}~{}{\in}~{}\mathbb{R}^{H\times W}$
and a ground truth range image
$\left(r_{i,j}\right)_{t}~{}{\in}~{}\mathbb{R}^{H\times W}$ is given by

$\displaystyle\mathcal{L}_{\text{R},t}$
$\displaystyle=\frac{1}{HW}\sum\limits_{i,j}\Delta
r_{i,j},\quad\textrm{with}\quad\Delta
r_{i,j}=\left\\{\begin{array}[]{cl}{\|\hat{r}_{i,j}-r_{i,j}\|}_{1}&,\text{if
}r_{i,j}>0\\\ 0&,\text{otherwise}\end{array}\right.,$ (3)

such that the range loss is only computed for valid ground truth points. We
train the mask output at time step $t$ with predicted probabilities
$\left(\hat{m}_{i,j}\right)_{t}~{}{\in}~{}\mathbb{R}^{H\times W}$ with a
binary cross-entropy loss

$\mathcal{L}_{\text{M},t}=\frac{1}{HW}\sum\limits_{i,j}-y_{i,j}\log(\hat{m}_{i,j})-(1-y_{i,j})\log(1-\hat{m}_{i,j}),$
(4)

where $y_{i,j}$ is 1 in case the ground truth point is valid and 0 otherwise.
Both losses only consider the predicted 2D range images and not the re-
projected point clouds to compute the loss faster. For comparison to 3D-based
methods like MoNet [9], we fine-tune our model on a 3D point-based loss. A
common metric for comparing 3D point clouds is the Chamfer distance [33]. At
time $t$, we re-project the masked range image into a 3D point cloud
$\hat{\mathcal{S}}_{t}=\\{\bm{\hat{p}}\in\mathbb{R}^{3}\\},|\hat{\mathcal{S}}_{t}|=N,$
and compare it with the ground truth point cloud
$\mathcal{S}_{t}=\\{\bm{p}\in\mathbb{R}^{3}\\},|\mathcal{S}_{t}|=M,$ by
computing

$\mathcal{L}_{\text{CD},t}=\frac{1}{N}\sum_{\bm{\hat{p}}\in\hat{\mathcal{S}}_{t}}{\min_{\bm{p}\in\mathcal{S}_{t}}{\|\bm{\hat{p}}-\bm{p}\|}_{2}^{2}}+\frac{1}{M}\sum_{\bm{p}\in\mathcal{S}_{t}}{\min_{\bm{\hat{p}}\in\hat{\mathcal{S}}_{t}}{\|\bm{\hat{p}}-\bm{p}\|}_{2}^{2}}.$
(5)

The computation of the Chamfer distance is in general slow due to the search
for nearest neighbors. A fast to compute image-based loss using
$\mathcal{L}_{\text{R},t}$ and $\mathcal{L}_{\text{M},t}$ is about 2.5 times
faster compared to a loss including the 3D Chamfer distance. We propose a
training scheme with only range and mask loss for pre-training. This results
in a good initialization for fine-tuning including a Chamfer distance loss. We
provide an experiment on the training scheme in Sec. 4.4. Given a current time
step $t~{}{=}~{}T$, our total loss function for $F$ future time steps is given
by

$\mathcal{L}=\sum_{t=T+1}^{T+F}\mathcal{L}_{\text{R,t}}+\alpha_{\text{M}}\mathcal{L}_{\text{M,t}}+\alpha_{\text{CD}}\mathcal{L}_{\text{CD,t}},$
(6)

with tunable weighting parameters $\alpha_{\text{M}}$ and
$\alpha_{\text{CD}}$. During the experiments, we experienced that setting
$\alpha_{\text{M}}~{}{=}~{}1.0$ and $\alpha_{\text{CD}}~{}{=}~{}0.0$ for pre-
training and $\alpha_{\text{M}}~{}{=}~{}1.0$ and
$\alpha_{\text{CD}}~{}{=}~{}1.0$ for fine-tuning worked best. We train our
network parameters with the Adam optimizer [34], a learning rate of $10^{-3}$,
and accumulate the gradients for 16 samples. We use an exponential decaying
learning rate schedule. All models require less than 50 epochs to converge.

###  4 Experimental Evaluation

In this work, we present an approach for self-supervised point cloud
prediction using spatio-temporal 3D convolutional neural networks. We present
our experiments to show the capabilities of our method and to support our key
claims, which are: (i) predicting a sequence of detailed future 3D point
clouds from a given input sequence by a fast joint spatio-temporal point cloud
processing using temporal 3D convolutional networks, (ii) outperforming state-
of-the-art point cloud prediction approaches, (iii) generalizing well to
unseen environments, and operate online faster than a typical rotating 3D
LiDAR sensor frame rate of 10 Hz. We follow the experimental setting of MoNet
[9] with the KITTI Odometry dataset [35] and predict the next 5 frames from
the past 5 LiDAR scans captured at 10 Hz. Please see Appendix B for further
details on the experimental setting.

####  4.1 Qualitative Results

The first experiment supports our claim that our model can predict a sequence
of detailed future 3D point clouds using a temporal convolutional network and
demonstrate the effectiveness of our proposed skip connections and circular
padding. The situation in Fig. 4 from test sequence 08 shows an ego-vehicle
finishing a left turn while passing a close vehicle. This is a challenging
scenario because of the different scales of close and far objects and the ego-
motion of the car. Our approach pre-trained on range and mask loss only is
capable of predicting the future point clouds across 5 time steps and
maintains details in the scene. Without skip connections, small objects are
lost due to the compression (green circle). With typical zero padding, the
points belonging to the car in the back of the ego vehicle (black circle) are
not consistently predicted since they span across the range image borders,
while our approach can solve this problem by using the proposed circular
padding.

Ours (Pre-trained)
---
No Skip Connections
---
No Circular Padding
---
T+5
---
T+1
---
T+1
---
T+1
---
T+5
---
T+5
---
Figure 4: Qualitative comparison of our pre-trained method and two ablated
models with predicted (blue) and ground truth (red) point clouds from test
sequence 08. Given a past sequence of 5 point clouds at time $T$, the upper
row shows the first predicted frame at $T{+}1$ and the lower row the last
predicted frame at $T{+}5$. The green circle encloses an area where details
are lost without skip connections and the points in the black circle
demonstrate that circular padding maintains spatial consistency at the range
image border. Best viewed in color.

####  4.2 Quantitative Results

To support the claim that our method outperforms state-of-the-art point cloud
prediction approaches, we compare the quantitative results to multiple
baselines. Note that Lu et al. [9] operate and also evaluate on sub-sampled
point clouds, whereas our approach can predict full-size point clouds. To
account for this, Tab. 1 (Left) shows the performance of our method on sampled
point clouds with two different sampling rates used in MoNet [9]. For both
samplings, our method achieves a better result for larger prediction steps.
However, it can be seen that sampling points from the original point clouds
during evaluation strongly affects the Chamfer distance of our full-sized
prediction. Thus, we provide an evaluation on full-scale point clouds to
demonstrate the true performance of our approach, see Tab. 1 (Right). We
compare our results to an Identity baseline that takes the last received scan
as a constant prediction for all future time steps. The Constant Velocity and
Ray Tracing baselines estimate the transformation between the last two
received scans with a SLAM approach [36] as a constant velocity prediction of
the sensor. At each prediction step, the Constant Velocity baseline transforms
the last received point cloud according to the predicted motion. The Ray
Tracing baseline aggregates all five transformed past point clouds to render a
more dense range image. Since a lot of points in the scene are static, these
are very strong baselines that work well for constant ego-motion. Note that in
contrast to our method, all three baselines do not consider points belonging
to moving objects. The results in Tab. 1 show that our method outperforms all
baselines at larger prediction steps and that it can reason about the sensor’s
ego-motion and the motion of moving points. The standard deviation of the
Chamfer distance indicates that the results of our approach are more
consistent throughout the test set whereas the baseline results suffer from
outliers in case of moving objects or non-linear ego-motion.

| Sampled Point Clouds | Full-scale Point Clouds
---|---|---
Prediction | PointLSTM | Scene Flow | MoNet | Ours | Ours | Identity | Const. Vel. | Ray Tracing | Ours
Step |  [26] |  [9] | [9] | 32,768 pts | 65,536 pts | Baseline | Baseline | Baseline |
1 | 0.332 | 0.503 | $\mathbf{0.278}$ | 0.367 | 0.288 | 0.271 (0.144) |  $\mathbf{0.105}$ (0.067) | 0.158 (0.095) | 0.254 (0.124)
2 | 0.561 | 0.854 | 0.409 | 0.446 | $\mathbf{0.352}$ | 0.719 (0.448) |  $\mathbf{0.217}$ (0.221) | 0.253 (0.246) | 0.310 (0.188)
3 | 0.810 | 1.323 | 0.549 | 0.546 | $\mathbf{0.428}$ | 1.216 (0.798) | 0.385 (0.443) | 0.388 (0.460) |  $\mathbf{0.378}$ (0.262)
4 | 1.054 | 1.913 | 0.692 | 0.638 | $\mathbf{0.509}$ | 1.727 (1.169) | 0.604 (0.711) | 0.556 (0.719) |  $\mathbf{0.448}$ (0.370)
5 | 1.299 | 2.610 | 0.842 | 0.763 | $\mathbf{0.615}$ | 2.240 (1.156) | 0.852 (0.968) | 0.751 (0.971) |  $\mathbf{0.547}$ (0.487)
Mean (Std) | 0.811 | 1.440 | 0.554 | 0.552 | $\mathbf{0.439}$ | 1.235 (1.192) | 0.433 (0.641) | 0.421 (0.627) |  $\mathbf{0.387}$ (0.331)
Table 1: Chamfer distance results on KITTI Odometry test sequences 08 to 10 in
$\left[\text{m}^{2}\right]$ with sampled (left) and full-scale (right) points
clouds. We first train our model on range and mask loss and then fine-tune
including a Chamfer distance loss for 10 epochs as described in Sec. 3.3. Best
mean results in bold for sampled and full-scale point clouds, respectively.
For full-scale point clouds, the standard deviation of the Chamfer distance is
given in parentheses. Please see Appendix E for a more detailed statistical
analysis.

####  4.3 Runtime

For the runtime, our approach takes on average 11 ms (90 Hz) for predicting 5
future point clouds with up to 131,072 points each on a system with an Intel
i7-6850K CPU and an Nvidia GeForce RTX 2080 Ti GPU, which is faster than the
sensor frame rate, i.e. 10 Hz, of a typical rotating 3D LiDAR sensor and
faster than existing approaches like MoNet [9] using the same GPU. Compared to
a runtime of 280 ms reported by Lu et al. [9] for MoNet for predicting 5
future point clouds with 65,536 points each, our approach is 25 times faster
for twice as many points.

####  4.4 Ablation Study

In this experiment, we provide a more detailed analysis of the effect of the
model architecture and the proposed training scheme. We average the metrics
over all validation samples and time steps as shown in Tab. 2. Note that we
always report the Chamfer distance but it is only used during training if
specified. First, we modify the size of the network by decreasing (Small) or
increasing (Large) the number of channels of our model by a factor of 2. The
larger model predicts slightly better masked range images but also requires 4
times more parameters. All losses increase without skip connections (No Skip
Connections). This is because a lot of objects are not predicted due to the
large feature compression. The effect of circular padding (No Circular
Padding) is not visible quantitatively because of the small number of affected
points and we refer to Sec. 4.1 for a qualitative evaluation. Finally, we show
that our fine-tuned model (Ours) achieves similar performance on the
validation set compared to a model optimized with the Chamfer distance loss
from the beginning (CD Loss From Start). Our proposed training scheme with
pre-training on fast-to-compute image-based losses gives a good pre-trained
model for fine-tuning with the 3D Chamfer distance loss while reducing
training time.

| Small | Large | No Skip | No Circular | CD Loss | Ours | Ours
---|---|---|---|---|---|---|---
Property |  |  | Connections | Padding | From Start | (Pre-trained) |
4 M Parameters | ✗ |  |  |  |  |  |
17 M Parameters |  |  | ✗ | ✗ | ✗ | ✗ | ✗
68 M Parameters |  | ✗ |  |  |  |  |
Skip Connections | ✗ | ✗ |  | ✗ | ✗ | ✗ | ✗
Circular Padding | ✗ | ✗ | ✗ |  | ✗ | ✗ | ✗
Loss Weight $\alpha_{\text{CD}}$ | 0.0 | 0.0 | 0.0 | 0.0 | 1.0 | 0.0 | 0.0 $\rightarrow$ 1.0
$\mathcal{L}_{\text{R}}$ $\left[\text{m}\right]$ | 0.867 | 0.741 | 1.258 | 0.801 | 0.805 | 0.798 | 0.858
$\mathcal{L}_{\text{Mask}}$ | 0.309 | 0.288 | 0.337 | 0.297 | 0.300 | 0.299 | 0.301
$\mathcal{L}_{\text{CD}}$ $\left[\text{m}^{2}\right]$ | 1.110 | 0.714 | 2.759 | 0.885 | 0.487 | 0.985 | 0.480
Table 2: Evaluation of models with varying architectures, padding, and
training schemes on the KITTI Odometry validation sequences 06 and 07. Each
variation (columns) is given by a combination of different design decisions
(rows) indicated by crosses (✗), e.g., variation ”Large” is a model with 68 M
parameters using skip connections and circular padding. The arrow
($\rightarrow$) indicates a pre-training and fine-tuning with two different
loss weights $\alpha_{\text{CD}}$ as discussed in Sec. 3.3.

####  4.5 Generalization

Figure 5: Predicted (blue) and ground truth (red) future point clouds on the
Apollo-SouthBay ColumbiaPark [37] test set at prediction step 5. Left: Model
only trained on KITTI odometry data. Right: After training one epoch on the
Apollo training set.

Finally, we aim at supporting our claim that our method generalizes well to
new, unseen environments. We first evaluate a model fine-tuned on KITTI
Odometry on the Apollo-SouthBay ColumbiaPark [37] test set. Note that this
data is from an unseen, _different_ environment (collected in the U.S.) with a
_different_ type of car setup, comparing to the training data (collected in
Germany). Our model achieves a mean Chamfer distance of $0.934$ $\text{m}^{2}$
on the Apollo test set. The performance gap compared to the results in Tab. 1
is due to the fact that the training data has a maximum range of 85 m but the
Apollo test data contains points up to a range of 110 m. Thus, we fine-tune
the model on the Apollo-SouthBay ColumbiaPark [37] training set for a single
epoch. This yields a mean Chamfer distance of $0.426$ $\text{m}^{2}$ and
demonstrates the generalization capability of our approach and the advantage
of self-supervised training. We provide an additional experiment on the effect
of using more training data in Appendix C. Fig. 5 shows that the fine-tuned
model reduces the Chamfer distance by predicting more accurate shapes and more
distant points. Note that compared to Fig. 4, the different sensor mount
causes a changed pattern of invalid points around the car. Our approach can
correctly predict the re-projection mask.

###  5 Conclusion

In this paper, we presented a novel approach to self-supervised point cloud
prediction with 3D convolutional neural networks. Our method exploits the
range image representation to jointly extract spatio-temporal features from
the input point clouds sequence. This allows us to successfully predict
detailed, full-scale, future 3D point clouds during online operation. We
implemented and evaluated our approach on different datasets and provided
comparisons to other existing techniques and supported all claims made in this
work. The experiments suggest that our method outperforms existing baselines
and generalizes well to different sensor mounts in unseen environments.

##### Acknowledgments

We thank Lu et al. [9] for providing implementation details and experimental
results.

This work has partially been funded by the Deutsche Forschungsgemeinschaft
(DFG, German Research Foundation) under Germany’s Excellence Strategy,
EXC-2070 – 390732324 – PhenoRob, by the European Union’s Horizon 2020 research
and innovation programme under grant agreement No 101017008 (Harmony), and by
the Federal Ministry of Food and Agriculture (BMEL) based on a decision of the
Parliament of the Federal Republic of Germany via the Federal Office for
Agriculture and Food (BLE) under the innovation support programme under
funding no 28DK108B20 (RegisTer).

### References

  * Chen et al. [2021] X. Chen, I. Vizzo, T. Läbe, J. Behley, and C. Stachniss.  Range Image-based LiDAR Localization for Autonomous Vehicles.  In _Proc. of the IEEE Intl. Conf. on Robotics & Automation (ICRA)_, 2021.
  * Yin et al. [2021] T. Yin, X. Zhou, and P. Krahenbuhl.  Center-Based 3D Object Detection and Tracking.  In _Proc. of the IEEE/CVF Conf. on Computer Vision and Pattern Recognition (CVPR)_ , 2021.
  * Meyer et al. [2021] G. P. Meyer, J. Charland, S. Pandey, A. Laddha, S. Gautam, C. Vallespi-Gonzalez, and C. K. Wellington.  LaserFlow: Efficient and Probabilistic Object Detection and Motion Forecasting.  In _Proc. of the IEEE Intl. Conf. on Robotics & Automation (ICRA)_, volume 6, pages 526–533, 2021.
  * Vizzo et al. [2021] I. Vizzo, X. Chen, N. Chebrolu, J. Behley, and C. Stachniss.  Poisson Surface Reconstruction for LiDAR Odometry and Mapping.  In _Proc. of the IEEE Intl. Conf. on Robotics & Automation (ICRA)_, 2021.
  * Chen et al. [2021] X. Chen, S. Li, B. Mersch, L. Wiesmann, J. Gall, J. Behley, and C. Stachniss.  Moving Object Segmentation in 3D LiDAR Data: A Learning-based Approach Exploiting Sequential Data.  _arXiv preprint_ , abs/2105.08971, 2021.
  * Qi et al. [2017] C. Qi, K. Yi, H. Su, and L. J. Guibas.  PointNet++: Deep Hierarchical Feature Learning on Point Sets in a Metric Space.  In _Proc. of the Conference on Neural Information Processing Systems (NeurIPS)_ , 2017.
  * Luo et al. [2018] W. Luo, B. Yang, and R. Urtasun.  Fast and Furious: Real Time End-to-End 3D Detection, Tracking and Motion Forecasting with a Single Convolutional Net.  In _Proc. of the IEEE Conf. on Computer Vision and Pattern Recognition (CVPR)_ , pages 3569–3577, 2018.
  * Weng et al. [2020] X. Weng, J. Wang, S. Levine, K. Kitani, and N. Rhinehart.  Inverting the Pose Forecasting Pipeline with SPF2: Sequential Pointcloud Forecasting for Sequential Pose Forecasting.  In _Proc. of the Conf. on Robot Learning (CoRL)_ , 2020.
  * Lu et al. [2020] F. Lu, G. Chen, Y. Liu, Z. Li, S. Qu, and T. Zou.  MoNet: Motion-based Point Cloud Prediction Network.  _arXiv preprint_ , abs/2011.10812, 2020.
  * Liu et al. [2019] X. Liu, C. R. Qi, and L. J. Guibas.  FlowNet3D: Learning Scene Flow in 3D Point Clouds.  In _Proc. of the IEEE/CVF Conf. on Computer Vision and Pattern Recognition (CVPR)_ , 2019.
  * Lu et al. [2021] F. Lu, G. Chen, S. Qu, Z. Li, Y. Liu, and A. Knoll.  PointINet: Point Cloud Frame Interpolation Network.  In _Proc. of the Conference on Advancements of Artificial Intelligence (AAAI)_ , 2021.
  * Ushani and Eustice [2018] A. K. Ushani and R. M. Eustice.  Feature Learning for Scene Flow Estimation from LIDAR.  In _Proc. of the Conf. on Robot Learning (CoRL)_ , pages 283–292, 2018.
  * Zhai et al. [2021] M. Zhai, X. Xuezhi, L. Ning, and X. Kong.  Optical Flow and Scene Flow Estimation: A Survey.  _Pattern Recognition_ , 114:107861, 02 2021.
  * Casas et al. [2018] S. Casas, W. Luo, and R. Urtasun.  IntentNet: Learning to Predict Intention from Raw Sensor Data.  In _Proc. of the Conf. on Robot Learning (CoRL)_ , pages 947–956, 2018.
  * Liang et al. [2020] M. Liang, B. Yang, W. Zeng, Y. Chen, R. Hu, S. Casas, and R. Urtasun.  PnPNet: End-to-End Perception and Prediction With Tracking in the Loop.  In _Proc. of the IEEE/CVF Conf. on Computer Vision and Pattern Recognition (CVPR)_ , 2020.
  * Laddha et al. [2020] A. Laddha, S. Gautam, G. P. Meyer, and C. Vallespi-Gonzalez.  RV-FuseNet: Range View based Fusion of Time-Series LiDAR Data for Joint 3D Object Detection and Motion Forecasting.  _arXiv preprint_ , abs/2005.10863, 2020.
  * Hoermann et al. [2018] S. Hoermann, M. Bach, and K. Dietmayer.  Dynamic occupancy grid prediction for urban autonomous driving: A deep learning approach with fully automatic labeling.  In _Proc. of the IEEE Intl. Conf. on Robotics & Automation (ICRA)_, 2018.
  * Song et al. [2019] Y. Song, Y. Tian, G. Wang, and M. Li.  2D LiDAR Map Prediction via Estimating Motion Flow with GRU.  In _Proc. of the IEEE Intl. Conf. on Robotics & Automation (ICRA)_, 2019.
  * Wu et al. [2020] P. Wu, S. Chen, and D. N. Metaxas.  MotionNet: Joint Perception and Motion Prediction for Autonomous Driving Based on Bird’s Eye View Maps.  In _Proc. of the IEEE/CVF Conf. on Computer Vision and Pattern Recognition (CVPR)_ , 2020.
  * Toyungyernsub et al. [2021] M. Toyungyernsub, M. Itkina, R. Senanayake, and M. J. Kochenderfer.  Double-Prong ConvLSTM for Spatiotemporal Occupancy Prediction in Dynamic Environments.  In _Proc. of the IEEE Intl. Conf. on Robotics & Automation (ICRA)_, 2021.
  * Srivastava et al. [2015] N. Srivastava, E. Mansimov, and R. Salakhutdinov.  Unsupervised Learning of Video Representations using LSTMs.  In _Proc. of the Int. Conf. on Machine Learning (ICML)_ , 2015.
  * Hochreiter and Schmidhuber [1997] S. Hochreiter and J. Schmidhuber.  Long Short-Term Memory.  _Neural Computation_ , 9(8):1735–1780, 1997.
  * Shi et al. [2015] X. Shi, Z. Chen, H. Wang, D.-Y. Yeung, W.-K. Wong, and W.-C. Woo.  Convolutional LSTM Network: A Machine Learning Approach for Precipitation Nowcasting.  In _Proc. of the Advances in Neural Information Processing Systems (NIPS)_ , 2015.
  * Wang et al. [2019] Y. Wang, L. Jiang, M.-H. Yang, L.-J. Li, M. Long, and L. Fei-Fei.  Eidetic 3D LSTM: A Model for Video Prediction and Beyond.  In _Proc. of the Int. Conf. on Learning Representations (ICLR)_ , 2019\.
  * Aigner and Korner [2018] S. Aigner and M. Korner.  FutureGAN: Anticipating the Future Frames of Video Sequences using Spatio-Temporal 3d Convolutions in Progressively Growing GANs.  In _Proc. of the IEEE/CVF Conf. on Computer Vision and Pattern Recognition (CVPR)_ , 2018.
  * Fan and Yang [2019] H. Fan and Y. Yang.  PointRNN: Point Recurrent Neural Network for Moving Point Cloud Processing.  _arXiv preprint_ , abs/1910.08287, 2019.
  * Deng and Zakhor [2020] D. Deng and A. Zakhor.  Temporal LiDAR Frame Prediction for Autonomous Driving.  In _Proc. of the Intl. Conf. on 3D Vision (3DV)_ , pages 829–837, 2020.
  * Caccia et al. [2019] L. Caccia, H. Van Hoof, A. Courville, and J. Pineau.  Deep Generative Modeling of Lidar Data.  In _Proc. of the IEEE/RSJ Intl. Conf. on Intelligent Robots and Systems (IROS)_ , 2019.
  * Bai et al. [2018] S. Bai, J. Kolter, and V. Koltun.  An Empirical Evaluation of Generic Convolutional and Recurrent Networks for Sequence Modeling.  _arXiv preprint_ , abs/1803.01271, 2018.
  * Ioffe and Szegedy [2015] S. Ioffe and C. Szegedy.  Batch normalization: Accelerating deep network training by reducing internal covariate shift.  _arXiv preprint_ , abs/1502.03167, 2015.
  * Maas et al. [2013] A. Maas, A. Hannun, and A. Ng.  Rectifier Nonlinearities Improve Neural Network Acoustic Models.  In _Proc. of the Int. Conf. on Machine Learning (ICML)_ , 2013.
  * Ronneberger et al. [2015] O. Ronneberger, P.Fischer, and T. Brox.  U-Net: Convolutional Networks for Biomedical Image Segmentation.  In _Medical Image Computing and Computer-Assisted Intervention_ , volume 9351 of _LNCS_ , pages 234–241. Springer, 2015.
  * Fan et al. [2017] H. Fan, H. Su, and L. Guibas.  A point set generation network for 3d object reconstruction from a single image.  In _Proc. of the IEEE/CVF Conf. on Computer Vision and Pattern Recognition (CVPR)_ , pages 605–613, 2017.
  * Kingma and Ba [2015] D. Kingma and J. Ba.  Adam: A Method for Stochastic Optimization.  In _Proc. of the Int. Conf. on Learning Representations (ICLR)_ , 2015\.
  * Geiger et al. [2012] A. Geiger, P. Lenz, and R. Urtasun.  Are we ready for Autonomous Driving? The KITTI Vision Benchmark Suite.  In _Proc. of the IEEE Conf. on Computer Vision and Pattern Recognition (CVPR)_ , pages 3354–3361, 2012.
  * Behley and Stachniss [2018] J. Behley and C. Stachniss.  Efficient Surfel-Based SLAM using 3D Laser Range Data in Urban Environments.  In _Proc. of Robotics: Science and Systems (RSS)_ , 2018.
  * Lu et al. [2019] W. Lu, Y. Zhou, G. Wan, S. Hou, and S. Song.  L3-Net: Towards Learning Based LiDAR Localization for Autonomous Driving.  In _Proc. of the IEEE/CVF Conf. on Computer Vision and Pattern Recognition (CVPR)_ , 2019.
  * Milioto et al. [2019] A. Milioto, I. Vizzo, J. Behley, and C. Stachniss.  RangeNet++: Fast and Accurate LiDAR Semantic Segmentation.  In _Proceedings of the IEEE/RSJ Int. Conf. on Intelligent Robots and Systems (IROS)_ , 2019.
  * Chen et al. [2019] X. Chen, A. Milioto, E. Palazzolo, P. Giguère, J. Behley, and C. Stachniss.  SuMa++: Efficient LiDAR-based Semantic SLAM.  In _Proc. of the IEEE/RSJ Intl. Conf. on Intelligent Robots and Systems (IROS)_ , 2019.

## Supplementary Material

###  Appendix A Point Cloud Sequence Representation

This section provides details on our range image representation. We convert
each 3D LiDAR point ${\mbox{\boldmath$p$}}=(x,y,z)$ to spherical coordinates
and map them to image coordinates $(u,v)$ resulting in a projection
$\Pi:\mathbb{R}^{3}\mapsto\mathbb{R}^{2}$

$\displaystyle(u,v)^{\top}=\left(\frac{1}{2}\left[1-\arctan(y,x)\,\pi^{-1}\right]W,\left[1-\left(\arcsin(z\,r^{-1})+\mathrm{f}_{\mathrm{up}}\right)\mathrm{f}^{-1}\right]H\right)^{\top},$
(7)

with $(H,W)$ as the height and width of the range image and
$\mathrm{f}~{}{=}~{}\mathrm{f}_{\mathrm{up}}~{}{+}~{}\mathrm{f}_{\mathrm{down}}$
is the vertical field-of-view of the sensor. Each range image pixel stores the
range value $r~{}{=}~{}{||{\mbox{\boldmath$p$}}||}_{2}$ of the projected point
or zero if there exist no projected point in the corresponding pixel. In case
multiple points project to the same pixel due to rounding, we keep the closer
one since closer points are more important than far-away points. A range image
pixel with coordinates $u$ and $v$ and range $r$ is re-projected by solving
Eq. (7) for pitch $\theta~{}{=}~{}\arctan(y,x)$ and yaw
$\gamma~{}{=}~{}\arcsin(z\,r^{-1})$ and computing

$\displaystyle(x,y,z)^{\top}=(r\cos(\theta)\cos(\gamma),r\cos(\theta)\sin(\gamma),r\sin(\theta))^{\top}.$
(8)

In contrast to the existing work using range images [38, 39], we found that
additionally using $x,y,z,$ and intensity values did not improve the
performance of point cloud prediction. We concatenate the range images along
the temporal dimension to obtain the input 3D tensor with a size of $(P,H,W)$
which is further processed by our encoder-decoder network.

###  Appendix B Experimental Setting

As explained in the experimental section of our main paper, we follow the
experimental setup of MoNet [9]. Thus, we train the models on sequences 00 to
05 of the KITTI Odometry dataset [35], validate on sequences 06 and 07, and
test on sequences 08 to 10. We project the point clouds from the Velodyne
HDL-64E laser scanner into range images of size
$H~{}{\times}~{}W~{}{=}~{}64~{}{\times}~{}2048$ according to Eq. (7). As
explained in Sec. 3.1, we map the normalized output values of the network to a
predicted range interval defined before training. Based on the minimum and
maximum range values of the training data, we set the minimum predicted range
to 1 m and the maximum range to 85 m for KITTI Odometry [35] and to 110 m for
Apollo [37].

###  Appendix C Amount of Training Data

This experiment demonstrates the advantage of a large amount of data for
training our point cloud prediction method. Since the approach is self-
supervised, data obtained from the LiDAR sensor can be directly used for
training without expensive labeling. Fig. 6 shows the total validation loss on
sequences 06 and 07 after 50 epochs with respect to the number of training
iterations. Note that the final number of training iterations within 50 epochs
depends on the total amount of samples in the training data. Fig. 6 shows the
validation loss with the previously explained experimental setting on the left
(yellow). If we train a model on a larger dataset containing the KITTI
sequences 00 to 05 plus 11 to 15, the validation loss decreases faster and
converges to a lower minimum (blue). For the largest dataset containing
sequences 00 to 05 plus 11 to 21 (green), the validation loss is further
reduced indicating that more data improves the model’s performance on the
unseen validation set. The individual loss components in Tab. 3 show that we
can improve all losses by using more data. This emphasizes that self-
supervised point cloud prediction profits from a large amount of sensor data
without the need for expensive labeling.

Figure 6: Total validation loss for a given number of training iterations with
different training data sizes. The solid line is a fitted exponential curve.

| 14,082 | 21,832 | 34,334
---|---|---|---
Loss | Samples | Samples | Samples
$\mathcal{L}_{\text{R}}$ [m] | 0.798 | 0.735 | $\mathbf{0.710}$
$\mathcal{L}_{\text{Mask}}$ | 0.299 | 0.290 | $\mathbf{0.285}$
$\mathcal{L}_{\text{CD}}$ $\left[\text{m}^{2}\right]$ | 0.985 | 0.795 | $\mathbf{0.768}$

Table 3: Final validation loss components for different numbers of training
dataset sizes. We train all methods for 50 epochs.

###  Appendix D Additional Qualitative Results

In this section, we provide additional analysis of our predicted masked range
images and the re-projected 3D point clouds. Our method can estimate the ego-
motion of the vehicle resulting in an accurate prediction of the environment
in Fig. 7. The moving cyclist in front of the ego vehicle is consistently
predicted in all five future frames which demonstrates that our predicted
range images like in Fig. 8 preserve details in the scene.

T+3
---
T
---
T+1
---
T+2
---
T+4
---
T+5
---
Figure 7: Last received point cloud at time T and the predicted next 5 future
point clouds. Ground truth points are shown in red and predicted points in
blue. The scene is taken from the KITTI Odometry sequence 08 that was not used
during training.  GT Range
---
Masked Prediction
---
Predicted Range
---
GT Mask
---
Predicted Mask
---
Figure 8: Ground truth and predictions of range image and re-projection mask
at the last prediction step. The range is color-coded with increasing range
from blue to yellow. Invalid points are shown in dark blue. For the mask,
yellow indicates a valid point whereas dark blue represents invalid points.
The masked prediction is the combination of predicted range and re-projection
mask.

###  Appendix E Statistical Analysis

For a more detailed quantitative analysis for full-scale point cloud
prediction, we compare box plots for each prediction step in Fig. 10 to Fig.
14 and across all prediction steps in Fig. 14. Each plot shows the median
(orange line), minimum (lower bar), maximum (upper bar), first quartile (lower
colorized rectangle), and third quartile (upper colorized rectangle) of the
Chamfer distances computed between the predicted and ground truth point clouds
for the test set. As shown in the main paper for the mean Chamfer distances,
the constant velocity baseline achieves the lowest median Chamfer distances
for short prediction horizons. Since all baselines ignore moving objects for
the prediction, our method outperforms them for larger time steps resulting in
a lower median and interquartile ranges. The statistical results support our
claims made in the paper.

Figure 9: Prediction step 1

Figure 10: Prediction step 2

Figure 11: Prediction step 3

Figure 12: Prediction step 4

Figure 13: Prediction step 5

Figure 14: All steps