##  Patches of developable surfaces bounded by NURBS curves

L<EMAIL_ADDRESS>
($\flat$) ETSI Navales,

Universidad Politécnica de Madrid

Avenida de la Memoria 4, 28040-Madrid, Spain.

###  1 Introduction

In this talk we review the problem of constructing a developable surface patch
bounded by two rational or NURBS (Non-Uniform Rational B-spline) curves.

NURBS curves are curves which are piecewise rational. That is, they are a
generalisation of spline curves, which are piecewise polynomial curves.
Similarly we define NURBS surfaces and solids. NURBS curves have been the
standard in Computer Aided Design [1] for a long time through IGES and STEP
specifications. They are described by a list of points called control polygon
and two list of numbers: the list of weights and the list of knots.

Developable surfaces are ruled surfaces with null Gaussian curvature [2]. This
implies that they can be constructed from planar surfaces by just cutting,
rolling and folding, so that metric properties such as lengths and angles
between curves and areas are preserved. These geometric properties are of
great interest for steel and textile industry, since these are pieces designed
in the plane and then combed into space.

For instance, in naval architecture sheets of steel are adapted to fit into
the hull of a ship [3, 4, 5]. If these sheets are combed just with a folding
machine, the costs are lower than if they require the use of heat. In textile
industry cloth is planar and is cut and sewn to produce garments and the
quality is improved if it is not stretched [6]. They have also been used for
designing facades in architecture [7] and in automobile industry [8].

In geometric design the standard relies on the use of rational B-spline curves
and surfaces (NURBS), which are described by control polygons or nets and
lists of weights and knots. In the case of ruled surfaces, the control net is
formed by just the control polygons of the bounding curves, since segments can
be described by just their endpoints, which form their control polygon.

This problem has been addressed in several ways [9], but the key drawback is
that when we require the developable surface to be NURBS and bounded by NURBS
curves, the possibilities are restricted [10, 11].

In some cases, the bounding curves are planar and lie on parallel planes [12,
13] and this simplifies the problem.

For instance, one can obtain general solutions for developable NURBS surfaces
bounded by NURBS curves [15, 16, 17, 18, 19].

Considering the dual space in projective geometry has been also profitable,
since developable surfaces may be seen as envelopes of lines of planes [20,
21].

For this reason our proposal is to consider developable surface patches which
are not NURBS, though bounded by NURBS curves [22]. In fact, we are able to
obtain every possible solution, good or bad, to this problem in our framework.

Our approach is based on performing a reparametrisation on one of the bounding
curves of the surface patch and imposing that the resulting surface be
developable [22]. One would expect a differential equation for the
reparametrisation function. However, the condition happens to be algebraic.
Another approach [23] reduces the null Gaussian curvature condition to
quadratic equations.

Moreover, if the bounding curves are (piecewise) polynomial or rational of
degree $n$, the developability condition is an algebraic equation of degree
$2n-2$. The degree may be lowered to $n-1$ if the curves are (piecewise)
polynomial and lie on parallel planes.

###  2 Methods

We start with a ruled surface parametrised by $b(t,v)$ and bounded by two
parametrised curves, $c(t)$, $d(t)$,

$\displaystyle b(t,v)=(1-v)c(t)+vd(t),\quad t,v\in[0,1].$

For given $c(t)$ and $d(t)$, this ruled surface will not be developable in
general, but, if it is developable, this feature shall not depend on the
chosen parametrisation.

The Gaussian curvature $K(t,v)$ is the quotient of the determinants of the
second, $B(t,v)$, and first, $G(t,v)$, fundamental forms of the surface patch
parametrised by $b(t,v)$,

$K(t,v)=\frac{\det B(t,v)}{\det G(t,v)}.$

The first fundamental form $G(t,v)$ is the usual scalar product restricted to
tangent vectors to the surface at the point $b(t,v)$ and is used, for
instance, for calculating areas,

$\mathrm{Area}(S)=\int_{D}\sqrt{\det G(t,v)}\,dtdv,$

but, since its determinant is positive and appears at the denominator of
$K(u,v)$, we need not pay attention to it, though due to Gauss’ Theorema
Egregium [2], the Gaussian curvature is an intrinsic property and may be
written in terms of just the first fundamental form.

On the other hand, the second fundamental form $B(t,v)$ is the extrinsic
curvature of the surface at the point $b(t,v)$ with $\nu$ as unitary normal to
the surface and can be constructed with the projections of the second
derivatives of the parametrisation along the unitary normal $\nu(t,v)$ to the
surface at $b(t,v)$,

$B(t,v)=\left(\begin{array}[]{cc}b_{tt}\cdot\nu&b_{tv}\cdot\nu\\\
b_{vt}\cdot\nu&b_{vv}\cdot\nu\\\ \end{array}\right)_{(t,v)}.$

In this sense, the second fundamental form is used to compute the normal
curvature at points on the surface.

For a ruled surface parametrised as in (2),

$\det
B(t,v)=\left|\begin{array}[]{cc}\left((1-v)c^{\prime\prime}(t)+vd^{\prime\prime}(t)\right)\cdot\nu&\left(d^{\prime}(t)-c^{\prime}(t)\right)\cdot\nu\\\
\left(d^{\prime}(t)-c^{\prime}(t)\right)\cdot\nu&0\\\
\end{array}\right|_{(t,v)}=-\left(\left(d^{\prime}(t)-c^{\prime}(t)\right)\cdot\nu\right)^{2},$

we see that the Gaussian curvature is always non-positive and hence there are
no elliptic points on a ruled surface, since $\det B(t,v)$ is always positive.

Moreover, vanishing Gaussian curvature is equivalent to vanishing
$\left(d^{\prime}(t)-c^{\prime}(t)\right)\cdot\nu$ on the points of the ruled
surface.

Since at a point $b(t,v)$ on the surface we have
$b_{t}(u,v)=(1-v)c^{\prime}(t)+d^{\prime}(t)$, $b_{v}(t,v)=d(t)-c(t)$ as two
tangent vectors to the surface, we may construct a normal vector to the
surface at $b(t,v)$ as $N=b_{t}\times b_{v}$, and then having a vanishing
Gaussian curvature is equivalent to having a vanishing triple product,

$c^{\prime}(t)\cdot
d^{\prime}(t)\times(d(t)-c(t))=\det\left(c^{\prime}(t),d^{\prime}(t),d(t)-c(t)\right)=0,\qquad
t\in[0,1]$ (1.1)

at all points of the surface patch.

Figure 1: Developability is equivalent to having the same tangent plane along
each ruling of a ruled surface

As we see in Fig 1, this is equivalent to requiring that the three vectors
$c^{\prime}(t)$, $d^{\prime}(t)$ and $d(t)-c(t)$ lie on the same plane for all
values of $t$. This also implies that the normal to the surface is the same
for all points on the segment (ruling) linking $c(t)$ and $d(t)$.

Our contribution to deal with this problem is based on reparametrisation of
one of the bounding curves by a function $T(t)$,

$\tilde{b}(t,v)=(1-v)c(t)+vd(T(t))$ (1.2)

and require $\tilde{b}(t,v)$ to satisfy the null Gaussian curvature condition.

###  3 Results

The developability condition (1.1) applied to parametrisations such as (1.2)
can be seen to be algebraic in $T(t)$, since the dependence on the derivative
$T^{\prime}(t)$ is factored out by the determinant,

$\det\left(c^{\prime}(t),\dot{d}(T),d(T)-c(t)\right)=0,$ (1.3)

where the dot stands for derivation with respect to $T$.

In the case of (piecewise) polynomial or rational curves $c(t)$, $d(t)$,
further consequences may be derived:

Theorem 1: Let $c(t)$, $d(T)$, $t,T\in[0,1]$ be rational curves of degree $n$.
The parameterized ruled surface,

$b(t,v)=(1-v)c(t)+vd(T(t)),\qquad t,v\in[0,1],$

is a developable surface if the reparameterization function $T(t)$ satisfies
the algebraic equation

$\det\left(c^{\prime}(t),\dot{d}(T),d(T)-c(t)\right)_{T=T(t)}=0,$

and is a real monotonically increasing function of $t$.

This equation is of degree $2n-2$ at most. If both curves are (piecewise)
polynomial and lie on parallel planes, the equation is of degree $n-1$ at
most. We may see an example in Fig 2.

Figure 2: Developable patch bounded by two cubic spline curves

The price to pay is that solutions of this algebraic equation will not be
rational or polynomial in general and $\tilde{b}(t,v)$ will no longer be
NURBS.

Since the condition on the reparametrisation is algebraic, the number of
possible solutions is finite, but not all of them are geometrically
acceptable.

Figure 3: Developable surface with a regression area

For being a reparametrisation, $T(t)$ must be a monotonically increasing
function. Otherwise, we would have unpleasant regression areas with more than
on ruling through some points of the curves (See Fig 3). This can be checked
with the help of

$T^{\prime}(t)=\left.\frac{\det\left(c^{\prime\prime}(t),\dot{d}(T),d(T)-c(t)\right)}{\det\left(\ddot{d}(T),c^{\prime}(t),d(T)-c(t)\right)}\right|_{T=T(t)},$

which we derive from the null Gaussian curvature condition.

This implies that monotonicity is granted if

$\mathrm{sgn}\left(c^{\prime\prime}(t)\cdot\nu(t)\right)=\mathrm{sgn}\left.\left(\ddot{d}(T)\cdot\nu(t)\right)\right|_{T=T(t)},$

where $\nu(t)$ is the unitary normal to the ruled surface along the segment at
$t$. This means that the normal curvatures of both curves must have the same
sign for each value of $t$.

Hence, acceptable solutions just appear if both curves are qualitatively
similar regarding their curvature.

Theorem 2: Let $c(t)$, $d(T)$, $t,T\in[0,1]$ be parameterized curves. Let
$T(t)$ be a reparameterization function so that

$b(t,v)=(1-v)c(t)+vd(T(t)),\qquad t,v\in[0,1],$

is a developable surface. $T(t)$ is a monotonically increasing function if and
only if for all $t$,

$\mathrm{sgn}\left(c^{\prime\prime}(t)\cdot\nu(t)\right)=\mathrm{sgn}\left.\left(\ddot{d}(T)\cdot\nu(t)\right)\right|_{T=T(t)},$

where $\nu(t)$ is the unitary normal to the surface along the ruling at $t$.

Or equivalently, for the normal curvatures $k_{n,c}$, $k_{n,d}$ of both curves

$\mathrm{sgn}\left(k_{n,c}(t)\right)=\mathrm{sgn}\left(k_{n,d}(T(t)\right),$

for all values of $t$.

In the case of parameterizations of class $C^{k}$ of differentiability, $T(t)$
is of class $C^{k-1}$.

###  4 Conclusions

We have produced a new approach for dealing with the problem of constructing a
developable surface patch between two parametrised curves $c(t)$ and $d(t)$.

This approach is grounded on performing a reparametrisation of one of the
curves and the developability condition is not a differential equation, but an
algebraic equation, and provides all possible solutions to the problem.

In the case of (piecewise) polynomial or rational curves of degree $n$, the
developability condition is an algebraic equation of degree $2n-2$. Since the
most usual degree in Computer Aided Design is three, this means we are dealing
with a fourth degree equation, which can be handled either by numerical or
analytical methods.

For (piecewise) polynomial curves of degree $n$, lying on parallel planes, the
algebraic equation is of degree $n-1$.

Requiring that the reparametrisation function be a monotonically increasing
function, in order to avoid regression areas on the developable surface patch,
is achieved if the sign of the normal curvatures of both bounding curves is
the same at the endpoints of each ruling.

With this approach, it is easy to control the final class of differentiability
of the surface in the piecewise case.

## References

  * [1] Farin, G, Curves and surfaces for CAGD: a practical guide, 5th Edition. San Francisco, CA, USA, Morgan Kaufmann Publishers Inc., 2002\.
  * [2] Struik, D.J., Lectures on classical differential geometry, 2nd Edition, New York, Dover Publications Inc., 1988.
  * [3] Kilgore, U., Developable hull surfaces, in: Fishing Boats of the World, Vol. 3, Fishing News Books Ltd., Surrey, 425–431, 1967.
  * [4] Chalfant, J.S., Maekawa, T. Design for manufacturing using B-spline developable surfaces., _J. Ship Research_ 42: 207–215, 1998.
  * [5] Pérez, F., Suárez, J, Quasi-developable B-spline surfaces in ship hull design, _Computer-Aided Design_ 39: 853–862, 2007.
  * [6] Rose, K., Sheffer, A., Wither, J., Cani, M.P., Thibert, B., Developable surfaces from arbitrary sketched boundaries, in: A. Belyaev, M. Garland (Eds.), Eurographics Symposium on Geometry Processing (2007), The Eurographics Association, 163–172, 2007.
  * [7] Pottmann, H., Asperl, A., Hofer, M., Kilian, A., Architectural geometry, Exton, Bentley Institute Press, 2007.
  * [8] Frey, W.H., Bindschadler, D., Computer aided design of a class of developable Bézier surfaces, Tech. rep., GM Research Publication R&D-8057, 1993.
  * [9] Pottmann, H., Wallner, J., Computational line geometry, Mathematics and Visualization, Berlin, Springer-Verlag, Berlin, 2001.
  * [10] Fernández-Jambrina, L, Bézier developable surfaces, _Comput. Aided Geom. Design_ 55: 15–28, 2017.
  * [11] Fernández-Jambrina, L, Characterisation of rational and NURBS developable surfaces in Computer Aided Design, _Journal of Computational Mathematics_ 39: 550–568, 2021.
  * [12] Aumann, G., Interpolation with developable Bézier patches, _Comput. Aided Geom. Design_ 8: 409–420, 1991.
  * [13] Maekawa, T., Design and tessellation of B-spline developable surfaces, _ASME Transactions Journal of Mechanical Design_ 120: 453–461, 1998.
  * [14] Chu, C.H., Séquin, C.H., Developable Bézier patches: properties and design, _Computer Aided Design_ 34: 511–527 (2002).
  * [15] Aumann, G., A simple algorithm for designing developable Bézier surfaces, _Comput. Aided Geom. Design_ 20: 601–619, 2003
  * [16] Aumann, G., Degree elevation and developable Bézier surfaces, _Comput. Aided Geom. Design_ 21: 661–670, 2004.
  * [17] Fernández-Jambrina, L, B-spline control nets for developable surfaces, _Comput. Aided Geom. Design_ 24: 189–199, 2007.
  * [18] Cantón, A., Fernández-Jambrina, L., Interpolation of a spline developable surface between a curve and two rulings, F _rontiers of Information Technology & Electronic Engineering_ 16:173–190 (2015).
  * [19] Cantón, A., Fernández-Jambrina, L, Non-degenerate developable triangular Bézier patches, in: J.-D. Boissonnat, P. Chenin, A. Cohen, C. Gout, T. Lyche, M.-L. Mazure, L. Schumaker (Eds.), Curves and Surfaces, Vol. 6920 of Lecture Notes in Computer Science, Springer Berlin / Heidelberg, 207–219, 2012.
  * [20] Bodduluri, R.M.C, Ravani, B., Design of developable surfaces using duality between plane and point geometries, _Computer Aided Design_ 25: 621–632, 1993.
  * [21] H. Pottmann, G. Farin, Developable rational Bézier and $B$-spline surfaces, Comput. Aided Geom. Design 12 (5) (1995) 513–531.
  * [22] Fernández-Jambrina, L, Pérez-Arribas, F, Developable surfaces bounded by NURBS curves, _Journal of Computational Mathematics_ 38: 693–709, 2020.
  * [23] Tang, C., Bo, Wallner, J., Pottmann, H., Interactive design of developable surfaces, _ACM Trans. Graph._ 35: 1–12, 2016.


#  Metal-semiconductor behavior along the line of stacking order change
in gated multilayer graphene

W. Jaskólski<EMAIL_ADDRESS>Institute of Physics, Faculty of Physics,
Astronomy and Informatics, Nicolaus Copernicus University in Torun,
Grudziadzka 5, 87-100 Toruń, Poland

###### Abstract

We investigate gated multilayer graphene with stacking order change along the
armchair direction. We consider some layers cracked to release shear strain at
the stacking domain wall. The energy cones of graphene overlap along the
corresponding direction in the $k$-space, so the topological gapless states
from different valleys also overlap. However, these states strongly interact
and split due to atomic-scale defects caused by the broken layers, yielding an
effective energy gap. We find that for some gate voltages, the gap states
cross and the metallic behavior along the stacking domain wall can be
restored. In particular cases, a flat band appears at the Fermi energy. We
show that for small variations of the gate voltage the charge occupying this
band oscillates between the outer layers.

multilayer graphene; topological states; defects in graphene

##  I Introduction

Stacking order domain walls in multilayer graphene are most usually considered
along the zigzag direction. They can be created when some layers are
stretched, corrugated, or delaminated along the armchair direction
Ju_Nature_2015 ; Pelc_2015 ; Lane_2018 ; Anderson_PRB_2022 . Stacking order
change causes the appearance of topological states. When the multilayer is
gated and the gap opens, the one-dimensional topological gapless states appear
at valleys $K$ and $K^{\prime}$ separated by $\frac{4}{3}\pi/a$, Vaezi_2013 ;
Zhang_Pnas_2013 ; Alden_2013 ; San_Jose_2014 ; Jaskolski_2020 ; Ju_Nature_2015
. The gapless states are then valley-protected as long as there is no valley
mixing caused, for example, by atomic-scale defects. It was however shown that
even under strong defects like bonds breaking in some layers at the domain
wall, some gap states that connect the valence band (VB) and conduction band
(CB) continua persist, and the systems preserve metallic character along the
stacking domain wall Jaskolski_2019 ; Jaskolski_2021 ; Jaskolski_2023 .

The situation is different if we allow stacking domain wall along the armchair
direction. Firstly, because along the armchair direction, the valleys appear
at the same $K=0$ Yin_NC_2016 . Therefore, the gapless states, which have
opposite slopes $E(k)$ at each valley, overlap and may no longer be valley-
protected. Secondly, because such stacking change causes strong shear strain
at the interface between different stackings, at least for narrow domain
walls.

Figure 1:  Schematic representations of the investigated systems. (a) BLG with
stacking order change AB/BA. (b) TLG with stacking orders ABC/BAC. (c) FLG
with stacking order change ABCA/BACB. (d) Lateral view of sequences of carbon
atoms in the consecutive layers along the armchair direction. Vertical broken
lines mark the interface between different stackings. The layers from bottom
to top are marked in colors: blue, red, green, and dark grey. The bottom
(blue) layer in (c) is covered and not seen.

In this work, we investigate the behavior of gapless states in such systems.
We consider multilayers in which the interface bonds are broken in some layers
to release any shear strain in the system. Fig. 1 shows how the unstrained
bilayer (BLG), trilayer (TLG), and four-layer graphene (FLG) with stacking
domain walls along the armchair direction can be formed. One can easily
imagine that such cracks of layers may be present in realistic multilayers
Cracks_Beitner_NL_2023 ; Jang_Nanoscale_2017 . The bonds broken in some layers
act as atomic-scale defects. Therefore, the gapless states strongly interact
and split yielding the appearance of the effective energy gap and thus
semiconducting properties of such multilayers. However, we show that for some
gate voltages, some gap states cross and the energy gap disappears. The
metallic behavior along the stacking domain wall can be restored. Transitions
from semiconductor to metal can be achieved by changing the gate voltage.
Moreover, for some gate voltages, the separated gap states in FLG connect at
the Fermi level ($E_{F}$) leading to a kind of flat band at $E_{F}$, with
charge occupying this band localized mainly in the outer layers.

##  II System and method

We investigate Bernal stacked bilayer and rhombohedral four-layer graphene
with stacking order domain walls. Fig. 1 demonstrates how the consecutive
layers are deposited. Panels (a) and (c) show the investigated BLG and FLG
systems, and panel (b) shows an intermediate system of the TLG. Only small
fragments of the multilayers are shown, while the investigated systems are
infinite in both in-plane directions. The bottom layer is always a uniform
graphene, as well as the top layer in the case of FLG. The nodes of the FLG
top layer have the same in-plane positions as the nodes of the bottom layer.
The remaining layers are broken, they consist of two nonconnected half-
infinite graphene layers shifted relative to each other by $a_{C-C}$ along the
armchair direction. This allows for the formation of a strainless interface
with different stacking orders of rhombohedral multilayers on both sides. Note
that due to the presence of the stacking domain wall, the systems studied are
periodic only in the direction defined by the domain wall, i.e., in the
armchair direction. Fig. 1 (d) shows the stacking order of layers on both
sides of the domain wall. For BLG they are simply Bernal stackings AB and BA,
respectively, while for FLG we have ABCA and CBAC stackings.

To investigate the electronic structure of the defined multilayers we
calculate the local density of states (LDOS) at the interface. The
calculations are performed using the Green function matching technique for a
system composed of a conductor connected to two semi-infinite leads
Nardelli_1999 . In our case, the conductor consists of the interface region,
as defined in Fig. 1 (a)-(c), connected on the left- and right-hand sides to
two semi-infinite multilayers. The local density of states is calculated as
${\rm{LDOS}}(E)=-\frac{1}{\pi}{\rm{Im[TrG_{C}}}(E)]$, where ${\rm{G_{C}}}$ is
the interface Green function determined by the Hamiltonians of the conductor
and the leads. The Hamiltonians are calculated using a tight-binding
approximation. The values of intra-layer and inter-layer hopping parameters
are $t_{i}=2.7$ eV and $t_{e}=0.27$ eV, respectively, following the
Castro_2007 ; Ohta_2006 ; Jaskolski_2016 ; Jaskolski_2D . Since the system is
periodic in the armchair direction, the local density of states is
$k$-dependent, i.e., LDOS($E$,$k$), where $k$ is the wave vector corresponding
to this periodicity.

Figure 2:  LDOS for BLG with stacking order domain wall along the armchair
direction, as visualized in Fig. 1 (a), for different gate voltages $V$
applied between the layers. (a) $V=0.025$ eV, (b) $V=0.075$ eV, (c) $V=0.15$
eV, (d) $V=0.4$ eV. The Fermi level lies in the middle of the gap defined by
the energy continua.  Figure 3:  LDOS for FLG with stacking order domain wall
along the armchair direction, as visualized in Fig. 1 (c), for different gate
voltages $V$. (a) $V=0.025$ eV, (b) $V=0.075$ eV, (c) $V=0.15$ eV, (d) $V=0.4$
eV. The Fermi level lies in the middle of the energy gap.

##  III Results and discussion

We perform calculations for defined graphene multilayers in a perpendicular
external electric field. Several values of the field are considered, namely,
E=7, 22, 44, and 120 mV/Å. Assuming the interlayer distance $d=3.35$ Å the
values of E translate into gate voltages between the neighbor layers
$V=0.025$, 0.075, 0.15, and 0.4 eV. The voltages we consider are available in
the experiment and are of the same order as typically used for multilayer
graphene systems Jaskolski_2D ; Li_NN_2016 ; Sena_Peeters_PRB_2014 ;
LC_Carbon_2024 .

In Fig. 2 we show LDOS calculated in the BLG system for four values of the
gate voltage $V$. Of the four gapless states that may appear at all at $K=0$
in the BLG with a stacking domain wall Yin_NC_2016 , there is little left in
the energy gap. For $V=0.4$ eV, four separated flat bands are visible as a
remnant of interacting topological gapless states (two of them are situated
very close to the band continua). This is caused by the atomic-scale crack in
the upper layer. For smaller $V$, fewer bands are visible in the energy gap.
Finally, for $V=0.025$ eV, the band that ”grows up” from the VB continuum
connects with the conduction band and the system takes on a metallic character
along the domain wall.

Figure 3 shows LDOS for FLG system. Now, the situation is more complex. The
stacking order change may, in principle, lead to even eight gapless states at
$K=0$, but here the number of bands in the energy gap depends strongly on $V$.
For very small $V=0.025$ eV the overall picture resembles two pairs of
intersecting topological states with opposite slopes. They slightly repeel at
the crossing points, so a small gap opens at the $E_{F}$. For higher values,
e.g., $V=0.075$ eV, this small gap closes, the band crossing points approach
the cone center and two new bands detach from the VB and CB continua.
Increasing farther the gate, the gap opens (for $V\approx 0.1$ eV), but closes
again for $V$ higher than 0.3 eV and the system restores again its metallic
character.

Especially interesting is the case when $V\approx 0.092$ eV. The corresponding
LDOS is presented in Fig. 4 (a). The four bands in the energy gap show flat
areas and a short flat (in the range of $k$) band is formed also at the Fermi
level where two bands touch. The increased density of this flat band at
$E=0.14$ eV is well seen in Fig. 4, where the projected density of states,
i.e., PDOS($E$)=$\int_{k\in IBZ}{\rm{LDOS}}(E,k)dk$ is shown. The distribution
of band densities in the layers is visualized in Fig. 4 (c) and (d). Of the
two bands that touch at the Fermi level, the band below the $E_{F}$ is
localized in two lower layers, while the band above the $E_{F}$ localizes in
the two upper layers. Slightly away from the cone center, i.e., for $k=0.1$,
the distribution over layers is uniform. However, for $k=0$, i.e., just at the
flat band, a fraction of the charge moves to the outer layers. The similar
effect of leakage of charge to the outer layers was recently observed also in
LC_Carbon_2024  for a system of double twisted trilayers.

In low temperatures the flat band at $E_{F}$ is only half-filled. Note, that
in graphene systems such a band can give rise to magnetic effects
Kusakabe_book ; Bouzerar_PRB_2023 ; Sawada_APL_2014 ; Pons_PRB_2020  and
superconductivity Cao_N_2018 ; Han_NN_2023 ; Yin_PRB_2020 ; Chittari_PRL_2019
. The width of the flat band at $E_{F}$ is about 5 meV. In Fig. 5 we show the
LDOS($E$,$k=0$) of this band resolved into the layers, for two values of the
gate voltage, $V=0.091$ eV and $V=0.093$ eV, for which the falt band is still
present at $E_{F}$. This result confirms that at $k=0$ the density is
localized mainly in the outer layers (blue and black). More importantly, we
find that a slight increase of $V$ by only a few mV, moves the charge
occupying the flat band at $E_{F}$, from the top to the bottom layer.

Figure 4:  (a) LDOS for FLG and $V=0.092$ eV. (b) PDOS, i.e., LDOS integrated
over $k$ in the IBZ. (c) and (d) LDOS distribution in the layers for $k=0.1$
and $k=0$, respectively. The density in consecutive layers, from bottom to
top, is shown in colors: blue, red, green, and black.  Figure 5:  The
LDOS($E$,$k=0$) distribution in the layers for FLG and two different gate
voltages $V$. Left - $V=0.091$ eV, right - $V=0.093$ eV. The density in
consecutive layers is shown in colors: blue, red, green, and black.

The calculations have been also performed for trilayer (Fig. 1 (b)), but since
the results do not differ substantially from those for BLG, they are not
presented. We have also checked that wider separation of the half-graphene
layers does not essentially change the results. The band crossing points and
the separation of bands can occur for different voltages, but the overall
picture of the gap states and their behavior with the changing $V$ is similar.

##  IV Conclusions

We have studied the electronic structure of gated bilayer and four-layer
graphene with stacking order domain walls created along the armchair
direction. Bernal and rhombohedral stackings have been considered for BLG and
FLG, respectively. Since the change of stacking across the armchair direction
causes shear strain at the domain wall, we allow some layers to crack, i.e.,
to break C-C bonds along the domain wall to release the strain.

We have investigated the local density of states along the domain wall, i.e.,
at the interface of the adjacent stacking orders. The graphene energy cones
$E(k)$ overlap in the corresponding direction in the $k$-space. Thus, all the
possible topological gapless states appear at the same valley, $K=0$. Since
the broken bonds act as atomic-scale defects and lead to the removal of valley
protection of gapless topological states, these states interact and split
yielding an effective energy gap. We have, however, found that for some gate
voltages, the metallic character of the multilayer along the domain wall is
restored. In the case of BLG, it happens for very small gates. In the case of
the four-layer graphene, the gap states cross or anticross, i.e., they close
or open the energy gap, depending on the gate.

In the case of four-layer graphene and for a specific value of the gate
voltage, some bands in the gap connect forming a kind of a short flat band at
the Fermi energy. We have demonstrated that the density of this band is moved
to the outer layers. Finally, we have also shown that for very small changes
in the gate voltage, the charge that fills the flat band at the $E_{F}$, can
change its localization between the bottom and the top layers. This gives us a
tool for changing the layer localization of one-dimensional currents that at
$E_{F}$ can flow along the stacking domain wall.

## References

  * [1] L. Ju, Z. Shi, N. Nair, Y. Lv, C. Jin, J. Velasco, Jr., C. Ojeda-Aristizabal, H. A. Bechtel, M. C. Martin, A. Zettl, J. Analytis, and F. Wang, Nature (London) 520, 650 (2015).
  * [2] M. Pelc, W. Jaskolski, A. Ayuela, and L. Chico, Phys. Rev. B 92, 085433 (2015).
  * [3] T. L. M. Lane, M. Andelkovic, J. R. Wallbank, L. Covaci, F. M. Peeters, and V.I. Falko, Phys. Rev. B 97, 045301 (2018).
  * [4] P. Anderson, Y. Huang, Y. Fan, S. Qubbaj, S. Coh, Q. Zhou, and C. Ojeda-Aristizabal, Phys. Rev. B 105, L081408 (2022).
  * [5] A. Vaezi, Y. Liang, D. H. Ngai, L. Yang, and E.-A. Kim, Phys. Rev. X 3, 021018 (2013).
  * [6] F. Zhang, A. H. McDonald, and E. J. Mele, Proc. Natl. Acad. Sci., 110 10546 (2013).
  * [7] J. S. Alden, A. W. Tsen, P. Y. Huang, R. Hovden, L. Brown, J. Park, D. A. Muller, and P. L McEuen, Proc. Natl. Acad. Sci. 110, 11256 (2013).
  * [8] P. San-Jose, R. V. Gorbachev, A. K. Geim, K. S. Novoselov, and F. Guinea, Nano Lett. 14, 2052 (2014).
  * [9] W. Jaskolski and G. Sarbicki, Phys. Rev. B 102, 035424 (2020).
  * [10] W. Jaskolski, Phys. Rev. B 100, 035436 (2019).
  * [11] W. Jaskolski, Mol. Phys. 120, e2013554 (2022).
  * [12] W. Jaskolski, Solid State Communication 360, 115043 (2023).
  * [13] L.-J. Yin, H. Jiang, J.-B. Qiao, and L. He, Nat. Commun. 7, 11760 (2016).
  * [14] D. Beitner, S. Amiaty, S. S. Atri, A. McEllistrim, T. Coen, V. I. Falko, S. Richter, M. B. Shalom, and H. Suchowski, Nano Lett. 23, 10758 (2023).
  * [15] B. Jang, B. Kim, J.-H. Kim, H.-J. Lee, T. Sumigawa, and T. Kitamura, Nanoscale 9, 17325 (2017).
  * [16] M. B. Nardelli, Phys. Rev. B 60, 7828 (1999).
  * [17] E. V. Castro, K. S. Novoselov, S. V. Morozov, N. M. R. Peres, J. M. B. L. dos Santos, J. Nilsson, F. Guinea, A. K. Geim, and A. H. C. Neto, Phys. Rev. Lett. 99, 216802 (2007).
  * [18] T. Ohta, A. Bostwick, T. Seyller, K. Horn, and E. Rotenberg, Science 313, 951 (2006).
  * [19] W. Jaskolski, M. Pelc, L. Chico, and A. Ayuela, Nanoscale 8, 6079 (2016).
  * [20] W. Jaskolski, M. Pelc, G. w. Bryant, L. Chicko, and A. Ayuela, 2D Materials 5 025006 (2018).
  * [21] J. Li, K. Wang, K. J. McFaul, Z. Zern, Y. Ren, K. Watanabe, T. Taniguchi, Z. Qiao, and J. Zhu, Nat. Nanotechnol. 11, 1060 (2016).
  * [22] S. H. R. de Sena, J. M. Pereira, F. M. Peeters, and G. A. Farias, Phys. Rev. B 89, 035420 (2014).
  * [23] F. Perez Riffo, S. Vizcaya, E. Menendez-Proupin, J. M. Florez, L. Chico, and E. Suarez Morell, Crbon 222 118952 (2024).
  * [24] K. Kusakabe, in Carbon-based magnetism, Ed. T. L. Makarova and F. Palacio, Elsevier, Amsterdam (2006).
  * [25] G. Bouzerar, Phys. Rev. B 107 184441 (2023).
  * [26] K. Sawada, F. Ishii, and M. Saito, App. Phys. Lett. 104 143111 (2014).
  * [27] R. Pons, A. Mielke, and T. Stauber, Phys. Rev. B 102 235101 (2020).
  * [28] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Nature 556 43 (2018).
  * [29] T. Han, Z. Lu, G. Scuri, J. Sung, J. Wang, T. Han, K. Watanabe, T. Taniguchi, H. Park, and L. Ju, Nat. Nanotechnol. 19 181 (2023).
  * [30] L.-J. Yin, L.-Z. Yang, L. Zhang, Q. Wu, X. Fu, L.-H. Tong, G. Yang, Y. Tian, L. Zhang, and Z. Qin, Phys. Rev. B 102 241403(R) (2020).
  * [31] B. L. Chittari, G. Chen, Y. Zhang, F. Wang, and J. Jung, Phys. Rev. Lett. 122, 016401 (2019).


Finite-volume energy shift of the three-pion ground state

Fabian Müller
Helmholtz-Institut für Strahlen- und Kernphysik (Theorie) and Bethe Center for Theoretical Physics, Universität Bonn, 53115 Bonn, Germany
Akaki Rusetsky
Helmholtz-Institut für Strahlen- und Kernphysik (Theorie) and Bethe Center for Theoretical Physics, Universität Bonn, 53115 Bonn, Germany
Tbilisi State  University,  0186 Tbilisi, Georgia
Tiansu Yu
Helmholtz-Institut für Strahlen- und Kernphysik (Theorie) and Bethe Center for Theoretical Physics, Universität Bonn, 53115 Bonn, Germany

Using the framework of nonrelativistic effective field theory, the finite-volume
ground-state energy shift is calculated up to and including $O(L^{-6})$
for the system of three pions in the channel with the total isospin $I=1$.
The relativistic corrections are included
perturbatively, up to the same order in the inverse of the box size $L$.
The obtained explicit expression, together with the known result for the system with maximal isospin $I=3$,
can be used for the extraction of two independent effective three-body couplings
from the measured ground-state spectrum of three pions.

§ INTRODUCTION

The extraction of various hadronic observables in the three particle sector is certainly a
challenging task for lattice QCD. Such an extraction has become truly realistic only recently,
owing to the rapid growth of computational resources available, as well as the progress
on the methods and algorithms [1, 2, 3, 4, 5, 6, 7]
(see also the closely related
papers [8, 9], where
the simulations were performed within the scalar $\varphi^4$-theory). There is, however,
more to this task than merely carrying out larger simulations which require more
resources. Namely, the extraction of the hadronic parameters, defined in the infinite
volume, from the observables, measured on a finite lattices, has proven to be significantly
more challenging than in the one- and two-particle sectors. The problem has been,
in principle, solved by the three-particle quantization condition, which has been derived
only recently in three different
settings [10, 11, 12, 13, 14, 15]. In its turn, this development has sparked activities in the field, as seen, e.g., from [16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39] (note also earlier
work on the related issues [40, 41, 42, 43]). The reader is referred to [44] for a recent review on the three-particle problem in a finite volume. An alternative approach
to the problem focuses on the perturbative shift of the free
three-particle energy levels. Albeit more restricted from
the beginning (e.g., it is clear that the resonances cannot be handled in this
manner), the approach is very popular at present. The reason for this is that
the perturbative approach is very transparent and convenient for the analysis of
lattice data. In case of identical particles,
the expansion parameter is given by $a/L$, where $a$ denotes the
two-body scattering length and $L$ is the spatial size of the cubic box. Up
to and including $O(L^{-5})$, only two-body threshold parameters – the scattering
length $a$ and the effective range $r$ – enter the expression for the energy shift, and the latter appears at this order only in the excited states. At order $L^{-6}$, a single additional
coupling $\eta_3$, characterizing the three-body nonderivative interaction,
appears in the perturbative energy shift, and hence the measured three-body energies at different values of $L$ can be used to extract the value of this coupling from data. Moreover, the perturbative approach can be directly generalized for the systems of four and more particles (which is presently
not the case for the exact quantization condition), and the measurements
in the many-particle sector can be used to restrict the fit of $\eta_3$, since the same coupling appears in all $N$-particle energies up to and including $O(L^{-6})$. An example of such a fit is given in the recent article [9].

Application of the perturbation theory to the calculation of the energy shifts
in a finite volume has a decades-long history [45, 46, 47, 48]. This development culminated in
Refs. [49, 50] where, using the nonrelativistic effective field theory, the $N$-particle ground-state
shift is derived up to and including order $L^{-6}$ and $L^{-7}$, respectively.
The relativistic corrections at $O(L^{-6})$ are calculated in the same approach
in Ref. [51] where, in addition, the inclusion of the electromagnetic interactions is considered. Further, carrying out the threshold expansion
in the nonrelativistic limit, the expression for the energy shift of the ground and first excited states is derived in Ref. [30] at $O(L^{-6})$.
On the other hand, the threshold
expansion of the relativistic quantization condition leads to the relativistic expression of the ground-state energy at the same order [17]. Finally, in Ref. [9], where the nonrelativistic
framework has been used,
the relativistic corrections in the ground and first excited states have
been calculated up to and including $O(L^{-6})$.

All above studies have one thing in common – they are carried out for the system of $N$ identical particles. Taking into account that the simplest three-particle system, which can be simulated in lattice QCD, is that of three pions, one may conclude that the expressions, obtained in the above papers, apply to the
case with maximal isospin only. It is,
however, interesting to simulate systems with a different total isospin, in order to extract
all low-energy
three-pion couplings. In addition, studying the three-pion decays on the lattice,
one does not encounter, in general, a three-pion system with maximal isospin in the final
A perturbative expression for the energy shift for $I\neq 3$
at the required accuracy is, however, unavailable in the literature (to the best of
our knowledge, a single study [25], which treats different
isospin channels, deals with an exact quantization condition and the threshold expansion is not performed there). Our present study aims at closing this gap. As already mentioned,
the results of this study will be certainly useful, for example,
in the analysis of lattice data on the three-pion decays, which are expected in the
near future.

The layout of the paper is as follows. In Sec. <ref> we write down
the nonrelativistic two- and three-body Lagrangians and perform the matching.
Further, in Sec. <ref> we construct the three-pion states with
different total isospin $I$ and calculate the potentials – the matrix elements
of the interaction Hamiltonian between the three-pion states. In Sec. <ref>, we derive the general expression for the perturbative energy shift up to and including $O(L^{-6})$ in the state with
the total isospin $I=1$. Finally, Sec. <ref> contains our conclusions.

§ NON-RELATIVISTIC EFFECTIVE LAGRANGIANS AND MATCHING

Our calculations closely follow the path of
Refs. [9, 49], where the identical particles have been considered.
Throughout this paper, we shall use the nonrelativistic effective field theory framework,
including relativistic effects perturbatively. We further assume isospin symmetry and
denote the common mass of pions by $M$. The Lagrangian, which will be used in the
calculations up to and including order $L^{-6}$, is given by

L=∑_i=±,0π_i^†(i∂_t-M+∇^2/2M+∇^4/8M^3)π^_i+L_2+L_3  .
Here, $\pi_\pm$ and $\pi^{}_0$ denote nonrelativistic pion fields, and ${\cal L}_2$ and ${\cal L}_3$ are the Lagrangians in the two- and three-particle sectors, respectively. Note that only leading relativistic
corrections in the kinetic term contribute up to the order $L^{-6}$.

Further, at the order we are working, it suffices to retain in ${\cal L}_2$ the terms up to
the two spatial
derivatives. Then, this Lagrangian takes the form

L_2    =     1/2  c_1 π_0^†π_0^†π^_0 π^_0
- 1/8  d_1( π_0^†π_0^†π^_0 π^_0+)

+    2c_2 ( π_+^†π_0^†π_
+ π^_0 + π_-^†π_0^†π^_- π^_0 )
- 1/2d_2 ( π_+^†π_0^†π^_+ π^_0
+ π_-^†π_0^†π^_- π^_0  +  )

+    c_3 ( π_+^†π_-^†π^_0 π^_0 + h.c. )
- 1/4 d_3 ( π_+^†π_-^†π^_0 π^_0 + π_0^†π_0^†π^_+ π^_- +  )

+    2c_4  π_+^†π_-^†π^_+ π^_-
- 1/2 d_4( π_+^†π_-^†π^_+ π^_-

+    1/2 c_5 ( π_+^†π_+^†π^_+ π^_+
+ π_-^†π_-^†π^_- π^_- )
- 1/8 d_5 ( π_+^†π_+^†π^_+ π^_+
+ π_-^†π_-^†π^_- π^_-  +)  .
Here, $\leftrightnabla = (\rightnabla - \leftnabla)/2$ is the Galilean-invariant derivative.
Note also that there are more terms with two derivatives, which describe the
center-of-mass motion of two-particle pairs. Since these terms do not contribute to the
energy shift of the ground state we are interested in [9],
we do not display them here explicitly. In addition, it can be checked that the P-wave
terms do not contribute to the ground-state energy at this order
(albeit they contribute to the excited states). These terms are omitted as well.

In a finite box, the Fourier transform of the field takes the form
π_i(x,t)=1/L^3 ∑_pe^-ip^0t+ipxa_i(p)  ,
where the creation and annihilation operators obey the commutation relations:
[a_i(p),a_j^†(q)]=L^3δ_ijδ_pq  .
The normalized one-particle states are given by
|π_i(p)⟩=1/L^3/2 a^†_i(p)|0⟩  .
If dimensional regularization is used, the couplings $c_i$ and $d_i$
can be related to the parameters of the effective-range expansion in the $\pi\pi$
scattering phase shift through the matching:
p δ_I(p) = -1/a_I + 1/2 r_I p^2 + ⋯,
I=0,2  .
The matching
condition takes the form
c_1    =    -4π/M 1/3(a_0 + 2a_2),            d_1    =    -4π/M 1/3(a_0^2 r̂^_0 + 2a_2^2 r̂^_2)  ,

c_2    =    -4π/M 1/2a_2,      d_2    =    -4π/M 1/2a_2^2 r̂^_2  ,

c_3    =    -4π/M 1/3(a_2 - a_0),      d_3    =    -4π/M 1/3(a_2^2 r̂^_2 - a_0^2 r̂^_0)  ,

c_4    =    -4π/M 1/6(2a_0 + a_2),      d_4    =    -4π/M 1/6(2a_0^2 r̂^_0 + a_2^2 r̂^_2)  ,

c_5    =    -4π/M a_2,      d_5    =    -4π/M a_2^2 r̂^_2   .

In the above equations, the quantities $\hat r^{}_I$ include relativistic corrections, which stem from the matching to the relativistic amplitude [9]:
r̂^_I=r^_I-1/a^_IM^2  .
Further, in the calculations up to and including $O(L^{-6})$, only the nonderivative
Lagrangians in the three-particle sector contribute. It is possible to construct two linearly
independent terms with this property and, hence, the Lagrangian can be written as
L_3 =
η_1 ( π_+^†π^_+ + π_0^†π^_0 + π_-^†π^_- )^3
+ η_2 ( π_+^†π^_+ + π_0^†π^_0 + π_-^†π^_- )
( 2 π_+^†π_-^†- π_0^†π_0^†)
( 2 π^_+ π^_- - π^_0π^_0 ) .
The low-energy couplings $\eta_1$ and $\eta_2$ are ultraviolet-divergent – these divergences
cancel the pertinent divergences, which arise in the perturbation theory at higher orders.
We shall see the examples of such cancellations below.
Note also that $\eta_1$ is a counterpart of the coupling $\eta_3$, introduced
in Ref. [49] for the case of three identical particles. Namely, one
has  $\eta_3=-6\eta_1$. On the contrary, the coupling $\eta_2$ does not emerge in the
sector with maximal isospin. In order to fix this coupling, one has to measure the
energy shift in the sector with $I=1$.

It is well known that the couplings
$\eta_1$ and $\eta_2$, via matching, can be traded for the “divergence free” threshold
amplitudes in the three-particle sector.
Such an amplitude can be defined in different ways. In the present context, the
definition, introduced in Ref. [9], is the most convenient,
and we stick to it. In order to match two constants, we consider two relativistic on-shell
amplitudes $3\pi^+\to 3\pi^+$ and $3\pi^0\to 3\pi^0$, denoted by ${\cal M}_+$ and
${\cal M}_0$, respectively. Next, we have to choose the
directions from which the threshold is approached. This can be done, e.g., in the
following way:
p_1=λe_y  ,    p_2=λ(√(3)/2 e_x-1/2 e_y)  ,    p_3=-λ(√(3)/2 e_x+1/2 e_y)  ,
and ${\bf p}_i'=-{\bf p}_i$ for $i=1,2,3$. Here, ${\bf p}_i,{\bf p}'_i$ are the initial
and final momenta of particles and ${\bf e}_{x,y,z}$ denote unit vectors in the direction
of pertinent axes.

The limit $\lambda\to 0$ in the above amplitudes is singular. In order to arrive at the
regular quantities, one has to subtract the one-particle irreducible piece (pole part) first.
An exact definition of this quantity can be found in Ref. [9].
Then, it can be proven that
=1/λ M_a^(-1)+lnλ/M M_a^(l)
+M_a^(0)+⋯  ,    a=+,0  .
Here, ellipses stand for the terms that vanish, as $\lambda\to 0$.

The couplings $\eta_1$ and $\eta_2$ can be related to the regular threshold amplitudes
${\cal M}_a^{(0)}$. This relation can be obtained exactly in the same way as
in Ref. [9], and we quote here only the final result:

36η_1^r   =   -288π^2/M a_2^3r̂_2-576π^2/M^3 a_2^2
+1/(2M)^3 M_+^(0)-192πa_2^4/M (3√(3)δ̅^(d)-4πδ̅^(e))  ,

36(η_1^r+η_2^r)   =   -32π^2 /M (a_0+2a_2)(a_0^2r̂_0+2a_2^2r̂_2)
+1/(2M)^3 M_0^(0)

-   64√(3)π/3M (a_0^4+4a_0^3a_2+9a_0^2a_2^2+2a_0a_2^3+11a_2^4)

+   128π^2/27M (2a_0^4+28a_0^3a_2+57a_0^2a_2^2+10a_0a_2^3+65a_2^4)δ̅^(e)  .
Here, $\eta_{1,2}^r$ denote the finite parts of the couplings $\eta_{1,2}$
(dimensional regularization with the minimal subtraction prescription has been used,
and $d$ denotes the number of space dimensions):
η_1   =   μ^2(d-3){-16πa_2^4/3M 1/d-3 (3√(3)-4π)+η_1^r}  ,

η_2   =   μ^2(d-3){
16π/243M 1/d-3 (a_0-a_2)[ 2π(2a_0^3 + 30 a_0^2 a_2 + 87 a_0 a_2^2 +97 a_2^3) - 9√(3)(a_0 + 2a_2)(a_0^2 + 3a_0a_2 + 8a_2^2) ]

+   η_2^r}  .
Furthermore, the renormalization scale was chosen at $\mu=M$ in the above

The quantities $\bar\delta^{(d)}$ and $\bar\delta^{(e)}$, which emerge in
Eq. (<ref>), are equal to
δ̅^(d,e)=δ^(d,e)-Γ'(1)-ln4π  ,    δ^(d)=-1.090…  ,    δ^(e)=3.926…  .
where $\delta^{(d)}$ and $\delta^{(e)}$ are the quantities defined in Appendix A of
Ref. [9].

To summarize, one can express two unknown couplings $\eta_1^r$ and $\eta_2^r$ through the
regular, real-valued threshold amplitudes ${\cal M}_+^{(0)}$ and ${\cal M}_0^{(0)}$, which
can be further extracted from the measured energy levels on the lattice. At this stage,
we do not see much advantage in
working in terms of the amplitudes rather than couplings,
since this amounts to the replacing of one set of
unknown constants by another one in the fit.
Note, however, that the amplitudes are observables,
whereas the couplings are not. For this reason, comparison with other approaches is
made easier, if the final result is rewritten in terms of the amplitudes.

§ THE FOCK STATE AND THE POTENTIALS

In a finite cubic box with size $L$, the Fock state of the theory consists of the vectors, containing freely moving particles $\pi_\pm$ and $\pi^{}_0$ with the discretized momenta ${\bf p}=2\pi{\bf n}/L$, where ${\bf n}\in\mathbb{Z}^3$. For instance, the closure relation in the three-particle sector
with a total charge +1 is written in the following form:

1   =   1/2! ∑_p_1,p_2,p_3

+   1/2! ∑_p_1,p_2,p_3
|π_0(p_1)π_0(p_2)π_+(p_3)⟩⟨π_0(p_1)π_0(p_2)π_+(p_3)|  .
Note that here the “odd” particle is always labeled by the index $3$. Instead, one could
sum up over three permutations of indices $(123)$, $(312)$, $(231)$, in each of the above sums,
dividing everything by $3!$ instead of $2!$.

Next, we wish to rewrite this sum as a sum over the three-pion states with a definite
total isospin $I$. It is a straightforward exercise to construct such states, using the table
of the Clebsh-Gordan coefficients. Below, we give maximal vectors $I_z=I$ in each
of the representation with $I=0,1,2,3$:
|3,3⟩   =    |π_+ π_+ π_+⟩  ,

|2,2⟩_A    =    1/√(2)
( |π_+π_0π_+⟩- |π_0π_+π_+⟩)  ,

|2,2⟩_S    =    1/√(6)
( 2|π_+π_+π_0⟩-  |π_0π_+π_+⟩-
|π_+π_0π_+⟩)  ,

|1,1⟩_A    =    1/2  ( |π_+π_0π_0⟩- |π_0π_+π_0⟩- |π_+π_-π_+⟩+ |π_-π_+π_+⟩)  ,

|1,1⟩_S_1    =    1/√(3)  ( |π_+π_-π_+⟩- |π_0π_0π_+⟩+ |π_-π_+π_+⟩)  ,

|1,1⟩_S_2    =     1/2√(15)
(  6|π_+π_+π_-⟩+  |π_+π_-π_+⟩+ |π_-π_+π_+⟩- 3|π_0π_+π_0⟩-3|π_+π_0π_0⟩+ 2|π_0π_0π_+⟩)  ,

|0,0⟩   =   1/√(6)  ( |π_+π_0π_-⟩- |π_+π_-π_0⟩- |π_0π_+π_-⟩+ |π_0π_-π_+⟩+ |π_-π_+π_0⟩- |π_-π_0π_+⟩)  .

Here, it is implicitly assumed that the particles in the ket-vectors carry momenta
${\bf p}_1,{\bf p}_2,{\bf p}_3$, respectively, and the subscripts $S/A$ stand for symmetric/antisymmetric wave functions with respect to $1\leftrightarrow 2$. Further,
the multiplicity of the irreducible representations, corresponding to $I=0,1,2,3$
is $1,3,2,1$, respectively, i.e., the unperturbed states in the sectors $I=1,2$ are,
generally, degenerate, and hence the perturbation theory for the degenerate states should be used in these sectors. Note also that this basis is the counterpart of the basis,
used in Ref. [25], Eqs. (2.56) and (2.57) (there, $I_z=0$
for all $I$ has been chosen).

Let us now focus on the unperturbed ground state with ${\bf p}_1={\bf p}_2={\bf p}_3=0$.
From Eq. (<ref>) it is clear that only the vectors
with $I=3$ and $I=1$ survive, whereas for $I=0,2$, the unperturbed ground-state
vector vanishes identically – these states will not be seen on the lattice.
Further, note that the solution for $I=3$ is
already known in the literature, whereas the expression for the shift in the $I=1$ state,
which will be given below, is new. Providing this expression,
our formalism that describes the ground state of three pions will become complete.

At the next step, we wish to rewrite the closure relation in terms of the states with
a definite total isospin. We shall work with the full basis of states with $I_z=+1$, where
the states with $I=1$ are also contained. In the Fock space, there are six basis vectors
with $I_z=+1$:
|v_1⟩   =   |π_+π_+π_-⟩  ,    |v_2⟩   =   |π_+π_-π_+⟩  ,    |v_3⟩   =   |π_-π_+π_+⟩  ,

|v_4⟩   =   |π_0π_0π_+⟩  ,    |v_5⟩   =   |π_0π_+π_0⟩  ,    |v_6⟩   =   |π_+π_0π_0⟩  .
We wish to express these vectors, as linear combinations of the vectors $|I,I_z\rangle$
with $I=1,2,3$ and $I_z=1$. The explicit expression of these vectors is given by:

|f_1⟩   =    1/√(15)  (  2|π_+π_+π_-⟩+2|π_+π_-π_+⟩+2|π_-π_+π_+⟩- |π_0π_0π_+⟩- |π_0π_+π_0⟩- |π_+π_0π_0⟩)  ,

|f_2⟩   =    1/2√(3)
(  2|π_+π_+π_-⟩-  |π_+π_-π_+⟩- |π_-π_+π_+⟩+ 2|π_0π_0π_+⟩- |π_0π_+π_0⟩-|π_+π_0π_0⟩)  ,

|f_3⟩   =    1/2  ( |π_+π_0π_0⟩- |π_0π_+π_0⟩- |π_+π_-π_+⟩+ |π_-π_+π_+⟩)  ,

|f_4⟩   =   1/2√(3) (
2|π_+π_+π_-⟩-|π_+π_-π_+⟩-|π_-π_+π_+⟩-2|π_0π_0π_+⟩+|π_0π_+π_0⟩+|π_+π_0π_0⟩)  ,

|f_5⟩   =   1/2 (
|π_+π_0π_0⟩-|π_0π_+π_0⟩+|π_+π_-π_+⟩-|π_-π_+π_+⟩)  ,

|f_6⟩   =   1/√(15) (
|π_+π_+π_-⟩+|π_+π_-π_+⟩+|π_-π_+π_+⟩+2|π_0π_0π_+⟩+2|π_0π_+π_0⟩+2|π_+π_0π_0⟩)  .
Note that the vectors $|f_1\rangle$ and $|f_2\rangle$ are linear combinations of the
vectors $|1,1\rangle_{S_1}$ and $|1,1\rangle_{S_2}$ (the same linear combinations, as
in Ref. [25], Eqs. (C.14) and (C.17), again for $I_z=0$).
The reason for such choice will become clear below.

It is now straightforwardly seen that the basis vectors of one basis can be expressed as linear
combinations of the basis vectors of other basis: $|v_a\rangle=C_{ab}|f_b\rangle$.
The orthogonal $6\times 6$ matrix $C$ is given by:
2/√(15)     1/√(3)     0     1/√(3)     0    1/√(15)

-2/√(15)     1/2√(3)    1/2      1/2√(3)    -1/2      -1/√(15)

-2/√(15)     1/2√(3)    -1/2      1/2√(3)    1/2      -1/√(15)

-1/√(15)     1/√(3)     0     -1/√(3)    0     2/√(15)

1/√(15)     1/2√(3)    1/2     -1/2√(3)    1/2    -2/√(15)

1/√(15)     1/2√(3)    -1/2     -1/2√(3)    -1/2    -2/√(15)
Using $CC^T=C^TC=1$, the closure relation (<ref>) can be finally rewritten in the isospin basis:
|f_a⟩⟨f_a|  .
The Lippmann-Schwinger equation, which can be obtained in the nonrelativistic effective theory, can be partially diagonalized in this basis. In order to derive this equation, we write down the canonical Hamiltonian, obtained from Eq. (<ref>), in the following form:
H=H_0+H_I=H_0+H_c+H_d+H_r+H_η  .
Here ${\bf H}_0$ is the free nonrelativistic Hamiltonian. The states, which were constructed above, are the eigenstates of ${\bf H}_0$:
H_0|π_i(p_1)π_j(p_2)π_k(p_3)⟩   =   1/2M (p_1^2+p_2^2+p_3^2)|π_i(p_1)π_j(p_2)π_k(p_3)⟩

≐    E_p|π_i(p_1)π_j(p_2)π_k(p_3)⟩  ,    i,j,k=±,0  .
Further, various terms in the interaction Lagrangian are given by
H_c   =   ∫d^3x{
-1/2  c_1 π_0^†π_0^†π^_0 π^_0
-2c_2 ( π_+^†π_0^†π^_+ π^_0
+ π_-^†π_0^†π^_- π^_0 )
- c_3 ( π_+^†π_-^†π^_0 π^_0 + h.c. )
- 2c_4  π_+^†π_-^†π^_+ π^_-

-   1/2 c_5 ( π_+^†π_+^†π^_+ π^_+
+ π_-^†π_-^†π^_- π^_- )}  ,

H_d   =   ∫d^3x{
1/8  d_1( π_0^†π_0^†π^_0 π^_0+)
+ 1/2d_2 ( π_+^†π_0^†π^_+ π^_0
+ π_-^†π_0^†π^_- π^_0  +  )

+    1/4 d_3 ( π_+^†π_-^†π_0 π^_0 + π_0^†π_0^†π^_+ π^_- +  )
+ 1/2 d_4( π_+^†π_-^†π^_+ π^_-

+    1/8 d_5 ( π_+^†π_+^†π^_+ π^_+
+ π_-^†π_-^†π^_- π^_-  +)}  ,

H_r   =   -1/8M^3∑_i=±,0
∫d^3x{π_i^†∇^4π^_i}  ,

H_η   =   ∫d^3x{
-η_1 ( π_+^†π^_+ + π_0^†π^_0 + π_-^†π^_- )^3
- η_2 ( π_+^†π^_+ + π_0^†π^_0 + π_-^†π^_- )
( 2 π_+^†π_-^†- π_0^†π_0^†)
( 2 π^_+ π^_- - π^_0π^_0 )} .

The integration in the above expressions is carried out in a cubic box with a size $L$.

In the operator form, the Lippmann-Schwinger equation can be written as
T(E)=H_I+H_IG_0(E)T(E)  ,    G_0(E)=1/E-H_0  .
Let us now sandwich this equation by the three-particle states with the definite total isospin
$I$, use closure relations and take into account the fact that ${\bf T}$, ${\bf H}_I$ and ${\bf G}_0$ are diagonal in $I$. In particular, we are primarily interested in the sector $I=1$,
where the equation simplifies to

⟨k|T_γβ(E)|q⟩  ,
where, the symbol $p$ denotes a set of momenta $({\bf p}_1,{\bf p}_2,{\bf p}_3)$, and so on. Further, the matrix elements are defined as
⟨f_α,p|H_I|f_β,q⟩  ,    ⟨p|T_αβ(E)|q⟩=1/3!
⟨f_α,p|T(E)|f_β,q⟩  .
Here, in difference with the indices $a$ and $b$, the indices $\alpha,\beta,\gamma$
run from 1 to 3, picking out the vectors $f_1,f_2,f_3$ from the sector $I=1$ only. Thus, the partial
of the Lippmann-Schwinger equation is achieved in the basis of the isospin
vectors (<ref>).

The matrix elements of the interaction Hamiltonian in the above equation can be calculated
straightforwardly. Below, we merely list different contributions to this matrix element
in the sector with $I=1$:

⟨f_1,p|H_c|f_1,q⟩   =    4π/M L^3 1/27 (5a_0 +4a_2)[ V_1 + V_2 + V_3 + V^+_1 + V^+_2 +  V^+_3]  ,

⟨f_1,p|H_c|f_2,q⟩   =    4π/M L^3 √(5)/54 (a_0 - a_2)[ 2(V_1 + V_2 -2V_3) - 3(V^-_1+V^-_2) - (V^+_1+V^+_2-2V^+_3) ]  ,

⟨f_2,p|H_c|f_2,q⟩   =    4π/M L^3 1/108 (4a_0 + 5a_2)[ (V_1 + V_2 + 4V_3) -(2V^+_1+2V^+_2-V^+_3) ]  ,

⟨f_1,p|H_c|f_3,q⟩   =    -4π/M L^3 √(5)/18√(3) (a_0 - a_2)[ 2(V_1-V_2) + (V^-_1-V^-_2-2V^-_3) - (V^+_1 - V^+_2)]  ,

⟨f_2,p|H_c|f_3,q⟩   =    -4π/M L^3 1/36√(3) (4a_0 + 5a_2)[ (V_1 - V_2) -(V^-_1 - V^-_2+V^-_3) + (V^+_1-V^+_2)]   ,

⟨f_3,p|H_c|f_3,q⟩   =    4π/M L^3 1/36 (4a_0+5a_2)
(V_1+V_2-V^+_3)  ,

V_1 = δ_𝐩_1,𝐤_1δ_𝐩_2+𝐩_3,𝐤_2+𝐤_3,     V_2 = δ_𝐩_2,𝐤_2δ_𝐩_1+𝐩_3,𝐤_1+𝐤_3,     V_3 = δ_𝐩_3,𝐤_3δ_𝐩_1+𝐩_2,𝐤_1+𝐤_2,

V^±_1    =    δ_𝐩_2,𝐤_3 δ_𝐩_1+𝐩_3,𝐤_1+𝐤_2 ±δ_𝐩_3,𝐤_2 δ_𝐩_1+𝐩_2,𝐤_1+𝐤_3  ,

V^±_2   =    δ_𝐩_1,𝐤_3 δ_𝐩_2+𝐩_3,𝐤_1+𝐤_2 ±δ_𝐩_3,𝐤_1 δ_𝐩_1+𝐩_2,𝐤_2+𝐤_3  ,

V^±_3    =    δ_𝐩_1,𝐤_2 δ_𝐩_2+𝐩_3,𝐤_1+𝐤_3 ±δ_𝐩_2,𝐤_1 δ_𝐩_1+𝐩_3,𝐤_2+𝐤_3  .
Further, the matrix elements of the operator ${\bf H}_d$ are also given  by
Eq. (<ref>), with the replacement
$a_I^{}\to a_I^{2}\hat r_I$ and $V^{}_i,V_i^\pm\to \hat V^{}_i,\hat V_i^\pm$. The latter are again given by
Eqs. (<ref>) and (<ref>), with the following replacements done:
δ_𝐩_i,𝐤_lδ_𝐩_j+𝐩_k,𝐤_m+𝐤_n →1/16  δ_𝐩_i,𝐤_lδ_𝐩_j+𝐩_k,𝐤_m+𝐤_n [ (𝐩_j - 𝐩_k )^2 + (𝐤_m - 𝐤_n )^2 ].
Next, the relativistic correction in the kinetic term has
the following matrix elements:
⟨f_1,p|H_r|f_1,q⟩   =    -1/96 M^3
∑_i=1^3 (𝐤^4_i + 𝐩^4_i) ( V_123 + V_132 + V_231  + V_213 + V_312 + V_321 )  ,

⟨f_1,p|H_r|f_2,q⟩   =    0  ,

⟨f_2,p|H_r|f_2,q⟩   =    -1/192M^3 ∑_i=1^3( 𝐤^4_i + 𝐩^4_i) ( 2 V_123 - V_132 + 2V_213  - V_231 - V_312 - V_321 )  ,

⟨f_1,p|H_r|f_3,q⟩   =    0  ,

⟨f_2,p|H_r|f_3,q⟩   =    1/64√(3)M^3 ∑_i=1^3( 𝐤^4_i + 𝐩^4_i) ( V_132  - V_231 + V_312 - V_321 )  ,

⟨f_3,p|H_r|f_3,q⟩   =    -1/192 M^3 ∑_i=1^3( 𝐤^4_i + 𝐩^4_i) ( 2 V_123 - 2V_213 + V_132  - V_231 - V_312 + V_321 )  ,
V_ijk = δ_𝐩_1,𝐤_i δ_𝐩_2,𝐤_j δ_𝐩_3,𝐤_k  .
Finally, the single non-zero matrix element of the three-particle Hamiltonian is given by:
⟨f_1,p|H_η|f_1,q⟩   =
-2/L^6  (3η_1+5η_2)  δ_𝐩_1+𝐩_2+𝐩_3,𝐤_1+𝐤_2+𝐤_3  .

§ PERTURBATION THEORY

The energy level shift in the nonrelativistic effective field theory can be calculated by using the
ordinary Rayleigh-Schrödinger perturbation theory. A small technical issue might arise here since,
as we have seen, there are multiple states with $I=1$, so, the application of the perturbation
theory for the degenerate states might be needed. This is, however, not the case for the
ground state – as seen from Eq. (<ref>), only the state $f_1$ has a nonvanishing
ground-state wave function. Note that this is different in the excited states – in this case,
the perturbation theory for the degenerate states should indeed be used.

In order to derive the perturbative expansion for the ground-state energy,
at the first step we single out the contribution from the
unperturbed ground-state level ${\bf k}_1={\bf k}_2={\bf k}_3=0$.
It can be seen that Eq. (<ref>) is equivalent to the
following system of linear equations:

⟨p|T_αβ(E)|q⟩   =   ⟨p|Ω_αβ(E)|q⟩+∑_γ=1^3⟨p|Ω_αγ(E)|0⟩1/E-E_0 ⟨0|T_γβ(E)|q⟩  ,

⟨p|Ω_αβ(E)|q⟩   =   ⟨p|V_αβ|q⟩+∑_k≠0∑_γ=1^3⟨p|V_αγ|k⟩1/E-E_k
⟨k|Ω_γβ(E)|q⟩  .
The shifted energy is given by the position of the pole of the scattering matrix.
Setting now external momenta to zero, $p=q=0$,
it is seen using Eq. (<ref>) that the shifted ground-state pole position obeys
the secular equation,

E-E_0-Ω̂(E)=0  ,
where $\hat\Omega(E) =\Omega_{11}(E)$ and $\Omega(E)$ is a $3\times 3$ matrix, whose components coincide
with $\langle 0|\boldsymbol{\Omega}_{\alpha\beta}(E)|0\rangle$ (we remind the reader that
all components of this matrix, except $\alpha=\beta=1$, vanish at $p=q=0$).

At the next step, we recall that the quantity
$\hat\Omega(E)$ does not contain the singular denominator $1/(E-E_0)$, which is excluded
in the sum with $k\neq 0$.
Consequently, one may expand $\hat\Omega(E)$ in Taylor series in the vicinity of the unperturbed ground state,
Ω̂(E)=Ω̂(E_0)+(E-E_0)Ω̂'(E_0)+1/2 (E-E_0)^2Ω̂”(E_0)
+O((E-E_0)^3)  .
The coefficients of this expansion can be further expanded in the inverse powers of $L$:
Ω̂(E_0)    =    Ω̂_3(E_0)/L^3 + Ω̂_4(E_0)/L^4 + Ω̂_5(E_0)/L^5 + Ω̂_6(E_0)/L^6 + O(L^-7)  ,

Ω̂'(E_0)    =    Ω̂'_2(E_0)/L^2 + Ω̂'_3(E_0)/L^3 + O(L^-4)  ,

Ω̂”(E_0)    =    Ω̂”_0(E_0) +O(L^-1)  .
Since $E-E_0=O(L^{-3})$, the terms, displayed above, are sufficient to calculate the energy
shift up to and including order $L^{-6}$. Using iterations, one straightforwardly obtains
E-E_0   =   E_3/L^3+E_4/L^4+E_5/L^5+E_6/L^6+O(L^-7)  ,

E_3   =   Ω̂_3(E_0)  ,

E_4   =   Ω̂_4(E_0)  ,

E_5   =   Ω̂_5(E_0)+Ω̂_2'(E_0)Ω̂_3(E_0)  ,

E_6   =   Ω̂_6(E_0)+Ω̂_3'(E_0)Ω̂_3(E_0)+Ω̂_2'(E_0)Ω̂_4(E_0)+1/2 Ω̂_0”(E_0)Ω̂_3(E_0)^2  .
The matrix $\Omega(E)$ can be expanded similarly.
The coefficients $\Omega_i(E_0),\Omega_i'(E_0),\Omega_i''(E_0)$ can be calculated by iterating
Eq. (<ref>). Defining a $3\times 3$ matrix $V_{pq}$, whose
elements are given by $\langle p|{\bf V}_{\alpha\beta}|q\rangle$, and attaching
superscripts “$c$,” “$d$,” “$r$,” “$\eta$,” to the individual contributions,
one gets a compact expression:

Ω_3(E_0)/L^3    =    V^c_00,

Ω_4(E_0)/L^4    =    -∑_p ≠0 V^c_0pV^c_p0/E_p,

Ω_5(E_0)/L^5    =    ∑_p,k ≠0 V^c_0pV^c_pkV^c_k0/E_p E_k,

Ω_6(E_0)/L^6    =    -∑_p,k,q ≠0 V^c_0pV^c_pqV^c_qkV^c_k0/E_p E_q E_k -∑_p ≠0 V^c_0pV^d_p0+V^d_0pV^c_p0/E_p + V^η_00 + ∑_p,k ≠0 V^c_0pV^r_pkV^c_k0/E_p E_k,

Ω'_2(E_0)/L^2    =    -∑_p ≠0 V^c_0pV^c_p0/E_p^2,

Ω'_3(E_0)/L^3    =    2∑_p,k ≠0 V^c_0pV^c_pkV^c_k0/E_p^2 E_k ,

Ω”_0(E_0)    =    2∑_p ≠0 V^c_0pV^c_p0/E_p^3.
Taking the component $\alpha=\beta=1$ from the above expressions, we finally arrive at the expressions for
$\hat\Omega_i(E_0)$, $\hat\Omega'_i(E_0)$, $\hat\Omega''_i(E_0)$, we are looking for.

Putting everything together, we obtain the final expression for the energy shift:

E_3   =   4 π/3 M (5a_0 + 4a_2),

E_4   =   -4/3 M (5a_0^2 + 4a_2^2) I,

E_5   =    4/3 M π [ (5a_0^3 + 4a_2^3) I^2
- 1/9 (55a_0^3 -120 a_0^2 a_2 -60 a_0 a_2^2 + 44 a_2^3) J  ]  ,

E_6   =   -4/3 M π^2 [ (5a_0^4 +4a_2^4) I^3 - 1/9 (185a_0^4 -160 a_0^3 a_2 - 80 a_0 a_2^3 + 136 a_2^4) IJ

-    1/27 (5a_0 + 4a_2) (115a_0^3 + 420 a_0^2 a_2 + 300 a_0 a_2^2 + 56 a_2^3) K

+    8/27  (10a_0^4 + 140 a_0^3 a_2 + 285 a_0^2 a_2^2 +  50 a_0 a_2^3 + a_2^4) Q^r

+    8/9  (5a_0^4 + 20 a_0^3 a_2 + 45 a_0^2 a_2^2 +  10 a_0 a_2^3 + a_2^4) R^r ]

+   8π^2/3M (5a_0^3 r_0 + 4a_2^3 r_2) - 4π^2/3M^3 (5a_0^2 + 4 a_2^2)

+    64 π/243 M[ 9√(3) (5a_0^4 + 20 a_0^3 a_2 + 45 a_0^2 a_2^2 +  10 a_0 a_2^3 + a_2^4)

-    2π(10a_0^4 + 140 a_0^3 a_2 + 285 a_0^2 a_2^2 +  50 a_0 a_2^3 + a_2^4) ]ln(μL)

-   6 η_1^r - 10 η_2^r   .
Finally, we give a list of all momentum sums that appear in the above equation. These
are the same quantities that enter the expression in case of the maximal isospin.
The quantities $I,J,K$ are finite in the dimensional regularization and
are given by
I   =   ∑_n≠01/n^2=-8.91363291781⋯  ,

J   =   ∑_n≠01/n^4=16.532315959⋯  ,

K   =   ∑_n≠01/n^6=8.401923974433⋯  .
The quantities $Q^r$ and $R^r$ represent the finite part of the double sums over momenta in
the minimal subtraction renormalization scheme. Introducing a formal notation
of a sum in $d$ dimensions,[The sums in $d$ dimensions can be interpreted as follows. One first performs the Poisson transform and singles out the infinite-volume part, which corresponds to an integral instead of a sum. The dimensional regularization and minimal subtraction is applied then to this integral. The remainder is ultraviolet-convergent and the limit $d\to 3$ can be performed there without a problem.] we get
=μ^2(d-3){1/48π^2 (ln(μL)-1/2(d-3)
+1/(2π)^6 Q^r}  ,

1/L^2d ∑_p≠01/p^4
=μ^2(d-3){-√(3)/32π^3 (ln(μL)-1/2(d-3)
+1/(2π)^6 R^r}  .
The finite parts are given by
Q^r   =   -102.1556055⋯  ,

R^r   =   19.186903⋯  .
Note that, in the above equations, ${\bf p}=2\pi{\bf n}/L$, ${\bf q}=2\pi{\bf m}/L$,
and the summation is carried out over ${\bf n},{\bf m}\in\mathbb{Z}^3$.

Equation (<ref>) represents our final result. It gives the expression for the
energy shift of the ground state with $I=1$ up to and including $O(L^{-6})$.
Together with the known expression for the
case $I=3$, it provides a full framework for analyzing the lattice data on ground states
in the three-pion system. Two independent nonderivative couplings $\eta_{1,2}$ can be
separated and extracted from data in a result of this analysis.

§ CONCLUSIONS

Using the nonrelativistic effective Lagrangian approach, we have evaluated the finite-volume
shift of the ground state of three pions in the state with total isospin $I=1$ up to and including $O(L^{-6})$. Since three pions
with zero momenta can have only $I=1$ or $I=3$, the expression, which is derived
in the present paper, provides the last remaining missing piece in the description of the three-pion
ground state, because the expression in the case of $I=3$ has been known in the
literature. With this result at hand, one can also study the three-pion decays into
all isospin channels.

An immediate application of the obtained expression could be the extraction of two
independent three-body couplings $\eta_{1,2}$ from the lattice data. Moreover, the
analysis of the three-particle states might put constraints on the parameters in the
two-particle sector – the scattering lengths $a_I$ and the effective radii $r_I$. According
to our past experience in the scalar $\varphi^4$-theory [8, 9],
these constraints are most effectively implemented when one performs a
simultaneous fit of
all datasets in the two- and three-particle sectors. The input from the three-particle sector might
be interesting, because the two-particle scattering parameters (especially those corresponding to the total isospin $I=0$) are not very well determined on the lattice at present
(for the recent developments in the two-particle sector, we refer the reader, e.g.,
to [52, 53, 54, 55]).

Albeit we are working in the nonrelativistic effective theory, our results contain a full
set of relativistic corrections up to and including $O(L^{-6})$. In general, such corrections
can be systematically included in the nonrelativistic approach up to a desired accuracy.

The above work is a rather straightforward generalization of the result obtained for
the maximal isospin. One has, however, to deal with subtle issues, for example, with the degeneracy of levels, which are absent in the former case. The level splitting does not occur
for the ground state, but one will have to face it, deriving expressions for the
shifts of the excited states. A framework for the degenerate states can be worked
out along the path
which is pretty similar to the one considered in the present paper.
The pertinent explicit formulae are, however,
not  needed for the study of the ground state, and are therefore not displayed.

§ ACKNOWLEDGMENTS

The authors would like to thank U.-G. Meißner, F. Romero-López and C. Urbach
for interesting discussions.
The work of F. M. and  A. R.  was  supported in part by the DFG (CRC 110 “Symmetries
and the Emergence of Structure in QCD”, Grant no. TRR110). A. R., in addition,
thanks Volkswagenstiftung
(Grant no. 93562) and the Chinese Academy of Sciences (CAS) President's
International Fellowship Initiative (PIFI) (Grant no. 2021VMB0007) for the partial
financial support.

[1]
S. R. Beane, W. Detmold, T. C. Luu, K. Orginos, M. J. Savage and A. Torok,

Phys. Rev. Lett. 100 (2008) 082004

[arXiv:0710.1827 [hep-lat]].

[2]
B. Hörz and A. Hanlon,

Phys. Rev. Lett. 123 (2019) 142002

[arXiv:1905.04277 [hep-lat]].

[3]
C. Culver, M. Mai, R. Brett, A. Alexandru and M. Döring,

Phys. Rev. D 101 (2020) 114507

[arXiv:1911.09047 [hep-lat]].

[4]
M. Fischer, B. Kostrzewa, L. Liu, F. Romero-López, M. Ueding and C. Urbach,

[arXiv:2008.03035 [hep-lat]].

[5]
T. D. Blanton, F. Romero-López and S. R. Sharpe,

Phys. Rev. Lett. 124 (2020) 032001

[arXiv:1909.02973 [hep-lat]].

[6]
M. T. Hansen, R. A. Briceño, R. G. Edwards, C. E. Thomas and D. J. Wilson,

Phys. Rev. Lett. 126 (2021) 012001
[arXiv:2009.04931 [hep-lat]].

[7]
A. Alexandru, R. Brett, C. Culver, M. Döring, D. Guo, F. X. Lee and M. Mai,

Phys. Rev. D 102 (2020) 114523
[arXiv:2009.12358 [hep-lat]].

[8]
F. Romero-López, A. Rusetsky and C. Urbach,

Eur. Phys. J. C 78 (2018) 846

[arXiv:1806.02367 [hep-lat]].

[9]
F. Romero-López, A. Rusetsky, N. Schlage and C. Urbach,

JHEP 60 (2021) 60
[arXiv:2010.11715 [hep-lat]].

[10]
M. T. Hansen and S. R. Sharpe,

Phys. Rev. D 90 (2014) 116003 [arXiv:1408.5933 [hep-lat]].

[11]
M. T. Hansen and S. R. Sharpe,

Phys. Rev. D 92 (2015) 114509 [arXiv:1504.04248 [hep-lat]].

[12]
K. Polejaeva and A. Rusetsky,

Eur. Phys. J. A 48 (2012) 67 [arXiv:1203.1241 [hep-lat]].

[13]
H. W. Hammer, J. Y. Pang and A. Rusetsky,

JHEP 09 (2017) 109

[arXiv:1706.07700 [hep-lat]].

[14]
H. W. Hammer, J. Y. Pang and A. Rusetsky,

JHEP 10 (2017) 115

[arXiv:1707.02176 [hep-lat]].

[15]
M. Mai and M. Döring,

Eur. Phys. J. A 53 (2017)  240

[arXiv:1709.08222 [hep-lat]].

[16]
M. T. Hansen and S. R. Sharpe,

Phys. Rev. D 93 (2016) 014506 [arXiv:1509.07929 [hep-lat]].

[17]
M. T. Hansen and S. R. Sharpe,

Phys. Rev. D 93 (2016) 096006 [arXiv:1602.00324 [hep-lat]].

[18]
R. A. Briceño, M. T. Hansen and S. R. Sharpe,

Phys. Rev. D 95 (2017)  074510

[arXiv:1701.07465 [hep-lat]].

[19]
S. R. Sharpe,

Phys. Rev. D 96 (2017)  054515

[arXiv:1707.04279 [hep-lat]].

[20]
R. A. Briceño, M. T. Hansen and S. R. Sharpe,

Phys. Rev. D 98 (2018) 014506

[arXiv:1803.04169 [hep-lat]].

[21]
R. A. Briceño, M. T. Hansen and S. R. Sharpe,

Phys. Rev. D 99 (2019) 014516

[arXiv:1810.01429 [hep-lat]].

[22]
T. D. Blanton, F. Romero-López and S. R. Sharpe,

JHEP 03 (2019) 106

[arXiv:1901.07095 [hep-lat]].

[23]
R. A. Briceño, M. T. Hansen, S. R. Sharpe and A. P. Szczepaniak,

Phys. Rev. D 100 (2019) 054508

[arXiv:1905.11188 [hep-lat]].

[24]
F. Romero-López, S. R. Sharpe, T. D. Blanton, R. A. Briceño and M. T. Hansen,

JHEP 10 (2019) 007

[arXiv:1908.02411 [hep-lat]].

[25]
M. T. Hansen, F. Romero-López and S. R. Sharpe,

JHEP 07 (2020) 047

[arXiv:2003.10974 [hep-lat]].

[26]
U.-G. Meißner, G.-Rios and A. Rusetsky,

Phys. Rev. Lett.  114 (2015) 091602,
Erratum: [Phys. Rev. Lett.  117 (2016) 069902]
[arXiv:1412.4969 [hep-lat]].

[27]
T. D. Blanton and S. R. Sharpe,

Phys. Rev. D 102 (2020) 054515

[arXiv:2007.16190 [hep-lat]].

[28]
T. D. Blanton and S. R. Sharpe,

Phys. Rev. D 102 (2020)  054520

[arXiv:2007.16188 [hep-lat]].

[29]
M. Döring, H. W. Hammer, M. Mai, J. Y. Pang, A. Rusetsky and J. Wu,

Phys. Rev. D 97 (2018) 114508

[arXiv:1802.03362 [hep-lat]].

[30]
J. Y. Pang, J. J. Wu, H. W. Hammer, U.-G. Meißner and A. Rusetsky,

Phys. Rev. D 99 (2019)  074513

[arXiv:1902.01111 [hep-lat]].

[31]
J. Y. Pang, J. J. Wu and L. S. Geng,

Phys. Rev. D 102 (2020) 114515
[arXiv:2008.13014 [hep-lat]].

[32]
A. W. Jackura, S. M. Dawid, C. Fernández-Ramírez, V. Mathieu, M. Mikhasenko, A. Pilloni, S. R. Sharpe and A. P. Szczepaniak,

Phys. Rev. D 100 (2019)  034508

[arXiv:1905.12007 [hep-ph]].

[33]
Y. Meng, C. Liu, U.-G. Meißner and A. Rusetsky,

Phys. Rev. D 98 (2018) 014508

[arXiv:1712.08464 [hep-lat]].

[34]
M. Mai and M. Döring,

Phys. Rev. Lett. 122 (2019) 062503

[arXiv:1807.04746 [hep-lat]].

[35]
M. Mai, M. Döring, C. Culver and A. Alexandru,

Phys. Rev. D 101 (2020)  054510

[arXiv:1909.05749 [hep-lat]].

[36]
S. König and D. Lee,

Phys. Lett. B 779 (2018) 9

[arXiv:1701.00279 [hep-lat]].

[37]
P. Klos, S. König, H. W. Hammer, J. E. Lynn and A. Schwenk,

Phys. Rev. C 98 (2018)  034004

[arXiv:1805.02029 [nucl-th]].

[38]
S. König,

Few Body Syst. 61 (2020)  20

[arXiv:2005.01478 [hep-lat]].

[39]
T. D. Blanton and S. R. Sharpe,

[arXiv:2011.05520 [hep-lat]].

[40]
S. Kreuzer and H.-W. Hammer,

Phys. Lett. B 694 (2011) 424 [arXiv:1008.4499 [hep-lat]].

[41]
S. Kreuzer and H.-W. Hammer,

Eur. Phys. J. A 43 (2010) 229 [arXiv:0910.2191 [nucl-th]].

[42]
S. Kreuzer and H.-W. Hammer,

Phys. Lett. B 673 (2009) 260 [arXiv:0811.0159 [nucl-th]].

[43]
S. Kreuzer and H.-W. Grießhammer,

Eur. Phys. J. A 48 (2012) 93 [arXiv:1205.0277 [nucl-th]].

[44]
M. T. Hansen and S. R. Sharpe,

Ann. Rev. Nucl. Part. Sci. 69 (2019) 65

[arXiv:1901.00483 [hep-lat]].

[45]
T. D. Lee, K. Huang and C. N. Yang,

Phys. Rev. 106 (1957) 1135.

[46]
K. Huang and C. N. Yang,

Phys. Rev.  105 (1957) 767.

[47]
T. T. Wu,

Phys. Rev. 115 (1959) 1390.

[48]
S. Tan,

Phys. Rev. A 78 (2008) 013636

[arXiv:0709.2530 [cond-mat.stat-mech]].

[49]
S. R. Beane, W. Detmold and M. J. Savage,

Phys. Rev. D 76 (2007) 074507

[arXiv:0707.1670 [hep-lat]].

[50]
W. Detmold and M. J. Savage,

Phys. Rev. D 77 (2008) 057502

[arXiv:0801.0763 [hep-lat]].

[51]
S. R. Beane, W. Detmold, R. Horsley, M. Illa, M. Jafry, D. J. Murphy, Y. Nakamura, H. Perlt, P. E. L. Rakow and G. Schierholz, et al.

[arXiv:2003.12130 [hep-lat]].

[52]
T. Yagi, S. Hashimoto, O. Morimatsu and M. Ohtani,

[arXiv:1108.2970 [hep-lat]].

[53]
Z. Fu,

Phys. Rev. D 87 (2013)  074501

[arXiv:1303.0517 [hep-lat]].

[54]
C. Helmes et al. [ETM],

JHEP 09 (2015) 109

[arXiv:1506.00408 [hep-lat]].

[55]
L. Liu, S. Bacchio, P. Dimopoulos, J. Finkenrath, R. Frezzotti, C. Helmes, C. Jost, B. Knippschild, B. Kostrzewa and H. Liu, et al.

Phys. Rev. D 96 (2017)  054516

[arXiv:1612.02061 [hep-lat]].


# Parallelized Instantaneous Velocity and Heading Estimation of Objects using
Single Imaging Radar

Nihal Singh1, Dibakar Sil2, and Ankit Sharma2  Steradian Semiconductors
Email<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

The development of high-resolution imaging radars introduce a plethora of
useful applications, particularly in the automotive sector. With increasing
attention on active transport safety and autonomous driving, these imaging
radars are set to form the core of an autonomous engine. One of the most
important tasks of such high-resolution radars is to estimate the
instantaneous velocities and heading angles of the detected objects (vehicles,
pedestrians, etc.). Feasible estimation methods should be fast enough in real-
time scenarios, bias-free and robust against micro-Dopplers, noise and other
systemic variations. This work proposes a parallel-computing scheme that
achieves a real-time and accurate implementation of vector velocity
determination using frequency modulated continuous wave (FMCW) radars. The
proposed scheme is tested against traffic data collected using an FMCW radar
at a center frequency of 78.6 GHz and a bandwidth of 4 GHz. Experiments show
that the parallel algorithm presented performs much faster than its
conventional counterparts without any loss in precision.

###### Index Terms:

Velocity estimation, heading angle, parallel-computing, Doppler radar, radar
signal processing, velocity profile, GPU kernel.

##  I Introduction

High-resolution imaging radars are touted as the eyes of an autonomous engine
by major automotive companies and organizations around the globe [1]. While
these radars are typically proposed for geoscience applications[2], traffic
monitoring [3], and military surveillance [4], increasing attention on active
transport safety has shifted the focus towards advanced driver assistance
systems (ADAS) and autonomous driving.

Frequency-Modulated Continuous Wave (FMCW) imaging radars carry out sensing
and monitoring of the environment by determining various attributes of
surrounding objects. This process of object detection entails measuring
differences in frequency and/or phase between the transmitted signals, and
received echoes. Initial physical layer radar data processing results in a
point cloud (PC) (refer Fig. 2) with each point including attributes such as
range, velocity, azimuth angle, elevation angle, etc., with respect to the
observing radar’s axes of reference.

Figure 1: Algorithm Chain for Radar Processing Stack up until Heading Angle
Estimation

Within the PC, an object appears as a dense collection of points, also termed
as a cluster. As represented in Fig. 1, the step after PC generation is the
clustering feature vector extraction. This process comprises selection of
attributes to construct a feature vector, to base the clustering process on.
Typically, the feature vector includes the Cartesian counterparts of the
detected ranges but could include other attributes of the detected points as
well. Subsequently, clustering involves grouping of points similar to each
other in the feature vector sense into respective clusters (refer Fig. 3).
This is usually carried out by employing spatial clustering algorithms such as
Nearest Neighbour [5], DBSCAN [6], Grid-DBSCAN [7] etc. After segregating the
PC into clusters, it is desired to determine the velocity profile and heading
of each such cluster.

Recent literature shows that separation of closely spaced clusters by standard
density-based methods gives unreliable categorizations of certain points.
Furthermore, these algorithms are not noise robust and require further
developments for accurate clustering in applications utilizing real radar data
[8]. Oscillatory motion of sub-parts of an object induces modulations in the
primary reflected signal, termed as micro-Doppler. A common example of this is
the signal produced by rotating wheels of vehicles in motion [9]. Accurate
estimation of the instantaneous velocity and heading angle is thus hindered by
the influence of these effects.

Figure 2: [Left] Camera Image depicting a subset of the captured Radar Frame’s
Field of View, [Right] Visualization of generated 3D Point Cloud (78.6 GHz
center freq., 4GHz BW)

Figure 3: [Left] 2D Clustering output along Longitudinal and Lateral axes.
Cluster points corresponding to Truck is represented with ( ), cluster
corresponding to motorcycle is respresented with ( ), and noise points are
shown with ( ), [Right] Result of inlier detection on cluster points
corresponding to the truck. (78.6 GHz center freq., 4GHz BW)

Typical velocity profile estimation by radar is achieved by using more than
one sensor [10] or restricting the detection to only linear motion [11].
Another commonly used approach is tracking [12], though tracking solutions
take multiple iterations of range and projected velocity determination to
converge to the object’s correct vector velocity, and therefore are not
instantaneous. In [13], Lidman _et al_. have proposed an algorithm to estimate
the instantaneous velocity by extracting the shape of the cluster. However,
this approach is only feasible when measurements from multiple Doppler radars
are combined for accurate shape formation.

This paper proposes and evaluates instantaneous velocity profile and heading
estimation of objects in real scenarios using a single Multiple-Input and
Multiple-Output (MIMO) FMCW imaging radar. The main contributions of this work
are the real-time implementations of:

  * •

Inlier determination algorithm for the purpose of real-time filtering of
random and systemic data variations from the radar detected Doppler velocity
and azimuth angle.

  * •

Least squares approximation method to simultaneously compute instantaneous
vector velocities of multiple objects.

Section $2$ describes the problem statement and motivates the need for a
better real-time solution to determine instantaneous vector velocity. Section
$3$ explains the proposed parallelized solution. Section $4$ highlights the
results obtained by employing the proposed method and Section $5$ provides the
concluding remarks.

##  II Problem Formulation

For practical radar applications in the automotive sector, accurate and real-
time determination of detected objects’ attributes is essential. Work
presented in this paper uses a single MIMO imaging radar module from Steradian
Semiconductors [14] to capture radar data in common automotive scenarios.

Range, Doppler speed and angle of arrival (AoA) from the said radar are
calculated using frames of data transmitted and echoed back to the radar. Each
frame is composed of a set of frequency modulated chirps, where the inter-
chirp time is $T_{c}$, and $N_{c}$ such chirps constitute a single frame.
Pulse compression with the received echo produces an intermediate frequency
(IF) signal that gives range information within the digital baseband samples
of a chirp, Doppler information across chirps’ samples and AoA information
from samples across antenna elements [15].

With the starting frequency of a chirp signal as $f_{0}$, and slope $S_{u}$,
instantaneous frequency of the chirp signal can be written as

$f(t)=f_{0}+S_{u}t,\quad 0<t<T_{c}$ (1)

where $T_{c}$ is the chirp ON time. Phase within the support of $t$ can be
represented by

$\displaystyle\phi(t)$ $\displaystyle=\int_{0}^{t}2\pi f(t)dt$
$\displaystyle=2\pi(f_{0}t+\frac{1}{2}S_{u}t^{2}).$ (2)

and the final transmitted chirp can be written as

$\displaystyle s(t)$ $\displaystyle=e^{\mathrm{j}\phi(t)}$
$\displaystyle=e^{\mathrm{j}2\pi(f_{0}t+\frac{1}{2}S_{u}t^{2})}.$ (3)

After the formation of clusters, the main task at hand is to estimate and
assign a velocity profile and heading to real life objects represented by said
clusters. Extensive testing on radar data had led to the observation that
standard estimation methods are not inherently immune to the effects of micro-
Dopplers, multipath reflections [16], clustering inaccuracies, etc. Robust
determination of object velocity and heading thus necessitates an inlier
estimation algorithm as a prerequisite for the purpose of noise elimination
[11, 10]. Furthermore, standard inlier subset detection algorithms like RAndom
SAmple Consensus (RANSAC) [17] are iterative in nature and scale with both the
number of points in each cluster, and the numbers of clusters itself [18, 19].

Radar data captured in urban scenarios has shown that the number of detected
clusters in the PC vary significantly with the infrastructure and prevalent
traffic conditions. Existing radar processing algorithms determine the
velocity profile and heading angle of a single object at a time, resulting in
drastic variation of total processing time based on the neighbourhood of the
ego-vehicle.

Hence, an end-to-end time-optimized algorithm is needed to increase the
throughput of the radar processing stack and improve crash avoidance and
collision warning for advanced driver assistance systems and autonomous
vehicles. This work discusses such an algorithm that leverages the power of
graphics processing units (GPUs) to run thousands of iterations in parallel by
exploiting simultaneous thread execution capabilities.

Figure 4: Starting, Intermediate, and Final Camera Images corresponding to
frames resulting in undesired values for LSQ only solution. Radar data was
captured by a vehicle accelerating towards a motorbike (78.6 GHz center freq.,
4GHz BW)

Figure 5: Comparison between LSQ only, and RANSAC plus LSQ fits on unfiltered
input data corresponding to the frames represented by the camera images above.
Performing RANSAC on the input dataset with a threshold of MAD results in
inlier categorization of points corresponding to the true projected Doppler of
the motorbike in the camera frames of Fig 4. Without this filtering pre-
requisite, LSQ alone is skewed easily by the presence of varying projected
Dopplers. (78.6 GHz center freq., 4GHz BW)

##  III Proposed Solution

The two proposed solutions for parallelized RANSAC inlier detection and the
subsequent vector velocity detection are detailed in the following
subsections.

###  III-A Parallelized RANSAC Algorithm

The rapid development of high-performance manycore GPUs has led to the
widespread access of computing platforms with increasing instruction
throughputs. These platforms have enabled parallel-computing based algorithms
to become an attractive option for practical radar processing applications.
Fig. 3 re-iterates the need for an inlier estimation algorithm by illustrating
how velocity and heading determination can be heavily skewed by the presence
of modulated frequency components in the input dataset.

The algorithm discussed in this section performs simultaneous categorization
of all points in the PC as either inliers or outliers of their respective
clusters. This concurrency is achieved by scheduling numerous thread
instances, with each thread corresponding to a linear fit generated by a pair
of seed points. While $N\\_clusts$ represents the number of clusters detected
in a frame, the $max\\_trials$ parameter determines the number of such fits
evaluated for each cluster, and is usually adjusted based on the hardware
specification of the user’s platform (4)

$\begin{split}thread\\_count=max\\_trials\times N\\_clusts\end{split}$ (4)

where $thread\\_count$ represents the number of threads that the GPU can
schedule in parallel without significant increase in execution time. This
reconstructed version of the standard RANSAC algorithm operates on the
projected velocity ($\mathbf{v_{r}}$) and azimuth angle ($\theta$) of all
points belonging to formed clusters.

Algorithm 1 RANSAC GPU Kernel

0: Random no.s spanning indices for every cluster

0: $tIdx\leq max\\_trials\times N\\_clusts$

$curr\\_cluster\leftarrow tIdx/N\\_clusts$

$seedpt1\leftarrow$ $normalize$($data[RND1]$)

$seedpt2\leftarrow$ $normalize$($data[RND2]$)

$m\leftarrow$ $Slope$($seedpt1,seedpt2$)

$c\leftarrow$ $Y\,intercept$($seedpt1,seedpt2$)

for all $index$ in $curr\\_cluster$ do

$test\\_pt\leftarrow$ norm ($data[index]$)

$aa=-m$, $bb=1$, $cc=-c$

$distn=aa\times test\\_pt[0]+bb\times test\\_pt[1]+cc$

$distd=\sqrt{aa^{2}+bb^{2}}$

$dist=$ abs$(distn/distd)$

if $dist\leq thresh[curr\\_cluster]$ then

$inliers[tIdx,index]=1$

$num[tIdx]=num[tIdx]+1$

end if

end for

The median absolute deviation (MAD) of normalized $\mathbf{v_{r}}$ values
determines the thresholding corridor around the linear fit (5). This
statistical parameter was selected after comparison with other known data
variability measures for both real and simulated datasets. Shown below is the
evaluation of this parameter, where $\mathbf{v_{r_{i}}}$ is the projected
velocity of each cluster point, and $med(x)$, denotes the median of the $x$
vector.

$\begin{split}MAD=\frac{1}{n}\sum_{i=1}^{n}[\mathbf{v_{r_{i}}}-med(\mathbf{v_{r}})]\end{split}$
(5)

Within each thread, while two seed points are used to generate the fit
hypothesis, the remaining cluster points are iterated through to calculate
respective perpendicular distances to the fit.

The inlier mask for a given cluster is chosen based on the thread index
corresponding to iteration that results in the largest number of post-
thresholding points. The combined inlier subset of the PC is determined based
on the logical OR operation on the parallely generated inlier masks for each
cluster. The result of the RANSAC filtering is much easily observed in the
straight line motion of a motorbike traveling ahead of the ego-vehicle in Fig
5.

Figure 6: [Left] Detected relative velocity plotted across frames shown in Fig
4. Comparison is drawn between the approach implemented RANSAC plus LSQ, and
LSQ only, [Right] Estimated Heading Angle is plotted across the same range of
frames. Similar comparison is plotted to observe the difference in variability
in the estimated angle (78.6 GHz center freq., 4GHz BW)

###  III-B Instantaneous velocity and Heading angle estimation by Least-
Squares Kernel

One of the important problems for ADAS is to estimate the velocity profile and
heading angle of the objects detected by the ego-vehicle. Keeping in mind the
need for real-time sensing and detection, a GPU kernel is proposed to
implement the least squares approximation method. This kernel instantiates a
thread in parallel for each cluster to ensure optimal processing time,
independent of the number of clusters. Let, $v_{x}$ and $v_{y}$ be the
velocities in the longitudinal and lateral directions, respectively (refer
Fig. 2 for axes representation). Also, let $\mathbf{v_{r_{i}}}$ and
$\mathbf{\theta_{i}}$ denote the vector of detected radial Doppler speeds and
the azimuth angles of all the points on the $i^{th}$ cluster. For any detected
point $j$ on cluster $i$, we have

${v_{r_{i,j}}}=v_{x_{i,j}}\text{cos}({\theta_{i,j}})+v_{y_{i,j}}\text{sin}({\theta_{i,j}}),$
(6)

where $j\in[1,N_{i}]$ with $N_{i}$ denoting the number of points making the
$i^{th}$ cluster. Also, $v_{x_{i,j}}$ and $v_{y_{i,j}}$ denote the vector
velocity components of a detected point on the cluster. Under the rigid body
assumption, each reflecting point on a cluster will have the same vector
velocity. Using vector notations, the velocity vector
$\mathbf{v_{i}}=(v_{x_{i}},v_{y_{i}})^{T}$ is fixed for all the points on
cluster $i$. Therefore, Eq. (6) can be rewritten as follows:

$\mathbf{v_{r_{i}}}=A_{i}\mathbf{v_{i}},\implies\mathbf{v_{i}}=A_{i}^{\dagger}\mathbf{v_{r_{i}}},$
(7)

where

$A_{i}=\begin{bmatrix}\text{cos}(\theta_{i1})&\text{sin}(\theta_{i1})\\\
\text{cos}(\theta_{i2})&\text{sin}(\theta_{i2})\\\ \vdots&\vdots\\\
\text{cos}(\theta_{iN_{i}})&\text{sin}(\theta_{iN_{i}})\end{bmatrix}.$

and $A_{i}^{\dagger}$ denotes the matrix inverse of $A_{i}$ [20]. Following
this, we compute the heading angle $(\phi_{i})$ which is defined as the arc
tangent of the ratio between $v_{y_{i}}$ and $v_{x_{i}}$ with respect to the
ADAS X-axis [21].

$\phi_{i}=\arctan(\frac{v_{y_{i}}}{v_{x_{i}}}).$ (8)

Algorithm 2 GPU kernel for component velocity $(v_{x},v_{y})$ and heading
angle $(\phi)$ estimation

0: $A_{i}$ and $v_{r_{i}}$ for every cluster (Eq. (6), Eq. (7))

0: $i\leq num\\_of\\_clusters$

$\mathbf{v_{i}}=A_{i}^{\dagger}\mathbf{v_{ri}}$

Set $\mathbf{v_{i}[1]}$ as the lateral velocity,
$v_{y_{i}}\leftarrow\mathbf{v_{i}[1]}$

Set $\mathbf{v_{i}[0]}$ as the longitudinal velocity,
$v_{x_{i}}\leftarrow\mathbf{v_{i}[0]}$

$\phi_{i}\leftarrow\arctan(\frac{v_{y_{i}}}{v_{x_{i}}})$

For each cluster $i$, coordinate pair $(v_{xi},v_{yi})$ denotes the vector
velocity of the cluster. Each of the cluster’s vector velocity gets computed
in parallel, subsequent to inlier detection using RANSAC.

##  IV Results

The implementation proposed in this work is compared against a standard LSQ
approach based on tests run on captured radar data. Both approaches were
evaluated on the basis of their ability to $(1)$ operate in real-time and
$(2)$ accurately determine instantaneous velocity and heading angle of
detected objects. As seen in Fig. 6, the presented approach has much better
overall performance with lesser erroneous spikes in the determined velocity,
and no output inconsistency in estimated heading angle. The depicted overall
improvement is a result of the real-time outlier eradication facilitated by
the proposed parallelized RANSAC algorithm and heading estimation.

TABLE I: RANSAC Execution Time Comparisons Number of | Points per | Parallelized | sklearn
---|---|---|---
Clusters | Cluster | RANSAC | RANSAC
8 | $100$ | $1.64ms$ | $70.50ms$
8 | $150$ | $1.77ms$ | $76.58ms$
16 | $100$ | $1.68ms$ | $141.15ms$
16 | $150$ | $1.82ms$ | $152.97ms$
32 | $100$ | $1.73ms$ | $281.95ms$
32 | $150$ | $1.89ms$ | $305.85ms$
64 | $100$ | $1.98ms$ | $564.48ms$
64 | $150$ | $2.23ms$ | $611.91ms$

This proposed approach was further compared against the commonly used scikit-
learn RANSACRegressor [22] for varying input conditions as shown in Table I.
An NVIDIA GeForce GTX 1060 3GB [23] was used to run the GPU kernels discussed
in Sections $3.1$ and $3.2$. The approach presented in this work was observed
to be much faster due to its ability to perform inlier [outlier] detection on
points using a twofold parallelization scheme, $(1)$ across multiple detected
points on a cluster and $(2)$ across multiple clusters themselves. Scikit-
learn RANSACRegressor could filter out noise from only one cluster at a time,
resulting in total time taken scaling linearly with the number of clusters.
With real-time applications as autonomous driving and when the number of
detected clusters scale up depending on traffic scenarios, it is evident that
such a twofold parallelization scheme is an extremely important feature.
Similar comparison is also shown between the proposed LSQ Kernel and a
standard LSQ approach using numpy.linalg.pinv [24] as shown in Table II. This
comparison also indicates the significant overall speed-up achieved by using
the parallelized pseudoinverse algorithm.

TABLE II: LSQ Execution Time Comparisons Number of | Points per | Parallelized | NumPy
---|---|---|---
Clusters | Cluster | pseudoinverse | pseudoinverse
8 | $100$ | $1.81ms$ | $2.08ms$
8 | $150$ | $1.83ms$ | $2.10ms$
16 | $100$ | $1.81ms$ | $4.18ms$
16 | $150$ | $1.83ms$ | $4.21ms$
32 | $100$ | $1.81ms$ | $8.74ms$
32 | $150$ | $1.83ms$ | $8.77ms$
64 | $100$ | $1.81ms$ | $16.84ms$
64 | $150$ | $1.83ms$ | $17.19ms$

##  V Conclusion

The algorithm presented in this paper showcases a significant speed-up of
existing estimation algorithms along with supporting quantitative metrics,
after extensive verification across various practical scenarios. This work
demonstrates the need for robust and optimized inlier detection methods for
time critical applications such as ADAS, traffic monitoring, etc. The proposed
algorithm is a significant step forward in the development of practically
realizable ADAS, and will scale directly with enhancements in GPU hardware and
programming interfaces as well.

## References

  * [1] ”DENSE - 24/7 Automotive Sensing System”. Available: $https://dense247.eu/home/$
  * [2] A. Brenner, and L. Roessing,“Radar imaging of urban areas by means of very high-resolution SAR and interferometric SAR,” IEEE Transactions on Geoscience and Remote Sensing, Vol. 46, No. 10, 2971–2982,(2008). IEEE.
  * [3] F. Bovolo, and C. Marin, and L. Bruzzone,“A hierarchical approach to change detection in very high resolution SAR images for surveillance applications,” IEEE Transactions on Geoscience and Remote Sensing, Vol. 51, No. 4, 2042–2054,(2012). IEEE.
  * [4] S. Sun, A. Petropulu and H. Poor,“MIMO radar for advanced driver-assistance systems and autonomous driving: Advantages and challenges,” IEEE Signal Processing Magazine, Vol. 37, No. 4, 98–117,(2020). IEEE.
  * [5] J. H. Friedman, F. Baskett and L. J. Shustek, ”An Algorithm for Finding Nearest Neighbors,” in IEEE Transactions on Computers, vol. C-24, no. 10, pp. 1000-1006, Oct. 1975, doi: 10.1109/T-C.1975.224110.
  * [6] M. Ester, H. Kriegel, J. Sander, and X. Xu and others, ”A density-based algorithm for discovering clusters in large spatial databases with noise”, 1996 Kdd, Vol. 96, No. 34, pp. 226–231.
  * [7] H. Darong, and W. Peng,“Grid-based DBSCAN algorithm with referential parameters,” Physics Procedia, Vol. 24, 1166–1170,(2012). Elsevier.
  * [8] S. Wang, C. Gao, Q. Zhang, H. Zeng, and J. Bai, ”The Latest Research on Clustering Algorithms Used for Radar Signal Sorting”, in Recent Trends in Intelligent Computing, Communication and Devices, pp. 799–805, 2020, Spinger.
  * [9] V. Chen, The micro-Doppler effect in radar, Artech House, 2019.
  * [10] D. Kellner, and M. Barjenbruch and J. Klappstein and J. Dickmann, and K. Dietmayer, ”Instantaneous full-motion estimation of arbitrary objects using dual Doppler radar”, in 2014 IEEE Intelligent Vehicles Symposium Proceedings, pp. 324–329, 2014, IEEE.
  * [11] D. Kellner, M. Barjenbruch, K. Dietmayer, J. Klappstein, and J. Dickmann, ”Instantaneous lateral velocity estimation of a vehicle using Doppler radar”, in Proceedings of the 16th International Conference on Information Fusion, pp. 877–884, 2013, IEEE.
  * [12] D. Kellner, M. Barjenbruch, J. Klappstein, and J. Dickmann,“Tracking of extended objects with high-resolution Doppler radar,” IEEE Transactions on Intelligent Transportation Systems, Vol. 17, No. 5, 1341–1353,(2015). IEEE.
  * [13] P. Lidman, and S. Luu, ”Clustering, shape extraction and velocity estimation applied to radar detections,” Masters Thesis. 2018
  * [14] G. Rajendran, A. Sivadas, “Method, Apparatus and Device for Doppler Compensation in a Time Switched MIMO Radar System ,” U.S. Patent US20200088867A1, March, 20, 2008.
  * [15] G. L. Charvat, A. J. Fenn and B. T. Perry, ”The MIT IAP radar course: Build a small radar system capable of sensing range, Doppler, and synthetic aperture (SAR) imaging,” 2012 IEEE Radar Conference, Atlanta, GA, 2012, pp. 0138-0144, doi: 10.1109/RADAR.2012.6212126.
  * [16] R. Feger, and C. Pfeffer, and W. Scheiblhofer, C. Schmid, M. Lang, and A. Stelzer,“A 77-GHz cooperative radar system based on multi-channel FMCW stations for local positioning applications,” IEEE transactions on microwave theory and techniques, Vol. 61, No. 1, 676–684,(2012). IEEE.
  * [17] K. Derpanis, “Overview of the RANSAC Algorithm,” Image Rochester NY, Vol. 4, No. 1, 2–3,(2010).
  * [18] O. Chum, and J. Matas, Vol. Optimal randomized RANSAC, No. IEEE Transactions on Pattern Analysis and Machine Intelligence, 30,(8). 1472–1482. 2008 IEEE,“
  * [19] R. Schnabel, R. Wahl, and R. Klein, ”Efficient RANSAC for point-cloud shape detection”, in Computer graphics forum, Vol. 26, No. 2, pp. 214–226, 2007, Wiley Online Library.
  * [20] S. Sawyer, ”Generalized Inverses: How to Invert a Non-Invertible Matrix” (2008). Available: https://www.math.wustl.edu/ sawyer/handouts/GenrlInv.pdf
  * [21] ”Coordinate Systems in Automated Driving Toolbox”, Get Started with Automated Driving Toolbox, Mathworks. Available: https://in.mathworks.com/help/driving/ug/coordinate-systems.html
  * [22] Pedregosa,“Scikit-learn: Machine Learning in Proc,” Journal of Machine Learning Research, Vol. 12, 2825–2830,(2011).
  * [23] D. Kirk. ”NVIDIA cuda software and gpu parallel computing architecture”. In Proceedings of the 6th international symposium on Memory management (ISMM ’07). Association for Computing Machinery, 2007, New York, NY, USA, 103–104. DOI:https://doi.org/10.1145/1296907.1296909
  * [24] Harris, C.R., Millman, K.J., van der Walt, S.J. et al. Array programming with NumPy. Nature 585, 357–362 (2020). DOI: 0.1038/s41586-020-2649-2.