# Application of coherent density fluctuation model to study the nuclear
matter properties of finite nuclei within relativistic mean-field formalism

Ankit Kumar1,2<EMAIL_ADDRESS>H. C. Das1,2 Manpreet Kaur1,2 M. Bhuyan3,4
<EMAIL_ADDRESS>S. K. Patra1,2 1Institute of Physics, Sachivalya Marg,
Bhubaneswar-751005, India 2Homi Bhabha National Institute, Training School
Complex, Anushakti Nagar, Mumbai 400094, India 3Department of Physics, Faculty
of Science, University of Malaya, Kuala Lumpur 50603, Malaysia 4Institute of
Research and Development, Duy Tan University, Da Nang 550000, Vietnam

###### Abstract

We obtained a density-dependent analytical expression of binding energy per
nucleon for different neutron-proton asymmetry of the nuclear matter (NM) with
a polynomial fitting, which manifests the results of effective field theory
motivated relativistic mean-field (E-RMF) model. This expression has the edge
over the Br$\ddot{u}$ckner energy density functional [Phys. Rev. 171, 1188
(1968)] since it resolves the Coster-Band problem. The NM parameters like
incompressibility, neutron pressure, symmetry energy, and its derivatives are
calculated using the acquired expression of energy per nucleon. Further, the
weight function calculated by E-RMF densities are folded with calculated NM
parameters within coherent density fluctuation model to find the properties of
closed/semi-closed-shell even-even 16O, 40Ca, 48Ca, 56Ni, 90Zr, 116Sn, and
208Pb nuclei. The values obtained for the neutron pressure $P^{A}$, symmetry
energy $S^{A}$ and its derivative $L_{sym}^{A}$ known as slope parameter, lie
within a narrow domain whereas there is a large variation in isoscalar
incompressibility $K^{A}$ and surface incompressibility $K_{sym}^{A}$ while
moving from light to heavy nuclei. The sizable variation in $K^{A}$ and
$K_{sym}^{A}$ for light and heavy nuclei depicts their structural dependence
due to the peculiar density distribution of each nucleus. A comparison of
surface quantities calculated in the present work has also been made with ones
obtained via Br$\ddot{u}$ckner energy density functional.

###### pacs:

##  I Introduction

The correlations among the nuclear matter (NM) and finite nuclei in terms of
symmetry energy and its coefficients play a crucial role not only in nuclear
physics but also in astrophysics. The isospin dependence of symmetry energy
imparts information about the isovector component of the nuclear interaction,
which is directly connected with the skin thickness of the nuclei. Eventually,
different studies such as the island of stability of exotic nuclei, the
dynamics of heavy-ion collisions, dipole polarizability, properties of neutron
star (NS), core-collapse of massive compact stars and the nucleosynthesis
process through neutrino convection at high-density region upon the symmetry
energy and its coefficients Lattimer and Prakash (2000, 2001); Steiner  _et
al._ (2005); Li  _et al._ (2008); Lynch  _et al._ (2009); Tsang  _et al._
(2012); Horowitz  _et al._ (2014); Antonov  _et al._ (2016); Bhuyan  _et al._
(2018); Kaur  _et al._ (2020a, b); Naik  _et al._ (2019). Therefore, it is
indispensable to determine the symmetry energy and its coefficients for finite
nuclei. Recently, many efforts have been made on theoretical as well as
experimental fronts to probe the isospin dependence of symmetry energy and its
coefficients, which is an ultimate bridge between finite nuclei and infinite
NM Steiner  _et al._ (2005); Tsang  _et al._ (2012); Horowitz  _et al._
(2014); Antonov  _et al._ (2016); Bhuyan  _et al._ (2018); Kaur  _et al._
(2020a, b); Naik  _et al._ (2019). Moreover, in some recent works, it is
established that the kink in the symmetry energy of finite nuclei over the
isotopic chain infers the appearance of shell/sub-shell closures Antonov  _et
al._ (2016); Bhuyan  _et al._ (2018); Kaur  _et al._ (2020a, b).

In recent works Bhuyan  _et al._ (2018); Quddus  _et al._ (2020), the
Br$\ddot{u}$ckner energy density functional Brueckner  _et al._ (1968, 1969)
has been used within the coherent density fluctuation model (CDFM) to
calculate the properties of nuclei. This functional had been fitted to the
kinetic and potential energy parts to get the analytical expression of binding
energy per particle E/A in the local density approximation (LDA) of the
Thomas-Fermi approach. It is relevant to point out that Br$\ddot{u}$ckner
energy functional does not respect the “Coester-Band” i.e. in Thomas-Fermi
approach NM saturates at $\rho\sim 0.2$ fm-3 instead of $\rho\sim 0.15$ fm-3
Coester  _et al._ (1970); Brockmann and Machleidt (1990). To have some
meaningful correlations while extrapolating to higher densities, the nuclear
equation of state (EOS) must satisfy the nuclear saturation properties. To
address this problem, we have fitted the NM saturation plots for different
values of asymmetry parameter, obtained using effective field theory motivated
relativistic mean-field (E-RMF) model Kumar  _et al._ (2018, 2017) with
standard NL3 and recently developed G3 parameter sets for the first time. The
different NM parameters such as incompressibility, symmetry energy, and its
derivatives are obtained employing the derived expression of E-RMF density
functional. Subsequently, theses NM parameters are used along with E-RMF
densities within the CDFM to find the corresponding quantities for magic/semi-
magic 16O, 40Ca, 48Ca, 56Ni, 90Zr, 116Sn, and 208Pb nuclei.

The paper is organized as follows: The E-RMF approach and the fitting
procedure to get the analytical polynomial expression of $E/A$ is discussed in
Sec. II. The CDFM is also discussed in this section. Sec. III is assigned to
the discussion of the results obtained from the calculations. A brief summary
and conclusions are presented in Sec. IV.

##  II Effective field theory motivated relativistic mean field (E-RMF) Model

In this section, we briefly describe the formalism of recently developed E-RMF
model. The E-RMF Lagrangian density is constructed by taking the interactions
of isoscalar (scalar $\sigma$, vector $\omega$) and isovector (scalar
$\delta$, vector $\rho$) mesons with nucleons and among themselves. The E-RMF
Lagrangian is discussed in Refs. Kumar  _et al._ (2018, 2017); Furnstahl  _et
al._ (1996, 1997). The E-RMF is considered to be one of the most successful
model to reproduce the ground state properties of not only $\beta-$stable
nuclei but also predicts quite reasonably the properties of drip-lines and
superheavy nuclei Kumar  _et al._ (2018, 2017). During last few decades, the
application of this formalism to nuclear astrophysics is at forefront. It
predicts the structure of NS and explains the tidal deformability
satisfactorily Malik  _et al._ (2018). The energy density functional for a
nucleon-meson interacting system is given as Kumar  _et al._ (2018):

$\displaystyle{\cal E}({r})$ $\displaystyle=$
$\displaystyle\sum_{\alpha=p,n}\varphi_{\alpha}^{\dagger}({r})\Bigg{\\{}-i\mbox{\boldmath$\alpha$}\\!\cdot\\!\mbox{\boldmath$\nabla$}+\beta\bigg{[}M-\Phi(r)-\tau_{3}D(r)\bigg{]}+W({r})+\frac{1}{2}\tau_{3}R({r})+\frac{1+\tau_{3}}{2}A({r})-\frac{i\beta\mbox{\boldmath$\alpha$}}{2M}\\!\cdot\\!\bigg{(}f_{\omega}\mbox{\boldmath$\nabla$}W({r})$
(1)
$\displaystyle+\frac{1}{2}f_{\rho}\tau_{3}\mbox{\boldmath$\nabla$}R({r})\bigg{)}\Bigg{\\}}\varphi_{\alpha}(r)+\left(\frac{1}{2}+\frac{\kappa_{3}}{3!}\frac{\Phi({r})}{M}+\frac{\kappa_{4}}{4!}\frac{\Phi^{2}({r})}{M^{2}}\right)\frac{m_{s}^{2}}{g_{s}^{2}}\Phi^{2}({r})-\frac{\zeta_{0}}{4!}\frac{1}{g_{\omega}^{2}}W^{4}({r})+\frac{1}{2g_{s}^{2}}\left(1+\alpha_{1}\frac{\Phi({r})}{M}\right)\bigg{(}\mbox{\boldmath$\nabla$}\Phi({r})\bigg{)}^{2}$
$\displaystyle-\frac{1}{2g_{\omega}^{2}}\left(1+\alpha_{2}\frac{\Phi({r})}{M}\right)\bigg{(}\mbox{\boldmath$\nabla$}W({r})\bigg{)}^{2}-\frac{1}{2}\left(1+\eta_{1}\frac{\Phi({r})}{M}+\frac{\eta_{2}}{2}\frac{\Phi^{2}({r})}{M^{2}}\right)\frac{m_{\omega}^{2}}{g_{\omega}^{2}}W^{2}({r})-\frac{1}{2e^{2}}\bigg{(}\mbox{\boldmath$\nabla$}A({r})\bigg{)}^{2}-\frac{1}{2g_{\rho}^{2}}\bigg{(}\mbox{\boldmath$\nabla$}R({r})\bigg{)}^{2}$
$\displaystyle-\frac{1}{2}\left(1+\eta_{\rho}\frac{\Phi({r})}{M}\right)\frac{m_{\rho}^{2}}{g_{\rho}^{2}}R^{2}({r})-\Lambda_{\omega}\bigg{(}R^{2}(r)\times
W^{2}(r)\bigg{)}+\frac{1}{2g_{\delta}^{2}}\left(\mbox{\boldmath$\nabla$}D({r})\right)^{2}+\frac{1}{2}\frac{{m_{\delta}}^{2}}{g_{\delta}^{2}}D^{2}(r)\;,$

Here, $\Phi$, $W$, $R$ and $D$ are the re-defined fields for $\sigma$,
$\omega$, $\rho$ and $\delta$ mesons given as $\Phi=g_{s}\sigma_{0}$,
$W=g_{\omega}\omega_{0}$, $R=g_{\rho}\vec{\rho_{0}}$ μ and
$D=g_{\delta}\delta_{0}$, respectively. $M$, $m_{\sigma}$, $m_{\omega}$,
$m_{\rho}$ and $m_{\delta}$ are the masses of nucleon, $\sigma$, $\omega$,
$\rho$ and $\delta$ mesons, respectively. From Eq. (1), we obtain the energy
density ${\cal{E}}_{nucl.}$ Kumar  _et al._ (2018, 2017) by considering that
the exchange of mesons create an uniform field, where the nucleon oscillates
in a simple harmonic motion. From the E-RMF energy density, the equation of
motions for the mesons and the nucleons are derived using the Euler-Lagrange
equation. A set of coupled differential equations are obtained and solved
self-consistently Kumar  _et al._ (2018). The scalar and vector densities,

$\displaystyle\rho_{s}(r)$ $\displaystyle=$
$\displaystyle\sum_{\alpha}\varphi_{\alpha}^{\dagger}({r})\beta\varphi_{\alpha},$
(2) $\displaystyle\rho_{v}(r)$ $\displaystyle=$
$\displaystyle\sum_{\alpha}\varphi_{\alpha}^{\dagger}({r})\tau_{3}\varphi_{\alpha},$
(3)

are evaluated from the converged solutions within spherical harmonics. The
vector density $\rho_{v}(r)$ is further used within CDFM to find out the
weight function $|F(x)|^{2}$, which is an important quantity to calculate the
incompressibility ($K^{A}$), symmetry energy ($S^{A}$), neutron pressure
($P^{A}$) and surface symmetry coefficient ($K_{sym}^{A}$) for the
closed/semi-closed-shell even-even nuclei.
The expression for energy density of infinite and isotropic NM are obtained
from the energy-momentum tensor:

$\displaystyle T_{\mu\nu}=$
$\displaystyle\sum_{i}\partial_{\nu}\phi_{i}\frac{\partial\cal
L}{\partial(\partial^{\mu}\phi_{i})}-g_{\mu\nu}{\cal L}.$ (4)

The zeroth component of the energy-momentum tensor $T_{00}$ gives the energy
density of the system as a function of baryon density as:

$\displaystyle{\cal{E}}(k)_{nucl.}$ $\displaystyle=$
$\displaystyle\frac{2}{(2\pi)^{3}}\int
d^{3}kE_{i}^{\ast}(k)+\frac{m_{s}^{2}\Phi^{2}}{g_{s}^{2}}\Bigg{(}\frac{1}{2}+\frac{\kappa_{3}}{3!}\frac{\Phi}{M}+\frac{\kappa_{4}}{4!}\frac{\Phi^{2}}{M^{2}}\Bigg{)}$
(5)
$\displaystyle+\rho_{b}W-\frac{1}{4!}\frac{\zeta_{0}W^{4}}{g_{\omega}^{2}}-\frac{1}{2}m_{\omega}^{2}\frac{W^{2}}{g_{\omega}^{2}}\Bigg{(}1+\eta_{1}\frac{\Phi}{M}+\frac{\eta_{2}}{2}\frac{\Phi^{2}}{M^{2}}\Bigg{)}$
$\displaystyle+\frac{1}{2}\rho_{3}R-\frac{1}{2}\Bigg{(}1+\frac{\eta_{\rho}\Phi}{M}\Bigg{)}\frac{m_{\rho}^{2}}{g_{\rho}^{2}}R^{2}-\Lambda_{\omega}(R^{2}\times
W^{2})$ $\displaystyle+\frac{1}{2}\frac{m_{\delta}^{2}}{g_{\delta}^{2}}D^{2},$

Figure 1: (Color online) The numerical fitting of nuclear matter E/A as a
function of baryon number density for NL3 parameter set using equation (6)
with a different number of terms. The $\delta$ represents the mean deviation
between the RMF and fitted data. Here, 5, 8, and 12 terms in the fitting stand
for the summation of $i$ from 3 to 7, 10, and 14 as given in Eq. (6). More
details can find in the text. Figure 2: (Color online) The nuclear matter E/A
as a function of baryon number density for different asymmetry
$\alpha=\frac{\rho_{n}-\rho_{p}}{\rho_{n}+\rho_{p}}$. For symmetric NM
$\alpha=0$ and $\alpha=1$ for pure neutron matter.

###  II.1 Fitting procedure

The important part of the present calculation is to convert the NM quantities
[Eq. (5)] from momentum, space to coordinate space, i.e. the reconstruction of
NM quantities at local density. The results of our calculations are shown in
Fig. 2 for NL3 and G3 parameter sets. The NL3 set gives a stiff equation of
state (EOS) as compared to the G3 force. This is because the NM
incompressibility $K$ at the saturation for NL3 is 271.76 MeV, and that of G3
is 243.96 MeV. We consider that the NM is composed of tiny spherical pieces
described by a local density function $\rho_{0}(x)=3A/4\pi x^{3}$. Using this
consideration, the fitted binding energy function (Fig. 2) of E-RMF is
embedded in the following equation,

$\displaystyle{\cal E}(x)$ $\displaystyle=$ $\displaystyle
C_{k}\rho_{0}^{2/3}(x)+\sum_{i=3}^{14}(b_{i}+a_{i}\alpha^{2})\rho_{0}^{i/3}(x),$
(6)

where $C_{k}$ is the kinetic energy coefficient given as
$C_{k}=37.53[(1+\alpha)^{5/3}+(1-\alpha)^{5/3}]$ within the Thomas-Fermi
approach. To find out the exact nature of the $E/A$ in position space, we use
a polynomial fitting which consists of several numbers of terms (Eq. (6)). The
fitted graph using the different number of terms for NL3 parameter set along
with the mean deviation coefficient $\delta$ is depicted in Fig. 1. The mean
deviation is calculated using the formula
$\delta=\Big{[}\sum_{j=1}^{N}(E/A)_{j,Fitted}-(E/A)_{j,RMF}\Big{]}/N$. Here,
$(E/A)_{j,Fitted}$ is the binding energy obtained from the polynomial fitting,
and $(E/A)_{j,RMF}$ is the binding energy per nucleon from the RMF functional
with N being the number of data points. Firstly, we fitted our data using 5
terms in the expansion of Eq. (6) (means $i$ runs from 3 to 7) and found that
the deviation of fitted data is around 14 percent from the actual value. .
Similarly, for 8, 10, and 12 terms, the deviation are 6.5%, 3.5%, and 0.2%,
respectively. We inspect that as we increase the number of terms in the
expansion of Eq. (6), the deviation of the fitting gets reduced considerably
and the fitting converges to more candidness. We also observe that around
10-12 terms are required in the expansion for an appropriate fitting of most
of the available RMF parameter sets. For the best fit, we take 12 number of
terms in the present calculation, and the coefficients are obtained from the
polynomial fitting, whose values are given in Table 1.

The NM parameters $K^{NM}$, $S^{NM}$, $L_{sym}^{NM}$ and $K_{sym}^{NM}$ are
obtained from the following standard relations Kumar  _et al._ (2018); Chen
and Piekarewicz (2014)

$\displaystyle K^{NM}$ $\displaystyle=$ $\displaystyle
9\rho^{2}\frac{\partial^{2}}{\partial\rho^{2}}\bigg{(}\frac{\cal
E}{\rho}\bigg{)}\Big{|}_{\rho=\rho_{0}},$ (7) $\displaystyle S^{NM}$
$\displaystyle=$ $\displaystyle\frac{1}{2}\frac{\partial^{2}({\cal
E}/\rho)}{\partial\alpha^{2}}\Big{|}_{\alpha=0},$ (8) $\displaystyle
L_{sym}^{NM}$ $\displaystyle=$ $\displaystyle 3\rho\frac{\partial
S(\rho)}{\partial\rho}\Big{|}_{\rho=\rho_{0}}=\frac{3P}{\rho_{0}},$ (9)
$\displaystyle K_{sym}^{NM}$ $\displaystyle=$ $\displaystyle
9\rho^{2}\frac{\partial^{2}S(\rho)}{\partial\rho^{2}}\Big{|}_{\rho=\rho_{0}}.$
(10)

which are given as follow using Eq. (6)

$\displaystyle K^{NM}$ $\displaystyle=$
$\displaystyle-150.12\,\rho_{0}^{2/3}(x)+\sum_{i=4}^{14}i\,(i-3)\,b_{i}\,\rho_{0}^{i/3}(x),$
(11) $\displaystyle S^{NM}$ $\displaystyle=$ $\displaystyle
41.7\,\rho_{0}^{2/3}(x)+\sum_{i=3}^{14}a_{i}\,\rho_{0}^{i/3}(x),$ (12)
$\displaystyle L_{sym}^{NM}$ $\displaystyle=$ $\displaystyle
83.4\,\rho_{0}^{2/3}(x)+\sum_{i=3}^{14}i\,a_{i}\,\rho_{0}^{i/3}(x),$ (13)
$\displaystyle K_{sym}^{NM}$ $\displaystyle=$
$\displaystyle-83.4\,\rho_{0}^{2/3}(x)+\sum_{i=4}^{14}i\,(i-3)\,a_{i}\,\rho_{0}^{i/3}(x),$
(14)

Table 1: The coefficients of the analytical expression for NM binding energy per particle as a function of density $\rho(x)$ and the asymmetric factor $\alpha=\frac{\rho_{n}-\rho_{p}}{\rho_{n}+\rho_{p}}$. The values are given for NL3 and G3 parameter sets. | NL3 | G3
---|---|---
b3 | -3449.92 | -490.15
b4 | 93386.65 | -465.80
b5 | -1233527.10 | 7107.17
b6 | 9041665.48 | -53960.91
b7 | -41166214.95 | 284155.27
b8 | 123164197.67 | -938303.73
b9 | -248225071.34 | 2066363.99
b10 | 338087637.81 | -3133853.65
b11 | -305682367.52 | 3246326.72
b12 | 174988863.34 | -2188895.63
b13 | -57095582.73 | 861872.30
b14 | 8030800.13 | -149719.19
a3 | -1098.99 | 391.87
a4 | 43110.63 | -5565.05
a5 | -636205.52 | 80413.45
a6 | 5283025.49 | -639847.05
a7 | -27638744.86 | 3139872.96
a8 | 96251325.61 | -10022304.90
a9 | -228719960.80 | 21277231.89
a10 | 372188746.60 | -30256280.52
a11 | -407653299.90 | 28503817.38
a12 | 286972591.80 | -17096240.48
a13 | -117150348.00 | 5921783.83
a14 | 21061682.62 | -903228.53

The densities of closed/semi-closed-shell spherical nuclei 16O, 40,48Ca, 56Ni,
90Zr, 116Sn and 208Pb are calculated using E-RMF formalism. These densities
are used as input in the CDFM (described in the following sub-section) to
calculate the weight function, which is a key quantity acting as a bridge
between NM parameters in $x-$space and finite nuclei in $r-$space (using LDA).
To match with the $r-$ and $x-$space together, we construct the total density
of the nucleus with the superposition of an infinite number of $Fluctons$,
following the approach of CDFM discussed below.

###  II.2 Coherent Density Fluctuation model

In the CDFM, we use the NM quantities $K^{NM}$, $S^{NM}$, $L_{sym}^{NM}$, and
$K_{sym}^{NM}$ from Eqs. (11–14) to extract their values for finite nuclei
Antonov  _et al._ (1980, 1982, 1985, 1994). Within CDFM, the one-body density
matrix (OBDM) $\rho$ (r, r′) of a finite nucleus is written as the coherent
superposition of OBDM $\rho_{x}$ (r, r′) for spherical pieces of NM termed as
Fluctons Bhuyan  _et al._ (2018); Gaidarov  _et al._ (2011),

$\rho_{x}({\bf r})=\rho_{0}(x)\,\Theta(x-|{\bf r}|),$ (15)

with $\rho_{o}(x)=\frac{3A}{4\pi x^{3}}$. The generator coordinate x is the
radius of a sphere consisting of Fermi gas having all the A nucleons
distributed uniformly within it. It is suitable to apply for such a system the
OBDM expressed as below Bhuyan  _et al._ (2018); Antonov  _et al._ (1982);
Gaidarov  _et al._ (2011, 2012),

$\rho({\bf r},{\bf r^{\prime}})=\int_{0}^{\infty}dx|F(x)|^{2}\rho_{x}({\bf
r},{\bf r^{\prime}}),$ (16)

where, $|F(x)|^{2}$ is the weight function (WF). The coherent superposition of
OBDM $\rho_{x}({\bf r},{\bf r^{\prime}})$ is given as:

$\displaystyle\rho_{x}({\bf r},{\bf r^{\prime}})$ $\displaystyle=$
$\displaystyle 3\rho_{0}(x)\frac{J_{1}\left(k_{f}(x)|{\bf r}-{\bf
r^{\prime}}|\right)}{\left(k_{f}(x)|{\bf r}-{\bf r^{\prime}}|\right)}$ (17)
$\displaystyle\times\Theta\left(x-\frac{|{\bf r}+{\bf
r^{\prime}}|}{2}\right),$

where J1 is the first order spherical Bessel function and $k_{f}$is the Fermi
momentum of nucleons inside the $Flucton$ having radius $x$ and
$k_{f}(x)=(3\pi^{2}/2\rho_{0}(x))^{1/3}=\gamma/x$, ($\gamma\approx
1.52A^{1/3}$). The Wigner distribution function of the OBDM of Eq. (17) is
given by,

$\displaystyle W({\bf r},{\bf k})=\int_{0}^{\infty}dx\,|F(x)|^{2}\,W_{x}({\bf
r},{\bf k}).$ (18)

Here, $W_{x}({\bf r},{\bf k})=\frac{4}{8\pi^{3}}\Theta(x-|{\bf
r}|)\Theta(k_{F}(x)-|{\bf k}|)$. The density $\rho$ (r) in terms of the WF
within the CDFM approach is:

$\displaystyle\rho(r)$ $\displaystyle=$ $\displaystyle\int d{\bf k}W({\bf
r},{\bf k})$ (19)
$\displaystyle=\int_{0}^{\infty}dx\,|F(x)|^{2}\,\frac{3A}{4\pi
x^{3}}\Theta(x-|{\bf r}|),$

which is normalized to A, i.e., $\int\rho({\bf r})d{\bf r}=A$. In the
$\delta$-function limit, the Hill-Wheeler integral equation, that is the
differential equation for the WF in the generator coordinate is obtained
Antonov  _et al._ (1980). The $|F(x)|^{2}$ for a given density $\rho$ (r) is
expressed as

$|F(x)|^{2}=-\left(\frac{1}{\rho_{0}(x)}\frac{d\rho(r)}{dr}\right)_{r=x},$
(20)

with $\int_{0}^{\infty}dx|F(x)|^{2}=1$. We refer Bhuyan  _et al._ (2018);
Antonov  _et al._ (1980, 1982); Gaidarov  _et al._ (2011, 2012) for a detailed
analytical derivation. The CDFM allows us to make a transition from the
properties of NM to those of finite nuclei. The finite nuclear
incompressibility $K^{A}$, symmetry energy $S^{A}$, neutron pressure $P^{A}$
and surface incompressibility $K_{sym}^{A}$ for a finite nucleus are
calculated by weighting the corresponding quantities for infinite NM within
the CDFM, as given below Antonov  _et al._ (1994); Gaidarov  _et al._ (2011,
2012); Fuchs  _et al._ (1995); Antonov  _et al._ (2017)

$\displaystyle K^{A}=\int_{0}^{\infty}dx\,|F(x)|^{2}\ K^{NM}(\rho(x)).$ (21)
$\displaystyle P^{A}=\int_{0}^{\infty}dx\,|F(x)|^{2}\,P^{NM}(\rho(x)),$ (22)
$\displaystyle S^{A}=\int_{0}^{\infty}dx\,|F(x)|^{2}\,S^{NM}(\rho(x)),$ (23)
$\displaystyle
L_{sym}^{A}=\int_{0}^{\infty}dx\,|F(x)|^{2}\,L_{sym}^{NM}(\rho(x)),$ (24)
$\displaystyle K_{sym}^{A}=\int_{0}^{\infty}dx\,|F(x)|^{2}\
K_{sym}^{NM}(\rho(x)),$ (25)

The $K^{A}$, $P^{A}$, $S^{A}$, $L_{sym}^{A}$, and $K_{sym}^{A}$ in Eqs.
($\ref{K0}-$25) are the surface weighted average of the corresponding NM
quantities in the LDA limit for finite nuclei.

##  III Results and discussion

In the present work, we have derived the relativistic density functional
expression from the NM calculations for different neutron-proton asymmetry
($\alpha$) using the NL3 and the recently developed G3 parameter sets. It is
quite important to note that the Br$\ddot{u}$ckner density functional
Brueckner  _et al._ (1968) used within a CDFM to calculate symmetry energy and
related observables in previous studies Antonov  _et al._ (2016); Bhuyan  _et
al._ (2018); Quddus  _et al._ (2020) is inadequate. The underlying reason is
that the NM saturation curves fitted by Br$\ddot{u}$ckner energy density
functional could not reproduce the empirical values of saturation density
($\rho\sim 0.15$ fm-3) and binding energy per nucleon (E $\sim$ 16 MeV) of NM
(see Fig.1 of Ref. Brueckner  _et al._ (1968)). In other words,
Br$\ddot{u}$ckner energy density functional could not sort out the Coster-Band
problem Coester  _et al._ (1970). This issue was rectified partially by the
inclusion of three-body force in the nucleon-nucleon potential Fujita and
Miyazawa (1957); Pieper  _et al._ (2001). To address this issue, we emphasize
to obtain relativistic energy density functional since the RMF Lagrangian with
non-linear terms mimics the three-body effect in the nuclear potential and
resolves the Coster-Band problem.

Figure 3: (color online) The densities (dotted lines) and weight functions (solid lines) of 16O and 208Pb for NL3 (upper panel) and G3 (lower panel) parameter sets. Table 2: The properties such as nuclear incompressibility $K^{A}$, symmetric energy $S^{A}$, neutron pressure $P^{A}$ = $\rho_{0}$ $L_{sym}^{A}$/3, slope $L_{sym}^{A}$ and curvature $K_{sym}^{A}$ of nuclei, with RMF density functional, for NL3 and G3 parameter sets. The dimension of all the parameters are in MeV. NL3 |  16O |  40Ca |  48Ca |  56Ni |  90Zr |  116Sn |  208Pb
---|---|---|---|---|---|---|---
$K^{A}$ | 618.05 | 584.11 | 564.86 | 627.33 | 450.87 | 476.46 | 411.03
$P^{A}$ | 9.70 | 8.30 | 7.84 | 8.68 | 6.90 | 7.24 | 6.68
$S^{A}$ | 40.30 | 39.40 | 38.83 | 41.66 | 36.97 | 38.38 | 37.28
$L_{sym}^{A}$ | 118.89 | 119.79 | 120.54 | 130.30 | 117.01 | 121.20 | 117.95
$K_{sym}^{A}$ | 32.42 | -10.92 | -3.18 | -4.37 | 33.30 | 31.34 | 47.99
G3 |  16O |  40Ca |  48Ca |  56Ni |  90Zr |  116Sn |  208Pb
$K^{A}$ | 258.87 | 262.09 | 270.26 | 279.33 | 253.87 | 249.74 | 238.85
$P^{A}$ | 3.39 | 3.07 | 3.08 | 3.12 | 2.75 | 2.68 | 2.52
$S^{A}$ | 30.12 | 30.43 | 31.15 | 31.98 | 30.88 | 30.87 | 30.53
$L_{sym}^{A}$ | 51.98 | 51.92 | 52.71 | 53.57 | 51.40 | 51.18 | 50.28
$K_{sym}^{A}$ | -103.23 | -89.43 | -89.52 | -90.00 | -91.55 | -93.07 | -96.39
Table 3: The properties such as symmetric energy $S^{A}$, neutron pressure $P^{A}$ and curvature $K_{sym}^{A}$ of nuclei for NL3 and G3 parameter sets using Br$\ddot{u}$ckner’s density functional within CDFM. The dimension of all the parameters are in MeV. NL3 |  16O |  40Ca |  48Ca |  56Ni |  90Zr |  116Sn |  208Pb
---|---|---|---|---|---|---|---
$S^{A}$ | 25.33 | 27.28 | 28.36 | 29.91 | 29.13 | 29.38 | 29.28
$P^{A}$ | -0.66 | 0.05 | -0.05 | 0.02 | 0.95 | 1.25 | 1.57
$K_{sym}^{A}$ | -281.34 | -275.28 | -296.82 | -324.29 | -267.48 | -257.88 | -236.89
G3 |  16O |  40Ca |  48Ca |  56Ni |  90Zr |  116Sn |  208Pb
$S^{A}$ | 25.39 | 27.28 | 28.42 | 29.40 | 29.05 | 29.19 | 29.12
$P^{A}$ | -0.41 | 0.31 | 0.32 | 0.53 | 1.23 | 1.48 | 1.61
$K_{sym}^{A}$ | -262.45 | -261.99 | -227.39 | -288.89 | -251.25 | -241.56 | -229.63

The NM binding energy per particle for variable $\alpha$ within the RMF models
(NL3 and G3 parameter sets) are shown in Fig. 2. These curves for different
values of $\alpha$ are fitted with a polynomial [Eq. (6)] where the first term
is the kinetic energy taken from Thomas-Fermi approximation and the second
term presents the potential energy. In the potential energy terms, the
obtained coefficients $a_{i}$ and $b_{i}$ are shown in Table 1. It is
important to note that the expression $E/A$ of NM in Eq. (5) is obtained as a
function of the density of spherical pieces of NM $\rho_{0}(x)$, called
Fluctons, as the expansion variable. The Eq. (6) is, in fact, equivalent
energy expression of Eq. (5) and therefore the subsequent expressions of
symmetry energy and other quantities (Eqs. (21–25)) obtained from Eq. (6) also
encompass the relativistic characteristics.

Figure 4: (color online) The proton (circle) and neutron (square) surface
diffusion parameter for the considered nuclei in Table 2 for NL3 (open symbol)
and G3 (solid symbol) parameter sets along with the experimental charge (open
triangle) diffuseness estimates De Vries  _et al._ (1987); Nadjakov  _et al._
(1994). See the text for more details.

Further, within CDFM, NM parameters together with weight functions will be
used to evaluate the corresponding quantities of finite closed/semi-closed-
shell 16O, 40,48Ca, 56Ni, 90Zr, 116Sn and 208Pb nuclei (see Eqs. (21–25)). The
self-consistently calculated RMF densities are the crucial ingredient to find
weight function of nuclei and both the density and weight function for 16O and
208Pb nuclei are shown in Fig. 3 as representative cases. The dashed lines
present the density of nuclei, and the solid lines depict the variation of
weight function. It is evident that although the density at the centre of 16O
and 208Pb nuclei is high, the value of weight function in that region is very
small. Moreover, the weight function exhibits a bell shape with a maxima,
which corresponds to the surface region where the density is significantly
reduced compared to the central region. Hence, in the surface region, the
weight function makes a substantial contribution to the calculation of
symmetry energy and other quantities. Due to this reason, the symmetry energy,
neutron pressure, etc. are labelled as surface properties. It is also
important to note that the total symmetry energy $S$ can be made a partition
into volume $S_{V}$ and surface $S_{S}$ components with an analogy of two
connected capacitors. A detailed discussion can be obtained in Refs. Kaur  _et
al._ (2020a); Danielewicz (2003, 2004, 2007); A. E. L.Dieperink and P. Van
Isacker (2007); Danielewicz and Lee (2009); Antonov  _et al._ (2018).

Figure 5: (color online) Left panel: The energy density ($\cal{E}$) vs density
($\rho$). Right panel: The density distribution of 16O and 208Pb for NL3 and
G3 parameter sets.

To determine the surface effect on the symmetry energy of finite nuclei, we
calculate the surface diffuseness parameters for all the nuclei by using the
nucleon density distributions. The equivalent nuclear surface diffuseness
parameter can be obtained using the relation
$a_{i}\approx-\rho_{i}/\frac{d\rho_{i}}{dr}$, where i stands for
proton-($a_{p}$), neutron-($a_{n}$) and charge-density ($a_{ch}$) density
distributions Bhuyan  _et al._ (2020) and reference therein. The results for
the surface diffuseness parameter are obtained for proton (circle) and neutron
(square) by using the respective densities from the RMF approach for NL3 (open
symbol) and G3 (solid symbol) parameter sets. The experimental diffuseness
parameter (open triangle) are estimated from the charge density distribution
De Vries  _et al._ (1987); Nadjakov  _et al._ (1994). All the results along
with the experimental estimates are given in Fig. 4. From the figure, one can
observe that there is a constant difference in the calculated diffuseness
parameter for proton and the experimental data for the charge. This is due to
the charge effect in the proton density in the experiment, whereas we use
proton number density in which the charge effect is absent. Comparing the Fig.
4 and Table 2, one can find a linear correlation among the neutron diffuseness
parameter and symmetry energy. In other words, the trend appearing for the
symmetry energy exactly reflects in the neutron surface diffuseness parameter.
In contrast, the diffuseness parameter for proton follows the linear fall,
which is consistent with the experimental data. Hence, the similar trends in
the symmetry energy and surface diffuseness parameter for neutron show the
implication of surface effects in terms of weight function for the study of NM
quantities for finite nuclei.

The surface properties such as symmetry energy, neutron pressure, slope and
curvature of symmetry energy, incompressibility of NM Eqs. (11–14) are
computed using relativistic density functional Eqs. (6). These are further
folded with weight function to calculate the corresponding quantities of
closed/semi-closed-shell even-even nuclei, which are shown in Table 2. From
the table, one can notice a wide range of finite nuclear incompressibility
$K^{A}$ = 618.05 to 411.03 MeV and $K^{A}$ = 258.87 to 238.85 MeV for 16O -
208Pb with NL3 and G3 parameter set, respectively. To understand this nature,
we have given the density distribution of nucleons for both NL3 and G3 in Fig.
5 for 16O and 208Pb as two representative cases. The energy density functional
($\cal{E}$) for both NL3 and G3 sets are also displayed in the left panel of
the same figure. It is very much clear that the densities are significantly
different for both the forces, mostly in the central region. In the case of
G3, a comparatively flat distribution of nucleons appears as compared to NL3,
predicting a small variation in the weight function for both the sets.
However, if one analyzes the variation in the energy density $\cal{E}$ (or
E/A) with density $\rho$ for both NL3 and G3 sets, a huge difference between
these two sets are visible. For the NL3 case, the variation in energy density
$\cal{E}$ with $\rho$ is substantially larger than G3. As a result, we get a
small $K^{A}$ for G3 than NL3, although the NM incompressibility at the
saturation is almost comparable for the two force parameters. From Eq. (7), we
know that the incompressibility is the second derivative of EOS (E/A) for the
density. Therefore, the incompressibility predicted by NL3 is more than the G3
for all nuclei. The variation in other quantities can also be explained taking
into the behavior of E/A for various force parameters. From these results, we
cannot conclude about the mass dependence of the finite nuclear
incompressibility, for example, $K^{A}=$ 618, 627, and 411 MeV for 16O, 56Ni
and 208Pb, respectively for NL3 case. Similar uncertainty in $K^{A}$ for
different nuclei is also seen in the G3 case. However, for $P^{A}$, $S^{A}$
and $L_{sym}^{A}$, the variation is in a narrow range from 16O to 208Pb, i.e.,
minimum $P^{A}$ is 2.52 MeV for 208Pb and maximum is 3.39 MeV for 16O with G3.
Furthermore, $K_{sym}^{A}$ varies a lot depending on both the force parameter
as well as the mass of the nucleus. The very different values of $K^{A}$ and
$K_{sym}^{A}$ for light and heavy nuclei indicate the structural dependence of
the finite nuclei since the density distribution varies from the nucleus to
nucleus.

For the sake of comparison, the values of symmetry energy and other quantities
of nuclei obtained using Br$\ddot{u}$ckner density functional within CDFM are
shown in Table 3 for NL3 and G3 densities. The values of symmetry energy
$S^{A}$ linearly increases with the mass number and the lies in the range
25-29 MeV. Similar conclusions can be drawn for $P^{A}$ and $K_{sym}^{A}$ for
these considered nuclei. We get $P^{A}$ is negative for 16O for both parameter
sets, which can be correlated with the charge effect of proton for N = Z
nucleus. We also get negative $P^{A}$ for 48Ca for the NL3 case, which may be
connected with the reduced charge radius of 48Ca in comparison to 40Ca, mainly
an unsolved issue in the theoretical models. In comparison, the results using
relativistic density functional along with mean-field densities for both the
parameter sets give positive $P^{A}$ value for all the nuclei (see Table 2).
We get a large difference and/or a wide range of values for symmetry energy
curvature $K_{sym}^{A}$ in Table 3. This difference is due to the small
difference in the resultant energy density for relativistic and
Br$\ddot{u}$kner energy density functional (see Table 2).

##  IV Summary and Conclusion

In brief, we have fitted the NM saturation curves for different values of
asymmetry parameter, employing the E-RMF density functional with two well-
known G3 and NL3 parameter sets, in the wake of the solution of the Coster-
Band problem by relativistic approach. The newly fitted expression of E/A is
used to find NM parameters such as incompressibility, neutron pressure,
symmetry energy. These NM parameters, together with E-RMF densities, are used
in the calculations of surface properties of magic nuclei ranging from light
to heavy mass region, within the coherent density fluctuation model. The
values of symmetry energy, neutron pressure, L-coefficient show small
variation for different mass nuclei. On the other hand, the incompressibility
and surface incompressibility show large variation while moving from light 16O
to heavy 208Pb nuclei for NL3 compared to the G3 case. This is due to the
variation of E/A with density or EOS is stiff for NL3 than in the G3 case and
as a result, the incompressibility (second derivative of E/A w.r.t. density)
predicted by NL3 is more than G3 for all nuclei. Besides, the correlation
between the diffuseness parameter and symmetry energy is discussed in the
context of surface effects. The results using relativistic energy density
functions are also compared with ones obtained using Br$\ddot{u}$kner energy
density functional.

In other words, we discuss here that the direct use of Br$\ddot{u}$ckner
energy functional to evaluate the effective nuclear surface properties is not
adequate in the context of the Coster-Band issue. The present work provides
the way to estimate the properties of $\beta$-stable and $\beta$-unstable
nuclei more accurately by employing relativistic density functional within
CDFM. These quantities are quite significant to probe the structural aspects
of finite nuclei and can also facilitate to constrain the EOS of NM. Further,
we also compare the RMF results with Br$\ddot{u}$ckner density functional,
which reflects the model dependant on the results. A more systematic study
using the present functional for the nuclei throughout the nuclear chart in
this direction is highly welcomed.

Acknowledgement: One of the authors (MB) is thankful to the Institute of
Physics, Bhubaneswar for the local hospitality during the scientific visit. MB
acknowledges the support from FAPESP Project No. 2017/05660-0, FOSTECT Project
No. FOSTECT.2019B.04, and the CNPq - Brasil.

## References

  * Lattimer and Prakash (2000) J. M. Lattimer and M. Prakash, Physics Reports 333-334, 121 (2000).
  * Lattimer and Prakash (2001) J. M. Lattimer and M. Prakash, The Astrophysical Journal 550, 426 (2001).
  * Steiner  _et al._ (2005) A. Steiner, M. Prakash, J. Lattimer, and P. Ellis, Physics Reports 411, 325 (2005).
  * Li  _et al._ (2008) B.-A. Li, L.-W. Chen, and C. M. Ko, Physics Reports 464, 113 (2008).
  * Lynch  _et al._ (2009) W. Lynch, M. Tsang, Y. Zhang, _et al._ , Progress in Particle and Nuclear Physics 62, 427 (2009).
  * Tsang  _et al._ (2012) M. B. Tsang, J. R. Stone, F. Camera, _et al._ , Phys. Rev. C 86, 015803 (2012).
  * Horowitz  _et al._ (2014) C. J. Horowitz, E. F. Brown, Y. Kim, _et al._ , Journal of Physics G: Nuclear and Particle Physics 41, 093001 (2014).
  * Antonov  _et al._ (2016) A. N. Antonov, M. K. Gaidarov, P. Sarriguren, and E. Moya de Guerra, Phys. Rev. C 94, 014319 (2016).
  * Bhuyan  _et al._ (2018) M. Bhuyan, B. V. Carlson, S. K. Patra, and S.-G. Zhou, Phys. Rev. C 97, 024322 (2018).
  * Kaur  _et al._ (2020a) M. Kaur, A. Quddus, A. Kumar, M. Bhuyan, and S. Patra, Nucl. Phys. A 1000, 121871 (2020a).
  * Kaur  _et al._ (2020b) M. Kaur, A. Quddus, A. Kumar, M. Bhuyan, and S. K. Patra, Journal of Physics G: Nuclear and Particle Physics 47, 105102 (2020b).
  * Naik  _et al._ (2019) K. C. Naik, M. Kaur, A. Kumar, and S. K. Patra, Int. J. of Mod. Phys. E 28, 1950100 (2019).
  * Quddus  _et al._ (2020) A. Quddus, M. Bhuyan, and S. K. Patra, Journal of Physics G: Nuclear and Particle Physics 47, 045105 (2020).
  * Brueckner  _et al._ (1968) K. A. Brueckner, J. R. Buchler, S. Jorna, and R. J. Lombard, Phys. Rev. 171, 1188 (1968).
  * Brueckner  _et al._ (1969) K. A. Brueckner, J. R. Buchler, R. C. Clark, and R. J. Lombard, Phys. Rev. 181, 1543 (1969).
  * Coester  _et al._ (1970) F. Coester, S. Cohen, B. Day, and C. M. Vincent, Phys. Rev. C 1, 769 (1970).
  * Brockmann and Machleidt (1990) R. Brockmann and R. Machleidt, Phys. Rev. C 42, 1965 (1990).
  * Kumar  _et al._ (2018) B. Kumar, S. K. Patra, and B. K. Agrawal, Phys. Rev. C 97, 045806 (2018).
  * Kumar  _et al._ (2017) B. Kumar, S. Singh, B. Agrawal, and S. Patra, Nuclear Physics A 966, 197 (2017).
  * Furnstahl  _et al._ (1996) R. Furnstahl, B. D. Serot, and H.-B. Tang, Nuclear Physics A 598, 539 (1996).
  * Furnstahl  _et al._ (1997) R. Furnstahl, B. D. Serot, and H.-B. Tang, Nuclear Physics A 615, 441 (1997).
  * Malik  _et al._ (2018) T. Malik, N. Alam, M. Fortin, _et al._ , Phys. Rev. C 98, 035804 (2018).
  * Chen and Piekarewicz (2014) W.-C. Chen and J. Piekarewicz, Phys. Rev. C 90, 044305 (2014).
  * Antonov  _et al._ (1980) A. N. Antonov, V. A. Nikolaev, and I. Z. Petkov, Z. Phys. A Atoms and Nuclei 297, 257 (1980).
  * Antonov  _et al._ (1982) A. N. Antonov, V. A. Nikolaev, and I. Z. Petkov, Z. Phys. A Atoms and Nuclei 304, 239 (1982).
  * Antonov  _et al._ (1985) A. N. Antonov, V. A. Nikolaev, and I. Z. Petkov, Nuovo Cimento A 86, 23 (1985).
  * Antonov  _et al._ (1994) A. N. Antonov, D. N. Kadrev, and P. E. Hodgson, Phys. Rev. C 50, 164 (1994).
  * Gaidarov  _et al._ (2011) M. K. Gaidarov, A. N. Antonov, P. Sarriguren, and E. Moya de Guerra, Phys. Rev. C 84, 034316 (2011).
  * Gaidarov  _et al._ (2012) M. K. Gaidarov, A. N. Antonov, P. Sarriguren, and E. M. de Guerra, Phys. Rev. C 85, 064319 (2012).
  * Fuchs  _et al._ (1995) C. Fuchs, H. Lenske, and H. H. Wolter, Phys. Rev. C 52, 3043 (1995).
  * Antonov  _et al._ (2017) A. N. Antonov, D. N. Kadrev, M. K. Gaidarov, P. Sarriguren, and E. M. de Guerra, Phys. Rev. C 95, 024314 (2017).
  * Fujita and Miyazawa (1957) J.-i. Fujita and H. Miyazawa, Progress of Theoretical Physics 17, 360 (1957).
  * Pieper  _et al._ (2001) S. C. Pieper, V. R. Pandharipande, R. B. Wiringa, and J. Carlson, Phys. Rev. C 64, 014001 (2001).
  * De Vries  _et al._ (1987) H. De Vries, C. W. De Jager, and C. De Vries, Atomic Data and Nucl. Data Tables 36, 495 (1987).
  * Nadjakov  _et al._ (1994) E. Nadjakov, K. Marinova, and Y. Gangrsky, Atomic Data and Nucl. Data Tables 56, 133 (1994).
  * Danielewicz (2003) P. Danielewicz, Nucl. Phys. A 727, 233 (2003).
  * Danielewicz (2004) P. Danielewicz, (2004), arXiv:0411115 .
  * Danielewicz (2007) P. Danielewicz, “Symmetry energy,” in  _Opportunities with Exotic Beams_ (2007) pp. 142–151.
  * A. E. L.Dieperink and P. Van Isacker (2007) A. E. L.Dieperink and P. Van Isacker, Eur. Phys. J. A 32, 11 (2007).
  * Danielewicz and Lee (2009) P. Danielewicz and J. Lee, Int. J. of Mod. Phys. E 18, 892 (2009).
  * Antonov  _et al._ (2018) A. N. Antonov, D. N. Kadrev, M. K. Gaidarov, P. Sarriguren, and E. Moya de Guerra, Phys. Rev. C 98, 054315 (2018).
  * Bhuyan  _et al._ (2020) M. Bhuyan, R. Kumar, S. Rana, D. Jain, S. K. Patra, and B. V. Carlson, Phys. Rev. C 101, 044603 (2020).


# Searching optimal conditions for quantum gates application with the new
3-body Förster resonances in Rb and Cs Rydberg atoms

K.-L. Pham, S. Lepoutre, P. Pillet and P. Cheinet<EMAIL_ADDRESS>Laboratoire Aimé Cotton, CNRS, Univ. Paris-Sud, Université Paris-Saclay, Bât.
505, 91405 Orsay, France  I. N. Ashkarin, I. I. Beterov, D. B. Tretyakov, V.
M. Entin, E. A. Yakshina and I. I. Ryabtsev<EMAIL_ADDRESS>Rzhanov
institute of semiconductor physics, Pr. Lavrentyeva 13, 630090 Novosibirsk,
Russia

###### Abstract

Abstract Three body resonant interactions between Rydberg atoms are considered
in order to perform few-body quantum gates. So far, the resonances found in
cesium or rubidium atoms relied on an adjacent two-body resonance which ceases
to exist for principal quantum numbers above $n\simeq 40$. We have proposed
recently a new class of 3-body interaction resonances in alkali-metal Rydberg
atoms [P. Cheinet et al., Quant. Elect. 50, 213 (2020)], which circumvienes
this limit. We investigate here the relative strength between this new class
of 3-body interaction resonance and quasi-forbidden 2-body interaction
resonances in rubidium and cesium Rydberg atoms. We then identify the best
case scenario for detecting and using this 3-body interaction.

##  I Introduction

Rydberg atoms Gallagher (1994) offer an ideal platform for quantum simulation
experiments Weimer  _et al._ (2010) or quantum computing Saffman  _et al._
(2010). The high control gained on ultra-cold atoms trapping Barredo  _et al._
(2016); Endres  _et al._ (2016) enabled new physics studies with Rydberg atoms
de Léséleuc  _et al._ (2019); Browaeys and Lahaye (2020); Madjarov  _et al._
(2020) involving the particularly large interactions displayed.

In this prospect, interaction resonances tuned with electric field Safinya
_et al._ (1981) have been extensively studied in alkali-metals Gallagher  _et
al._ (1982); Stoneman  _et al._ (1987); Anderson  _et al._ (1998); Mourachko
_et al._ (1998); Vogt  _et al._ (2006); Reinhard  _et al._ (2008); Ryabtsev
_et al._ (2010); Altiere  _et al._ (2011); Nipper  _et al._ (2012); Tretyakov
_et al._ (2014); Ravets  _et al._ (2015) with a recent attention to their
coherence properties Ravets  _et al._ (2014); Yakshina  _et al._ (2016).

Stark tuned Förtser resonances have recently been revisited towards few-body
resonances Gurian  _et al._ (2012); Faoro  _et al._ (2015); Tretyakov  _et
al._ (2017); Liu  _et al._ (2020) and their use for quantum computation is now
under investigation with theoretical studies on their coherence Ryabtsev  _et
al._ (2018) and on specific Toffoli gate protocol Beterov  _et al._ (2018).

In these investigations, the 3-body resonance presented below in equation (1)
benefits from the nearby 2-body Förster resonance presented below in equation
(2) to lead to a strong interaction. Considering rubidium or cesium atoms in a
given principal quantum number $n$, an orbital quantum number $\ell=1$ often
referred to as $p$ and a total momentum quantum number $j=3/2$, these
resonances can be written:

$\left|3\times
np_{3/2}\right\rangle\Leftrightarrow\left|ns_{1/2}+(n+1)s_{1/2}+np^{\prime}_{3/2}\right\rangle$
(1) $\left|2\times
np_{3/2}\right\rangle\Leftrightarrow\left|ns_{1/2}+(n+1)s_{1/2}\right\rangle$
(2)

with the $\left|np^{\prime}_{3/2}\right\rangle$ state corresponding to a
Rydberg atom in the same level but a different Zeeman sub-level. The 2-body
resonances of type (2) exist in rubidium and cesium atoms but are limited to
principal quantum numbers of around $n\simeq 40$ at most and, as a
consequence, the related 3-body resonances stop as well. But such interaction
resonances gain in strength with increasing $n$ and obtaining a 3-body
resonance for unlimited $n$ is thus of high interest.

One possibility to overcome the limit is to use atoms in different initial
states as proposed in Beterov  _et al._ (2018). But it then requires the use
of several phase-locked lasers at different frequencies which represents an
additional technical difficulty that may prevent a successfull implementation.
More recently, a new class of 3-body resonance was proposed Cheinet  _et al._
(2020) to overcome the $n$ limit while still using all atoms in a single
initial state:

$\left|3\times
np_{3/2}\right\rangle\Leftrightarrow\left|ns_{1/2}+(n+1)s_{1/2}+np_{1/2}\right\rangle.$
(3)

In this new scenario, the use of $\left|np_{1/2}\right\rangle$ state allows to
accomodate for larger intermediate state detunings in the 3-body resonance
such that it survives the disappearence of the 2-body resonance. But it also
implies a weaker 3-body interaction by typically an order of magnitude. In a
previous study Pelle  _et al._ (2016), we have demonstrated that quasi-
forbidden resonances can have non negligible strength for principal quantum
numbers of around $n\simeq 30$ compared to the 2-body interaction.

In this article, we thus study the relative strength between the new class of
3-body resonances and quasi-forbidden resonances nearby in rubidium or cesium
Rydberg atoms. Strong quasi-forbidden resonances will indeed perturb gate
protocols and may reduce significantly the resulting gate fidelity. This study
allows us to conclude on optimum scenario to perform a quantum gate from this
new resonant interaction scheme.

##  II  Theoretical methods

For this study, we follow the methodology already presented in Pelle  _et al._
(2016) that we remind briefly. We start by computing the Stark energies of all
states around the considered $\left|np_{3/2}\right\rangle$ starting state of a
given 3-body resonance of type (3). In this study, we restricted our
calculations to states $|m|=1/2$ to avoid searching through too numerous
resonances. It was demonstrated in Pelle  _et al._ (2016) that resonances
involving higher $|m|$ display similar strengths. We then obtain a resonance
map, as displayed from $\left|50p_{3/2}\right\rangle$ of cesium in figure 1:
we plot in black the Stark energies $E_{nlj}$ of all states in a given energy
interval above the starting state and in red we plot the symmetric energy
$2E_{50p_{3/2}}-E_{nlj}$. Any crossing between a black line and a red line
then corresponds to a possible interaction resonance. The strongest are the
ones for which at least one state follows the zero field selection rule
$\delta_{\ell}=\pm 1$ Pelle  _et al._ (2016). We have marqued in figure 1 the
first resonances implying states $\left|50s_{1/2}\right\rangle$ or
$\left|51s_{1/2}\right\rangle$. For the electric field range, we arbitrarily
choose to search up to around twice the 3-body resonance field.

Figure 1:  Stark Förster Resonance diagram from $\left|50p_{3/2}\right\rangle$
starting state of cesium at an energy of around -50.9 cm-1, limited to
$|m|=1/2$ states. (black) Stark energies $E_{nlj}$. (red) Symmetric energy
$2E_{50p3/2}-E_{nlj}$. The crossings between these curves are possible Förster
resonances. The blue circles show the resonances which follow the condition
$\delta_{\ell}=\pm 1$ for at least one state at zero field and for which we
calculate interaction strength.

In our search through different energy ranges in cesium, we find the following
resonances:

$\displaystyle\left|2\times np_{3/2}\right\rangle\Leftrightarrow$
$\displaystyle\left|ns_{1/2}+(n-3)l_{j}\right\rangle$ (4)
$\displaystyle\left|2\times np_{3/2}\right\rangle\Leftrightarrow$
$\displaystyle\left|(n+1)s_{1/2}+(n-4)l_{j}\right\rangle$ (5)
$\displaystyle\left|2\times np_{3/2}\right\rangle\Leftrightarrow$
$\displaystyle\left|(n+1)p_{1/2}+(n-2)d_{j}\right\rangle$ (6)
$\displaystyle\left|2\times np_{3/2}\right\rangle\Leftrightarrow$
$\displaystyle\left|(n-1)p_{1/2}+nd_{j}\right\rangle$ (7)
$\displaystyle\left|2\times np_{3/2}\right\rangle\Leftrightarrow$
$\displaystyle\left|(n+1)d_{j}+(n-3)d_{j}\right\rangle.$ (8)

When the resonances have been identified, we compute the corresponding dipole-
dipole interaction energy using the transition dipole moments between the two
starting $\left|np_{3/2}\right\rangle$ states and the two ending Stark states
Pelle  _et al._ (2016). They will depend, among other things, on the overlap
integral between the radial wavefunctions which is large for states close in
energy to the starting state.

The closest states leading to resonances and fulfilling the $\delta_{\ell}=\pm
1$ selection rule are the two states $\left|ns_{1/2}\right\rangle$ and
$\left|(n+1)s_{1/2}\right\rangle$. The first two resonances (4) and (5) are
thus the strongest as found previously Pelle  _et al._ (2016). We also note
that the resonance (8) is fully allowed regarding the selection rules but the
involved states are separated from the initial state by about 2 multiplicities
leading to small overlap integrals.

Due to significantly different quantum defects in rubidium Rydberg atoms, the
resonances found for this atom are of different type:

$\displaystyle\left|2\times np_{3/2}\right\rangle\Leftrightarrow$
$\displaystyle\left|(n-1)s_{1/2}+(n-1)l_{j}\right\rangle$ (9)
$\displaystyle\left|2\times np_{3/2}\right\rangle\Leftrightarrow$
$\displaystyle\left|(n+2)s_{1/2}+(n-4)l_{j}\right\rangle$ (10)
$\displaystyle\left|2\times np_{3/2}\right\rangle\Leftrightarrow$
$\displaystyle\left|(n-2)d_{5/2}+(n-2)l_{j}\right\rangle$ (11)
$\displaystyle\left|2\times np_{3/2}\right\rangle\Leftrightarrow$
$\displaystyle\left|(n-1)d_{5/2}+(n-3)l_{j}\right\rangle.$ (12)

We note that in rubidium, the strong quasi-forbidden resonances of type (4)
and (5) do not appear any more. They occur for fields of typically more than
twice the 3-body resonance field. The resonances involved are all of quasi-
forbidden type with large energy difference with the initial state, implying
small overlap integrals in the dipole-dipole interaction.

In the resonances (4-12) listed above, the generic label $l$ indicates that
the resonance occurs with states of the Stark multiplicity leading to a large
set of consecutive resonances. We do not compute the whole set but the first
ones in order to grasp their general behavior. The generic label $j$ is used
when the two possible values $j=l-1/2$ and $j=l+1/2$ lead to resonances. For
the sets of resonances with a multiplicity, we compute one out of the two as
they were found of similar strength Pelle  _et al._ (2016).

##  III  Results

Figure 2: Interaction strength at 10 $\mu m$ interatomic distance versus
Electric Field for various resonances in Cesium Rydberg atoms from initial
state $\left|np_{3/2}\right\rangle$. (a) Resonances at $n=50$. The 3-body
resonance is close to much stronger 2-body resonances. (b) resonances at
$n=70$. The 3-body resonance is still weaker than many close 2-body
resonances.

Förster resonances have already been studied on cesium Mourachko  _et al._
(1998); Vogt  _et al._ (2006); Nipper  _et al._ (2012) and rubidium Anderson
_et al._ (1998); Westermann  _et al._ (2006); Ryabtsev  _et al._ (2010) atoms
below $n\simeq$ 40 where the 2-body resonance (2) exists, together with the
3-body resonance (1). Therefore we concentrate on principal quantum numbers
above the limit and choose $n=50$ and $n=70$ to test the evolution with
increasing $n$. We present the results of the different resonances in cesium
and rubidium atoms in figures 2 and 3 which display their interaction strength
and compare to the corresponding 3-body resonance of type (3). As we consider
a potential application to quantum computation and the interaction strength is
highly dependent on distance, we compute the interaction strength for a
relatively large interatomic distance of 10 $\mu m$ which is likely to be well
controlled experimentally. All interaction strengths displayed in figures (2)
and (3) are also gathered in tables (1) and (2).

In cesium atoms, in figure (2), we can see that at $n=50$ the 3-body resonance
is embedded within numerous quasi-forbidden resonances which are up to a
hundred times stronger. It is therefore unlikely to achieve any reliable
quantum computation with this interaction scheme at this principal quantum
number $n$. Increasing $n$ will increase all interaction strengths but the
3-body one, which scales as $n^{11}$ will increase faster than the 2-body ones
which scale as $n^{4}$. We thus check at $n=70$ and see that the 3-body
interaction strength has indeed significantly increased comparatively to the
2-body ones. It also stands almost alone at its resonance field of 0.19 V/cm.
Nevertheless, it is still almost an order of magnitude weaker than the next
quasi-forbidden resonances of types (4) and (5). Increasing $n$ further will
keep strengthening the 3-body interaction compared to the quasi-forbidden
resonances. However the sensitivity to electric field increases as well, as
$n^{7}$, and will eventually become an issue. Another option is to use a
smaller interatomic distance. Indeed the 3-body interaction scales as
$1/R^{6}$ while the 2-body interactions scale as $1/R^{3}$ but the sensitivity
to the exact distance will then increase and might as well become an issue in
a gate protocol fidelity. It is thus interesting to check wether rubidium
atoms present better caracteristics.

Figure 3: Interaction strength at 10 $\mu m$ interatomic distance versus
Electric Field for various resonances in Rubidium Rydberg atoms from initial
state $\left|np_{3/2}\right\rangle$. (a) Resonances at $n=50$. The 3-body
resonance is embedded with numerous 2-body resonances of similar strengths.
(b) resonances at $n=70$. The 3-body resonance is now stronger and well
seperated from all surrounding 2-body resonances.

In rubidium Rydberg atoms, the energy difference between $|\,np_{3/2}\rangle$
and $|\,np_{1/2}\rangle$ is smaller than in cesium. This implies that the
consecutive 3-body resonance is stronger. We also explained that the
resonances (4) and (5) occur for significantly higher fields than the 3-body
resonance. Thus they should not impair a gate protocol. Instead, we found
resonances (9-12) surrounding the 3-body resonance with significant strength.
In figure 3, we present the strengths of these resonances together with the
3-body resonance, at $n=50$ and $n=70$. We can see that resonances of type
(12) are the strongest and are, at $n=50$, of similar strength and resonant
field as the 3-body resonance. Therefore it seems again unlikely to find
suitable experimental parameters at this principal quantum number to perform a
high fidelity 3-body quantum gate. Instead, at $n=70$, the 3-body resonance
has become the strongest interaction by almost an order of magnitude while,
like in cesium atoms, its resonance field becomes seperate from the resonance
field of the quasi-forbidden processes.

For increasing $n$, the relative strength of the 3-body resonance will keep
increasing compared to the spurious quasi-forbidden resonances and will
completely dominate them even at interatomic distances of up to 10 microns. We
thus conclude that rubidium atoms are indeed suited to implement 3-body
quantum gates based on this interaction resonance.

| n=50 | n=70
---|---|---
| F(V/cm) | I.S.(kHz) | F(V/cm) | I.S.(kHz)
ns+(n-3)l | 1.25 | 113.568 | 0.25 | 316.617
| 1.45 | 91.953 | 0.28 | 311.760
|  |  | 0.32 | 308.427
|  |  | 0.39 | 199.902
(n+1)s+(n-4)l | 1.62 | 79.728 | 0.30 | 234.933
| 2.05 | 50.453 | 0.36 | 239.299
|  |  | 0.43 | 213.391
|  |  |  |
(n+1)d5/2+(n-3)d5/2 | 1.48 | 8.250 |  |
(n+1)d3/2+(n-3)d5/2 | 1.65 | 3.970 | 0.05 | 20.874
(n+1)d5/2+(n-3)d3/2 | 1.66 | 2.890 | 0.11 | 19.118
(n+1)d3/2+(n-3)d3/2 | 1.84 | 1.580 | 0.10 | 6.990
(n-1)p1/2+nd5/2 | 0.73 | 2.380 |  |
(n-1)p1/2+nd3/2 | 1.04 | 1.080 |  |
(n-2)d3/2+(n+1)p1/2 | 0.60 | 2.860 |  |
ns+(n+1)s+np1/2 | 1.23 | 1.100 | 0.19 | 52.0
Table 1: Cesium computed resonance strengths. For each quasi-forbidden resonance of types (4 \- 8), the resonance electric field F and expected interaction strength (I.S.) is reported. It is compared at the end of the table with the 3-body resonance. | n=50 | n=70
---|---|---
| F(V/cm) | I.S.(kHz) | F(V/cm) | I.S.(kHz)
(n-1)s+(n-1)l | 0.90 | 0.130 | 0.109 | 0.106
| 1.30 | 0.120 | 0.113 | 0.088
|  |  | 0.117 | 0.081
|  |  | 0.122 | 0.077
(n+2)s+(n-4)l | 1.02 | 0.165 | 0.17 | 0.413
| 1.28 | 0.173 | 0.18 | 0.406
|  |  | 0.19 | 0.403
|  |  | 0.20 | 0.404
(n-2)d5/2+(n-2)l | 1.24 | 0.650 | 0.190 | 4.794
| 1.29 | 0.810 | 0.195 | 5.856
| 1.34 | 0.960 | 0.201 | 6.721
| 1.4 | 1.100 | 0.207 | 7.513
(n-1)d5/2+(n-3)l | 1.00 | 1.510 | 0.161 | 11.766
| 1.04 | 1.920 | 0.166 | 14.459
| 1.09 | 2.300 | 0.171 | 16.710
| 1.14 | 2.690 | 0.176 | 18.821
| 1.20 | 3.100 |  |
| 1.25 | 3.500 |  |
ns+(n+1)s+np1/2 | 1.22 | 3.493 | 0.14 | 172
Table 2: Rubidium computed resonance strengths. For each quasi-forbidden
resonance of types (9 \- 12), the resonance electric field F and expected
interaction strength (I.S.) is reported. It is compared at the end of the
table with the 3-body resonance.

##  IV Conclusion

In summary, we have investigated the quasi-forbidden resonances surrounding
the 3-body interaction resonance scheme (3) presented in Cheinet  _et al._
(2020) and evaluated their relative strengths both in cesium and rubidium
Rydberg atoms. This study points out that rubidium is better suited for
implementing quantum calculations with 3-body interactions than cesium,
provided the use of a high enough principal quantum number. Future studies
should concentrate on $n=70$ or higher in rubidium atoms. The next step
towards quantum calculations is to evaluate specific gate protocols
functionning with this interaction scheme. Expected key parameters to obtain
good fidelities are a high degree of control of atomic positions and of the
electric field to achieve high coherence times.

##  V Acknowledgements

The french team was supported by the EU H2020 FET Proactive project RySQ
(grant N. 640378). The russian team was supported by Novosibirsk State
University. This work was also supported by cooperation grant ECOMBI (CNRS
grant No PRC2312 and RFBR grant No 19-52-15010).

## References

  * Gallagher (1994) T. Gallagher,  _Rydberg atoms_ (Cambridge University Press, 1994).
  * Weimer  _et al._ (2010) H. Weimer, M. Müller, I. Lesanovsky, P. Zoller, and H. P. Büchler, Nature Physics 6, 382 (2010).
  * Saffman  _et al._ (2010) M. Saffman, T. G. Walker, and K. Mølmer, Rev. Mod. Phys. 82, 2313 (2010).
  * Barredo  _et al._ (2016) D. Barredo, S. De Léséleuc, V. Lienhard, T. Lahaye, and A. Browaeys, Science 354, 1021 (2016).
  * Endres  _et al._ (2016) M. Endres, H. Bernien, A. Keesling, H. Levine, E. R. Anschuetz, A. Krajenbrink, C. Senko, V. Vuletic, M. Greiner, and M. D. Lukin, Science 354, 1024 (2016).
  * de Léséleuc  _et al._ (2019) S. de Léséleuc, V. Lienhard, P. Scholl, D. Barredo, S. Weber, N. Lang, H. P. Büchler, T. Lahaye, and A. Browaeys,  365, 775 (2019).
  * Browaeys and Lahaye (2020) A. Browaeys and T. Lahaye, Nat. Phys. 16, 132 (2020).
  * Madjarov  _et al._ (2020) I. S. Madjarov, J. P. Covey, A. L. Shaw, J. Choi, A. Kale, A. Cooper, H. Pichler, V. Schkolnik, J. R. Williams, and M. Endres, Nature Physics 16, 857 (2020).
  * Safinya  _et al._ (1981) K. A. Safinya, J. F. Delpech, F. Gounand, W. Sandner, and T. F. Gallagher, Phys. Rev. Lett. 47, 405 (1981).
  * Gallagher  _et al._ (1982) T. F. Gallagher, K. A. Safinya, F. Gounand, J. F. Delpech, W. Sandner, and R. Kachru, Phys. Rev. A 25, 1905 (1982).
  * Stoneman  _et al._ (1987) R. C. Stoneman, M. D. Adams, and T. F. Gallagher, Phys. Rev. Lett. 58, 1324 (1987).
  * Anderson  _et al._ (1998) W. R. Anderson, J. R. Veale, and T. F. Gallagher, Phys. Rev. Lett. 80, 249 (1998).
  * Mourachko  _et al._ (1998) I. Mourachko, D. Comparat, F. de Tomasi, A. Fioretti, P. Nosbaum, V. M. Akulin, and P. Pillet, Phys. Rev. Lett. 80, 253 (1998).
  * Vogt  _et al._ (2006) T. Vogt, M. Viteau, J. Zhao, A. Chotia, D. Comparat, and P. Pillet, Phys. Rev. Lett. 97, 083003 (2006).
  * Reinhard  _et al._ (2008) A. Reinhard, T. Cubel Liebisch, K. C. Younge, P. R. Berman, and G. Raithel, Phys. Rev. Lett. 100, 123007 (2008).
  * Ryabtsev  _et al._ (2010) I. I. Ryabtsev, D. B. Tretyakov, I. I. Beterov, and V. M. Entin, Phys. Rev. Lett. 104, 073003 (2010).
  * Altiere  _et al._ (2011) E. Altiere, D. P. Fahey, M. W. Noel, R. J. Smith, and T. J. Carroll, Phys. Rev. A 84, 053431 (2011).
  * Nipper  _et al._ (2012) J. Nipper, J. B. Balewski, A. T. Krupp, B. Butscher, R. Löw, and T. Pfau, Phys. Rev. Lett. 108, 113001 (2012).
  * Tretyakov  _et al._ (2014) D. B. Tretyakov, V. M. Entin, E. A. Yakshina, I. I. Beterov, C. Andreeva, and I. I. Ryabtsev, Phys. Rev. A 90, 041403 (2014).
  * Ravets  _et al._ (2015) S. Ravets, H. Labuhn, D. Barredo, T. Lahaye, and A. Browaeys, Phys. Rev. A 92, 020701 (2015).
  * Ravets  _et al._ (2014) S. Ravets, H. Labuhn, D. Barredo, L. Béguin, T. Lahaye, and A. Browaeys, Nature Physics 10, 914 (2014).
  * Yakshina  _et al._ (2016) E. A. Yakshina, D. B. Tretyakov, I. I. Beterov, V. M. Entin, C. Andreeva, A. Cinins, A. Markovski, Z. Iftikhar, A. Ekers, and I. I. Ryabtsev, Phys. Rev. A 94, 043417 (2016).
  * Gurian  _et al._ (2012) J. H. Gurian, P. Cheinet, P. Huillery, A. Fioretti, J. Zhao, P. L. Gould, D. Comparat, and P. Pillet, Phys. Rev. Lett. 108, 023005 (2012).
  * Faoro  _et al._ (2015) R. Faoro, B. Pelle, A. Zuliani, P. Cheinet, E. Arimondo, and P. Pillet, Nature Communications 6, 8173 (2015).
  * Tretyakov  _et al._ (2017) D. B. Tretyakov, I. I. Beterov, E. A. Yakshina, V. M. Entin, I. I. Ryabtsev, P. Cheinet, and P. Pillet, Phys. Rev. Lett. 119, 173402 (2017).
  * Liu  _et al._ (2020) Z. C. Liu, N. P. Inman, T. J. Carroll, and M. W. Noel, Phys. Rev. Lett. 124, 133402 (2020).
  * Ryabtsev  _et al._ (2018) I. I. Ryabtsev, I. I. Beterov, D. B. Tretyakov, E. A. Yakshina, V. M. Entin, P. Cheinet, and P. Pillet, Phys. Rev. A 98, 052703 (2018).
  * Beterov  _et al._ (2018) I. I. Beterov, I. N. Ashkarin, E. A. Yakshina, D. B. Tretyakov, V. M. Entin, I. I. Ryabtsev, P. Cheinet, P. Pillet, and M. Saffman, Phys. Rev. A 98, 042704 (2018).
  * Cheinet  _et al._ (2020) P. Cheinet, K. L. Pham, P. Pillet, I. I. Beterov, I. N. Ashkarin, D. B. Tretyakov, E. A. Yakshina, V. M. Entin, and I. I. Ryabtsev, Quantum Electronics 50, 213 (2020).
  * Pelle  _et al._ (2016) B. Pelle, R. Faoro, J. Billy, E. Arimondo, P. Pillet, and P. Cheinet, Phys. Rev. A 93, 023417 (2016).
  * Westermann  _et al._ (2006) S. Westermann, T. Amthor, A. L. de Oliveira, J. Deiglmayr, M. Reetz-Lamour, and M. Weidemüller, The European Physical Journal D - Atomic, Molecular, Optical and Plasma Physics 40, 37 (2006).


# Phase Distribution in Probabilistic Movement Primitives, Representing Time
Variability for the Recognition and Reproduction of Human Movements

Vittorio Lippi, Raphael Deimel
Technische Universität Berlin, Fachgebiet Regelungssysteme, Sekretariat EN11,
Einsteinufer 17, -10587 Berlin, Germany
<EMAIL_ADDRESS><EMAIL_ADDRESS>

###### Abstract

Probabilistic Movement Primitives (ProMPs) are a widely used representation of
movements for human-robot interaction. They also facilitate the factorization
of temporal and spatial structure of movements. In this work we investigate a
method to temporally align observations so that when learning ProMPs,
information in the spatial structure of the observed motion is maximized while
maintaining a smooth phase velocity. We apply the method on recordings of hand
trajectories in a two-dimensional reaching task. A system for simultaneous
recognition of movement and phase is proposed and performance of movement
recognition and movement reproduction is discussed.

##  1 INTRODUCTION

### Overview.

Probabilistic movement primitives (ProMP) [Paraschos et al., 2013] are a
representation of movements used in robot control and human-robot interaction
(HRI) applications [Maeda et al., 2017] and provide several desirable
properties to model tasks for robot control and HRI[Paraschos et al., 2018].
One of these properties is temporal scaling: the movement trajectory is not a
direct function of time but function of a phase variable, $\phi(t)$, i.e. the
temporal evolution of $\phi$ determines the the velocity of the movement
independently of the spatial structure.

Figure 1: Features and Phase. In (a) the features $\Phi_{i}(\phi)$ are
represented as a function of the phase. In (b) a phase profile, based on the
beta function as a function of time. In (c) the features are represented as a
function of the time.

A ProMP represents a sample movement $y(t)$ as

$y(t)=\left[\begin{array}[pos]{c}q\\\
\dot{q}\end{array}\right]=\Phi(\phi(t))w+\epsilon$ (1)

where $q$ is the vector of variables describing the movement (usually joint
angles or hand effector position and orientation) and $\Phi$ is a $Nx2$ vector
computed with Gaussian functions and their derivatives:

$\Phi_{i}=\frac{e^{(\phi(t)-c_{i})^{2}/h_{i}}}{\sum^{N}_{k=1}e^{(\phi(t)-c_{k})^{2}/h_{k}}}$
(2)

where $c_{i}$ represents the center and $h_{i}$ expresses the spread of the
bell-shaped feature. The vector $w$ is drawn from a multinomial distribution
which is defined by its parameters $\theta$:

$p(w|\theta)=\mathcal{N}(\mu_{w},\Sigma_{w})$ (3)

In a general formulation, see for example [Colomé et al., 2014] the $\epsilon$
in Eq. 1 is characterized by a distribution such as
$p(\epsilon)=\mathcal{N}(0,\Sigma_{y})$. In this paper we assume that all the
variability observed on $y$ should be accounted by the distribution of $w$ and
that $\epsilon\approx 0$. In this way each observed trajectory $y_{i}$ is
represented by a single weight vector $w_{i}$. Vice versa, $w$ can be
interpreted as a compressed representation of the movement trajectory,
obtained from projecting trajectory values onto a lower dimensional subspace
using the Moore-Penrose pseudo-inverse of $\Phi$:

$w=\Phi^{+}y$ (4)

Which weight vector is most likely depends on the optimization criterium, e.g.
jerk minimization as regularization principle [Paraschos et al., 2018].
Alternative features for movement representation can be designed to guarantee
other properties such as minimal error at the final position and speed [Lippi
et al., 2012]. Due to the probabilistic nature of the ProMP framework
constrains on generated trajectories can be imposed by conditioning the
probability distribution[Paraschos et al., 2018]. Figure 1 illustrates the
features as either a function of time or phase. The time modulation performed
by the phase function also modulates the resolution of the low dimensional
representation during movement execution. In Figure 1 (c) features in the
middle of the movement are more dense than at the beginning or at the end, due
to the choice of a sigmoidal phase profile. We assume that part of the
variability observed among the repetitions of a movement is caused by shifted
and stretched phase signals. Time scaling can be performed by means of a
linear constant, i.e. $\phi=\alpha t$ [Maeda et al., 2017, Ewerton et al.,
2015]. In this work we estimate phase profile parameters for each observation
that minimize the variances over all observed $w$ (or, equivalently of the
response in phase domain), under the constraint of an assumed structure for
the phase profile. In particular we propose to restrict phase signals to the
beta function:

$\phi_{\Delta_{1},\Delta_{2}}(t)=\int_{0}^{t}\beta_{2,2}\left(\frac{\tau-\Delta_{1}}{\Delta_{1}-\Delta_{2}}\right)d\tau$
(5)

where $\beta_{2,2}$ is the incomplete beta function with parameters $a=b=2$.
This sigmoidal function, shown in 1 (b), is differentiable, monotonic and
saturates at $0$ and $1$. The parameters $\Delta_{1}$, $\Delta_{2}$ associated
with each trajectory depend on the distribution of the whole sample set.
Notice that weights $w$ are in general specific w.r.t. to the chosen phase
profile (see Eq. 4). The average of the sample set distribution is defined as:

$\overline{y}=\underset{\widetilde{y}}{argmin}\sum_{i=1}^{N}\int_{0}^{1}\left(y_{i}(\phi)-\widetilde{y}\left(\phi\right)\right)^{2}d\phi$
(6)

For the $i^{th}$ sample the phase is defined as:

$\phi_{i}=\underset{\widetilde{\phi}}{argmin}\int_{0}^{1}\left(y_{i}\left(\widetilde{\phi}\left(\tau(\phi)\right)\right)-\overline{y}(\phi)\right)^{2}d\phi$
(7)

The term $\tau(\phi)$ represents alignment between the average $\overline{y}$
and the $y_{i}$ sample, it express the fact that as $\widetilde{\phi}$ and
$\phi$ are both monotonically dependent on $t$, the former can be expressed as
a function of the latter, in particular:

$\tau(\phi)=arg_{t}\left(\phi(t)=\widetilde{\phi}(t)\right)$ (8)

As $\overline{y}$ depends on all phase profiles, Eq. 6 and Eq.7 need to be
optimized together to find the optimal phase profiles. The problem is solved
by iteratively optimizing $\overline{\phi}(t)$ and $\phi_{i}$, which yields
phase profile parameters for each sample movement. The phase profile
parameters of the average trajectory is fixed a priori to
$\Delta_{1}=\Delta_{2}=0$, as the time scaling when comparing two or more
movements (e.g. when computing the average) depends on the relationship
between the respective phase profiles and hence there is a degree of
redundancy. The described time scaling can also be performed on time domain
values. In order to take into account the movement primitive representation
the $y$ in Eq. 6 and Eq. 7 can be projected on the MP representation
$y_{mp}=\Phi w$. Once optimized time scaling parameters of the sample set are
obtained, we can obtain a probability distribution describing the movements
using empirical estimators from literature [Maeda et al., 2017, Paraschos et
al., 2018]. In the following paragraphs we demonstrate (a) how to identify a
model for movements, (b) how to recognize an observed movement given a set of
movement primitives, (c) how to estimate the current phase, (d) how to
integrate perception in the phase recognition process and (e) how to generate
a movement.

### Model identification.

The estimation of ProMP parameters from of a sample set can be based on
different principles, e.g. linear regression of each observation individually
[Paraschos et al., 2013] or maximizing the likelihood of the complete observed
data[Paraschos et al., 2018]. As introduced in the overview in this work we
assume $\epsilon\approx 0$ in Eq. 4 and that the phase profile can be chosen
to minimize the variability on $w$. This leads to a two step procedure: First,
each observed trajectory is assigned a set of parameters
$(w_{i},\Delta_{1i},\Delta_{2i})$. This is done on the basis of Eq. 6 and Eq.
7 to obtain the phase profile, and Eq. 4 to obtain $w$, then the distributions
for $w$ (Eq. 2) and $\phi(t)$ is obtained with an empirical estimator. An
example of model identification is shown in Fig.2. The training set is
generated artificially and consists of planar motion with random time profile
(with the form specified in Eq.5) and non-intersecting parabolic trajectories.
In Fig.2 (b) movements generated with the identified distribution are shown.
Notice that assuming a normal distribution of $w$ produces a distribution of
trajectories that is different from the more “uniform” training set.

### Movement recognition

An observed movement $y$ can be classified by comparing likelihoods between
ProMPs. The likelihood

$\begin{array}[h]{l}p\left(y(t_{i})|\mu_{k},\Sigma_{k}\right)\\\
=\int_{0}^{1}p\left(y(t_{i})\,|\,\Phi(\phi)\mu_{k},\Phi(\phi)\Sigma_{k}\Phi(\phi)^{T}\right)p(\phi|t_{i})d\phi\end{array}$
(9)

of observing the movement given model $k$ is computed and the most likely
model is selected, i.e. $L=\underset{k}{argmax}\,p(y|\theta_{k})$. Notice that
the probability of $y(t_{i})$ is considered independent of the previous
observations, i.e. $p(y(t_{i})|\theta_{k},y(t_{i-1}))=p(y(t_{i})|\theta_{k})$.
The probability of an observed sequence can be obtained by multiplying the
probabilities of the observed samples,
$p(y|\theta_{k})=\prod_{0}^{t_{f}}p(y_{t}|\theta_{k})$. In the described
examples we will include all the past observations of $y_{t}$, in the general
case not all points of the trajectory may be available, due to occlusion for
example. Applying the ProMP framework to such cases is demonstrated in [Maeda
et al., 2017].

### Phase recognition

Using the Bayes rule we can estimate the most likely phase of an observed
sample:

$p(\phi|y(t_{i}))=\frac{p(y(t_{i})|\Phi(\phi)\mu_{k},\Phi(\phi)\Sigma_{k}\Phi(\phi)^{T})p(\phi|t_{i})}{\sum_{k}^{M}p(\theta_{k})p(\phi|y(t),\theta_{k})}$
(10)

Notice that $p(\phi|t_{i})$ is assumed to be normally distributed and hence
parameterized by mean $\mu_{\phi}(t)$ and standard deviation
$\sigma_{\phi}(t)$ as functions over time. $\Delta_{1}$ and $\Delta_{2}$ are
not normally distributed due to the nonlinear relationship in Eq. 5. In most
prior work on the topic, $\phi$ is assumed to be a linear function over time
and that therefore $\phi(t)$ can be described by a normal distribution, e.g.
in [Ewerton et al., 2015]. Dynamic time warping (DTW) is another, nonlinear,
approach proposed recently [Ewerton et al., 2018].

### Perception.

Eq. 10 provides an estimate given a time-dependent but otherwise fixed
$p(\phi|t_{i})$ derived from the training set. Alternatively, we can use an
explicit phase estimate provided by some perceptual system. Examples of fusing
ProMP and perception have been proposed, e.g. in [Dermy et al., 2019] ProMPs
are used in the context of predicting human intentions. In general the sensor
input will have many more dimensions than $y$. In this work we will estimate
the phase from recent observations of the motion. We will consider a phase
estimator based on previous values of $y$ over a sliding time window to yield
an estimate $\phi^{*}(t_{i})$ of the current phase. Such estimator is used to
compute the probability of the observed $y(t_{i})$. Instead of integrating
over the distribution of probability of $\phi$ as shown in Eq. 9, only the
estimated value is used to compute $y(t_{i})$:

$\begin{array}[h]{l}p(y(t_{i})|\mu_{k},\Sigma_{k})\\\
=p\left(y(t_{i})\,|\,\Phi(\phi^{*}(t_{i}))\mu_{k},\Phi(\phi^{*}(t_{i}))\Sigma_{k}\Phi(\phi^{*}(t_{i}))^{T}\right)\end{array}$
(11)

### Movement Generation.

The model can be used to generate movements. Depending on the task the
movements can be generated deterministically using Eq. 1 or stochastically by
sampling $w$ from the distribution. In certain cases the term $\epsilon$ in
Eq. 1 may be nonzero due to actuation and external noise, but is not the case
in the presented example. Generating movements from ProMP can be used both for
robot control and to predict human movements in HRI, see, for example, [Oguz
et al., 2017]. Samples can obey constraints such as a specific target position
or velocity, by conditioning the distribution to the constraint prior to
sampling [Paraschos et al., 2013], which can be done parameterically (and
therefore efficiently) with normal distributions. A constraint can be
expressed as a target point $y^{*}(\phi_{i})$ to be reached at a given phase.
Exploiting the linear dependency between $w$ and $y$ the updated parameters
for $p(w|y_{\phi_{i}}=y^{*}(\phi_{i}))$ are

$\mu^{*}_{w}=\mu_{w}+\Sigma_{w}\Phi^{T}\\!(\phi)\left(\Phi(\phi)\Sigma_{w}\Phi^{T}\\!(\phi)\right)^{-1}\left(y^{*}-\Phi(\phi)\mu_{w}\right)$
(12)

$\Sigma^{*}_{w}=\Sigma_{w}+\Sigma_{w}\Phi^{T}\\!(\phi)\left(\Phi(\phi)\Sigma_{w}\Phi^{T}\\!(\phi)\right)^{-1}\Phi(\phi)\Sigma_{w}$
(13)

An example of generated movements is shown in Fig.2 (b,c,d). The conditioned
distributions in 2 (c) exhibit very small $\Sigma_{w}$ (i.e. all the
trajectories starting from a given point are indistinguishable in shape),
capturing the regularity of the training set. In Fig.2 (d) the distribution is
forced to pass through two via-points that are unlikely to be observed in the
same movement, i.e. it is conditioned on an outlier. As a consequence, the
resulting sampled trajectories do not resemble the observations presented in
the training set.

|
---|---
|
Figure 2: Parabolic lines synthetic data set. In (a) an artificial data-set of
trajectories with a parabolic shape and random duration is shown, with only 10
of the 100 produced samples shown but covering more intermediate positions.
The color gradient indicates the temporal evolution of the trajectory, the
black line the average trajectory in time. The resulting shape does not match
the presented examples well, especially in the end when sample count drops.
The data-set has been used as a training set for a proMP model. In (b)
trajectories generated with the proMP distribution are shown, the black line
represent the average (computed respect to phase). Notice that while the model
captured very accurately the shape of the movements, the gaussian distribution
of parameters results in a distribution of trajectory that is different from
the one presented in the training set. In (c) conditioned trajectories with
imposed via-points a time $t=0$ are shown, indicated by black stars. The
conditioned $\Sigma_{w}$ is so small that the trajectories generated for each
via-point are indistinguishable. In (d) the trajectory is conditioned to pass
through two via-points unlikely to appear in the same movement. The resulting
shape is dissimilar to those presented in the training set. For reference,
trajectories from the training set are shown in gray.

##  2 Dataset and Benchmark Problem

Figure 3: Reaching data set. A two-dimensional reaching task is performed by a
test participant on a computer screen using a mouse. The participant was asked
to perform the movement as fast as possible to reach circular targets that
were appearing in a position randomly decided among $4$. The resulting set is
composed by $25$ repetitions of $4$ distinct targets.

The Reaching data set consists in a reaching task: Participants have been
asked to move the mouse cursor from a starting point to targets appearing in 4
different positions. Each target identifies a different movement, associated
to a different ProMP model. The samples are shown in Fig. 3. Given the nature
of the task, consisting in reaching movements, it comes natural to perform the
segmentation on the basis of events, in particular: (a) a movement starts when
the mouse cursor is in the starting position and the target is shown on the
screen and (b) it stops when a target is reached. In general different
heuristics are possible for segmentation, notably the ProMP framework provides
the possibility to segment the signal exploiting an expectation-maximization
algorithm, treating segmentation as a latent variable to be optimized together
with the movement models parameters[Lioutikov et al., 2017]. The dataset
includes 100 movement samples used as training-set for the model
identification and $100$ samples used as test set for the movement recognition
and phase recognition. Movements were sampled at $100$ Hz. The number of
features was set to $N=9$, with parameters, in Eq. 2, $c_{k}=(k-1)/(N-1)$ and
$h_{k}=0.15$.

##  3 Results

Figure 4: Dataset ”mov 1” (see Fig. 3). (a) shows vertical and horizontal
coordinates over time, gray lines represent training data, black line the
average. In (b) the corresponding phase profiles are shown, in gray the ones
associated with the observations, in black the average one. In (c) the
resulting movements over phase are shown. Figure 5: Movement and phase
recognition accuracy evaluated on a test set of 100 samples (25 for each of
the 4 recognized movements). In (a) the relative error is shown as a function
of time (blue). As the different test samples have different durations. The
number of available movement samples compared at a given time is shown
(orange). The classification error decreases from random (75%) as more of the
sample movement is observed. In (b) the estimated phases are shown (gray).
Notice that the likelihood of every phase value is computed at each point in
time using Eq. 10 and the most likely phase is selected. The resulting phase
signal differs from the ones assumed to generate the training set (Eq. 5), and
in general is not guaranteed to be monotonic. Nevertheless, the average phase
(black) recovers the typical sigmoidal profile.

### Model identification, movement recognition and phase recognition.

The distribution of phase profiles obtained with Eq. 6, Eq. 7 is shown in Fig
4 (b). Each sample is transformed in an MP, i.e. the parameters $(w,\phi)$ are
extracted, $4$ ProMP models representing the $4$ movements are estimated. The
observed movement is identified as belonging to the most likely model to
produce it according to Eq. 9. The movement recognition is performed at each
time step. The classification accuracy increases with time as shown in Fig. 5
(a), and reaches $100\%$ at the end of the movement when movements are
unambiguously differentiated by their position. Phase recognition is also
performed for each time step using Eq. 10. The phase estimate is obtained
sample-by-sample and need not be monotonic as shown in 5 (b).

### Movement Generation.

|
---|---
|
Figure 6: Movements generated using the ProMP model for the observations set
”mov 1” (see Fig. 3). Trajectories in red are samples from the distribution
over $w$. Trajectories in blue are samples constrained to pass through a via-
point at phase $\phi^{*}=0.5$, indicated by a black star. Thick lines indicate
the average trajectories of the respective distributions. The red circle marks
the final target. In (a) the via-point is within the space of the presented
examples and hence the trajectories resemble thoses presented in the training
set. In (b) the via-point is an outlier to the positions presented in the
training set and hence the conditioned trajectories are not representative of
the observed movements. In (c) the conditioned distribution is used together
with classification and phase identification to produce a prediction of the
motion in progress (black dotted line). The most likely trajectory is plotted
with the blue thick line, and samples as thin blue lines. In (d) the
prediction is performed using a neural network to estimate the phase.

In order to generate concrete movements the parameters $w$ are sampled from
the distribution. Fig. 6 shows movements generated with and without
constraints. The constraint is expressed as a via-point to be reached at a
given phase. The constrained movements are generated from the conditioned
distribution with parameters from Eq. 12, Eq. 13. In Fig. 6(a) the imposed
via-point is within the space of the presented examples and hence trajectories
resemble those presented in the training set. In Fig. 6 (b) the via-point is
more distant from the positions presented in the training set and hence the
conditioned trajectories are less representative of the original distribution,
e.g. forming loops and missing the final target. In (c) the conditioned
distribution is used to predict a continuation of the current motion. The
prediction is also affected by the uncertainty about the current phase, which
is estimated according to Eq. 10.

### Perception.

Perceptual estimation of phase was implemented with a feedforward neural
network with three layers with 40, 20 and 10 neurons respectively, taking as
input the current time and a vector of 20 previous samples. The network has
been trained with the Levenberg-Marquardt backpropagation algorithm[Levenberg,
1944, Marquardt, 1963]. The results are shown in Fig. 7 where the
classification task (shown in Fig. 7) is repeated by exploiting the estimated
phase, according to Eq. 11. Classification performance using the neural
network is comparable to that obtained by integrating the phase distribution
(Eq. 9). In Fig. 6 (d) the phase estimated by the neural network is used to
perform a prediction on the basis of the observed trajectory. The result is
similar to that which is shown in Fig. 6 (c).

Figure 7: Classification and phase evaluation accuracy using a neural network
for the identification of phase. The training set and the test set are the
same as used in Fig. 5. In (a) the relative error is shown as a function of
time (blue). As the movements have different durations, the number of
available movement samples at a given time is indicated in orange. The
classification accuracy increases as more of the test movements are observed.
In (b) the identified phase profiles are shown in gray, where phase is
computed at each time step. The average phase (black) shows the typical
sigmoidal profile. Although in some cases the phase appears less accurate
using the neural network estimator, overall classification performance is
comparable to that of Fig. 5.

### A 7-DoF robot arm example.

In additon to previous low-dimensional examples, model identification is
applied to movements of a 7-DoF robotic arm.

Figure 8: The robot task configuration (A). The robot’s gripper is holding a
hook (B) which it uses to pick up the object (C).

The task consists in picking up an object using a hook and handing it to an
user in a single continuous movement. The setup is shown in Fig. 8. The
object’s loop represents a spatial constraint for the trajectory, since the
object is always picked up in the same position. Trajectory examples were
recorded by manually moving the robot, producing a total of $11$ examples. The
task was represented in seven dimensional joint space. The number of features
was set to $9$, hence the ProMP distribution has $63$ dimensions. Fig. 9.
shows the covariance matrices of the ProMP’s distributions learned on data
without and with optimized phase profiles. The use of phase modulation
decreases variation in adjacent feature weights in the covariance matrix
elements as temporal noise is externalized to the phase signal. A more compact
description of the concept is shown in Fig. 10, where the eigenvalues of the
two covariance matrices, i.e. the principal components, are compared. All the
non-zero eigenvalues are smaller with phase modulation.

|
---|---
Figure 9: Color representation of the covariance matrix for $w$, for the robot
arm movements set without phase modulation (left) and with phase scaling
(right). Most of the values in the covariance matrix are smaller. Figure 10:
Non-zero eigenvalues of the covariance matrices shown in Fig. 9. With the
phase modulation all eigenvalues are smaller, which indicates a more precise
reproduction of movements.

##  4 Discussion and Future Work

In this work we used ProMPs to model human movements and demonstrations to
robots in order to perform movement recognition and movement generation, with
particular emphasis on leveraging a parametric phase profile to modulate
movement execution speed. The approach showed to properly classify and
generalize the movements in a data-set of reaching movements. The proposed
parametric phase profile has a sigmoidal shape that is particularly suitable
to reaching tasks such as the one presented into the example and other ”real
world” tasks, such as pick and place operations where it is relevant to the
control of initial and final position. However, a sigmoidal phase profile
enforces zero velocity at the beginning and at the end of the movement.
Therefore, in tasks where final velocity must be controlled, other profiles
(e.g. a beta function with $a=2$, $b=1$) may be favorable. In general it
should be considered that time scaling of human movements can work differently
depending on the task and the working conditions[Zhang and Chaffin, 1999]. The
example reaching task is represented as a library of $4$ ProMPs describing
movements to $4$ possible targets. In general, the diversity of movements in a
task can be represented with a manifold in the space of hyper-parameters and
mixture distributions [Rueckert et al., 2015], or specifically with gaussian
mixtures describing the distribution of $w$ [Ewerton et al., 2015]. The
parameters can be constrained with linear relationships for dimensionality
reduction, as recently shown in [Colomé and Torras, 2018]. With human
movements, dependencies between target positions and the moving parameters can
be represented by a linear model[Avizzano and Lippi, 2011, Lippi et al.,
2012]. In all those cases phase can be estimated online while observing the
movement by modifying Eq. 10 to take into account the chosen distribution.
Further, a specific perceptual estimator dedicated to the identification of
phase using a feed-forward neural network was described and implemented. The
estimator led to a performance in terms of classification and prediction
accuracy that was comparable to that of the system exploiting the empirical
probability distribution of the phase directly. An advantage of using a neural
network is, that it avoids the potentially costly computation of the integral
in Eq.9 at each integration step. Future work on phase recognition will be
focus on the identification of state of dynamic system models tasked with
generating phases, e.g. phase-state machine models [Deimel, 2019]. Such models
could provide phase profiles and could be used in conjunction with a library
of ProMPs to control a task. Furthermore, current research on ProMPs includes
the design of feedback control systems considering physical interaction with
the environment and the users [Paraschos et al., 2018, Paraschos et al.,
2013]. Time scaling does not apply to the control of arm kinematics and the
external contact forces in the same way as for movements (this holds for each
kind of rescaling, also the linear one). A general framework that can scale
kinematics and contact instants in time should be integrated with a consistent
control of applied forces and torques.

## ACKNOWLEDGEMENTS

We gratefully acknowledge financial support for the project MTI-engAge
(16SV7109) by BMBF

## REFERENCES

  * Avizzano and Lippi, 2011 Avizzano, C. A. and Lippi, V. (2011).  A regression model for digital representation of juggling.  In BIO Web of Conferences, volume 1, page 00006. EDP Sciences.
  * Colomé et al., 2014 Colomé, A., Neumann, G., Peters, J., and Torras, C. (2014).  Dimensionality reduction for probabilistic movement primitives.  In Humanoid Robots (Humanoids), 2014 14th IEEE-RAS International Conference on, pages 794–800. IEEE.
  * Colomé and Torras, 2018 Colomé, A. and Torras, C. (2018).  Dimensionality reduction in learning gaussian mixture models of movement primitives for contextualized action selection and adaptation.  IEEE Robotics and Automation Letters, 3(4):3922–3929.
  * Deimel, 2019 Deimel, R. (2019).  A dynamical system for governing continuous, sequential and reactive behaviors.  In Proceedings of the Austrian Robotics Workshop, page (in press).
  * Dermy et al., 2019 Dermy, O., Charpillet, F., and Ivaldi, S. (2019).  Multi-modal intention prediction with probabilistic movement primitives.  In Human Friendly Robotics, pages 181–196. Springer.
  * Ewerton et al., 2015 Ewerton, M., Neumann, G., Lioutikov, R., Amor, H. B., Peters, J., and Maeda, G. (2015).  Modeling spatio-temporal variability in human-robot interaction with probabilistic movement primitives.  In Workshop on Machine Learning for Social Robotics, ICRA.
  * Ewerton et al., 2018 Ewerton, M., Rother, D., Weimar, J., Kollegger, G., Wiemeyer, J., Peters, J., and Maeda, G. (2018).  Assisting movement training and execution with visual and haptic feedback.  Frontiers in neurorobotics, 12:24.
  * Levenberg, 1944 Levenberg, K. (1944).  A method for the solution of certain non-linear problems in least squares.  Quarterly of applied mathematics, 2(2):164–168.
  * Lioutikov et al., 2017 Lioutikov, R., Neumann, G., Maeda, G., and Peters, J. (2017).  Learning movement primitive libraries through probabilistic segmentation.  The International Journal of Robotics Research, 36(8):879–894.
  * Lippi et al., 2012 Lippi, V., Avizzano, C. A., and Ruffaldi, E. (2012).  A method for digital representation of human movements.  In RO-MAN, 2012 IEEE, pages 404–410. IEEE.
  * Maeda et al., 2017 Maeda, G. J., Neumann, G., Ewerton, M., Lioutikov, R., Kroemer, O., and Peters, J. (2017).  Probabilistic movement primitives for coordination of multiple human–robot collaborative tasks.  Autonomous Robots, 41(3):593–612.
  * Marquardt, 1963 Marquardt, D. W. (1963).  An algorithm for least-squares estimation of nonlinear parameters.  Journal of the society for Industrial and Applied Mathematics, 11(2):431–441.
  * Oguz et al., 2017 Oguz, O. S., Sari, O. C., Dinh, K. H., and Wollherr, D. (2017).  Progressive stochastic motion planning for human-robot interaction.  In Robot and Human Interactive Communication (RO-MAN), 2017 26th IEEE International Symposium on, pages 1194–1201. IEEE.
  * Paraschos et al., 2018 Paraschos, A., Daniel, C., Peters, J., and Neumann, G. (2018).  Using probabilistic movement primitives in robotics.  Autonomous Robots, 42(3):529–551.
  * Paraschos et al., 2013 Paraschos, A., Daniel, C., Peters, J. R., and Neumann, G. (2013).  Probabilistic movement primitives.  In Advances in neural information processing systems, pages 2616–2624.
  * Rueckert et al., 2015 Rueckert, E., Mundo, J., Paraschos, A., Peters, J., and Neumann, G. (2015).  Extracting low-dimensional control variables for movement primitives.  In Robotics and Automation (ICRA), 2015 IEEE International Conference on, pages 1511–1518. IEEE.
  * Zhang and Chaffin, 1999 Zhang, X. and Chaffin, D. B. (1999).  The effects of speed variation on joint kinematics during multisegment reaching movements.  human movement science, 18(6):741–757.