# Entropy measurement of a strongly coupled quantum dot

Timothy Child<EMAIL_ADDRESS>Stewart Blusson Quantum Matter Institute,
University of British Columbia, Vancouver, British Columbia, V6T1Z4, Canada
Department of Physics and Astronomy, University of British Columbia,
Vancouver, British Columbia, V6T1Z1, Canada  Owen Sheekey  Stewart Blusson
Quantum Matter Institute, University of British Columbia, Vancouver, British
Columbia, V6T1Z4, Canada  Department of Physics and Astronomy, University of
British Columbia, Vancouver, British Columbia, V6T1Z1, Canada  Silvia Lüscher
Stewart Blusson Quantum Matter Institute, University of British Columbia,
Vancouver, British Columbia, V6T1Z4, Canada  Department of Physics and
Astronomy, University of British Columbia, Vancouver, British Columbia,
V6T1Z1, Canada  Saeed Fallahi  Department of Physics and Astronomy, Purdue
University, West Lafayette, Indiana, USA  Birck Nanotechnology Center, Purdue
University, West Lafayette, Indiana, USA  Geoffrey C. Gardner  Birck
Nanotechnology Center, Purdue University, West Lafayette, Indiana, USA  School
of Materials Engineering, Purdue University, West Lafayette, Indiana, USA
Michael Manfra  Department of Physics and Astronomy, Purdue University, West
Lafayette, Indiana, USA  Birck Nanotechnology Center, Purdue University, West
Lafayette, Indiana, USA  Elmore Family School of Electrical and Computer
Engineering, Purdue University, West Lafayette, Indiana, USA  School of
Materials Engineering, Purdue University, West Lafayette, Indiana, USA  Andrew
Mitchell  School of Physics, University College Dublin, Belfield, Dublin 4,
Ireland  Centre for Quantum Engineering, Science, and Technology, University
College Dublin, Dublin 4, Ireland  Eran Sela  Raymond and Beverly Sackler
School of Physics and Astronomy, Tel-Aviv University, IL-69978 Tel Aviv,
Israel  Yaakov Kleeorin  Center for the Physics of Evolving Systems,
University of Chicago, Chicago, IL, 60637, USA  Yigal Meir  Department of
Physics, Ben-Gurion University of the Negev, Beer Sheva 84105, Israel  The
Ilse Katz Institute for Nanoscale Science and Technology, Ben-Gurion
University of the Negev, Beer Sheva 84105, Israel  Joshua Folk
<EMAIL_ADDRESS>Stewart Blusson Quantum Matter Institute, University of
British Columbia, Vancouver, British Columbia, V6T1Z4, Canada  Department of
Physics and Astronomy, University of British Columbia, Vancouver, British
Columbia, V6T1Z1, Canada

###### Abstract

The spin 1/2 entropy of electrons trapped in a quantum dot has previously been
measured with great accuracy, but the protocol used for that measurement is
valid only within a restrictive set of conditions. Here, we demonstrate a
novel entropy measurement protocol that is universal for arbitrary mesoscopic
circuits and apply this new approach to measure the entropy of a quantum dot
hybridized with a reservoir, where Kondo correlations typically dominate spin
physics. The experimental results match closely to numerical renormalization
group (NRG) calculations for small and intermediate coupling. For the largest
couplings investigated in this work, NRG predicts a suppression of spin
entropy at the charge transition due to the formation of a Kondo singlet, but
that suppression is not observed in the experiment.

Entropy is a powerful tool for identifying exotic quantum states that may be
difficult to distinguish by more standard metrics, like conductance. For
example, bulk entropic signatures in twisted bilayer graphene indicate that
carriers in some phases with metallic conductivity retain their local moments,
as would normally be associated with a Mott insulator Saito  _et al._ (2021);
Rozen  _et al._ (2021); Lian (2021). Entropy has also been proposed as a tell-
tale characteristic of isolated non-abelian quasiparticles, whether Majorana
modes in a superconductor Smirnov (2015); Sela  _et al._ (2019) or excitations
of a fractional quantum Hall state Cooper and Stern (2009); Yang and Halperin
(2009); Ben-Shach  _et al._ (2013), distinguishing them from abelian analogs.

Quantifying the entropy of single quasiparticles is challenging due to the
small signal size, of order $k_{B}$, but first steps in this direction have
been made in recent years Hartman  _et al._ (2018); Kleeorin  _et al._ (2019).
Ref. Hartman  _et al._ , 2018 employed Maxwell relations to measure the
$k_{B}\ln{2}$ spin entropy of a single electron confined to a quantum dot (QD)
in GaAs via the temperature-induced shift of a Coulomb blockade charge
transition. That approach relied on the assumption of weak coupling between
the QD and the reservoirs, in order to fit based on the specific charging
lineshape known for that regime. In that weak-coupling regime, spin states are
pristine enough to serve as spin qubits Elzerman  _et al._ (2004); Petta  _et
al._ (2005); Hanson  _et al._ (2007); Barthel  _et al._ (2009); Bluhm  _et
al._ (2010); Nowack  _et al._ (2011); Shulman  _et al._ (2012) but the
underlying physics is very simple.

The weak-coupling approach of Ref. Hartman  _et al._ , 2018 is not applicable
to a broad class of mesoscopic devicesPyurbeeva and Mol (2021), which limits
its value in probing the complex Hamiltonians that may be implemented in such
systems. For example, a single-impurity Kondo effect is realized when the
localized spin is strongly coupled to a reservoir Goldhaber-Gordon  _et al._
(1998); Pustilnik and Glazman (2004). Recently, more complicated structures
including multiple dots have been engineered to host multi-channel Kondo
states Potok  _et al._ (2007); Keller  _et al._ (2015), a three-particle
simulation of the Hubbard model Dehollain  _et al._ (2020), and more. Entropy
measurements made on any of these systems would offer a significant advance in
their understanding, but are not yet possible.

Here, we develop a universal protocol for mesoscopic entropy measurement that
forgoes the simplifying assumptions of Ref. Hartman  _et al._ , 2018, then
apply it to investigate the entropy of the first electron as it enters a
quantum dot when strongly hybridized with a reservoir, the regime where Kondo
correlations are expected. The protocol is based on a Maxwell relation
expressed in integral form,

$\Delta S_{\mu_{1}\to\mu_{2}}=\int_{\mu_{1}}^{\mu_{2}}\frac{dN(\mu)}{dT}d\mu,$
(1)

that provides access to the entropy change, $\Delta S$, with chemical
potential $\mu$, based on measurements of the change in quantum dot
occupation, $N$, with temperature, $T$ Sela  _et al._ (2019). When the
coupling, $\Gamma$, between dot and reservoir is weak ($\Gamma\ll k_{B}T$),
the characteristic temperature of Kondo correlations is much less than the
experimental temperature, $T_{K}\ll T$, and the measurement matches well to
single-particle approximations. When $\Gamma\gtrsim k_{B}T$, the onset of
entropy as the electron enters the dot is strongly modified. The measurement
of this modified entropy signature is the primary result of this work.

Figure 1:  a) Scanning electron micrograph (colour added) of the device.
Electrostatic gates (gold) define the circuit in a 2D electron gas. Grey gates
were unused and remained grounded. Squares represent ohmic contacts to the
2DEG. The thermal electron reservoir (red) was alternated between base and
elevated temperatures. b) Current through the charge sensor for the $0\to 1$
charge transition in a weakly coupled regime, separated into the unheated
(100mK) and heated (130mK) parts of the interlaced measurement Child  _et al._
(2022), showing $N=0\rightarrow 1$ step height $I_{e}$. c,d) Change in
$I_{CS}$ from 100 to 130 mK, for weak (c) and strong (d) coupling between QD
and reservoir. c) includes fit to weakly-coupled theory, yielding $\Delta
S_{\rm{fit}}=(1.02\pm 0.01)k_{B}\ln 2$.

Measurements were performed on a circuit (Fig. 1a) defined by voltages applied
to electrostatic gates over a GaAs 2D electron gas sup , similar in
functionality to that in Ref. Hartman  _et al._ , 2018, with improvements made
in the design and entropy measurement protocol Child  _et al._ (2022). The
circuit includes the QD, a charge sensing quantum point contact, and an
electron reservoir that can be rapidly Joule-heated above the chip temperature
$T$ to an elevated $T+\Delta T$. Both $T$ and $T+\Delta T$ were calibrated by
fitting to thermally broadened charge transitions; except where noted, $T=100$
mK with $\Delta T\sim 30$ mK. Coupling between the QD and the thermal
reservoir is via a single tunnel barrier, with transmission controlled by
$V_{T}$.

Coarse adjustments in the QD energy level were made using gate voltage
$V_{P}$, with fine adjustments controlled by $V_{D}$. In practice, $V_{P}$ and
$V_{D}$ tune the chemical potential as it enters into Eq. 1 by adjusting the
difference in energy between $\mu$ in the reservoir and the energy needed to
add an electron to the dotPyurbeeva and Mol (2021). Throughout this paper we
report $V_{D}$ with respect to to the midpoint of the charge transition,
$\Delta V_{D}\equiv V_{D}-V_{D}(N=1/2)$.

$N$ was monitored via the current, $I_{CS}$, through the nearby charge sensor
Elzerman  _et al._ (2004), which was biased with a DC voltage between 50 and
300 µV. Fig. 1b illustrates weakly coupled $N=0\rightarrow 1$ transitions at
two different temperatures, with $I_{e}$ defined as the net drop in $I_{CS}$
across the transition due to the addition of $1e$ charge to the dot.
Measurements at the two temperatures were interlaced by alternated Joule
heating of the reservoir at 25Hz to reduce the impact of charge instability,
then averaged over several sweeps across the charge transition, see Ref. sup,
.

Figure 1c shows the difference in the detector current in the weakly coupled
regime between 100 and 130 mK (shown in Fig. 1b), $\Delta I_{CS}(V_{D})\equiv
I_{CS}(T+\Delta T,V_{D})-I_{CS}(T,V_{D})$. Note that $-\Delta I_{CS}$ is
plotted instead of $\Delta I_{CS}$ in order to connect visually with $\Delta
N$, which increases when $I_{CS}$ decreases. As in Ref. Hartman  _et al._ ,
2018, the lineshape of $\Delta I_{CS}(V_{D})$ in Fig. 1c may be fit to a non-
interacting theory for thermally-broadened charge transitions to extract the
change in entropy across the transition, $\Delta S_{\rm{fit}}$, not requiring
calibration factors or other parameters (see Ref. Hartman  _et al._ , 2018 for
details). For the data in Fig. 1c, this yields $\Delta S_{\rm{fit}}=(1.02\pm
0.01)k_{B}\ln{2}$, where the uncertainty reflects standard error among 5
consecutive measurements at slightly different $V_{T}$.

The limitation of this approach is illustrated by the very different lineshape
in Fig. 1d, reflecting the $0\to 1$ transition when the QD is strongly coupled
to the reservoir and Kondo correlations begin to play a role. Fitting this
data to thermally-broadened theory would yield a meaningless $\Delta
S_{\rm{fit}}>10k_{B}\ln{2}$ for the spin-1/2 electron, motivating the integral
approach from Eq. 1.

Before moving to the quantitative evaluation of entropy via Eq. 1, we note
that the different lineshapes of $\Delta I_{CS}(V_{D})$ in Figs. 1c and d—the
peak-dip structure in Fig. 1c contrasting with the simple peak in Fig. 1d—can
be understood as representing two temperature regimes for the Anderson
impurity model. Fig. 1c represents the high temperature limit, where $dN/dT$
is approximately a measure of the energy derivative of the density of states
in the QD, and thus exhibits positive and negative lobes. At sufficiently low
temperatures, the exact solution Tsvelick and Wiegmann (1983), and the
resulting Fermi liquid theory Mora  _et al._ (2015) predict a positive $dN/dT$
for all values of the chemical potential, from the empty level to the Kondo
regime through the mixed-valence regime, with a peak expected at a chemical
potential corresponding to $T_{K}(\mu)\sim T$, where the entropy is expected
to crossover from $S=0$ to $S=k_{B}\log 2$. Fig. 1d, corresponding to data
taken when $T\ll\Gamma$, demonstrate such all-positive $dN/dT$.

For a quantitative extraction of entropy beyond the weakly-coupled regime of
Fig. 1c, we follow the integral approach in Eq. 1 that makes no assumptions on
the nature of the quantum state. Evaluating Eq. 1 also provides a measurement
of $\Delta S(\mu)$ that is continuous across the charge transition, rather
than just comparing $N=0$ to $N=1$ values. To evaluate Eq. 1 from experimental
data, $dN(\mu)/dT$ is approximated by the ratio $\Delta N(V_{D})/\Delta
T=-\Delta I_{CS}(V_{D})/(I_{e}\Delta T)$. $\Delta T$ is expressed in units of
gate voltage using the corresponding lever armsup  so that the integral may be
evaluated over $V_{D}$, giving $\Delta S(V_{D})$.

Figure 2:  a) Change of $S$ in the QD across the $N=0\rightarrow 1$
transition, obtained by integrating $\Delta I_{CS}(V_{D})$ (Fig. 1c) following
Eq. 1. Dot occupation across the transition is shown in grey. Data obtained in
the weakly coupled limit, $V_{T}=-350$ mV corresponding to $\Gamma/k_{B}T\sim
1\times 10^{-4}$. $\Delta S_{0\rightarrow 1}=(0.99\pm 0.02)k_{B}\ln 2$ is the
net change $\Delta S$ across the complete transition. b) Comparison of $\Delta
S_{\rm{fit}}$, $\Delta S_{0\rightarrow 1}$, and $\Delta S_{max}$ (see text)
for $V_{T}$ covering approximately $10^{-5}<\Gamma/k_{B}T<10^{-1}$. c)
Variation of charge step $I_{e}$ and lever arm $\alpha$, measured
independently over the full range of $V_{T}$ explored in this experiment.
Dashed line: extrapolation of $\alpha$ into the strongly-coupled regime where
it cannot be measured directly.

We begin by confirming the integral approach in the weakly-coupled regime,
where the physics is simple. Fig. 2a shows the entropy change across the
$N=0\rightarrow 1$ charge transition for $\Gamma\ll k_{B}T$, calculated from
the data in Fig. 1b using Eq. 1. The resulting $\Delta S(\mu)$ indicates that
the change in dot entropy is non-monotonic as the first electron is added,
reaching a $k_{B}\ln{3}$ peak before settling to $k_{B}\ln{2}$. The
$k_{B}\ln{3}$ peak just above $\Delta V_{D}=0$ reflects a combination of
charge and spin degeneracy on the dot in the middle of the charge transition,
with three possible microstates
$\\{\ket{N=0},\ket{N=1,\uparrow},\ket{N=1,\downarrow}\\}$ all equally
probable. Charge degeneracy is gone after the transition, but spin degeneracy
remains, leaving two possible microstates
$\\{\ket{N=1,\uparrow},\ket{N=1,\downarrow}\\}$. The net change in entropy
from beginning to end, $\Delta S_{0\to 1}=(0.99\pm 0.02)k_{B}\ln{2}$, is
nearly identical to the $\Delta S_{\rm{fit}}=(1.02\pm 0.01)k_{B}\ln{2}$ from
Fig. 1c, despite very different sources of error for the two approaches.

Figure 2b compares the fit and integral approaches for weakly-coupled charge
transitions covering four orders of magnitude in $\Gamma$, with the coupling
tuned by $V_{T}$ (see Fig. 3b inset for calibration of $\Gamma$). The
consistency between $\Delta S_{0\to 1}$ and $\Delta S_{\rm{fit}}$ over the
full range of weakly-coupled $V_{T}$, in addition to the fact that $\Delta
S_{max}$ remains $k_{B}\ln$ 3 throughout this regime, confirms the accuracy of
the integral approach. Small deviations from $\Delta S_{0\to 1}=\Delta
S_{\rm{fit}}=k_{B}\ln{2}$, such as that seen around $V_{T}$= -330 mV, are
repeatable but sensitive to fine-tuning of all the dot gates; we believe they
are due to extrinsic degrees of freedom capacitively coupled to the dot
occupation, such as charge instability in shallow dopant levels in the GaAs
heterostructure.

After confirming the accuracy of Eq. 1 in the weakly coupled regime, we turn
to the regime $\Gamma\gtrsim k_{B}T$ ($V_{T}>-280$ mV), where the influence of
hybridization and Kondo correlations is expected to emerge. Fig. 3 shows the
crossover from $\Gamma\ll k_{B}T$ to $\Gamma\gg k_{B}T$, illustrating several
qualitative features. The $k_{B}\ln{3}$ peak in $\Delta S(\mu)$ decreases with
$\Gamma$, until no excess entropy is visible at the charge degeneracy point
for $\Gamma/k_{B}T\gtrsim 5$ (Figs. 3a). This suppression of the entropy
associated with charge degeneracy originates from the broadening by $\Gamma$
of the $N=1$ level due to hybridization with the continuous density of states
in the reservoir Sela  _et al._ (2019). At the same time, the total entropy
change $\Delta S_{0\rightarrow 1}$ remains $\sim k_{B}\ln{2}$ over the entire
range of $\Gamma$ explored in this experiment, reflecting the entropy of the
spin-1/2 electron trapped in the QD.

Figure 3:  Evolution of $S(\mu)$ from the weak (black) to strong (red)
coupling regimes, comparing data (panel a) to NRG calculations (panel c).
Measurements of occupation across the charge transition are fit to NRG (panel
b), leaving no free fit parameters for the $S(\mu)$ calculation. Panel b
inset: Coupling strength of the QD to the reservoir, $\Gamma/k_{B}T$,
extracted from fits, across the full range of $V_{T}$. Values
$\Gamma/k_{B}T\ll 1$ cannot be measured directly and are extrapolated (dashed
line).

To make quantitative comparison between theory and experiment, we employ NRG
simulations Legeza  _et al._ (2008); Tóth  _et al._ (2008), that yield $N$ as
a function of both $T$ and $\epsilon_{0}$, where $-\epsilon_{0}$ is the depth
of the dot level below the reservoir potential. From $N(T,\epsilon_{0})$,
$dN/dT$ and therefore $\Delta S$ may be extracted via Eq. 1. In order to make
a direct comparison with the experiment,
$\Delta\epsilon_{0}\equiv\epsilon_{0}-\epsilon_{0}(N=1/2)$ is defined like
$\Delta V_{D}$, centred with respect to the charge transition. NRG parameters
are calibrated to match those in the measurements by aligning the occupation
$N(\Delta\epsilon_{0}$) with the measured $N(\Delta V_{D}$)sup , from which
the appropriate $\Gamma/T$ calculation may be selected and the precise
connection between $\Delta\epsilon_{0}$ with $\Delta V_{D}$ is ensured. As
seen in Fig. 3b, the data/theory agreement is perfect in terms of dot
occupation, within the experimental resolution, giving confidence that
measured and calculated $\Delta S$ may be compared directly.

The calculations in Fig. 3c illustrate NRG predictions for $\Delta
S(\epsilon_{0})$ over the range of $\Gamma$ accessible in our measurements.
Matching the data, the peak in entropy due to charge degeneracy is suppressed
as $\Gamma$ increases above $\Gamma\gtrsim k_{B}T$, while the net entropy
change across the transition remains $k_{B}\ln 2$. At the same time, a clear
qualitative difference between data and NRG is the shift to the right seen in
NRG curves for higher $\Gamma$ (Fig. 3c), but not observed in the measurements
(Fig. 3a). This relative shift of NRG with respect to data is not explained by
an offset of $\Delta\epsilon_{0}$ with respect to $\Delta V_{D}$, as the two
are aligned by the occupation data (Fig. 3b).

At first glance, the shift of NRG curves to the right (to larger chemical
potential) with increasing $\Gamma$ may seem strange: why would the emergence
of $k_{B}\ln 2$ entropy in the dot be delayed until after the electron is
mostly localized in the dot? This behaviour is explained by the virtual
exchange interactions underlying the Kondo effect, which form a quasi-bound
singlet state between the localized spin and a cloud of delocalized spin for
temperatures below $T_{K}$. This results in a state with no magnetic moment
Kondo (1964) and, in the case of a single-electron QD, zero entropy. Thus, due
to the Kondo effect, we expect the entropy to remain zero for all chemical
potentials that obey $T<T_{K}(\mu)$. Since $T_{K}\propto
e^{-\pi(\epsilon_{0}-\mu)/\Gamma}$ in the (experimentally relevant) large-$U$
limit, where $U$ represents the QD charging energy, we expect the onset of
$k_{B}$ln 2 entropy to shift to larger values of $\mu$ as $\Gamma$ increases,
as seen in the NRG results. (Note that NRG curves are independent of the
values of bandwidth, $W$, and of $U$, as long as $U,W\gg T,\Gamma$ because, as
with the experimental plots, all curves are shifted such that $x=0$
corresponds to half filling.)

For the present device geometry, $\Delta I_{CS}$ collapses with increasing
$\Gamma$, limiting the strength of the coupling accessible in the measurement
to $\Gamma/k_{B}T\lesssim 25$ and preventing us from reaching the deep Kondo
regime far from the charge transition. $T\ll T_{K}$ may still be achieved
close to the charge transition, however, in the so-called mixed-valence regime
where $|\epsilon_{0}-\mu|\lesssim\Gamma$ and the electron is not yet strongly
localized in the dot. It remains a puzzle why the strong suppression of
entropy right at the charge transition, seen in NRG calculations for
$\Gamma/k_{B}T\geq 5$, is not observed in the data.

It has been suggested Silva and Levit (2003) that the charge measurement
itself can lead to dephasing of the Kondo singlet. Ref. Avinun-Kalish  _et
al._ , 2004 reported a suppression of the Kondo peak in the conductance due to
charge sensing. Those results implied stronger dephasing than had been
predicted theoretically, and led to new theoretical work proposing alternative
dephasing mechanismsKang (2005). In order to test for charge-sensor dephasing
in our measurement, the experiment was repeated at $V_{CS}$ down to 50 $\mu$V,
but no dependence on $V_{CS}$ was seen in the data sup . The discrepancy
between NRG and measurements therefore points to the relevance of other
degrees of freedom, such as other single-electron levels, that have not been
included in current theories. In the future, experiments that allow
simultaneous transport and entropy characterization of the Kondo state may
help to resolve this puzzle. We note that the entropy measurement presented
here is the first that could be sensitive to the dephasing of the Kondo state
itself, rather than dephasing of transport through the Kondo resonanceKang and
Khym (2007).

ACKNOWLEDGEMENTS: This project has received funding from European Research
Council (ERC) under the European Union’s Horizon 2020 research and innovation
program under grant agreement No 951541. Y. Meir acknowledges discussions with
A. Georges and support by the Israel Science Foundation (grant 3523/2020).
Experiments at UBC were undertaken with support from the Stewart Blusson
Quantum Matter Institute, the Natural Sciences and Engineering Research
Council of Canada, the Canada Foundation for Innovation, the Canadian
Institute for Advanced Research, and the Canada First Research Excellence
Fund, Quantum Materials and Future Technologies Program. S.F., G.C.G. and M.M.
were supported by the US DOE Office of Basic Energy Sciences, Division of
Materials Sciences and Engineering award DE-SC0006671 and QIS award DE-
SC0020138. AKM acknowledges funding from the Irish Research Council Laureate
Awards 2017/2018 through grant IRCLA/2017/169.

## References

  * Saito  _et al._ (2021) Yu Saito, Fangyuan Yang, Jingyuan Ge, Xiaoxue Liu, Takashi Taniguchi, Kenji Watanabe, J. I. A. Li, Erez Berg, and Andrea F. Young, “Isospin Pomeranchuk effect in twisted bilayer graphene,” Nature 592, 220–224 (2021).
  * Rozen  _et al._ (2021) Asaf Rozen, Jeong Min Park, Uri Zondiner, Yuan Cao, Daniel Rodan-Legrain, Takashi Taniguchi, Kenji Watanabe, Yuval Oreg, Ady Stern, Erez Berg, Pablo Jarillo-Herrero, and Shahal Ilani, “Entropic evidence for a Pomeranchuk effect in magic-angle graphene,” Nature 592, 214–219 (2021), 2009.01836 .
  * Lian (2021) Biao Lian, “Heating freezes electrons in twisted bilayer graphene,” Nature 592, 191–193 (2021).
  * Smirnov (2015) Sergey Smirnov, “Majorana tunneling entropy,” Physical Review B 92, 195312 (2015), 1507.00278 .
  * Sela  _et al._ (2019) Eran Sela, Yuval Oreg, Stephan Plugge, Nikolaus Hartman, Silvia Lüscher, and Joshua Folk, “Detecting the universal fractional entropy of Majorana zero modes,” Physical Review Letters (2019), 10.1103/physrevlett.123.147702, 1905.12237 .
  * Cooper and Stern (2009) NR Cooper and Ady Stern, “Observable bulk signatures of non-abelian quantum hall states,” Physical review letters 102, 176807 (2009).
  * Yang and Halperin (2009) Kun Yang and Bertrand I. Halperin, “Thermopower as a possible probe of non-abelian quasiparticle statistics in fractional quantum hall liquids,” Phys. Rev. B 79, 115317 (2009).
  * Ben-Shach  _et al._ (2013) G. Ben-Shach, C. R. Laumann, I. Neder, A. Yacoby, and B. I. Halperin, “Detecting Non-Abelian Anyons by Charging Spectroscopy,” Physical Review Letters 110, 106805 (2013), 1212.1163 .
  * Hartman  _et al._ (2018) Nikolaus Hartman, Christian Olsen, Silvia Lüscher, Mohammad Samani, Saeed Fallahi, Geoffrey C. Gardner, Michael Manfra, and Joshua Folk, “Direct entropy measurement in a mesoscopic quantum system,” Nature Physics 14, 1083–1086 (2018), 1905.12388 .
  * Kleeorin  _et al._ (2019) Yaakov Kleeorin, Holger Thierschmann, Hartmut Buhmann, Antoine Georges, Laurens W. Molenkamp, and Yigal Meir, “How to measure the entropy of a mesoscopic system via thermoelectric transport,” Nature Communications 10, 5801 (2019), 1904.08948 .
  * Elzerman  _et al._ (2004) J. M. Elzerman, R. Hanson, L. H. Willems van Beveren, B. Witkamp, L. M. K. Vandersypen, and L. P. Kouwenhoven, “Single-shot read-out of an individual electron spin in a quantum dot,” Nature 430, 431–435 (2004), cond-mat/0411232 .
  * Petta  _et al._ (2005) J. R. Petta, A. C. Johnson, J. M. Taylor, E. A. Laird, A. Yacoby, M. D. Lukin, C. M. Marcus, M. P. Hanson, and A. C. Gossard, “Coherent Manipulation of Coupled Electron Spins in Semiconductor Quantum Dots,” Science 309, 2180–2184 (2005).
  * Hanson  _et al._ (2007) R. Hanson, L. P. Kouwenhoven, J. R. Petta, S. Tarucha, and L. M. K. Vandersypen,  _Spins in few-electron quantum dots_ , Ph.D. thesis (2007).
  * Barthel  _et al._ (2009) C. Barthel, D. J. Reilly, C. M. Marcus, M. P. Hanson, and A. C. Gossard, “Rapid Single-Shot Measurement of a Singlet-Triplet Qubit,” Physical Review Letters 103, 160503 (2009), 0902.0227 .
  * Bluhm  _et al._ (2010) Hendrik Bluhm, Sandra Foletti, Diana Mahalu, Vladimir Umansky, and Amir Yacoby, “Enhancing the Coherence of a Spin Qubit by Operating it as a Feedback Loop That Controls its Nuclear Spin Bath,” Physical Review Letters 105, 216803 (2010), 1003.4031 .
  * Nowack  _et al._ (2011) K. C. Nowack, M. Shafiei, M. Laforest, G. E. D. K. Prawiroatmodjo, L. R. Schreiber, C. Reichl, W. Wegscheider, and L. M. K. Vandersypen, “Single-Shot Correlations and Two-Qubit Gate of Solid-State Spins,” Science 333, 1269–1272 (2011).
  * Shulman  _et al._ (2012) M. D. Shulman, O. E. Dial, S. P. Harvey, H. Bluhm, V. Umansky, and A. Yacoby, “Demonstration of Entanglement of Electrostatically Coupled Singlet-Triplet Qubits,” Science 336, 202–205 (2012), 1202.1828 .
  * Pyurbeeva and Mol (2021) Eugenia Pyurbeeva and Jan A Mol, “A Thermodynamic Approach to Measuring Entropy in a Few-Electron Nanodevice,” Entropy 23, 640 (2021), 2008.05747 .
  * Goldhaber-Gordon  _et al._ (1998) D. Goldhaber-Gordon, Hadas Shtrikman, D. Mahalu, David Abusch-Magder, U. Meirav, and M. A. Kastner, “Kondo effect in a single-electron transistor,” Nature 391, 156–159 (1998).
  * Pustilnik and Glazman (2004) Michael Pustilnik and Leonid Glazman, “Kondo effect in quantum dots,” Journal of Physics: Condensed Matter 16, R513 (2004), cond-mat/0401517 .
  * Potok  _et al._ (2007) R. M. Potok, I. G. Rau, Hadas Shtrikman, Yuval Oreg, and D. Goldhaber-Gordon, “Observation of the two-channel Kondo effect,” Nature 446, 167–171 (2007), cond-mat/0610721 .
  * Keller  _et al._ (2015) A. J. Keller, L. Peeters, C. P. Moca, I. Weymann, D. Mahalu, V. Umansky, G. Zaránd, and D. Goldhaber-Gordon, “Universal Fermi liquid crossover and quantum criticality in a mesoscopic system,” Nature 526, 237–240 (2015), 1504.07620 .
  * Dehollain  _et al._ (2020) J. P. Dehollain, U. Mukhopadhyay, V. P. Michal, Y. Wang, B. Wunsch, C. Reichl, W. Wegscheider, M. S. Rudner, E. Demler, and L. M. K. Vandersypen, “Nagaoka ferromagnetism observed in a quantum dot plaquette,” Nature 579, 528–533 (2020), 1904.05680 .
  * Child  _et al._ (2022) Timothy Child, Owen Sheekey, Silvia Lüscher, Saeed Fallahi, Geoffrey C. Gardner, Michael Manfra, and Joshua Folk, “A robust protocol for entropy measurement in mesoscopic circuits,” Entropy 24, 417 (2022).
  * (25) “See supplemental material at [url],” .
  * Tsvelick and Wiegmann (1983) AM Tsvelick and PB Wiegmann, “Exact results in the theory of magnetic alloys,” Advances in Physics 32, 453–713 (1983).
  * Mora  _et al._ (2015) Christophe Mora, Cătălin Paşcu Moca, Jan von Delft, and Gergely Zaránd, “Fermi-liquid theory for the single-impurity anderson model,” Physical Review B 92, 075120 (2015).
  * Legeza  _et al._ (2008) O Legeza, C P Moca, A I Toth, I Weymann, and G Zarand, “Manual for the Flexible DM-NRG code,” arXiv (2008), 0809.3143 .
  * Tóth  _et al._ (2008) A. I. Tóth, C. P. Moca, Ö. Legeza, and G. Zaránd, “Density matrix numerical renormalization group for non-Abelian symmetries,” Physical Review B 78, 245109 (2008), 0802.4332 .
  * Kondo (1964) Jun Kondo, “Resistance Minimum in Dilute Magnetic Alloys,” Progress of Theoretical Physics 32, 37–49 (1964).
  * Silva and Levit (2003) Alessandro Silva and Shimon Levit, “Peculiarities of the controlled dephasing of a quantum dot in the kondo regime,” EPL (Europhysics Letters) 62, 103 (2003).
  * Avinun-Kalish  _et al._ (2004) M Avinun-Kalish, Moty Heiblum, Alessandro Silva, Diana Mahalu, and V Umansky, “Controlled dephasing of a quantum dot in the kondo regime,” Physical review letters 92, 156801 (2004).
  * Kang (2005) Kicheon Kang, “Decoherence of the kondo singlet via a quantum point contact detector,” Phys. Rev. Lett. 95, 206808 (2005).
  * Kang and Khym (2007) Kicheon Kang and Gyong Luck Khym, “Entanglement, measurement, and conditional evolution of the kondo singlet interacting with a mesoscopic detector,” New Journal of Physics 9, 121 (2007).
  * Powell (1964) M. J. D. Powell, “An efficient method for finding the minimum of a function of several variables without calculating derivatives,” The Computer Journal 7, 155–162 (1964).

##  I Device Fabrication

The device was fabricated in a GaAs/AlGaAs heterostructure that hosts a 2D
Electron Gas (2DEG) 57nm below the surface and that had a 300 mK carrier
density of $2.42\times 10^{11}\text{cm}^{-2}$ with mobility $2.56\times
10^{6}\text{cm}^{2}/(\text{Vs})$. A UV laser writer was used to define the
mesas, followed by electron beam lithography to define NiAuGe ohmic contacts.
Additionally, 10 nm of HfO2 was deposited by atomic layer deposition to
improve gating stability. The electrostatic gates were fabricated with two
stages of electron beam lithography followed by electron beam evaporation: a
fine step for the inner parts, and a coarse step for the outer parts of the
gates. In the fine step, 2/12 nm of Pd/Au were deposited. In the coarse step,
10/150 nm of Ti/Au were deposited.

##  II Measurement Electronics

A custom built combined DAC/ADC unit was used to apply potentials to the gates
and heating QPCs as well as to record the voltage output of a Current to
Voltage Basel SP983c amplifier (https://www.baspi.ch/low-noise-high-stab-itov-
conv). The DAC/ADC unit is built from an Arduino Due and two Analog Devices
evaluations boards: the AD5764 DAC and AD7734 ADC. The Arduino is the
interface between the measurement PC and the DAC/ADC boards. The whole design
is based on the information provided at [http://opendacs.com/dac-adc-
homepage/] with some substantial modifications, particularly to the Arduino
code (https://github.com/folk-lab/FastDAC). The most significant modification
is to provide functionality to apply a synchronized square wave bias for
heating, whilst measuring continuously.

##  III Charge sensor mapping

The charge sensor is tuned to the most linear regime before each measurement
(Fig. S1a). In the limit of very strong coupling, however, the transition
becomes so broad that the non-linearity of the charge sensor may begin to play
a role. In our experiment, $\Delta I_{CS}$ is converted to $\Delta N$ assuming
a linear relation between the two, but when $I_{CS}$ is not linear in the
additional electrostatic potential provided by, e.g., cross capacitance with
$V_{D}$, this assumption is no longer valid.

To remove potential inaccuracy in the conversion between $\Delta I_{CS}$ and
$\Delta N$ due to non-linearity, $I_{CS}$ may be mapped back to an equivalent
charge sensor gate voltage using a measurement $I_{CS}(V_{QPC}$ (Fig. S1b).
Performing the entropy calculation using the equivalent gate voltage, instead
of $I_{CS}$, significantly reduces any impact of charge sensor non-linearity
in the measurement of $\Delta N$. In practice, however, no statistically
significant difference was observed in entropy calculations using the two
approaches, so charge sensor mapping was not used.

Figure S1:  a) A very wide sweep over a weakly coupled $0\to 1$ transition
using the plunger gate showing the response of the charge sensor, and the
usual alignment of the $0\to 1$ transition at the steepest, most linear part
of the charge sensor response. b) Data of charge sensor current, $I_{CS}$, vs
the charge sensor QPC tuning gate (see Fig. 1a) used for mapping $I_{CS}$to a
quantity in units of mV that is linear in electrostatic potential.

##  IV Determining $\Delta S$ by direct fitting of $dN/dT$ data

For systems where the addition of an electron to the charge sensed quantum dot
has a simple lineshape, as in the case of a QD weakly coupled to a thermal
reservoir, it is possible to extract the entropy change of the system by
fitting the $\Delta I_{CS}$ data directly as was done in Ref. Hartman  _et
al._ , 2018. This procedure is made possible through the application of the
Maxwell relation:

$\left(\frac{\partial\mu}{\partial T}\right)_{p,N}=-\left(\frac{\partial
S}{\partial N}\right)_{p,T}$

For a QD weakly coupled to a thermal reservoir, the charging lineshape takes
the form:

$N(V_{D},\Theta)=\tanh\left(\frac{V_{D}-V_{mid}(\Theta)}{2\Theta}\right)$

where $\Theta=\frac{k_{B}T}{\alpha e}$ and $V_{mid}(\Theta)$ is the plunger
gate voltage at $N=1/2$. $V_{mid}$ changes the quantum dot energy level–that
is, the energy required to add an electron to the dot–and therefore maps to
the chemical potential $\mu$ in equilibrium.

Differentiating $N$ with respect to $T$, one finds a lineshape that depends
explicitly on $\Delta S$:

$\partial N(V_{D},\Theta)\propto-\partial
T\left[\frac{V_{D}-V_{mid}(\Theta)}{2\Theta}-\frac{\Delta
S}{2k_{B}}\right]\times\cosh^{-2}\left(\frac{V_{D}-V_{mid}(\Theta)}{2\Theta}\right)$

after substituting

$\frac{\partial
V_{mid}}{\partial\Theta}=\frac{1}{k_{B}}\frac{\partial\mu}{\partial
T}=\frac{1}{k_{B}}\Delta S_{N-1\rightarrow N}$

.

Here $\partial N(V_{D},\Theta)$ is the difference in occupation of the QD for
the reservoir temperature changing from $T\rightarrow T+\Delta T$. Fitting
this equation to the $\Delta I_{CS}$ data, $\Delta S$ is obtained as a fit
parameter independent of the scaling of the data.

##  V Fitting NRG to Data

NRG calculations were carried out using the flexible DM-NRG codeLegeza  _et
al._ (2008); Tóth  _et al._ (2008) on the standard single impurity Anderson
model. In the calculation, we assume infinite interaction $U$, a constant
density of states in the reservoir with bandwidth $D=1$ exceeding all other
energy scales, and we keep 350 states per iteration with discretization
constant $\Lambda=2$. Results are given in arbitrary units of energy. Our
procedure for fitting the occupation, $N$, of NRG calculations to measured
data involves three steps: 1. Linearly interpolating over the 2D array of
calculations. 2. Adding terms (amplitude, constant, linear) to account for the
behaviour of the charge sensor in detecting the QD occupations. 3. Allowing
for an offset and scaling proportional to $\Theta$ ($T$ in units of gate
voltage) in the $\epsilon_{0}$ axis. We then use Powell’s method of
minimization Powell (1964) to find the best fitting parameters allowing all to
vary with the exception of $\Gamma$ and $\Theta$, for which only one is
allowed to vary. In the weakly coupled regime, it is reasonable to approximate
$\Gamma\sim 0$, and with that constraint, we are able to determine
$\Theta(V_{D})$. We find a linear relationship between $\Theta(V_{D})$ and
$V_{T}$ which implies a linearly changing lever arm, $\alpha$, as the system
temperature, $T$, is fixed. Note that the lever arm, $\alpha$, connects
$\Theta$ in units of gate voltage to temperature, $T$, in kelvin (Eq.S1),

$\alpha\Theta=k_{B}T$ (S1)

and is a measure of the strength of effect the plunger gate, $V_{D}$, has on
the QD. The linear change implies that as $V_{T}$ is varied, the strength of
effect of $V_{D}$ also varies. We attribute this to a change of shape of the
QD where it moves it further from $V_{D}$ for more positive $V_{T}$. For
measurements into the strongly coupled regime where $\Gamma\gg 0$, we force
the $\Theta$ parameter to follow the linear relationship found in the weakly
coupled regime, allowing $\Gamma$ to be a varying parameter. The fit
parameters found by comparing $N$ of NRG data to $I_{CS}$ of measured data can
then be used to directly compare between the NRG $dN/dT$ calculations and
$\Delta I_{CS}$ measurements.

##  VI Scaling from $\Delta I_{CS}$ to $dN/dT$

The complete procedure for scaling from $\Delta I_{CS}$ to $dN/dT$ is
comprised of two parts: Conversion of $\Delta I_{CS}$ to $\Delta N$, and
calculation of the corresponding $\Delta T$, expressed in equivalent mV on
$V_{D}$.

The procedure for scaling the $\Delta I_{CS}$ measurements to $dN/dT$ involves
scaling $\Delta I_{CS}\to\Delta N$, then dividing by $\Delta T$ as described
in the main text. The $\Delta I_{CS}\to\Delta N$ conversion is a
straightforward division by $I_{e}$, the net change of current through the
charge sensor for the addition of 1 full electron to the QD. Extracting
$I_{e}$ from the data is achieved by scaling NRG calculations of occupation
across the transition to the heated portion of $I_{CS}$ data, with additional
offset and linear terms added to the NRG to account for cross capacitance
between $V_{D}$ and the charge sensor.

$\Delta T$ is easily extracted in units of equivalent gate voltage ($V_{D}$)
for weakly coupled $V_{T}$ by fitting cold and hot occupation data to NRG. For
strongly coupled transitions, however, $\Delta T$ does not result in a
broadening of the transition lineshape, so it must be determined in another
way. The real temperature change of the reservoir does not depend on $V_{T}$,
of course, but the lever arm $\alpha$ does depend on $V_{T}$. We calculate
$\Delta T$($V_{T}$) in equivalent mV on $V_{D}$ by

  1. 1.

fitting hot and cold transitions for a range of weakly coupled $V_{T}$, to
determine both $\Delta T(V_{T})$ in equivalent $V_{D}$ and $\alpha(V_{T})$
through this range.

  2. 2.

$\alpha(V_{T})$ is observed to be linear in $V_{T}$, and extrapolated to
strongly-coupled $V_{T}$ (dashed line in Fig. 1c, main text).

  3. 3.

$\Delta T$ in equivalent $V_{D}$ is calculated for strongly coupled
transitions using $\alpha(V_{T})$ determined above.

##  VII Averaging Procedure

The data shown in the main text is the result of averaging measurements over
many sweeps over the transition. In the strongly coupled regime, as the
$\Delta I_{CS}$ signal becomes weaker, averaging data becomes particularly
important. It is often necessary to measure for 10’s of minutes or even hours
in order to obtain a reasonable signal to noise ratio, but the presence of
charge instability makes single slow measurements over the transition
unreliable. By repeatedly sweeping over the transition quickly, then aligning
each sweep based on a fit to the $I_{CS}$ data before averaging, we can
improve the signal to noise ratio of the corresponding $\Delta I_{CS}$ whilst
mitigating the effect of charge instability. This procedure of post-aligning
individual charge transition scans is followed for measurements in the range
$V_{T}<-230$ mV.

For more strongly coupled measurements ($V_{T}\geq-230$ mV), determining the
center of each individual scan through fitting is not reliable; as a result,
data is averaged without centring first. Charge instability (determined in the
weakly coupled regime) is on the order of 6 $\mu$eV without significant long-
term drift seen in the data. For the most strongly coupled measurements, where
the width of the transition is on the order of 400 $\mu$eV, the lack of
centring therefore is expected to have a negligible effect. Occasional larger
jumps in transition position ( 0.5 mV) do occur on a timescale of hours; care
is taken never to average data across such jumps.

Figure S2:  a) Markers illustrate a single measurement across the charge
transition, which takes 30 seconds to complete. No peak in $\Delta I_{CS}$ can
be seen in this raw data. After averaging 400 of such scans together (solid
line), however, a small peak in $\Delta I_{CS}$ is seen at $\Delta V_{D}=0$.
b) Raw $\Delta I_{CS}$ data (greyscale) for 400 scans as in panel a). Averaged
together, they yield the solid line in panel a).

##  VIII Lack of dependence on charge sensor bias

Figure S3:  The lineshape of $\Delta I_{CS}$, here plotted vs occupation
instead of $V_{D}$, shows no dependence on $V_{CS}$ within experimental noise,
though of course the magnitude of $I_{CS}$ and $\Delta I_{CS}$ scales linearly
with $V_{CS}$. The case of $\Gamma/k_{B}T=24$ is shown here. In particular,
$\Delta I_{CS}$ remains peaked at $N\sim 0.5$, in contrast to the NRG
calculation (solid line) in which the shifted peak reflects the screening of
spin entropy in the mixed valence regime due to the formation of the Kondo
singlet.


# Caging in Motion: Characterizing Robustness in Manipulation through Energy
Margin and Dynamic Caging Analysis

Yifei Dong†, Xianyi Cheng‡ and Florian T. Pokorny† ${\dagger}$: The authors
are with the division of Robotics, Perception and Learning, KTH Royal
Institute of Technology, 10044 Stockholm, Sweden, {yifeid<EMAIL_ADDRESS>${\ddagger}$: The author is with Carnegie Mellon University, Pittsburgh, PA,
15213, USA<EMAIL_ADDRESS>
###### Abstract

To develop robust manipulation policies, quantifying robustness is essential.
Evaluating robustness in general dexterous manipulation, nonetheless, poses
significant challenges due to complex hybrid dynamics, combinatorial explosion
of possible contact interactions, global geometry, etc. This paper introduces
“caging in motion”, an approach for analyzing manipulation robustness through
energy margins and caging-based analysis. Our method assesses manipulation
robustness by measuring the energy margin to failure and extends traditional
caging concepts for a global analysis of dynamic manipulation. This global
analysis is facilitated by a kinodynamic planning framework that naturally
integrates global geometry, contact changes, and robot compliance. We validate
the effectiveness of our approach in the simulation and real-world experiments
of multiple dynamic manipulation scenarios, highlighting its potential to
predict manipulation success and robustness.

##  I Introduction

Human manipulation is intriguing due to its dexterity, simplicity, and
remarkable robustness. Consider picking up a thin object from a table: this
involves complex interactions between the object, table, and hand, yet it
remains robust and fast. The robustness of human manipulation is often
attributed not to precise control over or knowledge of contact interactions
but to factors like hand kinematic configurations, environment geometry, joint
compliance, material softness, and so on. An empirical study in [1] discovers
similar robustness in in-hand manipulation with a soft robotic hand.

Nonetheless, there is a lack of methods in the current literature for
evaluating such general robustness. Existing research studying similar
robustness problems can generally be categorized into three areas: grasp
quality metrics, control stability analysis (including contraction analysis
[2], Lyapunov functions [3] and barrier functions [4]), and caging. Most
quality metrics in grasping [5] and non-prehensile manipulation [6] often
define robustness through object position or wrench space margins based on
first-order analysis. Control stability methods assess robustness by examining
system evolution from initial states, offering a local perspective that
struggles to account for the robustness under the combinatorial explosion of
possible contact interactions, the resulting hybrid dynamics, and global
geometric properties. In contrast, caging approaches, which confine an object
within a bounded space to prevent escape, offer a more global geometric
analysis of robustness.

This paper proposes the caging in motion approach to characterize robustness
in manipulation via energy margins. Our approach expands on prior quasi-static
cage [7] and soft fixture [8] analyses by adapting classical caging concepts
to a kinodynamic context, where dynamically feasible escape paths are
considered. Utilizing forward simulation enables natural consideration of
contact changes, robot joint stiffness, complex geometries, and system
kinodynamics. This methodology introduces a new lens by evaluating the
robustness of manipulation strategies through energy margins and broadens the
applicability of caging concepts to dynamic scenarios.

The main contributions of this paper are: 1\. We propose to characterize more
general robustness in manipulation with energy margins and caging analysis.
This approach evaluates manipulation strategies that consider global geometry,
kinodynamics, and complex contact events. 2\. We develop a sampling-based
kinodynamic planning framework for dynamic caging analysis and energy margin
calculation, designed to be inherently efficient in dynamic, contact-rich
environments. 3\. Through experimental validation in both real-world and
simulated settings, we demonstrate the practicality and effectiveness of our
approach, showcasing its benefits over a wrench-based baseline method.

##  II Related Work

###  II-A Caging

Classical caging has laid the groundwork for robotic nonprehensile
manipulation, focusing on preventing an object from escaping using minimal
constraints. This foundational concept, as introduced by Kuperberg [9] and
further developed by Rimon et al. [10], [11], has been pivotal in exploring
how robotic grasps can effectively cage objects, without necessitating force
closure or form closure. Rodriguez et al. [12] demonstrate utilizing caging to
achieve immobilizing grasps. Energy-bounded caging [13] relaxes the complete
caging in the presence of external forces such as gravity. Planar energy-
bounded cages [14] are synthesized by identifying configurations that maximize
the energy required for the object to escape. [15] proposes whole-hand
manipulation via caging by considering compliant underactuated hands as
obstacles and introducing an object energy field. Our previous work on soft
fixture [8] presents practical sampling-based methods to search for upper-
bound estimates of escape energy applicable in higher-dimensional
configuration spaces, such as simplified deformable objects. While classical
caging has offered foundational insights into object restraint, the “cage in
motion” approach we propose here represents a new perspective. It extends
traditional caging into dynamic settings, bridging the gap between static
caging theories and the dynamic realities of robotic manipulation. On the
other hand, our work first brings the caging analysis to characterize
manipulation robustness in motion.

###  II-B Robustness and Quality Metrics

Parallel to the evolution of caging, grasp quality measures that quantify the
goodness of a grasp have been a significant line of research. A comprehensive
review by Roa et al. [5] categorizes quality measures based on contact points
and hand configurations, providing a systematic framework for evaluating grasp
quality. Grasp Wrench Space (GWS) [16] has been a useful tool to determine the
quality of a grasp. Several metrics derived from the GWS can quantify the
“goodness” of a grasp, such as the volume of the GWS [17], the largest
perturbation wrench the grasp can resist in any direction [18], etc. More
recently, Lin et al. [19] introduce a task-specific grasp quality criterion
based on the distribution of task disturbance, exemplifying the importance of
task-specific optimization of grasps. Xu et al. [20] propose minimal work as a
grasp quality metric to plan grasps for 3D deformable hollow objects
considering wrench resistance. Machine learning approaches are also employed
to predict grasp success based on data labels derived from analytic quality
metrics [21]. This line of research thoroughly investigates wrench-based
quality metrics for prehensile grasping, especially force closure grasps. For
non-prehensile manipulation, a framework of quasi-static wrench-based analysis
is presented for balancing grasps in [22]. Makapunyo et al. [23] employ
sampling-based methods to determine a quality metric for partial caging. In
contrast to the previous literature, our work aims to provide a new paradigm
for robustness characterization that is independent of wrench information. We
focus instead on analyzing energy margin and dynamic caging within kinodynamic
settings. Thereby, we aim to encapsulate a broader spectrum of manipulation
scenarios, offering a more global and adaptable framework for assessing
manipulation robustness and success in dynamically evolving environments.

##  III Preliminaries

###  III-A Nomenclature

We denote by $\mathcal{X}\subset\mathrm{SE}(3)$ the configuration space
(C-space) of a 3D rigid object. The set of time derivatives of the object
configuration is denoted by $\dot{\mathcal{X}}\subset T\mathrm{SE}(3)$, where
$T\mathrm{SE}(3)$ is the tangent space associated with the Lie group
$\mathrm{SE}(3)$. The C-space of a robot end-effector and surrounding static
obstacles is denoted by
$\mathcal{Y}\subset\mathrm{SE}(3)\times\mathbb{R}^{n_{r}}$, where
$\mathrm{SE}(3)$ denotes the base pose, and $n_{r}$ refers to the number of
joints. Its set of time derivatives is denoted by $\dot{\mathcal{Y}}\subset
T\mathrm{SE}(3)\times\mathbb{R}^{n_{r}}$. $\mathcal{Y}$ might degenerate when,
for example, the base of the end-effector is fixed, and then
$\mathcal{Y}\subset\mathbb{R}^{n_{r}}$. For planar manipulation, $\mathcal{X}$
and $\mathcal{Y}$ degenerate to $\mathrm{SE}(2)$ and
$\mathrm{SE}(2)\times\mathbb{R}^{n_{r}}$, respectively. An element
$\boldsymbol{x}\in\mathcal{X}$, described as
$\boldsymbol{x}=(\boldsymbol{r}_{x},\boldsymbol{q}_{x})$, encompasses the
position $\boldsymbol{r}_{x}$ of the Center of Mass (CoM) and orientation (a
unit quaternion) $\boldsymbol{q}_{x}$. An element
$\dot{\boldsymbol{x}}\in\dot{\mathcal{X}}$ is described as
$\dot{\boldsymbol{x}}=(\dot{\boldsymbol{r}}_{x},\dot{\boldsymbol{q}}_{x})$,
which comprises the CoM linear velocity $\dot{\boldsymbol{r}}_{x}$ and the
derivative of quaternion $\dot{\boldsymbol{q}}_{x}$. Similarly for
${\boldsymbol{y}}\in\mathcal{Y}$ and
$\dot{\boldsymbol{y}}\in\dot{\mathcal{Y}}$, we have
$\boldsymbol{y}=(\boldsymbol{r}_{y},\boldsymbol{q}_{y},\boldsymbol{\alpha}_{y})$
and
$\dot{\boldsymbol{y}}=(\dot{\boldsymbol{r}}_{y},\dot{\boldsymbol{q}}_{y},\dot{\boldsymbol{\alpha}}_{y})$,
where $\boldsymbol{\alpha}_{y}$ and $\dot{\boldsymbol{\alpha}}_{y}$ refer to
the joint position and velocity of the end-effector, respectively.

Figure 1: Left: quasi-static soft fixture concepts of an escape path $\sigma$
that leads an object (orange) at $\boldsymbol{x}$ inside a bowl (green,
regarded as a static obstacle $\mathcal{X}_{\text{obs}}$) to a goal region
$\mathcal{X}_{\text{goal}}$. A capture set (blue, $\mathcal{X}_{\text{cap}}$)
is constructed by a bounded potential energy sublevel set with a minimum
energy threshold value. Right: “caging in motion” — an end-effector (with pose
$\boldsymbol{y}$, linear and angular velocity $\dot{\boldsymbol{r}}_{y}$,
$\boldsymbol{\omega}_{y}$) is employed to guide an object (orange) towards a
target area (top left corner) on a table. A capture set denotes the region
anticipated to be swept by the end-effector in the forthcoming time interval.

###  III-B Cage and Soft Fixture

Here we revisit the fundamental concept of caging. The free C-space
$\mathcal{X}_{\text{free}}$ indicates a set of configurations for which the
object does not penetrate any of the bodies (the end-effector or obstacles) in
the workspace.

######  Definition 1:

A cage occurs when the configuration of an object
$\boldsymbol{x}_{\text{init}}\in\mathcal{X}_{\text{free}}$ is situated in a
bounded path component of $\mathcal{X}_{\text{free}}$.

An object in a cage thus indicates its limited capacity to freely move beyond
a certain proximity to its starting configuration. Caging was relaxed as
partial caging by allowing escape paths through narrow passages in the free
space [23]. Partial caging further extends to energy-bounded caging [13] and
soft fixtures [8] that consider both geometric and potential energy
constraints, as introduced below.

######  Definition 2:

In scenarios where only conservative forces act within a quasi-static system,
an object at an initial configuration
$\boldsymbol{x}_{\text{init}}\in\mathcal{X}_{\text{free}}$ possesses potential
energy $E(\boldsymbol{x}_{\text{init}})$. Its _escape energy_
$\mathcal{E}(\boldsymbol{x}_{\text{init}})$ refers to the supremum of $e$ such
that the path component
$\mathcal{PC}_{\boldsymbol{x}_{\text{init}}}(\mathcal{X}_{e}(\boldsymbol{x}_{\text{init}}))$
of a sublevel set $\mathcal{X}_{e}(\boldsymbol{x}_{\text{init}})$ containing
$\boldsymbol{x}_{\text{init}}$ is bounded:

$\displaystyle\mathcal{E}(\boldsymbol{x}_{\text{init}})$
$\displaystyle=\sup\left\\{e\geq
0:\mathcal{PC}_{\boldsymbol{x}_{\text{init}}}(\mathcal{X}_{e}(\boldsymbol{x}_{\text{init}}))\text{
is bounded}\right\\},$ (1)

where

$\mathcal{X}_{e}(\boldsymbol{x}_{\text{init}})=\\{\boldsymbol{x}\in\mathcal{X}_{\text{free}}:E(\boldsymbol{x})\leq
E(\boldsymbol{x}_{\text{init}})+e\\}.$ (2)

The object is deemed to be in a _soft fixture_ if
$\mathcal{E}(\boldsymbol{x}_{\text{init}})>0$. Note that
$\mathcal{E}(\boldsymbol{x}_{\text{init}})$ is defined only if
$\mathcal{PC}_{\boldsymbol{x}_{\text{init}}}(\mathcal{X}_{0}(\boldsymbol{x}_{\text{init}}))$
is bounded.

In practice, sampling-based approaches can approximate a close upper bound of
$\mathcal{E}(\boldsymbol{x}_{\text{init}})$ by utilizing a goal region (Fig.
1, left) sufficiently far away from the obstacles with relatively lower energy
values and searching for geometrically feasible escape paths with low energy
cost leading to it [8].

###  III-C Assumptions

We introduce the assumptions that this paper relaxes compared to our prior
work [8].

1) System Motion: While the previous work only considers quasi-static
conditions, here we account for the movement of both the object and the end-
effector, while keeping other environmental elements static. We introduce
$\mathcal{Z}=\mathcal{X}\times\dot{\mathcal{X}}\times\mathcal{Y}\times\dot{\mathcal{Y}}$,
where $\mathcal{Z}_{\text{free}}$ represents penetration-free space. In this
context, the energy function $E$ is the system’s total mechanical energy (sum
of kinetic and potential energy).

2) Non-Conservative Forces: Given the existence of non-conservative forces,
such as friction, the previous definition of escape energy in Def. 1 is
problematic. Consider the example of wiping a table (Fig. 1, right), assuming
the tabletop with a rough surface is infinitely large. The escape energy
$\mathcal{E}(\boldsymbol{x}_{\text{init}})$ is always infinite in this
scenario because the object overcomes friction along any escape path to
arbitrarily far away. To utilize escape energy in this context, we define a
capture set
$\mathcal{Z}_{\text{cap}}(\boldsymbol{z}_{\text{init}})\subset\mathcal{Z}$
instead, considering the object as dynamically controllable by the end-
effector for achieving specific tasks. The effort of the object to escape from
$\mathcal{Z}_{\text{cap}}(\boldsymbol{z}_{\text{init}})$ thus indicates
manipulation robustness, implying an energy margin from failure.

3) Dynamic Escape Paths: We consider dynamically feasible escape paths that
lead the object out of the capture set
$\mathcal{Z}_{\text{cap}}(\boldsymbol{z}_{\text{init}})$. A dynamically
feasible path $\sigma:[0,1]\to\mathcal{Z}_{\text{free}}$, such that
$\sigma(0)=\boldsymbol{z}_{\text{init}}$ and
$\sigma(1)=\boldsymbol{z}_{\text{goal}}$, is generated by applying random
wrenches to the object to mimic perturbation or uncertainty during the
manipulation. Here, $\boldsymbol{z}_{\text{goal}}\in\mathcal{Z}_{\text{goal}}$
and $\mathcal{Z}_{\text{goal}}$ is a goal set of escape.

Besides the assumptions above, we only consider rigid objects under potential
energy fields, though our framework can potentially be extended to scenarios
with articulated or multiple objects.

###  III-D Problem Statement

In this work, we tackle the following problem: Given a specific manipulation
task objective and the state $\boldsymbol{z}_{\text{init}}$ of a rigid object
and an end-effector in motion, which is subject to the assumptions in Section
III-C, 1. Define metrics to quantify energy margins to failure for
characterizing robustness in dexterous manipulation (Section IV). 2. Compute
the energy margin metrics using kinodynamic motion planning and evaluate their
performance (Section V and VI).

##  IV Energy Margin for Caging in Motion

This section quantifies manipulation robustness with two methods based on
energy margin: the effort of escape and capture score. The former focuses on
the ease of failure, and the latter on the probability of staying safe.

###  IV-A Effort of Escape

We first measure the energy margin to failure by the minimum effort to escape
out of the capture set.

######  Definition 3:

The minimal _effort of escape_ is the minimal integral of the absolute power
$\lvert\dot{W}_{\text{ext}}(\sigma(t))\rvert$ of the external work
$W_{\text{ext}}(\sigma(t))$ required to be done to help the object escape from
the capture set $\mathcal{Z}_{\text{cap}}(\boldsymbol{z}_{\text{init}})$:

$\Omega_{\text{esc}}\left(\boldsymbol{z}_{\text{init}}\right)=\min_{\sigma\in\Sigma(\boldsymbol{z}_{\text{init}},\mathcal{Z}_{\text{goal}})}\int_{0}^{1}\lvert\dot{W}_{\text{ext}}(\sigma(t))\rvert
dt,$ (3)

where

$\dot{W}_{\text{ext}}(\sigma(t))=\frac{dE(\sigma(t))}{dt}-\dot{W}_{\text{fri}}(\sigma(t)).$
(4)

$\frac{dE(\sigma(t))}{dt}$ is the instantaneous rate of change of the
mechanical energy of the system, and $W_{\text{fri}}(\sigma(t))$ represents
the work done by friction at the instantaneous state and velocity of the
system. We denote
$\Sigma(\boldsymbol{z}_{\text{init}},\mathcal{Z}_{\text{goal}})$ as the set of
all escape paths starting from $\boldsymbol{z}_{\text{init}}$ terminating in
$\mathcal{Z}_{\text{goal}}$, where
$\mathcal{Z}_{\text{goal}}=\mathcal{Z}_{\text{cap}}^{c}(\boldsymbol{z}_{\text{init}})$
(the complement set of the capture set).

Eq. (4) follows the first law of thermodynamics and the thermal equilibrium
assumption. We employ the absolute value of the power
$\lvert\dot{W}_{\text{ext}}(\sigma(t))\rvert$ as a measure of the total
control input that we want to minimize because, for example, both speeding up
and slowing down a vehicle require fuel, and we might want to minimize the
total fuel used over a journey. Since paths in the path space
$\Sigma(\boldsymbol{z}_{\text{init}},\mathcal{Z}_{\text{goal}})$ can hardly be
enumerated or optimized analytically given the complexity introduced by
friction, we employ kinodynamic sampling-based planners to filter candidate
paths that upper bounds
$\Omega_{\text{esc}}\left(\boldsymbol{z}_{\text{init}}\right)$.

###  IV-B Capture Score

While the effort of escape assesses the minimal effort required to exit the
capture set, the capture score estimates the likelihood of remaining within
the capture set.

We first consider an energy cost field and its correlating probabilistic
distribution of the system. Specifically, a sequence of random wrenches
applied on the object at $\boldsymbol{z}_{\text{init}}$ lead it to
$\boldsymbol{z}^{0}$ with a cost of total effort $c^{0}$. By repeating the
process $M$ times from $\boldsymbol{z}_{\text{init}}$, the system terminates
in a list of states $\\{\boldsymbol{z}^{0},...,\boldsymbol{z}^{M}\\}$ with
corresponding costs $\\{c^{0},...,c^{M}\\}$. Data pairs in
$\\{(\boldsymbol{z}^{0},c^{0}),...,(\boldsymbol{z}^{M},c^{M})\\}$ constitute
an energy cost field that demonstrates the state space reachability in terms
of energy cost from $\boldsymbol{z}_{\text{init}}$. We can thus build a
probability mass function $L:\mathcal{Z}\to\mathbb{R}^{+}$ from the energy
cost field, approximating the likelihood of reaching a state
$\boldsymbol{z}^{m}\in\mathcal{Z}$ with a softmax function,

$L(\boldsymbol{z}^{m})=\frac{e^{-\lambda(c^{m}-c_{\text{min}})}}{\sum_{m=0}^{M}e^{-\lambda(c^{m}-c_{\text{min}})}},$
(5)

where $c_{\text{min}}=\min_{0\leq m\leq M}{c^{m}}$ and $\lambda$ is a hyper-
parameter. The approximated likelihood function indicates some states with
lower costs are probabilistically more reachable.

######  Definition 4:

The capture score is defined as the sum of likelihood values of samples in a
capture set $\mathcal{Z}_{\text{cap}}(\boldsymbol{z}_{\text{init}})$,

$\Omega_{\text{cap}}(\boldsymbol{z}_{\text{init}})=\sum_{m=0}^{M}\delta(\boldsymbol{z}^{m}\in\mathcal{Z}_{\text{cap}}(\boldsymbol{z}_{\text{init}}))\cdot
L(\boldsymbol{z}^{m}),$ (6)

where $\delta(\cdot)$ is an indicator function that equals $1$ if the
condition inside the brackets is satisfied and $0$ otherwise.

Similarly, we define a task success set
$\mathcal{Z}_{\text{suc}}\subset\mathcal{Z}$ independent of
$\boldsymbol{z}_{\text{init}}$, indicating the set of states symbolizing the
object accomplishes the manipulation objective. A capture score of success is
thus given by

$\Omega_{\text{suc}}(\boldsymbol{z}_{\text{init}})=\sum_{m=0}^{M}\delta(\boldsymbol{z}^{m}\in\mathcal{Z}_{\text{suc}})\cdot
L(\boldsymbol{z}^{m}),$ (7)

which is essentially a predictor for fulfilling a task-specific objective from
the state $\boldsymbol{z}_{\text{init}}$. In practice, we employ kinodynamic
motion planners rather than repetitive Monte Carlo rollouts from
$\boldsymbol{z}_{\text{init}}$. Sampling-based kinodynamic motion planners are
in general more time- and memory-efficient by using strategies such as biased
sampling and caching explored nodes.

##  V Computing Energy Margin Through Kinodynamic Planning

We develop kinodynamic motion planning algorithms to compute the energy
margins for characterizing robustness as defined in Section IV. We calculate
the effort of escape $\Omega_{\text{esc}}(\boldsymbol{z}_{\text{init}})$ by an
iterative tree search algorithm with an optimality guarantee (Section V-B).
The capture scores $\Omega_{\text{cap}}(\boldsymbol{z}_{\text{init}})$ and
$\Omega_{\text{suc}}(\boldsymbol{z}_{\text{init}})$ are obtained by growing an
expansive tree and approximating its correlating probabilistic distribution
(Section V-C).

###  V-A Planning Objectives

We consider a system with an object, an end-effector, and the environment. The
system state at time $k$ is

$\boldsymbol{z}(k)=\left[\begin{array}[]{llll}\boldsymbol{x}^{\top}(k)&\dot{\boldsymbol{x}}^{\top}(k)&\boldsymbol{y}^{\top}(k)&\dot{\boldsymbol{y}}^{\top}(k)\end{array}\right]^{\top}.$
(8)

A manipulation trajectory is denoted by
$\\{(k,\boldsymbol{z}(k)):k\in[0,K]\\}$, complying with the system dynamics
${\boldsymbol{z}}^{\prime}(k)=F\left(\boldsymbol{z}(k),\boldsymbol{\phi}(k)\right)$.
$\boldsymbol{\phi}$ denotes the robot control. The system dynamics is rolled
out through physics simulation. For each state $\boldsymbol{z}(k)$ along the
trajectory, Algo. 1 computes the minimal effort of escape
$\Omega_{\text{esc}}(\boldsymbol{z}(k))$, and Algo. 2 computes the capture
scores $\Omega_{\text{cap}}(\boldsymbol{z}(k))$ and
$\Omega_{\text{suc}}(\boldsymbol{z}(k))$. We assume that the end-effector is
not actuated in the short period when computing the metrics. Next, we detail
the two algorithms.

1 Initialize, $P$, $A$, $n$

2 $\sigma_{0}\leftarrow A(P_{\infty})$ $\triangleright$ Find a first escape
path

3 if  _$\sigma_{0}$ does not exists_ then

4 $P$ has no solution

5 Return $\infty$ $\triangleright$ Infinite cost to escape

6

7$c_{0}\leftarrow C(\sigma_{0})$ $\triangleright$ eq. (9). Initial cost

8 for  _$i=1,2,\ldots,n$_ do

9 $\sigma_{i}\leftarrow A(P_{c_{i-1}})$ $\triangleright$ Find a path of lower
cost

10 $c_{i}\leftarrow C(\sigma_{i})$ $\triangleright$ eq. (9). Update cost upper
bound

11

12$\Omega_{\text{esc}}(\boldsymbol{z}_{\text{init}})\leftarrow c_{n}$
$\triangleright$ eq. (3). Final cost

Return $\Omega_{\text{esc}}(\boldsymbol{z}_{\text{init}})$

Algorithm 1 Compute Effort of Escape

###  V-B Compute Effort of Escape

We frame finding the minimum effort to escape from the capture set
$\mathcal{Z}_{\text{cap}}(\boldsymbol{z}_{\text{init}})$ as an optimal
kinodynamic motion planning problem $P$ of which the goal set is the
complement set of the capture set
($\mathcal{Z}_{\text{goal}}=\mathcal{Z}_{\text{cap}}^{c}(\boldsymbol{z}_{\text{init}})$).
We employ the AO-$x$ meta algorithms [24], which are asymptotically optimal
motion planners that decompose $P$ into feasible motion planning problems in
the state-cost space. They iteratively find an upper bound approaching the
minimum effort of escape by running as subroutines feasible kinodynamic
planner $x$, such as Expansive Space Tree (EST) [25] or Rapidly-exploring
Random Tree (RRT) [26] algorithms.

The optimal motion planning problem
$P=(Q,\mathcal{Z}_{\text{free}},\mathcal{U},\boldsymbol{z}_{\text{init}},\mathcal{Z}_{\text{goal}},\mathcal{Z}_{b},\mathcal{U}_{b},G)$
produces a trajectory $\sigma(t):[0,1]\to\mathcal{Z}_{\text{free}}$ and
control $\boldsymbol{u}(t):[0,1]\to\mathcal{U}$ that minimizes the objective
functional:

$C(\sigma)=\int_{0}^{1}Q(\sigma(t),\boldsymbol{u}(t))dt=\int_{0}^{1}\lvert\dot{W}_{\text{ext}}(\sigma(t))\rvert
dt.$ (9)

Here, $Q$ denotes the incremental cost (terminal cost is $0$) following Eq.
(3). $\mathcal{U}_{b}$ is the set of control constraints and
$\boldsymbol{u}\in\mathcal{U}_{b}\subset\mathcal{U}$. $\mathcal{Z}_{b}$ refers
to the set of kinematic constraints and
$\boldsymbol{z}\in\mathcal{Z}_{b}\subset\mathcal{Z}_{\text{free}}$. The
dynamics is subjected to
$\boldsymbol{z}^{\prime}=G\left(\boldsymbol{z},\boldsymbol{u}\right)$. The
control input
$\boldsymbol{u}=\left[\begin{array}[]{llll}\boldsymbol{f}^{\top}&\boldsymbol{\tau}^{\top}\end{array}\right]^{\top}$
refers to a wrench (force $\boldsymbol{f}$ applied at its CoM and torque
$\boldsymbol{\tau}$) applied on the object, mimicking an external
perturbation.

Given the optimal motion planning problem $P$, a feasible motion planning
problem
$P_{\bar{c}}=(\mathcal{Z}_{\text{free}}\times\mathbb{R}^{+},\mathcal{U},(\boldsymbol{z}_{\text{init}},0),\hat{\mathcal{Z}}_{\text{goal}}^{\bar{c}},\mathcal{Z}_{b}\times\mathbb{R}^{+},\mathcal{U}_{b},\hat{G})$
is solved at each iteration. $P_{\bar{c}}$ augments the state in $P$ by an
auxiliary cost variable, which measures the accumulated cost from the root
$\boldsymbol{z}_{\text{init}}$, i.e. cost-to-come. Here,
$\hat{\mathcal{Z}}_{\text{goal}}^{\bar{c}}=\\{(\boldsymbol{z},c):\boldsymbol{z}\in\mathcal{Z}_{\text{goal}},c\in[0,\bar{c}]\\}$
denotes the augmented goal set with the range of cost space upper-bounded by
$\bar{c}$, and $\hat{G}$ refers to the augmented dynamics given by

$\hat{\boldsymbol{z}}^{\prime}=\left[\begin{array}[]{l}\boldsymbol{z}^{\prime}\\\
{c}^{\prime}\end{array}\right]=\left[\begin{array}[]{l}G(\boldsymbol{z},\boldsymbol{u})\\\
Q(\boldsymbol{z},\boldsymbol{u})\end{array}\right].$ (10)

We solve $P$ by converting it to a series of feasible motion planning
subproblems $P_{\bar{c}}$ and iteratively lowering the cost space upper bound
$\bar{c}$ utilizing a well-behaved, probabilistically complete feasible motion
planner $A$ (EST or RRT). A feasible planner is called well-behaved if there
is a nonzero chance that the planner does not produce the worst-possible path
in finite time (if one exists).

Specifically, we first run the subroutine algorithm $A(P_{\infty})$ with
infinite cost space upper bound $\bar{c}=\infty$ to find the first escape path
$\sigma_{0}$ by growing a tree rooted at $\boldsymbol{z}_{\text{init}}$ (Algo.
1, Line 2). The tree expands a new node $\boldsymbol{z}^{\prime}$ by enforcing
a randomized control $\boldsymbol{u}$ based on algorithm-specific strategies
from an existing node $\boldsymbol{z}$, complying with the system dynamics
$G$. The forward integration is performed in a physical simulator that allows
rebounding after a collision. A new extension is valid if
$\boldsymbol{z}^{\prime}$ abides by the augmented kinematic constraints
$(\boldsymbol{z}^{\prime},c(\boldsymbol{z}^{\prime}))\in\mathcal{Z}_{b}\times\mathbb{R}^{+}$.
If a path does not exist, it indicates the effort of escape is infinite
$\Omega_{\text{esc}}\left(\boldsymbol{z}_{\text{init}}\right)=\infty$ (Line
3-5). Once we found a path $\sigma_{0}$ with cost $c_{0}=C(\sigma_{0})$, we
solve a second subproblem $P_{c_{0}}$ (Line 8) that shrinks the cost space by
lowering the upper bound from $\infty$ to $c_{0}$. The tree in the previous
subproblem is cached and optionally pruned (discarding the nodes with costs
higher than $c_{0}$) such that we do not grow a new tree from scratch.
$A(P_{c_{0}})$ results in a new escape path $\sigma_{1}$ with cost $c_{1}$
(Line 9), and so on and so forth. Given the well-behavedness of $A$, the cost
$c_{i}$ approaches the optimal cost $c^{*}$ as the total number of iterations
approaches infinity. Therefore, we can approximately upper-bound the minimum
effort of escape
$\Omega_{\text{esc}}\left(\boldsymbol{z}_{\text{init}}\right)$ in finite ($n$)
loops.

In practice, we consider EST or RRT as the subroutine algorithms $A$. RRT
attempts to randomly sample a state in the state space and select an extension
from its nearest neighbor node already in the tree towards it. EST takes the
strategy of prioritizing extensions in the area with less density of existing
nodes in the tree.

1

2 Initialize, $P_{\infty}$,
$\mathcal{Z}_{\text{cap}}(\boldsymbol{z}_{\text{init}})$,
$\mathcal{Z}_{\text{suc}}$, $M$

3 $\triangleright$ Obtain energy cost field by growing a tree

4
$\\{(\boldsymbol{z}^{0},c^{0}),...,(\boldsymbol{z}^{M},c^{M})\\}\leftarrow\textsc{EST}(P_{\infty},M)$

5

6$\triangleright$ Approximate probabilistic distribution

7 $c_{\text{min}}=\min_{0\leq m\leq M}{c^{m}}$

8 $S\leftarrow{\sum_{m=0}^{M}e^{-\lambda(c^{m}-c_{\text{min}})}}$

9 for  _$m=1,2,\ldots,M$_ do

10
$L(\boldsymbol{z}^{m})\leftarrow\frac{1}{S}{e^{-\lambda(c^{m}-c_{\text{min}})}}$
$\triangleright$ eq. (5).

11

12$\triangleright$ Compute scores

13
$\Omega_{\text{cap}}(\boldsymbol{z}_{\text{init}})\leftarrow\sum_{m=0}^{M}\delta(\boldsymbol{z}^{m}\in\mathcal{Z}_{\text{cap}}(\boldsymbol{z}_{\text{init}}))\cdot
L(\boldsymbol{z}^{m})$ $\triangleright$ eq. (6).

14
$\Omega_{\text{suc}}(\boldsymbol{z}_{\text{init}})\leftarrow\sum_{m=0}^{M}\delta(\boldsymbol{z}^{m}\in\mathcal{Z}_{\text{suc}})\cdot
L(\boldsymbol{z}^{m})$ $\triangleright$ eq. (7).

Return $\Omega_{\text{cap}}(\boldsymbol{z}_{\text{init}})$,
$\Omega_{\text{suc}}(\boldsymbol{z}_{\text{init}})$

Algorithm 2 Compute Capture Scores

###  V-C Compute Capture Scores

To approximate the capture scores, we construct an energy cost field through
kinodynamic tree expansion. Achieving an evenly distributed sampling of state
space around the initial state $\boldsymbol{z}_{\text{init}}$ is critical for
this purpose. EST is thereby selected for its efficacy in promoting uniform
distribution across the state space through inverse density weighting, unlike
RRT, which may lead to uneven exploration. As detailed in Algo. 2, we consider
the cost-augmented unbounded feasible motion planning problem $P_{\infty}$
(Line 2-3) with the same settings as in Section V-B except an infinitely far
away goal set
$\mathcal{Z}_{\text{goal}}=\mathcal{Z}_{\infty}\subset\mathcal{Z}$ with
unconstrained kinematics $\mathcal{Z}_{b}=\mathcal{Z}_{\text{free}}$. EST
terminates after growing $M$ nodes in the tree. We thus obtain an energy cost
field of $M$ state-cost pairs,
$\\{(\boldsymbol{z}^{0},c^{0}),...,(\boldsymbol{z}^{M},c^{M})\\}$. We
thereafter approximate the probabilistic distribution (Line 4-8) and compute
the scores $\Omega_{\text{cap}}(\boldsymbol{z}_{\text{init}})$ and
$\Omega_{\text{suc}}(\boldsymbol{z}_{\text{init}})$ (Line 9-11) following the
procedure in Section IV-B.

##  VI Evaluation

In this section, we aim to validate the practicality and effectiveness of our
metrics and algorithms by conducting thorough simulation and real-world
experiments. Our primary assumption is that energy margins and caging analysis
can effectively characterize robustness in manipulation. To verify this, we
designed four simulated manipulation tasks and used our metrics to predict
their robustness and task success. More specifically, we generated multiple
trajectories of the tasks and recorded their system kinematic information and
ground-truth labels of robustness and success. We computed the metrics
offline, used them to make predictions and compared the results with the
ground-truth labels. Through the experiments, we validated the reliable
prediction capability of our metrics in environments with complex contact
events. Furthermore, we conducted ablation studies to verify the efficiency of
our kinodynamic planning algorithms and the robustness against different model
parameter errors.

Figure 2: Four examples of dynamic manipulation primitives under gravity — (i)
planar pushing, (ii) balanced transport, (iii) toppling and (iv) grasping from
a table, where the task objectives are all eventually achieved. Objects, end-
effectors and obstacles are shown in orange, green and black, respectively.

###  VI-A Task Description

Inspired by the manipulation primitive taxonomy in [27], [28], we provide
examples covering planar pushing, balanced transport, toppling, and grasping
from a table.

1) Planar Pushing: we consider the planar manipulation problem of pushing an
object on a horizontal plane towards a wall (Fig. 2-i). The task is considered
successful if a state $\boldsymbol{z}\in\mathcal{Z}_{\text{suc}}$ lies in the
task success set
$\mathcal{Z}_{\text{suc}}=\\{\boldsymbol{z}\in\mathcal{Z}:\boldsymbol{r}_{x}\in\mathcal{R}_{\text{suc}}\\}$,
and $\mathcal{R}_{\text{suc}}\subset\mathbb{R}^{2}$ refers to the set of
object positions where the longest edge of the convex hull of the object
aligns with the wall. The capture set
$\mathcal{Z}_{\text{cap}}(\boldsymbol{z}_{\text{init}})=\\{\boldsymbol{z}\in\mathcal{Z}:\boldsymbol{r}_{x}\in\mathcal{R}_{\text{cap}}(\boldsymbol{z}_{\text{init}})\\}$
includes state $\boldsymbol{z}\in\mathcal{Z}$ such that the CoM position of
the object $\boldsymbol{r}_{x}$ lies inside a subset of the workspace
$\mathcal{R}_{\text{cap}}(\boldsymbol{z}_{\text{init}})\subset\mathbb{R}^{2}$.
$\mathcal{R}_{\text{cap}}(\boldsymbol{z}_{\text{init}})$ indicates the
circular sector region in the workspace that the end-effector will sweep over
given the current instantaneous rotation center of the end-effector (Fig. 1,
right).

2) Balanced Transport: We aim to balance a cube on a steep slope with a
rectangular support surface so that it avoids the failure mode of falling off
and is transported to a target region (Fig. 2-ii). Similar to the selection of
the capture set $\mathcal{Z}_{\text{cap}}(\boldsymbol{z}_{\text{init}})$ and
the task success set $\mathcal{Z}_{\text{suc}}$ above, we add constraints on
the object CoM position $\boldsymbol{r}_{x}$ such that the capture condition
(maintaining on the support surface,
$\mathcal{R}_{\text{cap}}(\boldsymbol{z}_{\text{init}})$) and the task success
condition (reaching a goal region $\mathcal{R}_{\text{suc}}$) are satisfied.

3) Toppling: We consider the 2D task of toppling a box on the tabletop using a
spring-like fingertip manipulator (Fig. 2-iii). The task is accomplished when
the box rotates by $\pi/2$ rad pivoting on one of the corners in contact with
the table, i.e.
$\mathcal{Z}_{\text{suc}}=\\{\boldsymbol{z}\in\mathcal{Z}:p(\boldsymbol{q}_{x})=\pi/2\\}$,
where $p(\boldsymbol{q}_{x})$ is the $\mathrm{SO}(2)$ orientation of the box
in Euler angle. It often fails when sliding occurs on the table contact. The
capture set is defined as
$\mathcal{Z}_{\text{cap}}(\boldsymbol{z}_{\text{init}})=\\{\boldsymbol{z}\in\mathcal{Z}:\|\dot{\boldsymbol{r}}_{x,\text{pivot}}\|\leq
0.1\\}$, i.e. the magnitude of the sliding velocity of the table contact is
smaller than 0.1 m/s.

4) Grasping from a Table: The task of grasping a box lying on the desktop is
considered here, using a Robotiq three-finger adaptive robot gripper with
$n_{r}=12$ revolute joints in the fingers (Fig. 2-iv). The task is successful
when the box of edge length 10 cm is lifted above the desktop by 10 cm, while
it fails if it slips and falls from the gripper. The position-controlled
compliant gripper poses an external energy bound for the box to escape.
Constraints on the object CoM position $\boldsymbol{r}_{x}$ are considered
similarly as in the first two examples in selecting the sets
$\mathcal{Z}_{\text{cap}}(\boldsymbol{z}_{\text{init}})$ and
$\mathcal{Z}_{\text{suc}}$.

###  VI-B Baseline

We compare our approach with a manipulation quality score from [6], which is
designed to measure the robustness of contact modes (sticking, sliding,
disengaged, etc) against disturbance forces (e.g. avoid slipping for sticking
contacts). This characterizes manipulation robustness against failure modes
similar to the goal in this paper. However, our energy margin approach
quantifies manipulation robustness more globally and captures contact mode
transitions and loss of contact to a certain degree. To make a fair
comparison, we employ a hybrid force-based score $\Omega_{\text{force}}$ with
a weighted sum of the engaging score $\Omega_{\text{engage}}$, the sticking
score $\Omega_{\text{stick}}$ and a kinematic heuristic $\Omega_{\text{dist}}$
\- the minimum distance between any two points on the object and the end-
effector. Specifically, $\Omega_{\text{engage}}$ uses the magnitude of the
contact normal force $\lambda_{\perp}$ to evaluate how much normal disturbance
forces a contact can withstand, $\Omega_{\text{engage}}=\lambda_{\perp}$.
$\Omega_{\text{stick}}$ refers to the minimal amount of disturbance force
required to transition the sticking contact into sliding, which is the
distance between the contact force and the friction cone edges (i.e., a wrench
margin away from failure):

$\Omega_{\text{stick}}=\left(\mu\lambda_{\perp}-|\lambda_{\mathbin{\\!/\mkern-5.0mu/\\!}}|\right)\cos(\arctan\mu),$
(11)

where $\mu$ refers to the static friction coefficient.

Figure 3: Qualitative evaluation for the planar pushing example. Data in the
subfigures indicate $(\Omega_{\text{suc}},\Omega_{\text{cap}})$. (i) and (ii)
refer to a failed and a successful trajectory towards a task objective of
fixturing a rectangular box (yellow) on the wall (black), respectively. Three
screenshots are taken along each trajectory, for which we run Algo. 2 and
visualize the expansive tree after running 100 iterations. The nodes on the
tree are shown in colored dots (CoM position of the box) ranging from blue to
red, implying the energy cost field. More reddish dots indicate nodes with
higher cost-to-come $c$. Box and circular end-effector (green) configurations
of three random nodes on the tree are visualized in partial transparency. Note
the capture set $\mathcal{Z}_{\text{cap}}$ (light blue) and the task success
set $\mathcal{Z}_{\text{suc}}$ (gray).

###  VI-C Data Generation and Implementation

We collected trajectories and per-frame system states of the 4 tasks to
compute their energy margins offline and compare them with ground-truth
robustness and success labels.

1) Ground-truth Data Generation: In the simulation, we generated and recorded
50 trajectories for each of the four examples in Fig. 2 by randomizing the
initial states of the object and the end-effector, the friction coefficients,
etc. 10 frames ($K=10$) are evenly selected from each trajectory and the
system states across the trajectory $\boldsymbol{z}_{k},k\in\\{1,...,K\\}$ are
recorded. As an instance, two of such scripted trajectories with three states
each are illustrated in Fig. 3. The hybrid force-based score
$\Omega_{\text{force}}(\boldsymbol{z}_{k})$ is computed by recording the
friction coefficient $\mu$, the contact normal force $\lambda_{\perp}$ and
lateral force $\lambda_{\mathbin{\\!/\mkern-5.0mu/\\!}}$ in each frame. A
trajectory is labeled $1$ (successful, e.g. Fig. 3 bottom) if the task
objective is accomplished in the last recorded frame, i.e.
$\boldsymbol{z}_{K}\in\mathcal{Z}_{\text{suc}}$, otherwise $0$ (Fig. 3 top);
and a state in a frame $\boldsymbol{z}_{k}$ is labeled $1$ (captured) if the
state is contained in the capture set for the current and subsequent $\hat{k}$
recorded frames
$\boldsymbol{z}_{k^{\prime}}\in\mathcal{Z}_{\text{cap}}(\boldsymbol{z}_{{k}^{\prime}})$
for $k^{\prime}\in\\{k,...,k+\hat{k}\\}$.

2) Implementation Details: To analyze the effectiveness of our proposed
metrics in predicting the robustness and success of manipulation tasks, we
utilized the Area Under the receiver operating characteristic Curve (AUC) and
the Average Precision (AP) as evaluation tools. We run the algorithms in
Section V and obtain $\Omega_{\text{esc}}\left(\boldsymbol{z}_{k}\right)$,
$\Omega_{\text{cap}}(\boldsymbol{z}_{k})$ and
$\Omega_{\text{suc}}(\boldsymbol{z}_{k})$ for each recorded state
$\boldsymbol{z}_{k}$ in each trajectory of the examples. We thereafter compare
the scores with ground-truth labels using AUC and AP, where higher values
indicate better performance across all possible classification thresholds. In
Fig. 3, we demonstrated the expansive trees after running Algo. 2. The
distribution of the energy cost field with respect to the capture set
$\mathcal{Z}_{\text{cap}}$ and the task success set $\mathcal{Z}_{\text{suc}}$
directly indicates the corresponding quality score values and energy margins.
To demonstrate the capability of predicting task success using the capture
score of success $\Omega_{\text{suc}}(\boldsymbol{z}_{k})$, we designed a
trajectory-level score $\Omega_{\text{suc}}(\boldsymbol{z},{\bar{k}})$ which
is the weighted average of $\Omega_{\text{suc}}(\boldsymbol{z}_{k})$ of only
the first $\bar{k}$ states with an increasing weight from $k=1$ to
$k=\bar{k}$. We thereby have, for each scenario, a dataset comprising 500
state-level data points
$(\Omega_{\text{esc}}\left(\boldsymbol{z}_{k}\right),\Omega_{\text{cap}}(\boldsymbol{z}_{k}))$
and 50 trajectory-level data points
$\Omega_{\text{suc}}(\boldsymbol{z},{\bar{k}})$ with corresponding labels.

Task | $\Omega_{\text{cap}}$ |  $\Omega_{\text{esc}}$ (EST) |  $\Omega_{\text{esc}}$ (RRT) |  $\Omega_{\text{force}}$ (Baseline) | $\Omega_{\text{suc}}$
---|---|---|---|---|---
Pushing |  0.97 / 0.99 | 0.92 / 0.98 |  0.92 / 0.99 | 0.86 / 0.95 | 0.98 / 0.98
Balancing |  0.98 / 0.99 | 0.93 / 0.98 | 0.94 / 0.98 | 0.94 / 0.96 | 0.93 / 0.95
Toppling |  0.94 / 0.93 | 0.93 / 0.91 | 0.93 / 0.92 | 0.77 / 0.62 | 0.89 / 0.92
Grasping |  1.00 / 1.00 | 0.99 / 0.99 | 0.97 / 0.99 | 0.73 / 0.88 | 0.99 / 0.99
TABLE I: Energy margins evaluation in simulation (AUC/AP).

###  VI-D Quantitative Analysis

1) Overall Performance: Upon examination of the dataset for each example in
simulation (Fig. 2), AUC and AP were computed for both the robustness
predictions ($\Omega_{\text{cap}}$, $\Omega_{\text{esc}}$ and
$\Omega_{\text{force}}$, Table I left) and the success predictions
($\Omega_{\text{suc}}$, Table I rightmost column). The results indicate a high
level of predictive capability of our metrics. The AUC/AP values for the
robustness predictions using our methods ($\Omega_{\text{cap}}$,
$\Omega_{\text{esc}}$) consistently exceed 0.9 across all examples,
demonstrating strong discriminative power over the baseline
$\Omega_{\text{force}}$. Here, $\Omega_{\text{esc}}$ is computed by employing
both RRT and EST as subroutine algorithms $A$ for Algo. 1. The values for the
success predictions given $\bar{k}=5$ indicate that we can reliably predict
successful trajectories with high precision using the score
$\Omega_{\text{suc}}$ we proposed.

$\displaystyle{10^{1}}$$\displaystyle{10^{2}}$$\displaystyle{10^{3}}$$0.8$$0.82$$0.84$$0.86$$0.88$$0.9$$0.92$$0.94$#
IterationsAUC / AP of success
ratioAUCAPtime$\displaystyle{10^{1}}$$\displaystyle{10^{2}}$$\displaystyle{10^{3}}$$0$$0.5$$1$$1.5$Run-
time per data (s) Figure 4:  Time efficiency and performance (AUC/AP values
for $\Omega_{\text{suc}}(\boldsymbol{z},{\bar{k}})$) of Algo. 2 in a planar
pushing task.

2) Algorithm Efficiency: We examined the time efficiency of Algo. 2 with
respect to the number of iterations and the task success predictions (AUC/AP
values for $\Omega_{\text{suc}}(\boldsymbol{z},{\bar{k}})$) as an ablation
study. Fig. 4 illustrates that it takes about 1.75 sec to finish 1000
iterations in the simulated planar pushing task (Fig. 2-i). The AUC/AP values
increase given more iterations of the algorithm, which is expected because of
more nodes on the tree and better state space coverage. AUC/AP values reach
0.9 given only 100 iterations and about 0.2 sec of runtime, which confirms the
time efficiency of our algorithm.

3) Robustness against Different Model Parameter Errors: A study was conducted
(in Fig. 5) on the robustness of the caging-in-motion algorithms given
different unknown model parameter errors, such as the positions
${\boldsymbol{r}_{x}},{\boldsymbol{r}_{y}}$ and the velocities
$\dot{\boldsymbol{x}},\dot{\boldsymbol{y}}$ of the object and the end-
effector, the contact forces $\lambda_{\mathbin{\\!/\mkern-5.0mu/\\!}}$,
$\lambda_{\perp}$, the friction coefficient $\mu$, etc. For instance, we
simulated perturbed friction coefficient $\hat{\mu}$ by adding a perturbation
${\mu}_{e}$ randomly, uniformly sampled from the interval
$[0,{e}_{\text{max}}]$ on the true values $\mu$, where ${e}_{\text{max}}$ is
the maximal perturbation, i.e. $\hat{\mu}={\mu}_{e}+\mu$.

For the simulated planar pushing task, we randomized the perturbed values
$\hat{\mu}$ given increasing thresholds
${e}_{\text{max}}\in\\{0,{e}^{1}_{\text{max}},...,{e}^{s}_{\text{max}}\\}$ for
each trajectory $\boldsymbol{z}$. We then ran Algo. 2 and obtained the capture
score $\Omega_{\text{cap}}(\boldsymbol{z}_{k})$ and the AP values with respect
to scaled maximal thresholds
$\bar{e}_{\text{max}}\in\\{0,{e}^{1}_{\text{max}}/{e}^{s}_{\text{max}},...,1\\}$
given the ground-truth labels (Fig. 5). Similarly, we obtained the (simulated)
perturbed positions $\hat{\boldsymbol{r}}_{x},\hat{\boldsymbol{r}}_{y}$ and
the perturbed velocities
$\hat{\dot{\boldsymbol{x}}},\hat{\dot{\boldsymbol{y}}}$ for each state-level
data point ${\boldsymbol{z}}_{k}$ and the corresponding AP curves of velocity
and position. The perturbed contact forces
$\hat{\lambda}_{\mathbin{\\!/\mkern-5.0mu/\\!}}$, $\hat{\lambda}_{\perp}$ are
generated for the force-based score
$\Omega_{\text{force}}(\boldsymbol{z}_{k})$ with the AP curve of force. Note
that the maximal thresholds ${e}^{s}_{\text{max}}$ are of varying values for
different types of error but exhibit comparable intensity of disturbance
empirically so that the curves under the scaled thresholds
$\bar{e}_{\text{max}}$ are comparable. It can be seen that higher noises
affect the capability of positivity identification of our algorithm. However,
our methods (the friction, velocity, and position curves) demonstrate better
robustness against model parameter errors than the baseline (the force curve).

$0$$0.1$$0.2$$0.3$$0.4$$0.5$$0.6$$0.7$$0.8$$0.9$$1$$0.75$$0.8$$0.85$$0.9$$0.95$$1$scaled
max perturbation $\bar{e}_{\text{max}}$AP of robustness
metricsfrictionvelocitypositionforce (BL) Figure 5:  Robustness evaluation of
Algo. 2 under various modeling errors: This plot illustrates the algorithm
performance against different levels of estimation errors of various types,
including friction coefficient ($\mu$), velocity
($\dot{\boldsymbol{x}},\dot{\boldsymbol{y}}$), relative position
(${\boldsymbol{r}_{x}},{\boldsymbol{r}_{y}}$), and contact forces
($\lambda_{\mathbin{\\!/\mkern-5.0mu/\\!}},\lambda_{\perp}$). The graph
demonstrates how AP values for the capture score $\Omega_{\text{cap}}$ (blue,
orange and green) and the force-based baseline (BL) score
$\Omega_{\text{force}}$ (red) vary given increasing maximal error thresholds
($\bar{e}_{\text{max}}$).

###  VI-E Real-World Experiment

1) Setup and Data Collection: We conducted real-world experiments using an
Interbotix WidowX-200 robot arm (Fig. 6-i) for the planar pushing task of
various geometric shapes of 3D-printed objects (five shapes) and end-effectors
(two shapes) (Fig. 6-ii). The arm has a position error between 5-8 mm,
providing a testbed to evaluate our scores under system uncertainty as in
Section VI-D-3. We used a joystick to manually control the end-effector
positions to perform the planar pushing task in the horizontal desktop plane.
Visual markers were attached to the object and the end-effector to record
their poses over time using a Realsense D415 camera.

We have ten combinations of objects and end-effectors. For each combination,
we collected 12 trajectories, with 100 recorded states for each trajectory.
The velocities were solved given the poses and the frequency of the recording
(33 Hz). A human observer provided success labels of the task objective. Two
failed (Fig. 6-iv) and successful (Fig. 6-v) trajectories are also illustrated
as examples. Given the recorded states and system kinematic and dynamic
parameters (mass, moment of inertia, friction coefficients, geometric
information, etc.), the trajectories were replicated in the simulated
environment (Fig. 6-iii). We utilize the simulation for computing energy
margins offline.

Figure 6: Real-world experiments settings. (i) The planar pushing task using
an Interbotix WidowX-200 robot arm with a top-view camera. (ii) 5 3D-printed
objects and 2 end-effectors. (iii) Replication of the scene in (i) in
simulation for computing the energy margins. (iv,v) Screenshots of a
failed/successful trajectory from the camera view.

2) Geometry Effect Evaluation: Statistical analysis results of the AUC/AP
values for the scores $\Omega_{\text{cap}}(\boldsymbol{z}_{k})$ and
$\Omega_{\text{suc}}(\boldsymbol{z},\bar{k})$, are summarized in Table II
following the same procedure as in the simulation (Section VI-C) given
$\bar{k}=60,K=100$. The effort of escape $\Omega_{\text{esc}}$ is not
discussed here due to its inferior performance compared to the capture score
in Table I. The results presented mostly exceed 0.85, which supports the
discrimination and predictivity power of our scores and is in line with the
results in the simulation. The performance of the capture score of success
$\Omega_{\text{suc}}(\boldsymbol{z})$ of 3 of the 10 combinations of
geometries (the jaw end-effector with the triangle, concave and irregular
objects) is below the average though. It could be a reason for the failure
that the relatively complicated geometric shapes and the shape complementarity
in these three cases added uncertainty to the planar pushing. It is thereby
barely reasonable to predict whether the task objective will be achieved given
a partial observation of such a trajectory.

Object | End-effector |  Robust $\Omega_{\text{cap}}$ |  Success $\Omega_{\text{suc}}$
---|---|---|---
Rectangle | Jaw | 0.85 / 0.89 | 0.94 / 0.95
| Round | 0.89 / 0.98 | 1.00 / 1.00
Triangle | Jaw | 0.94 / 0.99 | 0.69 / 0.62
| Round | 0.93 / 0.99 | 0.94 / 0.93
Convex | Jaw | 0.91 / 0.92 | 1.00 / 1.00
| Round | 0.84 / 0.81 | 0.91 / 0.88
Concave | Jaw | 0.93 / 0.96 | 0.67 / 0.78
| Round | 0.97 / 0.99 | 0.94 / 0.93
Irregular | Jaw | 0.95 / 0.99 | 0.47 / 0.36
| Round | 0.94 / 0.99 | 0.93 / 0.87
TABLE II: Real-world experiments result (AUC/AP).

##  VII Conclusion

In this paper, we propose to quantify manipulation robustness through energy
margins to failure and compute them using kinodynamic motion planning
algorithms. We believe this approach is a step towards robustness
characterization of general dexterous manipulation. Our approach currently has
certain limitations, such as reliance on dynamic rollouts in the simulation,
which requires accurate system modeling and causes computational costs.
Additionally, our analysis is confined to single rigid objects and utilizes an
empirically defined capture set. Looking ahead, we aim to improve our method
by overcoming these limitations and applying our robustness characterization
in robust motion planning for general dexterous manipulations.

## Acknowledgements

This work was partially funded by the European Commission under the Horizon
Europe Framework Programme project SoftEnable, grant number 101070600,
https://softenable.eu/. We extend our gratitude to Robert Gieselmann and
Rafael Cabral for their discussions and insights on the idea presented in this
work.

## References

  * [1] A. Bhatt, A. Sieler, S. Puhlmann, and O. Brock, “Surprisingly robust in-hand manipulation: An empirical study,” _arXiv preprint arXiv:2201.11503_ , 2022.
  * [2] H. Tsukamoto, S.-J. Chung, and J.-J. E. Slotine, “Contraction theory for nonlinear stability analysis and learning-based control: A tutorial overview,” _Annual Reviews in Control_ , vol. 52, pp. 135–169, 2021.
  * [3] R. Tedrake _et al._ , “Lqr-trees: Feedback motion planning on sparse randomized trees.” in _Robotics: Science and Systems_ , vol. 2009, 2009.
  * [4] A. D. Ames, S. Coogan, M. Egerstedt, G. Notomista, K. Sreenath, and P. Tabuada, “Control barrier functions: Theory and applications,” in _European control conference_. IEEE, 2019, pp. 3420–3431.
  * [5] M. A. Roa and R. Suárez, “Grasp quality measures: review and performance,” _Autonomous robots_ , vol. 38, pp. 65–88, 2015.
  * [6] Y. Hou and M. T. Mason, “Criteria for maintaining desired contacts for quasi-static systems,” in _Proc. Int. Conf. Intell. Robot. Syst._ IEEE, 2019, pp. 6555–6561.
  * [7] S. Makita and W. Wan, “A survey of robotic caging and its applications,” _Advanced Robotics_ , vol. 31, no. 19-20, pp. 1071–1085, 2017.
  * [8] Y. Dong and F. T. Pokorny, “Quasi-static soft fixture analysis of rigid and deformable objects,” _Proc. Int. Conf. Robot. Automat._ , 2024.
  * [9] W. Kuperberg, “Problems on polytopes and convex sets,” in _DIMACS Workshop on polytopes_ , 1990, pp. 584–589.
  * [10] E. Rimon and A. Blake, “Caging 2d bodies by 1-parameter two-fingered gripping systems,” in _Proc. Int. Conf. Robot. Automat._ , vol. 2. IEEE, 1996, pp. 1458–1464.
  * [11] ——, “Caging planar bodies by one-parameter two-fingered gripping systems,” _IEEE Int. J. Robot. Res._ , vol. 18, no. 3, pp. 299–318, 1999.
  * [12] A. Rodriguez, M. T. Mason, and S. Ferry, “From caging to grasping,” _IEEE Int. J. Robot. Res._ , vol. 31, no. 7, pp. 886–900, 2012.
  * [13] J. Mahler, F. T. Pokorny, Z. McCarthy, A. F. van der Stappen, and K. Goldberg, “Energy-bounded caging: Formal definition and 2-d energy lower bound algorithm based on weighted alpha shapes,” _IEEE Robot. Automat. Lett._ , vol. 1, no. 1, pp. 508–515, 2016.
  * [14] J. Mahler, F. T. Pokorny, S. Niyaz, and K. Goldberg, “Synthesis of energy-bounded planar caging grasps using persistent homology,” _IEEE Trans. Automat. Sci. Eng._ , vol. 15, no. 3, pp. 908–918, 2018.
  * [15] R. R. Ma, W. G. Bircher, and A. M. Dollar, “Modeling and evaluation of robust whole-hand caging manipulation,” _IEEE Trans. Robot._ , vol. 35, no. 3, pp. 549–563, 2019.
  * [16] N. S. Pollard, “Synthesizing grasps from generalized prototypes,” in _Proc. Int. Conf. Robot. Automat._ , vol. 3. IEEE, 1996, pp. 2124–2130.
  * [17] A. T. Miller and P. K. Allen, “Examples of 3d grasp quality computations,” in _Proc. Int. Conf. Robot. Automat._ , vol. 2. IEEE, 1999, pp. 1240–1246.
  * [18] C. Ferrari, J. F. Canny, _et al._ , “Planning optimal grasps.” in _Proc. Int. Conf. Robot. Automat._ , vol. 3, no. 4, 1992, p. 6.
  * [19] Y. Lin and Y. Sun, “Task-based grasp quality measures for grasp synthesis,” in _Proc. Int. Conf. Intell. Robot. Syst._ IEEE, 2015, pp. 485–490.
  * [20] J. Xu, M. Danielczuk, J. Ichnowski, J. Mahler, E. Steinbach, and K. Goldberg, “Minimal work: A grasp quality metric for deformable hollow objects,” in _Proc. Int. Conf. Robot. Automat._ IEEE, 2020, pp. 1546–1552.
  * [21] A. Saxena, L. L. Wong, and A. Y. Ng, “Learning grasp strategies with partial shape information.” in _AAAI_ , vol. 3, no. 2, 2008, pp. 1491–1494.
  * [22] R. Krug, Y. Bekiroglu, D. Kragic, and M. A. Roa, “Evaluating the quality of non-prehensile balancing grasps,” in _Proc. Int. Conf. Robot. Automat._ IEEE, 2018, pp. 4215–4220.
  * [23] T. Makapunyo, T. Phoka, P. Pipattanasomporn, N. Niparnan, and A. Sudsang, “Measurement framework of partial cage quality based on probabilistic motion planning,” in _Proc. Int. Conf. Robot. Automat._ IEEE, 2013, pp. 1574–1579.
  * [24] K. Hauser and Y. Zhou, “Asymptotically optimal planning by feasible kinodynamic planning in a state–cost space,” _IEEE Trans. Robot._ , vol. 32, no. 6, pp. 1431–1443, 2016.
  * [25] D. Hsu, J.-C. Latombe, and R. Motwani, “Path planning in expansive configuration spaces,” in _Proc. Int. Conf. Robot. Automat._ , vol. 3. IEEE, 1997, pp. 2719–2726.
  * [26] S. M. LaValle, “Rapidly-exploring random trees: A new tool for path planning,” _Research Report 9811_ , 1998.
  * [27] J. Liu, F. Feng, Y. C. Nakamura, and N. S. Pollard, “A taxonomy of everyday grasps in action,” in _Proc. Int. Conf. Humanoid Robot._ IEEE, 2014, pp. 573–580.
  * [28] F. Ruggiero, V. Lippiello, and B. Siciliano, “Nonprehensile dynamic manipulation: A survey,” _IEEE Robot. Automat. Lett._ , vol. 3, no. 3, pp. 1711–1718, 2018.