#  $\bm{D^{*}}$ and $\bm{D^{*}_{s}}$ distribution amplitudes from Bethe-
Salpeter wave functions

Fernando E. Serna  LFTC, Universidade Cidade de São Paulo, Rua Galvão Bueno
868, São Paulo, SP 01506-000, Brazil  Departamento de Física, Universidad de
Sucre, Carrera 28 No. 5-267, Barrio Puerta Roja, Sincelejo, Colombia  Roberto
Correa da Silveira  LFTC, Universidade Cidade de São Paulo, Rua Galvão Bueno
868, São Paulo, SP 01506-000, Brazil  Bruno El-Bennich  LFTC, Universidade
Cidade de São Paulo, Rua Galvão Bueno 868, São Paulo, SP 01506-000, Brazil

###### Abstract

We report on the first calculation of the longitudinal and transverse light
front distribution amplitudes of the $D^{*}$ and $D^{*}_{s}$ mesons and their
first four moments. As a byproduct, we also obtain these distribution
amplitudes for the $\rho$, $\phi$, $K^{*}$ and $J/\Psi$ mesons and confirm a
prediction of lattice QCD for the vector kaon: while the longitudinal
distribution amplitude is almost symmetric, the transverse one is oblique
implying that the strange quark carries more momentum.

##  I Motivation

In relativistic quantum field theory the infinite degrees of freedom do not
allow for a straightforward definition of a particle’s wave function as in
quantum mechanics. In particular, in Quantum Chromodynamics (QCD) the
fundamental quark and gluon fields are not even observable. On the other hand,
the bound states of _valence_ quark-antiquark pairs can be described by a
Bethe-Salpeter wave function, the closest relative to a wave function in
quantum mechanics. Still, in the instant-form of QCD dynamics these wave
functions are defined in an infinite-body field theory in which particles
interact and their number is not conserved.

One could overcome this difficulty if the hadron’s light-front wave function
was known exactly, though realistic calculations of hadronic bounds states in
the front form are a challenging task. A different path to a sensible
definition of a wave function in quantum field theory is drawn by projecting
the Bethe-Salpeter wave functions on the light front. This yields the hadron’s
light-front distribution amplitude (LFDA) which describes the longitudinal
momentum distribution of valence quarks in the limit of negligible transverse
momentum. While they are non-measurable objects, they are widely being applied
in hadron and flavor physics.

For instance, the asymptotic LFDA of the pion,
$\phi(x,\mu)\stackrel{{\scriptstyle\mu\to\infty}}{{=}}6x(1-x)$, enters in the
expression of its elastic electromagnetic form factor at very large momentum
transfers Efremov:1978rn ; Lepage:1980fj . Since the LFDAs are scale-dependent
and become broader at smaller momenta, they directly influence the momentum
dependence of the elastic form factors in momentum regions accessible in
collider experiments Chang:2013nia ; Raya:2015gva ; Raya:2019dnh . Weak $B$
decays into two light(er) mesons are frequently treated as hard exclusive
processes in which the decay amplitude is factorized into perturbative short-
distance contributions and a nonperturbative transition amplitude. Here too,
the LFDAs enter both, the hard-scattering integrals and the heavy-to-light
transition amplitudes Beneke:2002jn ; Bauer:2005wb ; El-Bennich:2006rcn ; El-
Bennich:2009gqk ; Leitner:2010fq . More recently, the exclusive electroweak
production of $D_{s}^{(*)}$ mesons on an unpolarized nucleon was investigated
in the framework of collinear QCD factorization which also involves the heavy
meson’s LFDA Pire:2015iza ; Pire:2017lfj ; Pire:2017yge ; Pire:2021dad .

Beyond its numerous applications in hard exclusive processes, these one-
dimensional distributions provide a practical probability interpretation of
partons, as in this frame the particle number is conserved. Namely, the
distributions $\phi(x,\mu)$ express the light-front fraction of the hadron’s
momentum that a valence quark carries. Another compelling feature is that one
can observe the qualitative and quantitative impact of dynamical chiral
symmetry breaking (DCSB) on the LFDA at a given scale $\mu$. For instance, the
distribution amplitude $\phi_{\pi}(x,\mu)$ of the pion is a concave function
which clearly evolves from its asymptotic $\mu\to\infty$ form to a much
broader distribution Chang:2013pq . Similarly, the kaon’s distribution
amplitude, $\phi_{K}(x,\mu)$, is not symmetric about the midpoint $x=1/2$,
which expresses nothing but SU(3) flavor symmetry breaking, and that asymmetry
is exacerbated with increasing mass difference of the quarks Shi:2015esa ;
Serna:2020txe .

The question arises of how DCSB impacts antiquark-quark states in other
$J^{PC}$ channels and an extension to the vector mesons is natural. Moreover,
the LFDA of vector mesons arises in the collinear factorization of weak
$B$-decay amplitudes Beneke:2000wa  and in diffractive vector-meson production
Forshaw:2010py ; Forshaw:2012im . Within the combined framework of the Dyson-
Schwinger equation (DSE) and the Bethe-Salpeter equation (BSE) Bashir:2012fs
the LFDAs of the $\rho$ and $\phi$ mesons were calculated in Ref. Gao:2014bca
and later the LFDAs of heavy quarkonia were obtained in Ref. Ding:2015rkn . In
here, using a kindred DSE and BSE framework, we extend earlier work on $D$ and
$D_{s}$ distribution amplitudes Serna:2020txe  to those of their vector
partners and make predictions for the twist-2 LFDA of the $D^{*}$ and
$D_{s}^{*}$ mesons considering the two-quark Fock-state of their light front
wave function. Along the way, we compute the LFDA of the $\rho$, $K^{*}$
$\phi$ and $J/\Psi$ mesons and compare them with the distribution amplitudes
of other approaches Forshaw:2010py ; Forshaw:2012im ; Ball:2007zt ;
Boyle:2008nj ; Hua:2020gnw ; Ding:2015rkn .

##  II Twist-Two Distribution Amplitudes

A vector meson with total momentum $P$ and mass $m_{V}$, $P^{2}=-m_{V}^{2}$,
made of a quark and an antiquark of flavors $f$ and $g$ is described by four
twist-two distribution amplitudes, though only two of them are independent at
leading twist as a consequence of a Wandzura-Wilczek type of relation
Ball:2007zt . The two LFDAs we consider, $\phi^{\parallel}_{V}(x;\mu)$ and
$\phi^{\perp}_{V}(x;\mu)$, describe the fraction of total momentum on the
light front, $x=k^{+}/P^{+}=(k_{0}+k_{z})/(P_{0}+P_{z})$, carried by the quark
in longitudinally and transversely polarized mesons, respectively. They can be
extracted from the Bethe-Salpeter wave function, $\chi^{fg}_{V\nu}(k;P)$, with
the following projections onto the light front Gao:2014bca ; Ding:2015rkn :

$\displaystyle f_{V}\phi^{\parallel}_{V}(x;\mu)=\frac{m_{V}N_{c}{\cal
Z}_{2}}{\sqrt{2}\ n\cdot P}$ $\displaystyle\
\mathrm{Tr}_{D}\\!\int^{\Lambda}\\!\\!\frac{d^{4}k}{(2\pi)^{4}}\,\delta(n\cdot
k_{\eta}-x\,n\cdot P)$ $\displaystyle\times\,\gamma\cdotp
n\,n_{\nu}\,\chi^{fg}_{V\nu}(k,P)\ ,$ (1) $\displaystyle
f^{\perp}_{V}\phi^{\perp}_{V}(x;\mu)=-\frac{N_{c}{\cal Z}_{T}}{2\sqrt{2}}$
$\displaystyle\
\mathrm{Tr}_{D}\\!\int^{\Lambda}\\!\\!\frac{d^{4}k}{(2\pi)^{4}}\,\delta(n\cdot
k_{\eta}-x\,n\cdot P)$
$\displaystyle\times\,n_{\mu}\sigma_{\mu\rho}\,\mathcal{O}^{\perp}_{\rho\nu}\,\chi^{fg}_{V\nu}(k,P)\
,$ (2)

where $N_{c}=3$, $n=(0,0,1,i)$ is a light-like vector and
$\bar{n}=\tfrac{1}{2}(0,0,-1,i)$ its conjugate with $n^{2}=\bar{n}^{2}=0$,
$n\cdot P=-m_{V}$, $\bar{n}\cdot P=-m_{V}/2$ and $n\cdot\bar{n}=-1$.111We use
Euclidean metric with the Dirac algebra:
$\\{\gamma_{\mu},\gamma_{\nu}\\}=2\delta_{\mu\nu}$;
$\gamma_{\mu}^{\dagger}=\gamma_{\mu}$;
$\gamma_{5}=\gamma_{4}\gamma_{1}\gamma_{2}\gamma_{3}$,
tr$[\gamma_{4}\gamma_{\mu}\gamma_{\nu}\gamma_{\rho}\gamma_{\sigma}]=-4\,\epsilon_{\mu\nu\rho\sigma}$;
$\sigma_{\mu\nu}=(i/2)[\gamma_{\mu},\gamma_{\nu}]$; $a\cdot
b=\sum_{i=1}^{4}a_{i}b_{i}$. A time-like vector $P_{\mu}$ satisfies $P^{2}<0$.
In Eq. (2) the Dirac commutator $\sigma_{\mu\nu}$ is contracted with the
tensor Lu:2021sgg ,

$\mathcal{O}^{\perp}_{\rho\nu}=\delta_{\rho\nu}+n_{\rho}\bar{n}_{\nu}+\bar{n}_{\rho}n_{\nu}\
.$ (3)

In Eqs. (1) and (2), $\chi^{fg}_{V\nu}(k;P)=S_{f}(k_{\eta})$
$\Gamma^{fg}_{V\nu}(k,P)S_{g}(k_{\bar{\eta}})$ is the projected wave function,
where $\Gamma^{fg}_{V\nu}(k;P)$ denotes the Bethe-Salpeter amplitude (BSA) and
$S_{f}(k_{\eta})$ and $S_{g}(k_{\bar{\eta}})$ are respectively the quark and
antiquark propagators with momenta $k_{\eta}=k+\eta P$ and
$k_{\bar{\eta}}=k-\bar{\eta}P$. The details of their calculation, solving
numerically the DSE for the quarks of a given flavor and the BSE for a vector
meson, in particular the $D$ and $D^{*}$ mesons, are provided elsewhere
Serna:2020txe ; Mojica:2017tvh ; El-Bennich:2021ldv . The parameters
$\eta+\bar{\eta}=1$ define momentum fractions and $\Lambda$ is an ultraviolet
regularization mass-scale; no observables can depend on $\eta$, $\bar{\eta}$
and $\Lambda$ owing to Poincaré covariance. Furthermore, ${\cal
Z}_{2}(\mu,\Lambda)$ is the wave-function renormalization constant and ${\cal
Z}_{T}(\mu,\Lambda)$ is the tensor-vertex renormalisation constant of the
quark. Both constants as well as $f^{\perp}_{V}$ depend on the renormalization
scale $\mu$, whereas $f_{V}$ is renormalization-point independent and measures
the strength of the $\rho^{0}\to e^{+}e^{-}$ decay amplitude.

The expressions for $\phi^{\parallel}_{V}(x;\mu)$ and
$\phi^{\perp}_{V}(x;\mu)$ in Eqs. (1) and (2) are not amenable to
straightforward numerical integration. Instead, one computes Mellin moments
Chang:2013pq ,

$\displaystyle\langle x^{m}\rangle_{\parallel}$
$\displaystyle=\int^{1}_{0}x^{m}\,\phi^{\parallel}_{V}(x,\mu)\,dx\ ,$ (4)
$\displaystyle\langle x^{m}\rangle_{\perp}$
$\displaystyle=\int^{1}_{0}x^{m}\,\phi^{\perp}_{V}(x,\mu)\,dx\ ,$ (5)

from which one can reconstruct the distribution amplitudes on the domain
$x\in[0,1]$. The BSA normalization ensures that $\langle
x^{0}\rangle_{\parallel}=\langle x^{0}\rangle_{\perp}=1$ which in turn defines
the vector and tensor decay constants, $f_{V}^{\parallel}$ and
$f^{\perp}_{V}$.

Integrating both sides of Eqs. (1) and (2) and applying the Dirac-function
property $\int_{0}^{1}x^{m}\delta(a-xb)dx=\frac{a^{m}}{b^{m+1}}\,\theta(b-a)$,
leads to the expressions,

$\displaystyle\langle x^{m}\rangle_{\parallel}=\frac{m_{V}N_{c}{\cal
Z}_{2}}{\sqrt{2}\,f_{V}}$
$\displaystyle\,\mathrm{Tr}_{D}\\!\int^{\Lambda}\\!\frac{d^{4}k}{(2\pi)^{4}}\frac{(n\cdot
k_{\eta})^{m}}{(n\cdot P)^{m+2}}$ $\displaystyle\times\,\gamma\cdot
n\,n_{\nu}\,\chi^{fg}_{V\nu}(k,P)\,,$ (6) $\displaystyle\langle
x^{m}\rangle_{\perp}=-\frac{N_{c}{\cal Z}_{T}}{2\sqrt{2}\,f^{\perp}_{V}}$
$\displaystyle\,\mathrm{Tr}_{D}\\!\int^{\Lambda}\\!\frac{d^{4}k}{(2\pi)^{4}}\frac{(n\cdot
k_{\eta})^{m}}{(n\cdot P)^{m+1}}$
$\displaystyle\times\,n_{\mu}\sigma_{\mu\rho}\,\mathcal{O}^{\perp}_{\rho\nu}\,\chi^{fg}_{V\nu}(k,P)\,.\hskip
8.53581pt$ (7)

With this, we are in principle able to compute Mellin moments to arbitrary
order $m$. We do so by employing the numerical solutions of the quark
propagators for complex momenta defined by the parabolas,
$k_{\eta}^{2}=k^{2}-\eta^{2}m^{2}_{V}+2i\eta\,m_{V}|k|z_{k}$ and
$k_{\bar{\eta}}^{2}=k^{2}-{\bar{\eta}}^{2}m^{2}_{V}-2i\bar{\eta}\,m_{V}|k|z_{k}$,
where $z_{k}=k\cdot P/|k||P|$, $-1\leq z\leq+1$, and of the BSA of the vector
mesons El-Bennich:2021ldv . That is, other than in Ref. Serna:2020txe , we do
not rely on complex-conjugate pole parametrizations of the propagators nor on
Nakanishi representations of the BSA, as the latter introduce ambiguities when
fitted to numerical solutions. However, direct integration comes at the price
that we can only access moments up to $m_{\mathrm{max}}=4-6$, as the numerical
error of the integral becomes significant for larger moments. These moments,
though, are sufficient to reconstruct the desired LFDA.

| $\langle x\rangle_{\parallel,\perp}$ | $\langle x^{2}\rangle_{\parallel,\perp}$ | $\langle x^{3}\rangle_{\parallel,\perp}$ | $\langle x^{4}\rangle_{\parallel,\perp}$ | $a_{1}^{\parallel,\perp}$ | $a_{2}^{\parallel,\perp}$ | $\alpha^{\parallel,\perp}$
---|---|---|---|---|---|---|---
$\rho_{\parallel}$ | 0.500 | 0.312 | 0.226 | 0.161 | 0.0 | $0.003\pm 0.038$ | $0.908\pm 0.023$
$\rho_{\perp}$ | 0.500 | 0.312 | 0.218 | 0.160 | 0.0 | $-0.136\pm 0.007$ | $0.799\pm 0.006$
$\phi_{\parallel}$ | 0.500 | 0.296 | 0.195 | 0.134 | 0.0 | $-0.372\pm 0.010$ | $0.864\pm 0.010$
$\phi_{\perp}$ | 0.500 | 0.296 | 0.193 | 0.134 | 0.0 | $-0.386\pm 0.002$ | $0.870\pm 0.002$
$K^{*}_{\parallel}$ | 0.509 | 0.323 | 0.236 | 0.179 | $0.041\pm 0.027$ | $-0.191\pm 0.048$ | $0.643\pm 0.031$
$K^{*}_{\perp}$ | 0.528 | 0.351 | 0.262 | 0.204 | $0.119\pm 0.003$ | $0.122\pm 0.015$ | $0.840\pm 0.019$
Table 1: The first four Mellin moments, $\langle x^{m}\rangle_{\parallel}$ and $\langle x^{m}\rangle_{\perp}$, of the light vector mesons and the coefficients of their reconstructed Gegenbauer expansion (8). The errors on $a_{1}$, $a_{2}$ and $\alpha$ stem from the minimization.  $\rho_{\parallel,\perp}$ | $\langle\xi^{2}\rangle_{\parallel}$ | $\langle\xi^{2}\rangle_{\perp}$ | $\langle\xi^{4}\rangle_{\parallel}$ | $\langle\xi^{4}\rangle_{\perp}$ | $\langle\xi^{6}\rangle_{\parallel}$ | $\langle\xi^{6}\rangle_{\perp}$ | $\langle\xi^{8}\rangle_{\parallel}$ | $\langle\xi^{8}\rangle_{\perp}$
---|---|---|---|---|---|---|---|---
Herein | 0.263 | 0.250 | 0.136 | 0.127 | 0.090 | 0.081 | 0.062 | 0.044
DSE Gao:2014bca  | 0.231 | 0.252 | 0.109 | 0.126 | 0.065 | 0.079 | 0.044 | 0.056
QCDSR Ball:2007zt  | 0.234 | 0.238 | 0.109 | 0.111 | 0.063 | 0.065 | 0.042 | 0.043
HERA Forshaw:2010py ; Forshaw:2012im  | 0.227 | 0.260 | 0.105 | 0.130 | 0.062 | 0.079 | 0.041 | 0.054
LQCD Boyle:2008nj  | 0.240(40) |  |  |  |  |  |  |
$\phi_{\parallel,\perp}$ | $\langle\xi^{2}\rangle_{\parallel}$ | $\langle\xi^{2}\rangle_{\perp}$ | $\langle\xi^{4}\rangle_{\parallel}$ | $\langle\xi^{4}\rangle_{\perp}$ | $\langle\xi^{6}\rangle_{\parallel}$ | $\langle\xi^{6}\rangle_{\perp}$ | $\langle\xi^{8}\rangle_{\parallel}$ | $\langle\xi^{8}\rangle_{\perp}$
Herein | 0.186 | 0.182 | 0.077 | 0.073 | 0.042 | 0.039 | 0.026 | 0.024
DSE Gao:2014bca  | 0.233 | 0.253 | 0.111 | 0.127 | 0.067 | 0.080 | 0.046 | 0.056
QCDSR Ball:2007zt  | 0.245 | 0.238 | 0.115 | 0.111 | 0.068 | 0.065 | 0.045 | 0.043
LQCD Hua:2020gnw  | 0.212 | 0.250 | 0.097 | 0.127 | 0.057 | 0.081 | 0.039 | 0.058
$K^{*}_{\parallel,\perp}$ | $\langle\xi^{2}\rangle_{\parallel}$ | $\langle\xi^{2}\rangle_{\perp}$ | $\langle\xi^{4}\rangle_{\parallel}$ | $\langle\xi^{4}\rangle_{\perp}$ | $\langle\xi^{6}\rangle_{\parallel}$ | $\langle\xi^{6}\rangle_{\perp}$ | $\langle\xi^{8}\rangle_{\parallel}$ | $\langle\xi^{8}\rangle_{\perp}$
Herein | 0.272 | 0.298 | 0.146 | 0.164 | 0.097 | 0.109 | 0.072 | 0.080
QCDSR Ball:2007zt  | 0.227 | 0.227 | 0.104 | 0.104 | 0.060 | 0.060 | 0.039 | 0.039
LQCD Hua:2020gnw  | 0.200 | 0.292 | 0.088 | 0.162 | 0.050 | 0.111 | 0.032 | 0.084
Table 2: Comparison of $\langle\xi^{2m}\rangle_{\parallel,\perp}$ moments for
the $\rho$, $\phi$ and $K^{*}$ mesons. The QCDSR values are obtained with Eqs.
(4) and (5) employing the Gegenbauer expansion (8) with $\alpha=3/2$ and the
value for $a_{2}$ in Ref. Ball:2007zt . Similarly, we fit the tabulated values
of $\phi^{\parallel}_{V}(x;\mu)$ and $\phi^{\perp}_{V}(x;\mu)$ provided in
Ref. Hua:2020gnw  with the same Gegenbauer expansion and use them to calculate
the moments.

We proceed as in Refs. Chang:2013pq ; Shi:2015esa ; Serna:2020txe ;
Gao:2014bca ; Ding:2015rkn  and in the case of the light vector mesons we use
an expansion in terms of Gegenbauer moments $C_{n}^{\alpha}(2x-1)$, which form
a complete orthonormal set on $x\in[0,1]$ with respect to the measure
$[x(1-x)]^{\alpha-1/2}$, in order to reconstruct their two independent twist-
two LFDAs ($\bar{x}=1-x$):

$\phi^{\parallel,\perp}_{V\mathrm{rec.}}(x,\mu)=\mathcal{N}(\alpha)[x\bar{x}]^{\alpha-\tfrac{1}{2}}\Bigg{[}1+\sum^{N}_{n=1}a_{n}C^{\alpha}_{n}(2x-1)\Bigg{]}.$
(8)

This expansion is employed for neutral mesons as well as for flavored mesons,
which are not $C$-parity eigenstates. In case of the former, the odd
components $a_{n}$ vanish. In fitting the calculated moments in Eqs. (6) and
(7), we consider, besides the coefficients $a_{n}$, the power $\alpha$ itself
a parameter rather than projecting on the $\alpha=3/2$ basis. This allows to
limit the expansion to $N=2$ and considerably simplifies the fits discussed
below Chang:2013pq . The normalization is obtained as,

$\mathcal{N}(\alpha)=\frac{\Gamma(2\alpha+1)}{[\Gamma(\alpha+1/2)]^{2}}\ .$
(9)

The heavy vector mesons, i.e. the $D^{*}$, $D^{*}_{s}$ and $J/\psi$, are
parametrized with a different expression:

$\phi_{V\mathrm{rec.}}^{\parallel,\perp}(x,\mu)=\mathcal{N}(\alpha,\beta)\,4x\bar{x}\,e^{4\,\alpha
x\bar{x}+\beta(x-\bar{x})}\ .$ (10)

This functional form is more appropriate for a distribution amplitude with a
convex-concave-convex functional behavior that tends to a $\delta$-function in
the infinite heavy quark limit, as the use of an expansion, such as in Eq.
(8), leaves one no choice but to retain a large number of Gegenbauer moments.
A very similar functional expression is also found when the Nakanishi weight
function is extracted from the quarkonia’s Bethe-Salpeter wave function
Gao:2016jka . The normalization is given by Serna:2021xnr ,

$\displaystyle\mathcal{N}(\alpha,\beta)$ $\displaystyle=\
16\,\alpha^{5/2}\Bigg{[}4\sqrt{\alpha}\,\left(\beta\sinh\beta+2\alpha\cosh\beta\right)$
$\displaystyle+$ $\displaystyle\
\sqrt{\pi}\,e^{\alpha+\frac{\beta^{2}}{4\alpha}}\left(-2\alpha+4\alpha^{2}-\beta^{2}\right)$
$\displaystyle\times$
$\displaystyle\,\Bigg{\\{}\operatorname{Erf}\left(\frac{2\alpha-\beta}{2\sqrt{\alpha}}\right)+\operatorname{Erf}\left(\frac{2\alpha+\beta}{2\sqrt{\alpha}}\right)\Bigg{\\}}\Bigg{]}^{-1},$
(11)

in which the error function is defined as:
$\operatorname{Erf}(x)=\frac{2}{\sqrt{\pi}}\int^{x}_{0}dt\,e^{t^{2}}$.

We thus reconstruct the vector LFDAs by minimizing the sum,

$\epsilon_{\parallel,\perp}=\sum^{m_{\textrm{max}}}_{m=1}\left|\frac{\langle
x^{m}\rangle^{\textrm{rec.}}_{\parallel,\perp}}{\langle
x^{m}\rangle_{\parallel,\perp}}-1\right|\ ,$ (12)

where the moments $\langle x^{m}\rangle^{\textrm{rec.}}_{\parallel,\perp}$ are
calculated using Eqs. (4) and (5) and the expansion in either Eq. (8) or Eq.
(10), whereas $\langle x^{m}\rangle_{\parallel,\perp}$ denotes the moments in
Eqs. (6) and (7). It is useful to contrast our predictions for the
longitudinal and transverse LFDAs with those obtained using other approaches,
namely with lattice QCD (LQCD) Boyle:2008nj ; Hua:2020gnw , QCD sum rules
(QCDSR) Ball:2007zt  and with earlier calculations in the DSE-BSE framework
(DSE) Gao:2014bca . In order to do so we also compute the moments,

$\left\langle\xi^{2m}\right\rangle_{\parallel,\perp}=\int_{0}^{1}\,\xi^{2m}\,\phi^{\parallel,\perp}_{V}(x,\mu)dx\
,$ (13)

in terms of the difference of momentum fractions, $\xi=x-(1-x)=2x-1$.

Figure 1: Comparison of the longitudinal and transverse LFDAs for the $\rho$
(top panel) and $\phi$ (bottom panel) mesons with those of Refs. Gao:2014bca
(DSE) and Hua:2020gnw  (LQCD) at $\mu=2$ GeV. Error bands reflect the
uncertainties of the fit parameters in Table 1. The intervals $0\leq x<0.1$,
$1\geq x>0.9$ are shaded, as LQCD does not provide data for these momentum
fractions due to systematic errors. For comparison, we plot the asymptotic
LFDA $\phi(x,\mu)\stackrel{{\scriptstyle\mu\to\infty}}{{=}}6x\bar{x}$.

Figure 2:  _Top panel_ : longitudinal and transverse distribution amplitudes,
$\phi^{\parallel}_{K^{*}}(x,\mu)$ and $\phi^{\perp}_{K^{*}}(x,\mu)$ for
$\mu=2$ GeV. _Bottom panel_ : Comparison of our predictions for the $K^{*}$
with those of QCDSR Ball:2007zt  and LQCD Hua:2020gnw , where we replaced
$x\rightarrow 1-x$ in Eq. (12). The shaded areas and error bands are as in
Figure 1.

##  III Results

We begin with the light vector mesons and determine the coefficients
$a_{n}^{\parallel,\perp}$ of their Gegenbauer expansion via a least-square fit
of $\epsilon_{\parallel,\perp}$ (12) with the four moments $\langle
x\rangle,\langle x^{2}\rangle,\langle x^{3}\rangle,\langle x^{4}\rangle$. We
report their values and those of the corresponding $a_{n}^{\parallel,\perp}$
of the $\rho$, $\phi$ and $K^{*}$ mesons in Table 1 and compare them to other
results in Table 2 with the moments $\langle\xi^{2m}\rangle_{\parallel,\perp}$
(13).

The LFDAs for the $\rho$ and $\phi$ mesons are compared in Figure 1 with the
prediction of a DSE-based calculation and the LDFA reconstructed with moments
from LQCD, respectively. We infer that the distributions follow the expected
pattern: both LFDAs are symmetric about the midpoint, $x=1/2$. However, the
$\phi_{\rho}^{\parallel,\perp}(x,\mu)$ distributions are broad while
$\phi_{\phi}^{\parallel,\perp}(x,\mu)$ tend to the asymptotic form
$\phi(x)\stackrel{{\scriptstyle\mu\to\infty}}{{=}}6x\bar{x}$. In addition, we
observe that $\phi^{\parallel}_{\rho}(x,\mu)$ is slightly broader than
$\phi^{\perp}_{\rho}(x,\mu)$, the origin of which are the different values of
$a^{\parallel}_{2}$ and $a^{\perp}_{2}$ in Table 1. It appears from Table 2
that our calculated $\langle\xi^{2m}\rangle_{\parallel}$ moments for the
$\rho$ meson are overall about 11% larger, whereas the values for
$\langle\xi^{2m}\rangle_{\perp}$ are in very good agreement with those of Ref.
Gao:2014bca  and the HERA fit Forshaw:2010py ; Forshaw:2012im .

In the case of the $\phi$-meson, we note that
$\phi^{\parallel}_{\phi}(x,\mu)\approx\phi^{\perp}_{\phi}(x,\mu)$ since
$a^{\parallel}_{2}\approx a^{\perp}_{2}$ and
$\alpha^{\parallel}\approx\alpha^{\perp}$. We remark that our results for
$\phi_{\phi}^{\parallel,\perp}(x,\mu)$ differ from those in Ref. Gao:2014bca
as can be inferred from Fig. 1. The reason for this, despite a like-minded BSE
approach, is that we use a larger strange-quark mass, $m_{s}=166$ MeV at
$\mu=2$ GeV. With a lower value of $m_{s}\approx 100$ MeV we find similar
distributions as in Ref. Gao:2014bca . However, we prefer to renormalize the
DSE with a larger strange mass as it results in a more consistent description
of the $K$, $K^{*}$ and $\phi$ mesons.

| $\langle x\rangle_{\parallel,\perp}$ | $\langle x^{2}\rangle_{\parallel,\perp}$ | $\langle x^{3}\rangle_{\parallel,\perp}$ | $\langle x^{4}\rangle_{\parallel,\perp}$ | $\alpha_{\parallel,\perp}$ | $\beta_{\parallel,\perp}$
---|---|---|---|---|---|---
$J/\Psi_{\parallel}$ | 0.500 | 0.274 | 0.159 | 0.097 | $4.549\pm 0.411$ | $0.081\pm 0.051$
$J/\Psi_{\perp}$ | 0.500 | 0.259 | 0.139 | 0.076 | $12.703\pm 1.931$ | $0.004\pm 0.710$
$D^{*}_{\parallel}$ | 0.694 | 0.511 | 0.396 | 0.315 | $0.531\pm 0.207$ | $2.460\pm 0.131$
$D^{*}_{\perp}$ | 0.742 | 0.589 | 0.471 | 0.389 | $0.094\pm 0.001$ | $3.073\pm 0.001$
$D^{*}_{s\parallel}$ | 0.627 | 0.418 | 0.294 | 0.217 | $2.582\pm 0.651$ | $2.263\pm 0.296$
$D^{*}_{s\perp}$ | 0.655 | 0.465 | 0.346 | 0.272 | $0.448\pm 0.305$ | $1.832\pm 0.136$
Table 3: Mellin moments $\langle x^{m}\rangle_{\parallel}$ and $\langle
x^{m}\rangle_{\perp}$ of the $J/\Psi$, $D^{*}$ and $D^{*}_{s}$ mesons. Fitting
these moments with their definitions in Eqs. (4) and (5) and the corresponding
$\phi_{V}^{\parallel,\perp}(x,\mu)$ parametrization (10) yields
$\alpha_{\parallel,\perp}$ and $\beta_{\parallel,\perp}$; the fit errors arise
in the minimization process.  Figure 3: Longitudinal and transverse
distributions, $\phi^{\parallel}_{J/\Psi}(x,\mu)$ and
$\phi^{\perp}_{J/\Psi}(x,\mu)$, reconstructed from the moments in Table 3 with
Eqs. (4), (5) and (10). The error bands reflect the uncertainties in the fit
parameters $\alpha_{\parallel,\perp}$ and $\beta_{\parallel,\perp}$ in Table
3.

We now turn our attention to the $K^{*}$ and present the longitudinal and
transverse LFDAs in Figure 2 where we juxtapose them with predictions from
LQCD and QCDSR. Notably, the longitudinal distribution is a concave, nearly
symmetric function of $x$, much broader than the asymptotic form, which is a
consequence of the smallness of the $a^{\parallel}_{1}$ coefficient. The
transverse LFDA, on the other hand, is asymmetric around the midpoint and its
maximum is located at $x=0.78$, which clearly indicates SU(3) flavor symmetry
breaking and that the strange valence quark carries a larger amount of meson
momentum. The asymmetric shape is due to the similarity of the Gegenbauer
coefficients, $a^{\perp}_{1}\approx a^{\perp}_{2}$ whereas
$a^{\parallel}_{1}\ll a^{\parallel}_{2}$, see Table 1. This is in agreement
with a recent calculation in LQCD, though in that study
$\phi^{\parallel}_{K^{*}}(x,\mu)$ tends toward the asymptotic distribution
Hua:2020gnw . In contrast to these findings, QCDSR predicts
$\phi^{\parallel}_{K^{*}}(x,\mu)\approx\phi^{\perp}_{K^{*}}(x,\mu)$
Ball:2007zt .

As we noted earlier, the heavier vector charmonium and charmed mesons require
a modified description of their LFDA (10) to fit the moments. We report these
moments, $\langle x^{m}\rangle_{\parallel}$ and $\langle
x^{m}\rangle_{\perp}$, for the $J/\Psi$, $D^{*}$ and $D^{*}_{s}$ in Table 3.
The distributions $\phi^{\parallel}_{J/\Psi}(x,\mu)$ and
$\phi^{\perp}_{J/\Psi}(x,\mu)$ we then reconstruct are plotted in Figure 3.
They are reminiscent of their pseudoscalar counterpart, i.e. the LFDA of the
$\eta_{c}$, which exhibits the same _convex-concave-convex_ functional
behavior and is more sharply peaked than the asymptotic LFDA Serna:2020txe .
It turns out that the longitudinal distribution is broader and less localized
as a function of $x$ than the transverse distribution, an observation also
made in Ref. Ding:2015rkn .

Figure 4: Longitudinal and transverse LFDAs of the $D^{*}$ and $D_{s}^{*}$
mesons at $\mu=2$ GeV; error bands as in Figure 3.

We conclude this section with a first prediction of the $D^{*}$ and
$D_{s}^{*}$ meson distribution amplitudes which we compute with the
projections in Eqs. (1) and (2) of the Bethe-Salpeter wave functions. The
latter are taken from Ref. El-Bennich:2021ldv ; see Table 4 therein for the
corresponding masses and weak decay constants. The distributions we
reconstruct from $\langle x^{m}\rangle_{\parallel}$ and $\langle
x^{m}\rangle_{\perp}$ listed in Table 3 are shown in Fig. 4. Clearly, in both
cases the LFDAs are asymmetric and $\phi^{\parallel}_{D^{*}}(x,\mu)$ and
$\phi^{\perp}_{D^{*}}(x,\mu)$ peak at about $x\approx 0.8-0.85$, while
$\phi^{\parallel}_{D^{*}_{s}}(x,\mu)$ and $\phi^{\perp}_{D^{*}_{s}}(x,\mu)$
reach their maximum in the range $x\approx 0.65-0.8$.

This is readily interpreted as the charm quark carrying most of the light-
front momentum in the $D^{*}$ meson, but less so in the $D^{*}_{s}$ meson.
Interestingly, the transverse distributions are more asymmetric and the charm
seems to carry a larger fraction of the meson momentum than in the
longitudinal distribution. Arguably, this observation generalizes our results
for the $K^{*}$, where the much smaller mass difference between the strange
and up quarks leads to an almost symmetric form of
$\phi^{\parallel}_{K^{*}}(x,\mu)$ and to a broad yet asymmetric function
$\phi^{\perp}_{K^{*}}(x,\mu)$.

##  IV Conclusion

We extracted the LFDAs of the $\rho$, $\phi$, $K^{*}$, $J/\Psi$, $D^{*}$ and
$D^{*}_{s}$ mesons from their Bethe-Salpeter wave functions, which we
calculated in Refs. Mojica:2017tvh ; El-Bennich:2021ldv , with two projections
onto the light front given by Eqs. (1) and (2). The transverse LFDA of the
$\rho$ meson is in very good agreement with that obtained in a similar DSE-BSE
approach Gao:2014bca  and with the HERA fit Forshaw:2010py ; Forshaw:2012im ,
while our longitudinal moments, $\langle\xi^{m}\rangle_{\parallel}$, are
generally about 11% larger than those in the literature.

We then presented the first calculation of the
$\phi^{\parallel}_{K^{*}}(x,\mu)$ and $\phi^{\perp}_{K^{*}}(x,\mu)$ within the
DSE-BSE framework and confirm the functional form found with LQCD simulations
Hua:2020gnw : while the longitudinal distribution of the $K^{*}$ is almost
symmetric about the midpoint $x=1/2$, the transverse distribution is broad and
slanted, which we interpret as the strange quark carrying the larger fraction
of the meson’s momentum. In the heavy meson sector, both LFDAs of the $J/\Psi$
are alike with that of the $\eta_{c}$, i.e. they are symmetric and narrow, yet
not merely concave distributions.

Last not least, we extended our studies in Ref. Serna:2020txe  to the
longitudinal and transverse LFDAs of the $D^{*}$ and $D_{s}^{*}$ mesons, a
first calculation of these distributions to our knowledge. Our findings are in
line with observations for the pseudoscalar $D$ and $D_{s}$ mesons
Serna:2020txe : the distributions are asymmetric and reach their maximum at
large momentum fractions, namely $x\approx 0.65-0.85$. In other words, the
charm quark is most likely to carry the largest fraction of the
$D_{(s)}^{*}$-momentum, and this is even more so the case for the transverse
distribution.

We remind that we provided all the analytic parametrizations of the LFDAs
discussed in this work and the parameters are found in Tables (1) and (3).
Therefore, the LFDAs of the $J/\Psi$ and $D^{*}_{(s)}$ mesons can readily be
used in diffractive vector-meson production and are are of interest to the
experimental program of the Electron-Ion Collider.

## Acknowledgments

We acknowledge helpful discussions with Peter Tandy and Minghui Ding. B.E. and
F.E.S. participate in the Brazilian network project _INCT-Física Nuclear e
Aplicações_ , no. 464898/2014-5. This work was supported by the São Paulo
Research Foundation (FAPESP), grant no. 2018/20218-4, and by the National
Council for Scientific and Technological Development (CNPq), grant no.
428003/2018-4. F.E.S. is a CAPES-PNPD postdoctoral fellow financed by grant
no. 88882.314890/2013-01.

## References

  * (1) A. V. Efremov and A. V. Radyushkin, Theor. Math. Phys. 42 (1980), 97-110 doi:10.1007/BF01032111
  * (2) G. P. Lepage and S. J. Brodsky, Phys. Rev. D 22 (1980), 2157 doi:10.1103/PhysRevD.22.2157
  * (3) L. Chang, I. C. Cloët, C. D. Roberts, S. M. Schmidt and P. C. Tandy, Phys. Rev. Lett. 111 (2013) no.14, 141802 doi:10.1103/PhysRevLett.111.141802
  * (4) K. Raya, L. Chang, A. Bashir, J. J. Cobos-Martínez, L. X. Gutiérrez-Guerrero, C. D. Roberts and P. C. Tandy, Phys. Rev. D 93 (2016) no.7, 074017 doi:10.1103/PhysRevD.93.074017
  * (5) K. Raya, A. Bashir and P. Roig, Phys. Rev. D 101 (2020) no.7, 074021 doi:10.1103/PhysRevD.101.074021
  * (6) M. Beneke and M. Neubert, Nucl. Phys. B 651 (2003), 225-248 doi:10.1016/S0550-3213(02)01091-X
  * (7) C. W. Bauer, D. Pirjol, I. Z. Rothstein and I. W. Stewart, Phys. Rev. D 72 (2005), 098502 doi:10.1103/PhysRevD.72.098502
  * (8) B. El-Bennich, A. Furman, R. Kaminski, L. Lesniak and B. Loiseau, Phys. Rev. D 74 (2006), 114009 doi:10.1103/PhysRevD.74.114009
  * (9) B. El-Bennich, A. Furman, R. Kaminski, L. Lesniak, B. Loiseau and B. Moussallam, Phys. Rev. D 79 (2009), 094005 [erratum: Phys. Rev. D 83 (2011), 039903] doi:10.1103/PhysRevD.83.039903
  * (10) O. Leitner, J. P. Dedonder, B. Loiseau and B. El-Bennich, Phys. Rev. D 82 (2010), 076006 doi:10.1103/PhysRevD.82.076006
  * (11) B. Pire and L. Szymanowski, Phys. Rev. Lett. 115 (2015) no.9, 092001 doi:10.1103/PhysRevLett.115.092001
  * (12) B. Pire and L. Szymanowski, Phys. Rev. D 96 (2017) no.11, 114008 doi:10.1103/PhysRevD.96.114008
  * (13) B. Pire, L. Szymanowski and J. Wagner, Phys. Rev. D 95 (2017) no.9, 094001 doi:10.1103/PhysRevD.95.094001
  * (14) B. Pire, L. Szymanowski and J. Wagner, Phys. Rev. D 104 (2021) no.9, 094002 doi:10.1103/PhysRevD.104.094002
  * (15) L. Chang, I. C. Cloët, J. J. Cobos-Martínez, C. D. Roberts, S. M. Schmidt and P. C. Tandy, Phys. Rev. Lett. 110 (2013) no.13, 132001 doi:10.1103/PhysRevLett.110.132001
  * (16) C. Shi, C. Chen, L. Chang, C. D. Roberts, S. M. Schmidt and H. S. Zong, Phys. Rev. D 92 (2015), 014035 doi:10.1103/PhysRevD.92.014035
  * (17) F. E. Serna, R. C. da Silveira, J. J. Cobos-Martínez, B. El-Bennich and E. Rojas, Eur. Phys. J. C 80 (2020) no.10, 955 doi:10.1140/epjc/s10052-020-08517-3
  * (18) M. Beneke and T. Feldmann, Nucl. Phys. B 592 (2001), 3-34 doi:10.1016/S0550-3213(00)00585-X
  * (19) J. R. Forshaw and R. Sandapen, JHEP 11 (2010), 037 doi:10.1007/JHEP11(2010)037
  * (20) J. R. Forshaw and R. Sandapen, Phys. Rev. Lett. 109 (2012), 081601 doi:10.1103/PhysRevLett.109.081601
  * (21) A. Bashir, L. Chang, I. C. Cloet, B. El-Bennich, Y. X. Liu, C. D. Roberts and P. C. Tandy, Commun. Theor. Phys. 58 (2012), 79-134 doi:10.1088/0253-6102/58/1/16
  * (22) F. Gao, L. Chang, Y. X. Liu, C. D. Roberts and S. M. Schmidt, Phys. Rev. D 90 (2014) no.1, 014011 doi:10.1103/PhysRevD.90.014011
  * (23) M. Ding, F. Gao, L. Chang, Y. X. Liu and C. D. Roberts, Phys. Lett. B 753 (2016), 330-335 doi:10.1016/j.physletb.2015.11.075
  * (24) P. Ball, V. M. Braun and A. Lenz, JHEP 08 (2007), 090 doi:10.1088/1126-6708/2007/08/090
  * (25) P. A. Boyle et al. [RBC and UKQCD], PoS LATTICE2008 (2008), 165 doi:10.22323/1.066.0165
  * (26) J. Hua et al. [Lattice Parton], Phys. Rev. Lett. 127 (2021) no.6, 062002 doi:10.1103/PhysRevLett.127.062002
  * (27) Y. Lu, D. Binosi, M. Ding, C. D. Roberts, H. Y. Xing and C. Xu, Eur. Phys. J. A 57 (2021) no.4, 115 doi:10.1140/epja/s10050-021-00427-6
  * (28) F. F. Mojica, C. E. Vera, E. Rojas and B. El-Bennich, Phys. Rev. D 96 (2017) no.1, 014012 doi:10.1103/PhysRevD.96.014012
  * (29) B. El-Bennich and F. E. Serna, PoS CHARM2020 (2021), 025 doi:10.22323/1.385.0025
  * (30) F. Gao, L. Chang and Y. x. Liu, Phys. Lett. B 770 (2017), 551-555 doi:10.1016/j.physletb.2017.04.077
  * (31) F. E. Serna and B. El-Bennich, PoS CHARM2020 (2021), 047 doi:10.22323/1.385.0047


# Improving Kinodynamic Planners for Vehicular Navigation
with Learned Goal-Reaching Controllers

Aravind Sivaramakrishnan, Edgar Granados, Seth Karten, Troy McMahon, Kostas E.
Bekris The authors are with the Dept. of Computer Science at Rutgers, NJ, USA.
Email:  {aravind.siva, eg585, sak295, tm799<EMAIL_ADDRESS>This work was
supported by NSF HDR TRIPODS award 1934924. The results do not reflect the
sponsor’s positions.

###### Abstract

This paper aims to improve the path quality and computational efficiency of
sampling-based kinodynamic planners for vehicular navigation. It proposes a
learning framework for identifying promising controls during the expansion
process of sampling-based planners. Given a dynamics model, a reinforcement
learning process is trained offline to return a low-cost control that reaches
a local goal state (i.e., a waypoint) in the absence of obstacles. By focusing
on the system’s dynamics and not knowing the environment, this process is
data-efficient and takes place once for a robotic system. In this way, it can
be reused in different environments. The planner generates online local goal
states for the learned controller in an informed manner to bias towards the
goal and consecutively in an exploratory, random manner. For the informed
expansion, local goal states are generated either via (a) medial axis
information in environments with obstacles, or (b) wavefront information for
setups with traversability costs. The learning process and the resulting
planning framework are evaluated for a first and second-order differential
drive system, as well as a physically simulated Segway robot. The results show
that the proposed integration of learning and planning can produce higher
quality paths than sampling-based kinodynamic planning with random controls in
fewer iterations and computation time.

##  I Introduction

This work focuses on improving the efficiency of sampling-based motion
planning for vehicular systems that exhibit significant dynamics using machine
learning. Planning motions for such systems can be challenging when there is
no access to a local planner and the basic primitive to explore the state
space is forward propagation of controls. In this context, tree sampling-based
kinodynamic planners have been developed [1], some of which achieve asymptotic
optimality (AO) by propagating randomly sampled controls during each iteration
[2, 3, 4, 5]. While random controls are desirable for analysis purposes, in
practice they result in low quality (i.e., high cost) trajectories.

Learning candidate controls for edge propagation from previous planning
experience is a promising avenue for improving the practical efficiency of
sampling-based kinodynamic planners. Prior work on the integration of machine
learning for planning has primarily focused on geometric motion planning
problems that don’t deal with dynamics [6, 7], or requires access to a two-
point boundary value problem (BVP) solver to construct a supervised learning
dataset [8]. The latter may not be easily available for some vehicular models,
especially for higher-fidelity physics-based models.

Recent work has shown that deep reinforcement learning (DRL) can be used to
learn a local controller (or a policy) that steers the robot between two
collision-free states. The learned policy can be then used as an edge
connection strategy for a sampling-based motion planner, such as a
Probabilistic Roadmap (PRM) [9] or a Rapidly-exploring Random Tree (RRT) [10].
This approach, however, requires a large amount of training data and manual
specification of reward parameters.

Figure 1: Visualizations of solution trajectories found by the proposed method
for the first-order (top-left, bottom-left over a cost-map) and second-order
(top-right) differential drive systems, and a Segway robot (bottom-right)
simulated using the Bullet physics engine.

A representation of the planning environment, either in the form of binary
occupancy grids, or cost maps that represent the cost of navigating across the
environment, has been used as input to neural network models, which can then
bias the planner towards shortest paths [11, 12, 13, 14]. Although it is
possible to use these learned sampling strategies as input to a learned local
planner, some of these samples may correspond to dynamically infeasible states
or inevitable collision states (ICS). Recent work has proposed learning a
sampling strategy for local goals that can be reached while obeying
kinodynamic constraints and avoiding collisions [15], but assumes a steering
function that can move the robot from one state to another.

This work proposes a learning process for a controller and integrates it with
a framework for sampling-based kinodynamic planning to improve performance.
The controller is trained offline using DRL in an empty environment to
generate high quality controls towards desirable local goals. During planning,
a node expansion procedure generates local goals biasing towards the global
goal that are passed to the learned controller, which outputs the best
control. The benefit of this two-step approach is that the learning process
has to only deal with the system’s dynamics and needs to be trained once per
system. In this way it does not require a separate training process for each
new environment that a robot is being deployed. The local goal selection and
the planning framework are responsible for dealing with different problem
instances and selecting a direction of progress towards the global goal. For
the first expansion out of a node of the sampling-based tree, an informed
approach is used where the local goal is selected so as to bias expansion
towards the ultimate planning goal. Consecutive expansions follow a more
exploratory strategy and select randomly local goals in the vicinity of the
tree node.

The evaluation uses a first-order and second-order differential drive system
as well as a physically simulated Segway robot across different environments.
The testing environments are different from the obstacle-free environment used
for training the controller. It shows that the proposed framework results in
very low cost solutions that are achieved faster relative to traditional
random propagation approaches.

##  II Problem Setup and Notation

Consider a system with state space $\mathbb{X}$ \- divided into collision-free
($\mathbb{X}_{f}$) and obstacle ($\mathbb{X}_{obs}$) subsets - and control
space $\mathbb{U}$, governed by the dynamics $\dot{x}=f(x,u)$ (where
$x\in\mathbb{X},u\in\mathbb{U}$). The process $f$ can be an analytical
ordinary differential equation (ODE) or modeled via a physics engine like
Bullet [16]. A plan is defined by a sequence of piecewise-constant controls
executed in order. A plan of length $T$ induces a trajectory
$\tau\in\mathcal{T}$, where $\tau:[0,T]\mapsto\mathbb{X}_{f}$. For a start
state $x_{0}\in\mathbb{X}_{f}$ and a goal set $X_{G}\subset\mathbb{X}_{f}$, a
feasible motion planning problem admits a solution trajectory of the form
$\tau(0)=x_{0},\tau(T)\in X_{G}$. Each solution trajectory has a cost $g$
according to a function $\texttt{cost}:\mathcal{T}\mapsto\mathbb{R}^{+}$,
which the motion planner must minimize. In this setup, $X_{G}$ is defined as
the set of all states that reach a specific state $x_{G}\in\mathbb{X}_{f}$
with some predefined tolerance.

##  III Sampling-Based Motion Planning Framework

We focus on sampling-based motion planners that only have access to a forward
propagation model of the system’s dynamics. This class of algorithms explore
$\mathbb{X}_{f}$ by incrementally constructing a tree by sampling
$\mathbb{U}$. Algorithm 1 outlines the high-level operation of a sampling-
based motion planner (Tree-SBMP) that builds a tree of states reachable from
$x_{0}$.

Each iteration of Tree-SBMP starts with the selection of an existing tree node
to expand that corresponds to a robot state (Line 3). A candidate control is
considered for the selected node, and an end state is obtained by forward
propagation of the system dynamics (Lines 4-5). The resulting edge is added to
the tree if it is not in collision (Lines 6-7). By varying how these
operations are implemented, different algorithms can be obtained. After $N$
iterations (or a preset time limit), if the tree has states in $X_{G}$, the
best found solution trajectory according to $g$ can be returned.

1 $T\leftarrow\\{x_{0}\\}$;

2 while  _termination condition is not met_ do

3 $x_{select}\leftarrow\texttt{SELECT-NODE}(T)$;

4 $u\leftarrow\texttt{EXPAND-NODE}(x_{select})$;

5 $x_{new}\leftarrow\texttt{PROPAGATE}(x_{select},u)$;

6 if  _$(x_{select}\rightarrow x_{new})\in\mathbb{X}_{f}$_ then

7 EXTEND-TREE($T,x_{select}\rightarrow x_{new}$);

8 end if

9

10 end while

Algorithm 1 Tree-SBMP($\mathbb{X}_{f},\mathbb{X}_{o},\mathbb{U},x_{0},X_{G}$)

This work presents a version of EXPAND-NODE (Fig. 2) that selects a local goal
in the environment and uses a learned controller to predict a control that
makes progress to the local goal. The proposed node expansion framework is
integrated into the DIRT algorithm [4], which also achieves AO properties.
DIRT is selected due to its informed implementation of SELECT-NODE that biases
selection towards high quality nodes close to the unexplored parts of the
state space.

Figure 2: Illustration of the proposed expansion function.

DIRT also evaluates multiple candidate controls as part of EXPAND-NODE. This
Blossom expansion procedure ranks the resulting trajectories and processes
them in order so as to prioritize expansions that bring the system closer to
the goal according to a heuristic cost-to-go.

##  IV Proposed Method

The proposed pipeline has two stages. Offline: For each system, a controller
is trained using DRL to reach a goal state in an obstacle-free environment.
For every planning environment, a representation in the form of an obstacle
map or of a cost map is assumed available. Online:  Upon planner
initialization, a vector field is constructed for the given planning problem
to bias the planner. At every iteration, a local goal state is generated as
input to the learned controller. The controller then outputs a candidate
control to be propagated towards the local goal state. The local goal is
generated either in an exploitative manner using the goal-biasing vector field
or in an exploratory manner by using random sampling.

###  IV-A Learning the controller

Given the current state of the system $x\in\mathbb{X}$, the aim is to train a
controller $u_{RLC}=\pi(x,x_{G})$ to reach the goal set $X_{G}$ in the absence
of obstacles. Let $\rho_{\pi}$ be the distribution of trajectories induced by
the controller $\pi$, and $x_{0}$ the initial state sampled from an initial
state distribution $\rho(x_{0})$. The controller can be optimized using a DRL
algorithm that maximizes the expected return

$\mathcal{L}(\pi)=\mathbb{E}_{\tau\sim\rho_{\pi}}[\sum_{t=0}^{\infty}\gamma^{t}r(x_{t},x_{G})],\vspace{-0.05in}$

where $\gamma\in[0,1]$ is a discount factor and
$r:\mathbb{X}\times\mathbb{X}\mapsto\\{-1,0\\}$ is a bounded reward function
defined as $r(x_{t},x_{G})=0$ if $x_{t}\in X_{G}$ and $r(x_{t},x_{G})=-1$
otherwise. Alternative DRL algorithms can be used to maximize
$\mathcal{L}(\pi)$. This work employs Hindsight Experience Replay (HER) [17]
to facilitate learning from sparse rewards.

Figure 3: Off-policy DRL algorithm with online data collection. The replay
buffer $\mathcal{D}$ is implemented using Hindsight Experience Replay (HER).
The controller update step depends on the choice of the DRL algorithm used to
maximize the expected return.

As the focus is mainly on vehicular systems, $x_{G}$ is defined relative to
the robot’s current state. For the state variables of the system that are
transformation-invariant (e.g., current position and steering angle), the
initial state of the system is fixed at the origin. The values of the initial
state variables that are not transformation-invariant (e.g., velocities of the
wheels) are sampled uniformly within bounds. This re-parameterization reduces
data requirements while learning the controller.

Comparison alternative: A Supervised Learned Controller (SLC),
$u_{SLC}=\pi(\delta x)$ is also evaluated. It learns the inverse of a
discretized set of forward dynamics. SLC acts as a locally optimal surrogate
for the system’s dynamics. The control-duration space is discretized to form a
discrete set of plans $\mathcal{U}_{LC}$ and all plans in $\mathcal{U}_{LC}$
are forward propagated from a fixed initial state. The resulting end states
form a reachability region $\mathcal{X}$. A dataset is obtained by using a
nearest neighbor data structure to store the best plan that reaches a
neighborhood in $\mathcal{X}$ for every initial state. A neural network is
then trained to minimize the mean squared error (MSE) between the plan in the
dataset and the plan output by SLC given $\delta x$, which is defined as the
difference in states between an initial state $x_{0}\in\mathbb{X}$ and a local
goal state $x_{G}\in\mathcal{X}$.

###  IV-B Generating informed local goal states among obstacles

A two-dim. workspace $\mathbb{T}$ can be divided into two disjoint subsets:
the collision-free subset $\mathbb{T}_{f}$ and the obstacle subset
$\mathbb{T}_{o}$ containing all $n$ obstacles. The boundary of the environment
is also considered as an obstacle. Define points $p$ as projections of states
$x$ from $\mathbb{X}$ to $\mathbb{T}$. Then, let $p_{f}\in\mathbb{T}_{f}$ be
the projection in $\mathbb{T}$ of each tree state/node $x_{select}$ selected
for propagation. Then, the objective is to generate a local goal point
$p_{lg}\in\mathbb{T}_{f}$ so as to perform the propagation with the learned
controller. Given the global goal $p_{G}$, $p_{lg}$ is generated by a local
goal function $\phi(p_{f},p_{G})$ so that: i) $h^{*}(p_{lg})<h^{*}(p_{f})$
where $h^{*}$ is a cost-to-go function given $p_{G}$; and ii) $p_{lg}$
maximizes clearance from obstacles.

Figure 4: An example $\mathcal{MA}$ highlighted in red.

The local goal function $\phi(p_{f},p_{G})$ is implemented using a medial axis
$\mathcal{MA}\subset\mathbb{T}_{f}$. Given an obstacle surface
$\mathcal{O}_{i}\subset\mathbb{T}_{o}$ (e.g., an edge of a polygonal obstacle
in 2D) there exists a subset $P_{\mathcal{O}_{i}}\subset\mathbb{T}_{f}$ of all
free points, such that $\mathcal{O}_{i}$ is the closest obstacle surface. The
medial axis is the intersection of all such subsets:
$\mathcal{MA}=\bigcap_{i=0}^{n}P_{\mathcal{O}_{i}}$. To compute the medial
axis, the environment is discretized. The computation is performed once for
every planning environment (Fig 4).

Figure 5: (Left) Closeup of $\mathcal{MA}$ from Fig. 4. Given
$p_{f}\in\mathbb{T}_{f}$, the vectors $U_{rep}(p_{f})$ and $U_{att}(p_{f})$
give rise to $U(p_{f})$. (Right) The normalized vector field for a given goal
(green circle).

For a given $p_{f}\in\mathbb{T}_{f}$, 3 vectors are computed (Fig 5 left).
First, a repulsive vector $U_{rep}(p_{f})$ that points $p_{f}$ away from the
closest obstacle towards the nearest point on $\mathcal{MA}$. An attractive
vector $U_{att}(p_{f})$ that points to $p_{MA}^{*}\in(\mathcal{MA}\cup p_{G})$
such that: (a) $p_{f}$ has line of sight (with clearance) to $p_{MA}^{*}$, and
(b) $h^{*}(p_{MA}^{*})$ is minimum, i.e., makes the most progress to the goal
among all visible points from $p_{f}$. Both vectors are used to compute
$U(p_{f})=w\cdot U_{rep}(p_{f})+(1-w)\cdot U_{att}(p_{f})$, where
$w=|U_{rep}(p_{f})|/(|U_{rep}(p_{f})|+|U_{att}(p_{f})|))$. The magnitude of
$U$ is set so it reaches $\mathcal{MA}$ when applied at $p_{f}$. $U$ is the
final output of $\phi(p_{f},p_{G})$. The normalized version of $U$ is shown in
Fig 5 right. It is computed once offline given a discretization of the
environment.

###  IV-C Generating informed local goal states given a costmap

To use the learned controller $\pi$ in the traversability cost environment,
the environment is discretized on a grid in $\mathbb{T}$. For every
$p\in\mathbb{T},$ the traversability cost $c_{tv}(p)\in[0,1]$ is defined,
where $c_{tv}(p)<1$ is traversable (Fig. 6, left). The higher the cost, the
longer it takes for the robot to traverse. Equation (1) computes the path cost
of a trajectory $\tau$ of length $T$ as:

$\texttt{cost}(\tau)=\int_{0}^{T}\exp(Kc_{tv}(\tau(t)))dt\vspace{-.2in}$ (1)

Figure 6: (Left) Cost map of a planning environment where blue/red regions are
easy/hard to traverse. (Right) The $U$ vector field for a given goal state
(green circle).

Offline, once per planning problem, a modified version of the any-angle path
planning algorithm $\text{Theta}^{*}$ [18] is computed over the grid to create
a vector field (Fig. 6, right). Instead of checking for an obstacle-free line
of sight when smoothing the solution path, the modified version checks for a
line of decreasing cost of a minimum length. For a point $p\in\mathbb{T}$, the
vector $U(p)$ points to the next node $p^{*}$ along the lowest cost path to
the goal $p_{G}$ according to $\text{Theta}^{*}$. Online, the end point of $U$
is returned as the output of $\phi(p,p_{G}).$

##  V Results

The systems considered in our evaluation are an analytically simulated first-
order ($\texttt{dim}(\mathbb{X})=3,\texttt{dim}(\mathbb{U})=2$) and second-
order differential drive vehicles
($\texttt{dim}(\mathbb{X})=5,\texttt{dim}(\mathbb{U})=2$) and a Bullet-
simulated Segway robot
($\texttt{dim}(\mathbb{X})=3,\texttt{dim}(\mathbb{U})=2$).

###  V-A Controller training evaluation

The following algorithms were tested to learn $u_{RLC}=\pi(x,x_{G})$: Deep
Q-Networks (DQN) [19], Twin-Delayed Deep Deterministic Policy Gradient (TD3)
[20] and Soft Actor-Critic (SAC) [21]. For DQN, the discretized action space
$\mathcal{U}$ consists of $3^{\texttt{dim}(\mathbb{U})}$ controls. For each
dimension of the control space, the control can take either the minimum, mid
or maximum value. The training pipeline was implemented on a Google Colab CPU
instance. The DRL algorithm implementations were derived from Gym [22] and
Stable Baselines 3 [23].

Figure 7: Training performance of different controllers for the first-order
and second-order differential drive systems. The baseline is chosen to be the
highest possible cost that a goal-reaching trajectory can incur in an
obstacle-free environment.

Each controller is trained for 200 thousand forward propagation steps, and
their performance is evaluated using the average time to reach the goal
(trajectory cost). Fig. 7 shows the results for the analytical systems. The
costs are averaged over 30 trials, where each trial corresponds to 100
randomly sampled goals. DQN does not converge to a controller that finds low
cost goal-reaching trajectories due to the restrictive nature of the control
space. Although SAC and TD3 exhibit nearly similar performance after parameter
tuning, the best-performing SAC controller is selected for the proposed
expansion framework since its performance is closest to the baseline.
Subsequently, SAC+HER is also used to train the controller for the Segway
system.

(a) Results for first-order differential drive system in the city
environments.

(b) Results for second-order differential drive system in the city
environments.

(c) Results for first-order differential drive system in the cost map
environments.

(d) Results for Bullet-simulated Segway system.

Figure 8: For left and middle-left graphs, higher success rate at earlier
time/iteration is better. For right and middle-right graphs, lower path cost
is better. Each problem instance is run 30 times to account for different
random seeds.

###  V-B Planner integration evaluation

The integrated planning and learning solution is evaluated with three
different expansion functions using a Blossom expansion of five controls each.
(a) Random uses a blossom expansion of random controls in $\mathbb{U}$. (b)
SLC uses the supervised learning process with a local goal provided by $\phi$
as input for the first control; for subsequent controls, a random local goal
in $\mathcal{X}$ is provided as input to the supervised learned controller.
(c) RLC uses SAC with local goal provided by $\phi$ as input for the first
control; for subsequent controls, a random local goal in $\mathbb{X}$ is
provided as input to SAC. For the Segway system, since each forward
propagation is expensive, a Blossom number of 1 is used.

Across the experiments, the following metrics are measured for every planner:
(1) Average normalized cost of solutions found over time/number of planner
iterations, and (2) Ratio of experiments for which solution was found over
time/number of planner iterations. To account for the difficulty of different
planning problems, path costs are normalized by dividing by the best path cost
found for a problem across any planner. For the non-cost map problems, the
path cost is defined to be the duration of the solution plan. For the cost map
problems, the path cost is computed according to Equation 1 with $K=4.68$.

For the first and second-order systems, the planner performance is evaluated
on three maps from the city benchmarks [24]. Planner performance with the
first-order system is also evaluated on three cost maps from [25]. For each
map, ten different problem instances are generated by sampling collision-free
start and goal states that have a minimum separation of $M$ pixels. For the
Segway robot, three environments are defined and planner performance is
evaluated across all. Example solutions found by RLC are given in Fig 1. For
the city environments, the search trees are visualized.

In the city environments, RLC finds the lowest cost solutions very quickly for
both first (Fig 8(a)) and second-order (Fig 8(b)) systems. For the first-order
system, RLC finds solutions in fewer iterations than Random or SLC. For the
second-order system, SLC finds solutions in fewer iterations than RLC, but RLC
finds better quality solutions. In the cost map environments (Fig 8(c)), RLC
finds the lowest cost solutions quickly. For the same time budget, RLC and
Random find solutions to more problems compared to SLC, while RLC does so in
fewer iterations than Random. For the Segway (Fig 8(d)), RLC finds lower cost
solutions than both Random and SLC, and also finds more solutions in fewer
iterations / lesser time. Across experiments, both RLC and SLC achieve a high
success rate in fewer iterations than Random. The solutions obtained by RLC,
however, are of significantly lower cost than SLC and Random.

(a)

(b)

(c) Figure 9: Results of the ablation experiments for the second-order system.

###  V-C Ablation Studies

These experiments evaluate various choices for the planning process in the
case of the second-order vehicle: (a) should the first propagation be informed
given medial axis information, (b) should consecutive propagations use random
controls (Random), or use a random local goal and supervised learning (SLC) or
use a random local goal and SAC (Rand Goal), and (c) what is a good blossom
number. The planner using Medial Axis Goal finds lower cost solutions in fewer
iterations than the one using Random Goal (Fig. 9(a)). The planner using SLC
finds solutions in fewer iterations than using Rand Goal but RandGoal finds
higher quality solutions overall (Fig. 9(b)). RandGoal also remains
competitive in terms of the number of solutions found compared to Random. A
smaller Blossom number finds solutions faster (Fig 9(c)). The quality of
solutions found, however, improves considerably when a higher Blossom number
is used. The best-performing choices across these experiments are used for RLC
evaluated in Section V-B.

##  VI Conclusion

This paper proposes a learning-based node expansion method that improves the
path quality and computational efficiency of sampling-based kinodynamic
planners for vehicular systems. It consists of a learned controller that can
be trained with small data requirements for different robotic systems, and a
selection process of local goal states that is computed once for every
planning problem.

The inference time of the learned controller can further be improved by
storing it upon training to an integrated circuit, such as an FPGA.
Furthermore, the local goal state selection can also be learned in conjunction
with the controller. Similarly, the proposed method can also be integrated
with other learning primitives for planning, such as a learned selection
procedure. Next steps include simulations with more complex systems and
experiments on real vehicular systems.

## References

  * [1] S. M. LaValle and J. J. Kuffner Jr, “Randomized kinodynamic planning,” _IJRR_ , vol. 20, no. 5, pp. 378–400, 2001.
  * [2] Y. Li, Z. Littlefield, and K. E. Bekris, “Asymptotically optimal sampling-based kinodynamic planning,” _IJRR_ , vol. 35, no. 5, pp. 528–564, 2016.
  * [3] K. Hauser and Y. Zhou, “Asymptotically optimal planning by feasible kinodynamic planning in a state–cost space,” _T-RO_ , vol. 32, no. 6, pp. 1431–1443, 2016.
  * [4] Z. Littlefield and K. E. Bekris, “Efficient and asymptotically optimal kinodynamic motion planning via dominance-informed regions,” in _IROS_ , 2018\.
  * [5] M. Kleinbort, E. Granados, K. Solovey, R. Bonalli, K. E. Bekris, and D. Halperin, “Refined analysis of asymptotically-optimal kinodynamic planning in the state-cost space,” in _ICRA_ , 2020.
  * [6] A. H. Qureshi, Y. Miao, A. Simeonov, and M. C. Yip, “Motion planning networks: Bridging the gap between learning-based and classical motion planners,” _T-RO_ , vol. 37, pp. 48–66, 2021.
  * [7] T. Jurgenson and A. Tamar, “Harnessing reinforcement learning for neural motion planning.” in _R: SS_ , 2019.
  * [8] Y. Li, R. Cui, Z. Li, and D. Xu, “Neural network approximation based near-optimal motion planning with kinodynamic constraints using RRT,” _IEEE Transactions on Industrial Electronics_ , vol. 65, pp. 8718–8729, 2018\.
  * [9] A. Faust, O. Ramírez, M. Fiser, K. Oslund, A. Francis, J. O. Davidson, and L. Tapia, “PRM-RL: Long-range robotic navigation tasks by combining reinforcement learning and sampling-based planning,” _ICRA_ , 2018.
  * [10] H.-T. L. Chiang, J. Hsu, M. Fiser, L. Tapia, and A. Faust, “RL-RRT: Kinodynamic motion planning via learning reachability estimators from RL policies,” _RA-L_ , vol. 4, no. 4, pp. 4298–4305, 2019.
  * [11] B. Ichter, J. Harrison, and M. Pavone, “Learning sampling distributions for robot motion planning,” in _ICRA_ , 2018.
  * [12] R. Kumar, A. Mandalika, S. Choudhury, and S. Srinivasa, “LEGO: Leveraging experience in roadmap generation for sampling-based planning,” _IROS_ , pp. 1488–1495, 2019.
  * [13] J. Wang, W. Chi, C. Li, C. Wang, and M. Q.-H. Meng, “Neural RRT*: Learning-based optimal path planning,” _IEEE T-ASE_ , vol. 17, no. 4, pp. 1748–1758, 2020.
  * [14] B. Ichter, E. Schmerling, T. W. E. Lee, and A. Faust, “Learned critical probabilistic roadmaps for robotic motion planning,” in _ICRA_ , 2020, pp. 9535–9541.
  * [15] J. J. Johnson, L. Li, F. Liu, A. H. Qureshi, and M. C. Yip, “Dynamically constrained motion planning networks for non-holonomic robots,” in _IROS_ , 2020, pp. 6937–6943.
  * [16] E. Coumans and Y. Bai, “PyBullet, a python module for physics simulation for games, robotics and machine learning,” http://pybullet.org, 2016–2019.
  * [17] M. Andrychowicz, D. Crow, A. Ray, J. Schneider, R. H. Fong, P. Welinder, B. McGrew, J. Tobin, P. Abbeel, and W. Zaremba, “Hindsight experience replay,” in _NeurIPS_ , 2017.
  * [18] A. Nash, K. Daniel, S. Koenig, and A. Feiner, “Theta*: Any-angle path planning on grids,” in _AAAI_ , 2007.
  * [19] V. Mnih, K. Kavukcuoglu, D. Silver, A. A. Rusu, J. Veness, M. G. Bellemare, A. Graves, M. Riedmiller, A. K. Fidjeland, G. Ostrovski, _et al._ , “Human-level control through deep reinforcement learning,” _Nature_ , vol. 518, no. 7540, pp. 529–533, 2015.
  * [20] S. Fujimoto, H. Hoof, and D. Meger, “Addressing function approximation error in actor-critic methods,” in _ICML_ , 2018.
  * [21] T. Haarnoja, A. Zhou, P. Abbeel, and S. Levine, “Soft actor-critic: Off-policy maximum entropy deep reinforcement learning with a stochastic actor,” in _ICML_ , 2018.
  * [22] G. Brockman, V. Cheung, L. Pettersson, J. Schneider, J. Schulman, J. Tang, and W. Zaremba, “OpenAI gym,” 2016.
  * [23] A. Raffin, A. Hill, M. Ernestus, A. Gleave, A. Kanervisto, and N. Dormann, “Stable baselines3,” https://github.com/DLR-RM/stable-baselines3, 2019\.
  * [24] N. Sturtevant, “Benchmarks for grid-based pathfinding,” _Transactions on Computational Intelligence and AI in Games_ , vol. 4, no. 2, pp. 144 – 148, 2012\. [Online]. Available: http://web.cs.du.edu/ sturtevant/papers/benchmarks.pdf
  * [25] L. Jaillet, J. Cortés, and T. Siméon, “Transition-based rrt for path planning in continuous cost spaces,” in _IROS_. IEEE, 2008, pp. 2145–2150.