††thanks: These authors contributed equally††thanks: These authors contributed
equally††thanks: These authors contributed equally

# Supplementary Material for:
Tunable spin and valley excitations of correlated insulators in
$\Gamma$-valley moiré bands

Benjamin A. Foutty  Geballe Laboratory for Advanced Materials, Stanford, CA
94305, USA  Department of Physics, Stanford University, Stanford, CA 94305,
USA  Jiachen Yu  Geballe Laboratory for Advanced Materials, Stanford, CA
94305, USA  Department of Applied Physics, Stanford University, Stanford, CA
94305, USA  Trithep Devakul  Department of Physics, Massachusetts Institute of
Technology, Cambridge, Massachusetts 02139, USA  Carlos R. Kometter  Geballe
Laboratory for Advanced Materials, Stanford, CA 94305, USA  Department of
Physics, Stanford University, Stanford, CA 94305, USA  Yang Zhang  Department
of Physics, Massachusetts Institute of Technology, Cambridge, Massachusetts
02139, USA  Kenji Watanabe  Research Center for Functional Materials, National
Institute for Material Science, 1-1 Namiki, Tsukuba 305-0044, Japan  Takashi
Taniguchi  International Center for Materials Nanoarchitectonics, National
Institute for Material Science, 1-1 Namiki, Tsukuba 305-0044, Japan  Liang Fu
Department of Physics, Massachusetts Institute of Technology, Cambridge,
Massachusetts 02139, USA  Benjamin E. Feldman<EMAIL_ADDRESS>Geballe
Laboratory for Advanced Materials, Stanford, CA 94305, USA  Department of
Physics, Stanford University, Stanford, CA 94305, USA  Stanford Institute for
Materials and Energy Sciences, SLAC National Accelerator Laboratory, Menlo
Park, CA 94025, USA

##  I 1\. Spatial variation measured with scanning SET

In the main text, we focus on a single spatial location in each sample. In
Fig. S1, we present spatial dependence of inverse electronic compressibility
$d\mu/dn$ taken in each sample at magnetic field $B=11$ T, where the
incompressible states are strongest. Here, $\mu$ is the chemical potential and
$n$ is the carrier density, with the convention that increasing hole density
is negative. Both devices are highly homogeneous over micron length scales,
with little twist angle disorder in the regions examined. We estimate the
total twist angle disorder by fitting the positions of the measured
incompressible states and comparing their spatial dependence. In Sample S1, we
compare the difference in gate voltage between the peak at filling factor
$\nu=-1$, which occurs around back gate voltage $V_{\textrm{BG}}=-7$ V and the
$\nu=-\frac{1}{3}$ peak, which occurs around $V_{\textrm{BG}}=-3.2$ V. This
yields a twist angle that varies between $\theta=3.38^{\circ}$ and
$\theta=3.45^{\circ}$ in Area 1, where the measurements in the main text were
conducted (Fig. S1a-b, Fig. S2a). In a second region (Area 2) in Sample S1
that is about 4 $\mu$m away, the twist angle varies between
$\theta=3.42^{\circ}$ and $\theta=3.61^{\circ}$ (Fig. S1c-d, Fig. S2a). In
Sample S2, we fit the difference between the incompressible $\nu=-1$ gap at
$V_{\textrm{BG}}=-3.95$ V and the incompressible peak at $V_{\textrm{BG}}=-3$
V, and estimate a twist angle that varies between $\theta=2.55^{\circ}$ and
$\theta=2.65^{\circ}$ (Fig. S1e). We note that this peak may not be a charge-
ordered state, but appears to occur at a roughly constant filling factor
position and therefore can be used as a benchmark to estimate the twist angle
disorder that behaves better than the band edge, which can vary due to
underlying potential disorder in the sample. If we assume a constant band edge
and just treat the position of the $\nu=-1$ gap in gate voltage as the only
parameter, this leads to an even smaller estimate of twist angle disorder.
Across both devices, we measure significantly less twist angle inhomogeneity
than is typically observed across magic-angle twisted bilayer graphene devices
Uri  _et al._ (2020); Zondiner  _et al._ (2020); Yu  _et al._ (2021). This
difference likely results from the larger overall twist angles in these
devices which limit excessive relaxation/reconstruction and/or increased
stiffness due to the presence of more layers.

Figure S1: Spatial dependence at $B=11\textrm{ T}$. a, Measured locations in
Area 1 of Sample S1. $x$ and $y$ label relative spatial coordinates. b,
Inverse electronic compressibility d$\mu/$d$n$ as a function of back voltage
$V_{\textrm{BG}}$ at the positions shown in a, labelled by index. c, Locations
of measurements in a far-separated region of Sample S1. $X$ and $Y$ label
spatial coordinates within this second region. d, d$\mu/$d$n$ at positions
shown in c. e, Spatial variation of d$\mu$/d$n$ along a continuous linecut in
Sample S2. Figure S2: Optical images of the samples. a, Optical image of
Sample S1. The measurements in the main text are all taken in the vicinity of
the purple dot, while similar behavior was observed in the vicinity of the red
dot (Fig. S1c-d) and measurements of the adjacent Bernal bilayer region were
taken at the orange dot (Fig. S10). b, Optical image of Sample C2. Scale bars
in a and b are both 5 $\mu$m.

##  II 2\. Gap extraction

The gaps $\Delta_{\nu}$ presented in the main text correspond to the step in
the chemical potential at filling factor $\nu$, found by taking the difference
between the local maximum and local minimum in $\mu$ on either side of the
gap. In the scanning single-electron transistor (SET) measurements, we use the
DC measurement of $\mu(n)$ to extract charge gap sizes, rather than
integrating the AC measurement of $d\mu/dn$ because the high resistance of the
contacts and/or sample artificially enhances the a.c. signal of the
incompressible peaks.

In Sample S1, the $\nu=-1/3$ gap appears as $\mu(n)$ begins to turn up near
the mobility edge even in the d.c. measurement scheme. To accurately extract
the gap, we fit a quadratic polynomial curve to the surrounding data and
subtract this background before evaluating the gap size. Variation of the
bounds of this fitted background do not affect the measured gap size nor its
quantitative trend with $B$ (i.e., the effective $g$-factor). An example
showing data at 5 T is plotted in Fig. S3. This is the only gap for which we
perform any background subtraction; other gaps are well separated from the
mobility edge.

In the global capacitance measurements, we only have direct access to
$d\mu/dn$, rather than the chemical potential itself. To determine the charge
gaps, we convert the capacitance signal to $d\mu/dn$ (see Supplementary Sec.
13 below) and then integrate over density to find $\mu(n)$. The charge gaps
are then taken as the resulting step on this $\mu(n)$ curve.

Figure S3: Background fitting near $\nu=-1/3$. Background subtraction around
the $\nu=-1/3$ gap in Sample S1, measured at $B=5$ T. Black dashed lines show
local maximum and minimum values of $\mu_{+}$ and $\mu_{-}$ around the gap,
the difference of which is taken to be $\Delta_{-1/3}$

##  III 3\. Continuum modeling

We derive effective continuum models for twisted double bilayer WSe2 (tdWSe2)
from first principles calculations at high symmetry stacking configurations.

The low energy physics of hole-doped monolayer WSe2 consists of (local)
quadratic band maxima at the $\Gamma$ and $\pm K$ points of the Brillouin
zone. The $\Gamma$ valley band maximum is spin-degenerate, while the $\pm K$
band maxima are spin polarized with spin $s_{z}=\pm\frac{1}{2}$ due to spin-
orbit coupling. The structure of four-layer WSe2 is such that each layer is
rotated $180^{\circ}$ relative to its adjacent layers, so that interlayer
tunneling is strongly suppressed in the $\pm K$ valley bands.

Three different high-symmetry stacking configurations, MM′, MX′, and XX′, as
illustrated in the main text, are considered. First principles density
functional theory (DFT) calculations for four-layer WSe2, shown in Fig. 1e of
the main text, reveal that the valence band maximum is at $\Gamma$ for all
three stacking configurations. The energies of the top four valence bands at
$\Gamma$, and the top two at $K$, are listed in Table 1.

Band |  $E_{\Gamma}$ (MM′) |  $E_{\Gamma}$ (MX′) |  $E_{\Gamma}$ (XX′) |  $E_{K}$ (MM′) |  $E_{K}$ (MX′) |  $E_{K}$ (XX′)
---|---|---|---|---|---|---
1 | -4.5164 | -4.7022 | -4.6223 | -4.6523 | -4.8350 | -4.7064
2 | -4.6412 | -4.8144 | -4.7011 | -4.6602 | -4.8569 | -4.7168
3 | -4.9506 | -5.1330 | -5.0938 | - | - | -
4 | -5.3199 | -5.5013 | -5.3622 | - | - | -
Table 1: Energy (in eV) obtained from DFT at the $\Gamma$ and $K$ points, for
each of the three high symmetry stacking regions MM′, MX′, and XX′. Note that
each band is doubly degenerate due to inversion and time reversal symmetry.

We first consider an effective single-particle continuum model for the twisted
system of the form

$H_{\Gamma}=-\frac{|\mathbf{k}|^{2}}{2m_{\Gamma}}+V_{\Gamma}(\mathbf{r})$ (1)

where $m_{\Gamma}$ is the $\Gamma$ valley effective mass, and $V(\mathbf{r})$
describes a periodic potential with the moiré periodicity. The spin degeneracy
is implicit. The potential $V(\mathbf{r})$ can be obtained by taking the first
harmonic approximation,
$V_{\Gamma}(\mathbf{r})=2v\sum_{i=1,3,5}\cos(\mathbf{g}_{i}\cdot\mathbf{r}+\phi)+c$,
from which the three parameters, $v,c,\phi$, are chosen to match the energy of
the first $\Gamma$ band in Table 1 at the three positions
$\mathbf{R}_{\mathrm{MM}^{\prime}},\mathbf{R}_{\mathrm{MX}^{\prime}}$, and
$\mathbf{R}_{\mathrm{XX}^{\prime}}$. Here, $\mathbf{g}_{i}$ are moiré
reciprocal lattice vectors defined
$\mathbf{g}_{i}=\frac{4\pi}{\sqrt{3}a_{M}}(-\sin\frac{2\pi(i-1)}{6},\cos\frac{2\pi(i-1)}{6})$,
and $\mathbf{R}_{\mathrm{MM}^{\prime}}=0$,
$\mathbf{R}_{\mathrm{MX}^{\prime}}=\frac{a_{M}}{\sqrt{3}}(0,-1)$, and
$\mathbf{R}_{\mathrm{XX}^{\prime}}=\frac{a_{M}}{\sqrt{3}}(0,1)$. With this
definition, $V(\mathbf{R}_{\mathrm{MM}^{\prime}})=6v\cos(\phi)+c$,
$V(\mathbf{R}_{\mathrm{MX}^{\prime}})=6v\cos(\phi+2\pi/3)+c$, and
$V(\mathbf{R}_{\mathrm{XX}^{\prime}})=6v\cos(\phi-2\pi/3)+c$. We find the
parameters $v=17.9$ meV, $\phi=0.44$ (radians), and $c=-4613.6$ meV. We use
the moiré period $a_{M}=3.32\ \mathrm{\AA}/\sin(3.4^{\circ})\approx 56\
\mathrm{\AA}$, and effective mass $m_{\Gamma}=1.2m_{e}$, where $m_{e}$ is the
electron mass Movva  _et al._ (2018). The hole density obtained from filling
the first moiré band of this continuum model is shown in Fig. 1f of the main
text.

To incorporate the effect of a vertical displacement field, we now construct a
model which includes all four layer degrees of freedom. We assume the form for
the model Hamiltonian at zero displacement field
$H_{\Gamma}=-\frac{|\mathbf{k}|^{2}}{2m_{\Gamma}}\mathbb{I}+W_{\Gamma}(\mathbf{r})$
where $\mathbb{I}$ is the identity matrix and

$W_{\Gamma}(\mathbf{r})=\begin{pmatrix}V_{1\Gamma}(\mathbf{r})&T_{1\Gamma}(\mathbf{r})&0&0\\\
T_{1\Gamma}(\mathbf{r})&V_{2\Gamma}(\mathbf{r})&T_{2\Gamma}(\mathbf{r})&0\\\
0&T_{2\Gamma}(\mathbf{r})&V_{2\Gamma}(\mathbf{r})&T_{1\Gamma}(\mathbf{r})\\\
0&0&T_{1\Gamma}(\mathbf{r})&V_{1\Gamma}(\mathbf{r})\end{pmatrix},$ (2)

in which all $T$, $V$, are periodic functions with the moiré periodicity, and
are chosen such that the eigenvalues of ${W}_{\Gamma}(\mathbf{R}_{\alpha})$
with $\alpha=MM^{\prime},MX^{\prime},XX^{\prime}$, match with
$E_{\Gamma}(\alpha)$ in Table 1. To uniquely specify $T$, we have taken it to
be real and positive. At each $\mathbf{R}_{\alpha}$, there are four parameters
$\\{V_{1\Gamma}(\mathbf{R}_{\alpha}),V_{2\Gamma}(\mathbf{R}_{\alpha}),T_{1\Gamma}(\mathbf{R}_{\alpha}),T_{2\Gamma}(\mathbf{R}_{\alpha})\\}$,
which can be chosen to fit the four eigenvalues $E_{\Gamma}$. From these, we
can use the first harmonic approximation to obtain

$\begin{split}V_{1\Gamma}(\mathbf{r})/\mathrm{meV}&=-4799.5+16.4\sum_{i=1,3,5}2\cos(\mathbf{g}_{i}\cdot\mathbf{r}+0.62)\\\
V_{2\Gamma}(\mathbf{r})/\mathrm{meV}&=-5093.6+18.5\sum_{i=1,3,5}2\cos(\mathbf{g}_{i}\cdot\mathbf{r}+0.47)\\\
T_{1\Gamma}(\mathbf{r})/\mathrm{meV}&=218.6+1.2\sum_{i=1,3,5}2\cos(\mathbf{g}_{i}\cdot\mathbf{r}+1.84)\\\
T_{2\Gamma}(\mathbf{r})/\mathrm{meV}&=220.3-7.8\sum_{i=1,3,5}2\cos(\mathbf{g}_{i}\cdot\mathbf{r}+2.18)\end{split}$
(3)

which specifies the continuum model.

A similar construction can be made for the $K$ or $K^{\prime}$ valley bands.
We focus on the $K$ valley, which is related to the $K^{\prime}$ valley by
time reversal. The main difference from the $\Gamma$ bands is that interlayer
tunneling is suppressed, and hence
$H_{K}=-\frac{|\mathbf{k}|^{2}}{2m_{K}}+\mathrm{diag}(V_{1K}(\mathbf{r}),V_{2K}(\mathbf{r}),V_{2K}(\mathbf{r}),V_{1K}(\mathbf{r}))$,
and $\mathbf{k}$ is now measured relative to the $K$ point. The first two
energies $E_{K}$ at each stacking region specify the potentials $V_{1K}$ and
$V_{2K}$. However, there is now an additional ambiguity as the eigenvalues are
invariant under interchanging $V_{1K}\leftrightarrow V_{2K}$. Additional
information is needed to determine whether the first pair of bands should be
associated with $V_{1K}$ (the outer two layers) or $V_{2K}$ (the inner two
layers). From direct inspection of the DFT wave function, we find that the
first band has larger component on the outer two layers at all stacking
regions. Hence, $V_{1K}>V_{2K}$ and the potentials are fully specified in the
first harmonic approximation:

$\begin{split}V_{1K}(\mathbf{r})/\mathrm{meV}&=-4731.2+18.1\sum_{i=1,3,5}2\cos(\mathbf{g}_{i}\cdot\mathbf{r}+0.75)\\\
V_{2K}(\mathbf{r})/\mathrm{meV}&=-4744.6+19.5\sum_{i=1,3,5}2\cos(\mathbf{g}_{i}\cdot\mathbf{r}+0.76)\end{split}$
(4)

We use the effective mass $m_{K}=0.5m_{e}$ Movva  _et al._ (2018). The moiré
band structure for the $\Gamma$ and $K$ valley bands is shown in Fig S4.

Figure S4: Moiré bands from the continuum model. The $\Gamma$ (left) and $K$
(right) valley moiré bands obtained from the four layer continuum model.

To address the competition between $\Gamma$ and $K$ valley bands in a
displacement field, we model the displacement field by a layer-dependent
potential difference $\Delta$ via the term
$H_{\Delta}=-\Delta\;\mathrm{diag}(0,1,2,3)$. At $\nu=0$, the layer potential
is related to the displacement field $D$ by $\Delta=dD/\epsilon$, where
$\epsilon$ is the relative out-of-plane dielectric constant and $d$ is the
interlayer spacing.

We remark that the competition between $\Gamma$ and $K$ band edges from first
principles is often delicate, and can depend sensitively on relaxed lattice
structure and van der Waals density functional (see, for instance, Refs. Wu
_et al._ (2019) and Devakul  _et al._ (2021), in which different functionals
result in different $\Gamma$,$K$ band orderings in both twisted homobilayer
WSe2 and MoTe2). In addition, interaction effects may also result in a
renormalization of the band gap measured experimentally. To compare
quantitatively with experiments, we thus take as an additional tuning
parameter a constant energy offset $\Delta_{K}$ between $\Gamma$ and $K$ band
edges, modifying the Hamiltonian as $H_{\Gamma}\rightarrow
H_{\Gamma}-\Delta_{\Gamma K}/2$ and $H_{K}\rightarrow H_{K}+\Delta_{\Gamma
K}/2$. The experimentally measured gap is maximum at $D=0$, and decreases
quickly with non-zero $D$: this is consistent with the scenario in which the
state at $\nu=-1$ corresponds to filling a spin-polarized $\Gamma$ moiré band
and the gap being to the first $K$ or $K^{\prime}$ moiré band. We choose
$\Delta_{\Gamma K}=0.1$ eV such that the gap between the first $\Gamma$ and
$K$ bands roughly match the experimentally measured value at $D=0$. Once
$\Delta_{\Gamma K}$ is fixed, the dependence of the gap on $\Delta$ follows
from the continuum model, as shown in Fig. 3f of the main text.

Figure S5: $K$-valley band density distribution. Density distribution in the
first $K$ valley band without (left) and with (right) the electrostatic
potential from the filled $\Gamma$ band at $\nu=-1$, calculated from the
continuum model. Without interactions, the charges are localized at MM′,
whereas including the electrostatic potential shifts them to XX′.

To conclude this section, we comment on the nature of the lowest energy $K$
valley excitations at $\nu=-1$. In the absence of interactions, it is clear
from the moiré potential that holes prefer to localize on the MM′ sites for
both $\Gamma$ and $K$ valley. However, holes filling into the $K$ valley
beyond $\nu=-1$, will feel an effective potential resulting from electrostatic
repulsion with the filled $\Gamma$ band. This effect can be captured by a
Hartree potential in the continuum model, $V_{H}(\mathbf{r})=\int
d\mathbf{r}V_{C}(\mathbf{r}-\mathbf{r}^{\prime})\langle
n^{(\Gamma)}(\mathbf{r}^{\prime})\rangle$ where $\langle
n^{(\Gamma)}(\mathbf{r}^{\prime})\rangle$ is the density distribution of the
filled first $\Gamma$ band (obtained from the non-interacting continuum
model), and $V_{C}(\mathbf{r})$ is the Coulomb potential. The continuum model
description of the first $K$ band now includes $V_{H}(\mathbf{r})$ in addition
to the moiré potential. We further simplify by taking $\mathbf{r}$ to be
purely in the 2D plane and use a vertically displaced Coulomb potential
$V_{C}(\mathbf{r})=e^{2}/(4\pi\epsilon\sqrt{\mathbf{r}^{2}+d^{2}})$, where
$\epsilon=3.5\epsilon_{0}$ and $d=7$ Å. Figure S5 shows the hole density of
the first $K$ valley band with and without the Hartree term. As can be seen,
holes in the $K$ valley prefer to localize in the XX′ regions once
electrostatic repulsion from the filled $\Gamma$ band is taken into account,
resulting in an effective honeycomb lattice description near $\nu=-1$.

##  IV 4\. Alternative excitations at $\nu=-\frac{1}{3}$

In the main text, we mention an alternative explanation of the saturating gap
at $\nu=-\frac{1}{3}$, namely that below $B=7$ T, exchange coupling causes the
lowest energy additional hole to be anti-aligned with the spin-polarized
ground state, until Zeeman coupling overcomes this barrier and the additional
hole is aligned with the rest of the spins. In this section, we expand on that
picture and discuss why it is not supported by our data and calculations. In
this picture, the change in behavior around $B=7$ T occurs on the excitations
on the higher hole doping side of the gap (i.e., $\mu_{-}$ in the main text).
If the changing spin character of the additional hole excitation were driving
the saturating gap, then we would expect the discontinuity in the slope of the
chemical potential on the band edge to occur at $\mu_{-}$, rather than as we
observe at $\mu_{+}$. Additionally, the energy scale of $J$ that our data
would predict in this case is not necessarily consistent with other
microscopic parameters of the problem. To lowest order, we could estimate
exchange coupling $J$ as $\frac{1}{3}g^{*}\mu_{B}B^{*}$ where $g^{*}$ is the
effective $g$-factor that we measure for the gap, $\mu_{B}$ is the Bohr
magneton, $B^{*}=7$ T is where the behavior of the gap changes, and the factor
of $\frac{1}{3}$ is due to the fact that an added hole couples to three
nearest neighbors. This yields $J\approx 0.6$ meV, which is quite large
compared with other parameters in the problem. For example, if we take $J\sim
4t^{2}/U$ where $t$ is the hopping and $U$ is the onsite potential, then we
would estimate $J\approx 0.01$ meV, as $t\approx 0.2$ meV from the bandwidth
of the lowest moiré band from the continuum model while $U\approx 100$ meV due
to the strength of Coulomb interactions Zhang  _et al._ (2020).

##  V 5\. Extraction of chemical potential $\mu_{+}$ and $\mu_{-}$

The SET tip directly measures changes in chemical potential of the sample by
tracking a particular Coulomb blockade oscillation of the SET while the sample
density is tuned by a gate. In order to evaluate whether the discontinuity in
the slope of the gap at $\nu=-\frac{1}{3}$ comes from a change on the particle
or hole side of the gap, we need to evaluate the behavior of the chemical
potential on either side of the gap as a function of magnetic field, which is
the “slow” axis in our measurement. One approach is to pick a particular point
(say the band edge, or some location far from the gap) and choose this as a
reference point at which $\mu$ is set to zero at every magnetic field.
However, this can give different behavior depending on the choice of reference
point, especially when the overall behavior of $\mu(n)$ changes smoothly with
magnetic field. For example, in Fig. 4 in the main text, it is apparent that
the chemical potential at the band edge is controlled by the magnitude of
negative compressibility; if the minimum point is chosen to be the “reference”
point where $\mu=0$ for each value of $B$, then the trends in both $\mu_{+}$
and $\mu_{-}$ will mostly be determined by the change in the negative
compressibility, rather than identifying the behavior in the vicinity of the
gap itself.

To avoid such artifacts, we take advantage of the fact that the gap at
$\nu=-\frac{1}{3}$ is small compared to the background chemical potential. We
can fit a smooth polynomial background to the chemical potential on either
side of the gap at $\nu=-\frac{1}{3}$ and subtract this background, as in Fig.
S3. Then the value of the chemical potential on either side of the gap is
measured relative to the background chemical potential of filling the lowest
energy moiré band between 0 and 1. The precise values of $\mu_{+}$ and
$\mu_{-}$ and can vary slightly with details of the fit, but the fact that the
$\textrm{d}\mu_{+}/\textrm{d}B$ changes at $B=7$ T while
$\textrm{d}\mu_{-}/\textrm{d}B$ does not is consistent. We note that this
essentially is a statement of the evolution of the chemical potential relative
to the background filling of the lowest moiré band.

At $\nu=-1$, the behavior of $\mu(n)$ in both samples is dramatically
different on either side of the gap, preventing an equivalent type of
subtraction. Focusing on Sample S1, in Fig. S6 we show extracted chemical
potential on either side of the $\nu=-1$ gap with two different choices of
reference chemical potential. In Fig. S6a, we set the reference chemical
potential to $\nu=-1.25$. Practically, this means that we take the average
measured chemical potential in a certain region (red shaded box in Fig. S6c)
and subtract it from the measurement at each field such that it becomes the
“zero” point of the chemical potential, before then reading off $\mu_{+}$ and
$\mu_{-}$ on either side of the gap. In Fig. S6b, we perform the same analysis
but with the reference point taken near the middle of the first moiré band
(blue shaded box in Fig. S6c). In both of these cases, $\mu_{-}$ does not
change very much; most of the gap dependence thus relates to the change of
$\mu_{+}$. Given that we choose partially filled moiré bands as references,
this suggests that the dependence of $\mu_{-}$ on magnetic field is quite
close to the overall change in chemical potential of the bands themselves,
while $\mu_{+}$ varies strongly relative to this. This qualitatively agrees
with the theoretical picture described in the main text, as we would expect
that at all experimentally accessible fields, the lowest-energy excitation at
$\mu_{-}$ would be adding another spin-aligned hole, whose dependence on $B$
would be identical to filling the spin-polarized bands. The excitation at
$\mu_{+}$ on the other hand, would be due to the formation of spin-polarons,
which should have a different (stronger) dependence on magnetic field than the
surrounding moiré bands (Fig. 2i in the main text). Therefore, our data also
suggest spin polaron excitation are relevant at $\nu=-1$, but the smoking gun
signature of a change in gap size and chemical potential evolution with
magnetic field are not accessible in our measurement setup.

Figure S6: Behavior of chemical potential around $\nu=-1$. a, Chemical
potential on the left ($\mu_{-}$) and right ($\mu_{+}$) sides of $\nu=-1$
taking the reference $\mu$ at high densities around $\nu=-1.25$. b, Same as a,
but with the reference $\mu$ set to intermediate densities around $\nu=-0.5$.
c, $\mu(n)$ measured at $B=0$. The average $\mu$ within the red shaded box is
set to zero for the data in a, while the average $\mu$ within the blue shaded
box is set to zero for the data in b.

##  VI 6\. Spin $g$-factors

To quantitatively compare the experimental gaps and their increase in a
magnetic field to our theoretical model, we must consider the Zeeman
contribution to the effective $g$-factors of the various excitations. Our
measurement of a gap at filling $\nu$, $\Delta_{\nu}\equiv\Delta\mu$ is
equivalent to the particle-hole excitation energy, given by
$E[N_{\nu}+1]+E[N_{\nu}-1]-2E[N_{\nu}]$, where $E[N]$ is the ground state
energy at particle number $N$, and $N_{\nu}$ denotes the particle number at
filling $\nu$. First, we consider the case at high magnetic fields (i.e., the
right side of Fig. 2i in the main text), where the ground state is a spin-
polarized insulator with charge transfer energy gap $\Delta_{\textrm{CT}}$,
the high-density excitation is an additional interstitial hole with aligned
spin, and the low-density excitation is the absence of a hole. Then, we can
write down the energies as follows, which results in a gap that is independent
of magnetic field:

$E[N_{\nu}]=-\frac{1}{2}g_{s}N_{\nu}B,\quad
E[N_{\nu}+1]=\Delta_{\textrm{CT}}-\frac{1}{2}g_{s}(N_{\nu}+1)B,\quad
E[N_{\nu}-1]=-\frac{1}{2}g_{s}(N_{\nu}-1)B$
$\Delta_{\nu}=E[N_{\nu}+1]+E[N_{\nu}-1]-2E[N_{\nu}]=\Delta_{\textrm{CT}}-\frac{1}{2}g_{s}B+\frac{1}{2}g_{s}B=\Delta_{\textrm{CT}}$

At intermediate fields (i.e., the left side of Fig. 2i in the main text),
where the only change is that the hole excitation is the formation of a spin
polaron, we instead have

$E[N_{\nu}]=-\frac{1}{2}g_{s}N_{\nu}B,\quad
E[N_{\nu}+1]=\Delta_{\textrm{CT}}-E_{b}-\frac{1}{2}g_{s}(N_{\nu}+1)B,\quad
E[N_{\nu}-1]=-\frac{1}{2}g_{s}(N_{\nu}-2)B+\frac{1}{2}g_{s}B$
$\Delta_{\nu}=E[N_{\nu}+1]+E[N_{\nu}-1]-2E[N_{\nu}]=\Delta_{\textrm{CT}}-E_{b}-\frac{1}{2}g_{s}B+\frac{3}{2}g_{s}B=\Delta_{\textrm{CT}}-E_{b}+g_{s}B$

where $E_{b}>0$ is the spin polaron binding energy. Given that the bare spin
$g$-factor $g_{s}=2$, the measured effective $g$-factor of our gapped states
would be 2 in the absence of orbital/Berry curvature and exchange effects.

As mentioned in the main text, orbital and Berry curvature differences between
different moiré bands (at $K$ and $\Gamma$) can lead to additional
contributions to $g^{*}$ at $\nu=-1$. The momentum-resolved orbital
magnetization $m(k)$ can vary widely across the the moiré Brillouin zone (mBZ)
Grover  _et al._ (2022). Theoretical calculation of $m(k)$ requires
calculation of the Berry curvature on the single-particle moiré bands and can
depend sensitively on shape of the moiré bands, and therefore the twist angle
Grover  _et al._ (2022); He  _et al._ (2020); Devakul  _et al._ (2021).
Because our probe is only sensitive to the lowest energy excitation, the
measured gap will depend sensitively on the extremal value of $m(k)$ across
the mBZ, which can lead to different effective $g$-factors between samples
with different twist angles, and may contribute to the slightly different
$g$-factors measured at $\nu=-1$ in Samples S1 and S2.

##  VII 7\. Capacitance data from Sample C2

In Fig. S7, we present d$\mu$/d$n$ of Sample C2 as a function of $D$ and
$\nu$. The data are extracted from “top-gate capacitance” in which we apply a
10 mV a.c. excitation at 337 Hz to the top gate and measure the response on
the HEMT gate, which is electrically connected to the sample contacts. The
data were taken at $T=$ 4.2 K. The range that we can access on the lower
(upper) left side are limited by top (bottom) gate dielectrics.

Figure S7: Displacement field dependence in a second device. Measured
d$\mu$/d$n$ in Sample C2 as a function of displacement field $D$ and moiré
filling factor $\nu$.

##  VIII 8\. Estimating $D$ in SET measurements

In general, the effective layer potential difference in a singly-gated two-
dimensional (2D) heterostructure device is a difficult problem that depends in
part on details of the material being gated Kou  _et al._ (2014); McCann
(2006). Even the usual definition of an applied electric displacement field,
which takes the average of the displacement field on either side of the gated
2D material, is an approximation and not necessarily directly convertible to
microscopic parameters such as layer potential differences. Comparing our
scanning SET measurements, in which the samples have no top gate, to
measurements with a controlled displacement field is therefore nontrivial.
However, a lowest order estimate can be made by approximating the set-up as a
dual gated device where the SET tip is treated as an effective “top gate” held
at the same voltage as the sample. This voltage constraint is accurate because
the relative sample-tip voltage is held by d.c. feedback in our experiment to
minimize tip-induced doping. Additionally, the tip is held $\sim 50-100$ nm
away in vacuum, such that the effective “top-gate” capacitance will be quite
small compared to that of the back gate, and any small offsets in voltage will
not provide very much doping. We note that this is distinct from what has been
reported in recent STM measurements, in which the close tip-sample distance
requires the displacement field to be approximated by the voltage difference
between tip and back gate Kerelsky  _et al._ (2021). Because the SET tip is
far-separated and senses a spatial area underneath the tip that convolves both
tip size and distance, any voltage separation between tip and sample leads to
(small) non-uniform doping of the sample, resulting in broadening features
rather than reliably doping the sample, which limits our ability to change the
effective displacement field.

Using this approximation, the effective displacement field can be written as
$D_{\textrm{eff}}=\frac{1}{2\epsilon_{0}}(0-C_{b}(V_{b}-V_{0}))=-\frac{1}{2\epsilon_{0}}en$.
Here, $C_{b}$ is the back gate capacitance, $V_{b}$ is the back gate voltage,
$V_{0}$ is offset of the density due to intrinsic sample doping, and we’ve
used that $C_{b}(V_{b}-V_{0})/e=n$. Applying this expression to relevant
densities in our samples, in Sample S1 at $\nu=-1$, $D_{\textrm{eff}}=0.35$
V/nm, and in sample S2, $D_{\textrm{eff}}=0.18$ V/nm; this nonzero
displacement field explains the reduction in gap size relative to the $D=0$
value we measure in dual gated devices. We note that for Sample S1, this
estimated $D$ is near and slightly above the critical displacement field at
which we stop seeing insulating behavior at $\nu=-1$ in dual gated capacitance
devices, though the precise value varies slightly across the two devices
measured. However, the fact that we measure a gap at all at $\nu=-1$ in the
singly-gated devices is strongly suggestive that the lowest energy moiré bands
being filled are still at $\Gamma$, and suggest that the estimated
$D_{\textrm{eff}}$ should be taken as a rough approximation. Additionally,
global measurements in our dual gated devices inevitably average over
disorder, which may lead to lower overall gap sizes relative to local sensing
with the scanning SET. This will artificially decrease the critical
displacement field at which the $\nu=-1$ gap vanishes in the dual gated
devices.

##  IX 9\. Negative compressibility in Sample S1

Though we focus on Sample S2 in the discussion of negative
$\textrm{d}\mu/\textrm{d}n$ in the main text, we also observe negative
compressibility over large regions of hole doping in Sample S1, before an
upturn at low hole density (Fig. S8). The magnitude of negative
compressibility, as estimated from the average $|\textrm{d}\mu/\textrm{d}n|$
between $\nu=-.6$ and $\nu=-0.95$ also grows in a magnetic field, though no
clear saturation is observed.

Figure S8: Negative d$\mu/$d$n$ in Sample S1. a, $\mu(n)$ as measured in S1
across the entire density range. b, Average value of
$|\textrm{d}\mu/\textrm{d}n|$ between $n=-3.8\times 10^{12}$ cm-2 and
$n=-2.5\times 10^{12}$ cm-2.

##  X 10\. Fitting negative compressibility $\mu(n)$

Here we expand on the details of the fit to $\mu(n)$ in Sample S2 shown in
Fig. 4 in the main text. While the order of magnitude of the coefficient $A$
is robust to variation of fitting parameters/schemes, as is its increase by a
factor of $\sqrt{2}$ from $B=0$ to high fields, the precise value of $A$ is
somewhat sensitive to details of the fitting procedure. In Fig. S9, we present
the extracted $A$ from 3 different fits, changing the range of density
included in the fit. The blue data is that which is presented in the main
text, because fitting over the widest possible density range seems to provide
the best fit to the data. However, given the complicated details of the moiré
bands, it is not clear which of these is most physically justifiable. Fitting
across such a large density range means that we include filling over multiple
moiré bands, each of which will have a varying effective density of states as
they are filled. Additionally, fitting across the charge-ordered gap at
$\nu=-1$, around which the effective density $n$ may be better approximated by
$|n-n_{\nu}|$, where $n_{\nu}$ is the charge density at moiré band filling
$\nu$. It is apparent in these fits that between $n=-2\times 10^{12}\
\textrm{cm}^{-2}$ and $n=-2.5\times 10^{12}\ \textrm{cm}^{-2}$, just beyond
$\nu=-1$, the behavior of $\mu(n)$ varies sharply compared to the background.
For the red line, which fits to just the lowest moiré band between $\nu=0$ and
$\nu=-1$, the fit has no linear component. In the blue and green fits in Fig.
S9b, the size of the linear term is $\frac{\pi\hbar^{2}}{m^{*}}\approx
1.25\times 10^{-11}\ \textrm{meV/cm}^{-2}$ and
$\frac{\pi\hbar^{2}}{m^{*}}\approx 1.75\times 10^{-11}\ \textrm{meV/cm}^{-2}$,
respectively. Both of these terms would suggest effective masses much larger
than expected for moiré bands ($m^{*}=0.1-0.15m_{e}$), which is inconsistent
with the size of the Landau level gaps we observe simultaneously (Fig. S10).
Because these fits are over the entirety of multiple moiré bands, the
effective mass is not necessarily expected to remain constant across the fit.

Figure S9: Dependence of coefficient $A$ on fit parameters. a, Extracted
coefficient $A$ from a fit to $\mu(n)=\frac{\pi\hbar^{2}}{m^{*}}n-A\sqrt{n}$
where $m^{*}$ and $A$ are free parameters. Blue, green, and red data
correspond to fitting the entire curve out to $n=-5\times 10^{12}$ cm-2,
$n=-4\times 10^{12}$ cm-2, and $n=-2\times 10^{12}$ cm-2 ($\nu=-1$),
respectively. b, Measured $\mu(n)$ at $B=5T$ (black) overlaid with fitted
lines corresponding to those described in panel a.

##  XI 11\. Band edge determination

In order to quote accurate densities, we need to determine the location of the
band edge in gate voltage, which is set by both local (global) intrinsic
doping in the sample for SET (dual gated) measurements and by work function
differences between gates and sample Gustafsson  _et al._ (2018). For both
samples measured with the SET, the gate voltages at which we observe Landau
level gaps allows for unambiguous determination of the band edge, which also
matches well to the mobility edge (Fig. S10). In Sample S1, we don’t observe
Landau levels in the tdWSe2 region of the device (Fig. 2a), which may reflect
the large effective mass of the moiré bands and/or the proximity to charge
ordered states at the low fillings we measure. In comparison, Landau levels
are present in an adjacent Bernal bilayer WSe2 region of Sample S1, from one
of the “halves” of the flake that makes up the device. Assuming the entire
flake has similar intrinsic doping, this sets the gate voltage offset of the
band edge in the twisted double bilayer region of the device as well. In
Sample S2, we observe two-fold degenerate landau levels between $\nu=-1$ and
$\nu=-3$ near the noise floor of our measurement, with gap sizes of order
$100\ \mu$V at $B=11$ T. This corresponds to a filling factor range which we
were unable to explore in Sample S1.

In Samples C1/C2, we did not conduct measurements at nonzero magnetic field
and thus must rely on the capacitance measurements themselves to estimate the
band edge. In Fig. S11, we show a linetrace of raw capacitance signal taken at
a much higher contact gate voltage. This has the effect of improving the
sample contact resistance without affecting the device region measured via
capacitance to the top gate, as the top gate screens out any stray field from
the contact gates which are above it. In most of the measured parameter space,
the signal was indeed unaffected. However, it did lead to a nonzero signal
extending to lower and lower densities. We take the point where this signal
becomes zero within the signal to noise as zero density, and note that this
also agrees with the spacing between the $\nu=-1$ and $\nu=-3$ features in
Fig. 4 in the main text (i.e., the gate voltage range between $\nu=-1$ and
$\nu=-3$ is twice that between $\nu=-1$ and the band edge, $\nu=0$).

Figure S10: Landau levels and band edge determination. a, Landau levels
observed in a Bernal bilayer region in Sample S1. Dotted lines are filling
factor $\nu=3,4,\ldots,24,25$. No Landau levels were observed in the tdWSe2
region of the device. b, Landau levels in tdWSe2 beyond $\nu=-1$ in Sample S2.
These show a two-fold degeneracy; dotted lines are filling factor
$\nu=14,16,\ldots,24$. Figure S11: Capacitance signal measured at higher
contact gate voltages. A 9 mV a.c. signal is applied to the top gate at a
frequency of 911 Hz, and top and back gate voltages are swept to maintain
$D=0$ V/nm.

##  XII 12\. Temperature dependence of capacitance data

The temperature dependence of Sample C1 is shown up to $T=15$ K in Fig. S12.
The magnitude of the incompressible peak at $\nu=-1$ weakens with increasing
temperature. In Fig. S12a-b, $X$-channel (capacitive) data is presented in
units of the geometrical top-gate capacitance $c_{t}$. As the temperature is
increased, the charging of the tdWSe2 sample at low densities improves. This
is particularly notable at the maroon dashed line ($\nu=-\frac{2}{3}$) where
at $T=4.2$ K the capacitance is at $\approx 0.7c_{t}$ but it quickly rises as
the mobility edge extends rightwards. A more modest version of this effect is
also seen at $\nu=-1$ (black), whereas in the negative compressibility region
(pink), the temperature dependence is nonmonotonic, first increasing
(presumably) due to better charging before decreasing at the highest
temperatures. In Fig. S12c, we show line cuts at constant $D=-0.2$ V/nm at a
few temperatures. The incompressible peak at $\nu=-1$ is weaker at 12 K
relative to 4.2 K, but still clearly evident. In Fig. S13 we show the behavior
of the $\nu=-1$ gap as a function of $D$ at a few distinct temperatures. As in
previous reports, the gaps show a linear decrease in gap size as the
temperature increases, likely due to thermal broadening Li  _et al._ (2021).
The slope of gap size $\Delta_{-1}$ with $D$ is similar across various
temperatures.

Figure S12: Temperature dependence in Sample C1. a, Temperature dependence of
measured capacitance $C_{T}$ in units of geometric capacitance $c_{t}$. An 8
mV excitaiton at a frequency of 739.1 Hz and a constant d.c. voltage $-1.6$ V
is applied to the top gate. b, Line cuts as a function of temperature from
panel a at $\nu=-1$, $\nu=-0.875$, $\nu=-\frac{2}{3}$. c Measured d$\mu$/d$n$
at $D=-0.2$ V/nm at various temperatures. Figure S13: Temperature dependence
of the $\nu=-1$ gap.

##  XIII 13\. Estimation of quantum capacitance from global capacitance

As discussed in the Methods, we convert the capacitance signal into units of
the geometric capacitance $c_{t}$ by assuming $C_{T}=0$ in the gap between the
WSe2 valence and conduction bands and $C_{T}=c_{t}$ at high hole dopings
$>5\times 10^{12}\textrm{cm}^{-2}$. Using a lumped circuit capacitance model
Li  _et al._ (2021); Shi  _et al._ (2020), the measured top gate capacitance
$C_{T}$ is given by

$C_{T}=\frac{c_{t}c_{q}}{c_{t}+c_{b}+c_{q}}$ (5)

where $c_{t}$ and $c_{b}$ are the geometric top and bottom gate capacitances
and $c_{q}\equiv e^{2}(d\mu/dn)^{-1}$ is the quantum capacitance. As discussed
in the main text, this approximation will overestimate $d\mu/dn$ in the case
of high sample resistance or high frequency. Some of this enhancement can be
accounted for by including the simultaneous measurement of the out-of-phase
quadrature of the measured capacitance, which is equivalent to the dissipation
signal Goodall  _et al._ (1985); Shi  _et al._ (2022). In this case, the total
measured signal is given by

$V^{\textrm{complex}}\propto
C_{T}/c_{t}=\frac{c_{q}}{c_{t}+c_{b}+c_{q}}\frac{\tanh(\alpha)}{\alpha}$ (6)
$\alpha=\sqrt{i\omega\frac{c_{q}c_{t}}{c_{t}+c_{b}+c_{q}}R_{s}}$ (7)

where $\omega$ is the measurement frequency and $R_{s}$ is the sample
resistance. Note that in the small $\alpha$ limit, the imaginary (capacitive)
part of this signal agrees with Eq. 5 and the real (dissipative) part is zero.
This system can be solved for $c_{q}$ and $R_{s}$ using both the $x$ and $y$
quadratures, and the resulting $c_{q}$ yields $d\mu/dn$ as before. An example
of $x$ and $y$ quadrature capacitance signal and the extracted $d\mu/dn$ by
both of these methods are shown in Fig. S14. All data presented from
capacitance measurements in the main text show $d\mu/dn$ calculated using this
correction. The main effect of including the dissipation on the extracted
$d\mu/dn$ occurs in the regions of parameter space which exhibit negative
compressibility; there is also a consistent small decrease of the broad peak
at $\nu=-1$ (Fig. S14).

In Fig. S15 we show a comparison of the measured charge gap at $\nu=-1$ as a
function of applied displacement field $D$ at various frequencies. At 4.7 kHz,
the measured gap is clearly enhanced above the lower frequency curves. Below
1.7 kHz and 739 Hz, there is little measured change. This lack of change
suggests that the extraction of $d\mu/dn$ from Eqs. 6 and 7 performs
adequately at these measurement frequencies.

Figure S14: Example capacitance measurements on Sample C1. a, Line cuts
showing the capacitance signal in the $X$ (capacitance) and $Y$ (dissipation)
channels. b, Extracted $d\mu/dn$ signal using each respective method discussed
in the text. Data were measured at a temperature of 4.2 K with an 8 mV
excitaiton applied to the top gate at frequency 739.1 Hz and a constant d.c.
voltage $-1.6$ V on the top gate. Figure S15: Frequency dependence. Measured
gap $\Delta_{-1}$ in Sample C1 at various frequencies.

## References

  * Uri  _et al._ (2020) A. Uri, S. Grover, Y. Cao, J. A. Crosse, K. Bagani, D. Rodan-Legrain, Y. Myasoedov, K. Watanabe, T. Taniguchi, P. Moon, M. Koshino, P. Jarillo-Herrero, and E. Zeldov, Nature 581, 47 (2020), number: 7806 Publisher: Nature Publishing Group.
  * Zondiner  _et al._ (2020) U. Zondiner, A. Rozen, D. Rodan-Legrain, Y. Cao, R. Queiroz, T. Taniguchi, K. Watanabe, Y. Oreg, F. von Oppen, A. Stern, E. Berg, P. Jarillo-Herrero, and S. Ilani, Nature 582, 203 (2020), number: 7811 Publisher: Nature Publishing Group.
  * Yu  _et al._ (2021) J. Yu, B. A. Foutty, Z. Han, M. E. Barber, Y. Schattner, K. Watanabe, T. Taniguchi, P. Phillips, Z.-X. Shen, S. A. Kivelson, and B. E. Feldman, arXiv:2108.00009 [cond-mat] (2021), arXiv: 2108.00009.
  * Movva  _et al._ (2018) H. C. Movva, T. Lovorn, B. Fallahazad, S. Larentis, K. Kim, T. Taniguchi, K. Watanabe, S. K. Banerjee, A. H. MacDonald, and E. Tutuc, Physical Review Letters 120, 107703 (2018), publisher: American Physical Society.
  * Wu  _et al._ (2019) F. Wu, T. Lovorn, E. Tutuc, I. Martin, and A. MacDonald, Physical Review Letters 122, 086402 (2019), publisher: American Physical Society.
  * Devakul  _et al._ (2021) T. Devakul, V. Crépel, Y. Zhang, and L. Fu, Nature Communications 12, 6730 (2021), number: 1 Publisher: Nature Publishing Group.
  * Zhang  _et al._ (2020) Y. Zhang, N. F. Q. Yuan, and L. Fu, Physical Review B 102, 201115 (2020), publisher: American Physical Society.
  * Grover  _et al._ (2022) S. Grover, M. Bocarsly, A. Uri, P. Stepanov, G. Di Battista, I. Roy, J. Xiao, A. Y. Meltzer, Y. Myasoedov, K. Pareek, K. Watanabe, T. Taniguchi, B. Yan, A. Stern, E. Berg, D. K. Efetov, and E. Zeldov, arXiv:2203.01010 [cond-mat] (2022), 10.48550/arXiv.2201.06901, arXiv:2201.06901 [cond-mat] type: article.
  * He  _et al._ (2020) W.-Y. He, D. Goldhaber-Gordon, and K. T. Law, Nature Communications 11, 1650 (2020), number: 1 Publisher: Nature Publishing Group.
  * Kou  _et al._ (2014) A. Kou, B. E. Feldman, A. J. Levin, B. I. Halperin, K. Watanabe, T. Taniguchi, and A. Yacoby, Science 345, 55 (2014), publisher: American Association for the Advancement of Science.
  * McCann (2006) E. McCann, Physical Review B 74, 161403 (2006), publisher: American Physical Society.
  * Kerelsky  _et al._ (2021) A. Kerelsky, C. Rubio-Verdú, L. Xian, D. M. Kennes, D. Halbertal, N. Finney, L. Song, S. Turkel, L. Wang, K. Watanabe, T. Taniguchi, J. Hone, C. Dean, D. N. Basov, A. Rubio, and A. N. Pasupathy, Proceedings of the National Academy of Sciences 118, e2017366118 (2021), publisher: Proceedings of the National Academy of Sciences.
  * Gustafsson  _et al._ (2018) M. V. Gustafsson, M. Yankowitz, C. Forsythe, D. Rhodes, K. Watanabe, T. Taniguchi, J. Hone, X. Zhu, and C. R. Dean, Nature Materials 17, 411 (2018).
  * Li  _et al._ (2021) T. Li, J. Zhu, Y. Tang, K. Watanabe, T. Taniguchi, V. Elser, J. Shan, and K. F. Mak, Nature Nanotechnology 16, 1068 (2021).
  * Shi  _et al._ (2020) Q. Shi, E.-M. Shih, M. V. Gustafsson, D. A. Rhodes, B. Kim, K. Watanabe, T. Taniguchi, Z. Papić, J. Hone, and C. R. Dean, Nature Nanotechnology 15, 569 (2020), number: 7 Publisher: Nature Publishing Group.
  * Goodall  _et al._ (1985) R. K. Goodall, R. J. Higgins, and J. P. Harrang, Physical Review B 31, 6597 (1985), publisher: American Physical Society.
  * Shi  _et al._ (2022) Q. Shi, E.-M. Shih, D. Rhodes, B. Kim, K. Barmak, K. Watanabe, T. Taniguchi, Z. Papić, D. A. Abanin, J. Hone, and C. R. Dean, Nature Nanotechnology , 1 (2022), publisher: Nature Publishing Group.


# Task-Based Design and Policy Co-Optimization
for Tendon-driven Underactuated Kinematic Chains

Sharfin Islam12, Zhanpeng He13, Matei Ciocarlie2
https://roamlab.github.io/tentamorph/ 1Equal contributions. 2Dept. of
Mechanical Engineering, 3Dept. of Computer Science, Columbia University, New
York, USA

###### Abstract

Underactuated manipulators reduce the number of bulky motors, thereby enabling
compact and mechanically robust designs. However, fewer actuators than joints
means that the manipulator can only access a specific manifold within the
joint space, which is particular to a given hardware configuration and can be
low-dimensional and/or discontinuous. Determining an appropriate set of
hardware parameters for this class of mechanisms, therefore, is difficult -
even for traditional task-based co-optimization methods. In this paper, our
goal is to implement a task-based design and policy co-optimization method for
underactuated, tendon-driven manipulators. We first formulate a general model
for an underactuated, tendon-driven transmission. We then use this model to
co-optimize a three-link, two-actuator kinematic chain using reinforcement
learning. We demonstrate that our optimized tendon transmission and control
policy can be transferred reliably to physical hardware with real-world
reaching experiments.

##  I Introduction

Underactuated manipulators require a carefully designed transmission, often
tendon-driven, to take advantage of a reduced number of actuators in the
robot. Such designs range from planar serial chains with relatively few links
to complex, hyper-redundant continuum robots [1, 2, 3, 4, 5, 6, 7]. In all of
these cases, being able to reduce the number of actuators means that we can
build smaller and more lightweight designs, place actuation at more proximal
locations in the chain, and take advantage of passive compliance in the un-
actuated degrees of freedom.

However, the compromise of such designs is that, with fewer actuators than
degrees of freedom, underactuated manipulators can directly access only a
certain manifold within the overall state space. This manifold, which contains
the set of obtainable states for a particular hardware configuration, can be
low dimensional and/or discontinuous. These limitations affect our ability to
plan controllers that smoothly transition between various states to accomplish
a desired task. The process of tuning the hardware parameters in order to
ensure the set of accessible states matches the desired task can be slow and
cumbersome, particularly if, as is often the case, changes in the design
parameters of the underactuated transmission have a counter-intuitive effect
the overall behavior of the robot.

Figure 1: We optimize an underactuated, tendon-driven transmission for
kinematic chains. We formulate and parameterize a general model for N links
and M actuators. We apply our model to a three-link, two-actuator tentacle
that we co-optimize using reinforcement learning (top row, left). We then
validate our results on physical hardware (right).

Co-optimization, or the process of simultaneously searching the space of both
hardware and control, is a possible solution to the problem of ensuring that a
given hardware design is capable of a desired task. Such methods have been
attempted in the context of underactuated kinematic chains, but often
restricted to simulation [8, 9, 10], limited validation of the simulated
controller on real hardware [11], or only implemented on single-actuator
systems [12]. The fundamental difficulty of such approaches lies in
formulating a co-design problem that a) enables the use of non-trivial
controllers or policies, b) can be solved to an acceptable optimum point, c)
guarantees that the final result can be physically realized in real hardware,
and finally d) ensures that the optimized control policy also transfers to
real hardware without substantial loss of performance. This is a difficult set
of goals to achieve simultaneously, and, to the best of our knowledge, no
current method has done so for underactuated, tendon-driven transmissions with
multiple actuators.

In this paper, we start with an underactuated tendon transmission model for
general planar kinematic chains, which can capture a diverse set of hardware
configurations. We use a forward model that captures the behavior of the
system as a function of design parameters and control inputs. We then show
that this model enables end-to-end co-optimization of design and control
policies for a specified task. To solve the co-optimization problem, we adopt
MORPH [13], a recently introduced end-to-end co-optimization framework that
uses a proxy model to mimic the effect of hardware, and generate gradients of
hardware parameters w.r.t. tasks performance. This allows us to learn both
control and hardware parameters to accomplish desired tasks. We then show that
it is possible to build a real robotic manipulator with this optimized
transmission, and validate that our optimized control policy can be
transferred reliably to real hardware.

An additional feature of our hardware design is that some design parameters
can be modified without requiring re-assembly of the manipulator, and by using
the same set of fabricated components. This allows for different levels of
flexibility in robot behaviors. If the task requires explorations of diverse
behaviors, we can optimize all hardware parameters to explore a large solution
space. In the case that the task is simpler, we can optimize only the easily
adaptable parameters to avoid reassembly. While in this paper we focus on a
single kinematic chain optimized for simple reaching tasks, our directional
goal to enable underactuated, tendon transmission models that can be co-
optimized along with their control policies, which can in turn transfer to the
real world. Overall, the main contributions of this paper include:

  * •

We formulate a model of underactuated, tendon-driven, passively-compliant
planar kinematic chains that enables efficient end-to-end task-based
optimization of the design and control parameters.

  * •

We show that the results of the co-optimization process can be transferred to
real hardware implementing the optimized design parameters. The resulting
robots can then run the optimized control policy which also achieves sim2real
transfer in task performance.

  * •

To the best of our knowledge, this is the first time that task-based, policy
and design co-optimization methods have been demonstrated for underactuated
manipulators with multi-dimensional manifolds.

##  II Related Work

An underactuated manipulator is a mechanical system composed of links and
joints that has fewer actuators than degrees of freedom [14]. The transmission
of these manipulators are often tendon-driven, as a single actuated tendon can
be routed to control multiple joints. Reducing the number of bulky actuators
enables designers to build more compact and lightweight manipulators.
Moreover, tendon-transmission allow designers to dislocate the motors from the
joints. Dislocating the actuators reduces the inertia of the links and makes
it easier to make the manipulator robust to water, dust, and other difficult
environmental conditions [15]. With so many practical benefits, a myriad of
diverse underactuated, tendon-driven manipulators have been proposed over the
past several decades.

There is an extensive history of underactuated manipulator design. The first
underactuated manipulator, Hirose’s soft gripper introduced in 1978, was
designed to softly capture a diverse range of objects with uniform pressure
[1]. Over time, these manipulator designs became much more advanced and their
applications now extend to more robust and intricate grasping behaviors [16,
17, 18, 19, 20]. Currently, there is extensive research in the design and
control of underactuated manipulators for tasks even more dexterous than
grasping, such as in-hand manipulation [21, 22, 23].To realize these more
advanced capabilities, the tendon-transmission of these manipulators had to be
carefully designed, optimized, and further hand-tuned. This process is
cumbersome and time-consuming, but necessary. Co-optimization methods are a
possible tool to help design these complicated mechanisms. With recent
advancements in reinforcement learning, co-optimization is now more powerful
than ever.

Recently, researchers have considered reinforcement learning for task-driven
design and control co-optimization [24, 25, 26, 27]. The key of this line of
work is to derive policy gradient w.r.t both the control and design
parameters. Chen and He et al. [12] propose to integrate a differentiable
model of the hardware with a control policy to adapt hardware design via
policy gradients. A key limitation of this approach is the requirement of
differentiable physics simulation. In cases when the forward transition cannot
be modeled in a differentiable manner, researchers propose methods to learn
design parameters in the input or output space of a policy. For instance, Luck
et al. [28] propose to learn an expressive latent space to represent the
design parameters and condition a policy with a latent design representation.
Transform2Act [29] propose to have a transform stage in their policy that
estimates actions to modify a robot’s design and a control stage that computes
control sequences given a specific design. Our hardware design is not
differentiable. Hence, in this work, we apply MORPH [13], a method that co-
optimizes design and control in parameter space that does not require
differentiable physics.

##  III Method

We formalize our problem as follows. Our goal is to build a kinematic chain
(i.e., tentacle/trunk) robot that can achieve flexible behaviors, such as
reaching desired parts of the workspace, while maintaining the practical
design benefits of underactuation. For this class of mechanisms, the design
parameters and controller for a given task are innately coupled. Therefore, we
take a task-based hardware optimization approach to search for hardware
parameters that yield manifolds such that we can smoothly transition between
desired states.

Our method has three main components: 1) an under-actuated, tendon-driven
transmission design for compliant kinematic chains; 2) a model that captures
the forward dynamics of our designed robot; and 3) an end-to-end co-
optimization pipeline that can directly learn the parameters of the proposed
design along with a control policy for a given reaching task.

###  III-A Transmission design

Figure 2: Our flexible tendon-transmission design for compliant, underactuated
kinematic chains. In this design, $N$ links are driven by $M$ actuated flexion
tendons. We also implement a passive extension mechanism, which we can
precisely pre-tension. The parameters of our transmission are all flexion
tendon radii, extension tendon radii, and elastic tendon pre-elongations.

Our design, shown in Figs. 1 and 2, assume a kinematic chain with $N$ links
driven by $M$ motors. All motors are located inside the base and actuate
winches that are connected to the joints via tendons. For each motor $i$, we
denote the radius of the actuator as $R^{A}_{i}$. We assume that each motor is
driving a tendon that traverses the entire length of the chain, thereby
helping flex every joint. We denote the length of each link $j$ as $L_{j}$.
For a the corresponding joint $j$, each motor $i$ will wrap around a pulley.
We use ${R_{ij}}^{f}$ to denote the radius of the flexion pulley for this
joint and motor pair.

We assume all actuators provide flexion forces, and the transmission uses
purely passive extension mechanisms. At each joint, the mechanism features a
passive elastic tendon that stretches over a pulley of constant radius
(denoted by $R^{e}$) to provide a restoring extension torque. Each elastic
tendon can be pretensioned individually; we use ${l_{j}}^{pre}$ to denote the
pretensioning elongation of the passive tendon at joint $j$.

This mechanism has a number of desirable characteristics. Underactuation leads
to a small number of motors, and placing all motors inside the base frees the
links of the kinematic chain from any motors or electronics. However, the
movement of the robot is non-trivial to define or control. Critically, robot
movement is determined not only by the actuators but also by a number of
design parameters. These include all flexion tendon radii, extension tendon
radii, and elastic tendon pre-elongations. Furthermore, the ability of a robot
to reach specific points in its workspace is clearly also determined by the
lengths of the links. Our goal is to devise a procedure capable of optimizing
all these design parameters while providing a policy for controlling motor
movements in order to achieve a specific task.

In addition, our proposed design makes some of the design parameters easier to
change than others. In particular, we mount each elastic tendon on a linear
slider mechanism. The position of this slider is set by rotating a lead-screw,
thereby controlling the pre-tension elongations of the elastic extension
tendons to less than 1mm of precision. This means that some of our design
parameters (pulley radii, link lengths) are more difficult to change, as they
require a full reassembly, while others (elastic tendon elongations) are
easier to change depending on the task. We want to leverage this ability in
our co-optimization procedure.

###  III-B Forward actuation model for our transmission design

In order to enable a co-optimization routine for our transmission design, we
first need a forward actuation model that relates the hardware parameters and
the actuator inputs to the resulting state of the manipulator. The input to
this model consists of the servo angle of the winches that control flexion
tendons, which we write as a column vector $\theta_{A}$ of size
$\mathbb{R}^{M\times 1}$. The next state is defined by the set of joint angles
$\theta_{J}$ of size $\mathbb{R}^{N\times 1}$. In other words, the actuation
model must predict joint angles $\theta_{J}$ as a function of the actuator
commands $\theta_{A}$, as well as all design parameters described in the
previous section. We formulate this actuation model by searching for
$\theta_{J}$ that minimizes the total stored energy, denoted by $U$, while
still meeting the constraints imposed by the rigid flexion tendons.

We begin our model by first formulating this constraint: since the flexion
tendons cannot extend, we know that the tendons will either be in tension or
accumulate slack. This slack is a difference between the collective change in
length of the tendon due to the motion of the joints versus the change in
length due to servo-driven winch. The slack, therefore, is a function of the
flexion pulley radii, actuated winch radii, and joint angles.

The flexion radii can be composed into a matrix of size $\mathbb{R}^{M\times
N}$ as follows:

$\displaystyle\mathbf{R^{F}}=\begin{bmatrix}R_{11}^{F}&R_{12}^{F}&...&R_{1N}^{F}\\\
R_{21}^{F}&R_{22}^{F}&...&R_{2N}^{F}\\\ ...&...&...&...\\\
R_{M1}^{F}&R_{M2}^{F}&...&R_{MN}^{F}\\\ \end{bmatrix}$ (1)

The radii of the actuated winches can be similarly composed into a diagonal
matrix of size $\mathbb{R}^{M\times M}$.

$\displaystyle\mathbf{R^{A}}=\begin{bmatrix}R_{1}^{A}&0&...&0\\\
0&R_{2}^{A}&...&0\\\ ...&...&...&...\\\ 0&0&...&R_{M}^{A}\\\ \end{bmatrix}$
(2)

With these matrices, we can now define a column vector of the slack collected
in the tendons $(w)$.

$\displaystyle w=\mathbf{R^{F}}\theta_{J}-\mathbf{R^{A}}\theta_{A}$ (3)

The slack is important as it limits the set of obtainable states for any given
action. We know that the robot will settle at a configuration in this set of
states that will minimize the overall stored elastic energy $(U)$ [14]. This
stored energy comes from the preloaded elongation $l_{pre}$ of the elastic
tendons, but also the elongation due to joint flexion. The total elongation
for any j-th link $(l_{j})$ can be written as follows:

$\displaystyle\Delta l_{j}=R^{e}\theta_{j}+l^{pre}_{j}$ (4)

We can now pose the forward actuation model to solve for the next state as
$\theta_{J}^{\prime}=f(\theta_{J},\theta_{A})$ relating the state to the
action as a numerical optimization problem: given $\theta_{A}$,$\theta_{J}$,

$\displaystyle\operatorname*{argmin}_{\theta_{J}}$ $\displaystyle
U=\frac{1}{2}\sum_{j=1}^{N}kl_{j}^{2}$ (5) subject to
$\displaystyle\text{\hskip 10.00002pt }w\geq 0$ (6) $\displaystyle\text{\hskip
10.00002pt }det(diag(w))=0$ (7)

Since the flexion tendons are in-extensible, the change in length of the
tendon at the joints must always be greater than the change in length of the
tendon due to the actuator. The slack, therefore, must always be greater than
zero, as shown in Eq. 6. Additionally, if the joint angles are non-zero, the
tendon must always be in tension, and therefore, at least one element in
$\vec{w}$ must be zero. This constraint is shown in Eq. 7. We enforce these
constraints While minimizing the total stored energy given in Eq.5,

Actions and hardware parameters do not change the global energy landscape, but
instead, the section of the landscape that satisfies our constraints (see
Fig.3). Hence, optimizing control $\theta_{A}$ and hardware parameters (e.g.,
$\mathbf{R_{F}}$ and $R_{e}$) is finding appropriate energy manifolds whose
minimum energy states are ones we want to visit for task completion. There
exists some combination of hardware parameters that ensures these manifolds
are continuous are have a clear energy minima. Arriving at this specific set
of parameters relies on observing how the changes in parameters affect the
most suitable control strategy for solving a given task. Finally, in this
work, we focus on optimizing the following set of hardware parameters:
$\phi=(\mathbf{R^{F}},l_{pre},L)$.

Figure 3: Illustrations of our optimization-based forward actuation model. (a)
and (b) show the global energy landscape (z-axis). The yellow regions are
manifolds that satisfy our constraints (Eq.6 and Eq.7). The red arrows in (a)
represent optimization steps (Eq.5) that find the energy minimum inside the
manifold. Changing design parameters and control actions result in a change in
the manifold.

###  III-C Task-aware co-optimization of design and control

Equipped with the forward actuation model described in the previous section,
we can now focus on the problem of co-optimizing a set of design parameters
and a motor control policy for a given task. Using a standard reinforcement
learning (RL) formalism, we model each task as a Markov Decision Process
(MDP), or a tuple $(\mathcal{S},\mathcal{A},\mathcal{F},\mathcal{R})$, where
$\mathcal{S}$ is the state space, $\mathcal{A}$ is the action space,
$\mathcal{R}(\bm{s},\bm{a})$ is the reward function, and
$\bm{s}^{\prime}=\mathcal{F}(\bm{s},\bm{a})$ is the forward transition model,
where $\bm{s},\bm{s}^{\prime}\in\mathcal{S}$, and $\bm{a}\in\mathcal{A}$.

In this work, we use a state vector $s$ containing the joint angles
$\theta_{J}$, motor positions $\theta_{A}$, and 2D end-effector positions
$\mathrm{x}_{e.e.}$. Actions for a control policy are motor commands
$\theta_{A}$. The reward function $\mathcal{R}(\bm{s},\bm{a})$ encodes the
desired performance; for example, if we want the end-effector to reach a
specific point, the reward will comprise a negative distance to the goal. For
each task, we want to find the parameters of a control policy $\pi_{\theta}$
as well as the design parameters $\phi$ that optimize task performance, which
is evaluated by its expected returns:
$\mathbb{E}[\sum_{t=0}^{T}\mathcal{R}(\bm{s}_{t},\bm{a}_{t})]$.

The state transition model $\mathcal{F}$ requires additional considerations.
As described in Section III-B, it consists of an optimization-based model $f$
whose behaviors are determined by some of the hardware parameters from $\phi$.
However, this model is non-differentiable and, therefore, can not be used in a
standard policy gradient optimization routine.

To co-optimize the hardware design with control, we apply MORPH [13], an end-
to-end co-optimization method that uses a neural network proxy model $h^{nn}$
to approximate the forward transition $\mathcal{F}$. The proxy model and a
control policy are co-optimized with task performance while asking $h^{nn}$ to
be close to $\mathcal{F}$. In this work, instead of mimicking just the forward
transition, we consider the effect of hardware design parameters in task
space. Hence, we ask the hardware proxy to approximate both the forward
transition and forward kinematics: $q_{\phi}=f\circ T_{FK}$. Note that our
hardware parameters $\phi$ are encapsulated in different parts of $q$. For
example, link lengths only affect forward kinematics, while pulley radii and
preloaded tension govern the behaviors of $f$. Hence, the overall optimization
objective is:

$\displaystyle\max_{\theta,\psi}\mathbb{E}_{\pi,h^{nn}}[\underbrace{\sum_{t=0}^{T}\mathcal{R}(s_{t},\theta_{A})}_{\text{task
performance}}+\underbrace{\alpha||h^{nn}(\theta_{J},\theta_{A})-q(\theta_{J},\theta_{A})||^{2}}_{\text{hardware
constraint}}]$ (8)

where $\alpha$ is a constant number.

To derive explicit hardware design parameters, for every $N$ epochs, we also
use CMA-ES [30] to search $\phi$ to match the updated hardware proxy $h^{nn}$:

$\displaystyle\min_{\phi}||h^{nn}(\theta_{J},\theta_{A})-q_{\phi}(\theta_{J},\theta_{A})||^{2}$
(9)

The final algorithm is an iterative process: We first co-optimize the hardware
proxy and the control policy to improve task performance, then extract
explicit hardware parameters that match the current version of the hardware
proxy.

##  IV Experimental Set-up

###  IV-A Design and control co-optimization

To evaluate our method, we optimize the aforementioned design with several
goal-reaching tasks, i.e., reaching a goal location in 2D with its end-
effector. We have three sets of experiments demonstrating three ways to
utilize our co-optimization pipeline for flexible motor skills:

  1. 1.

Goal reaching via re-fabrication: We adapt all design and control parameters
to different goals but only train each policy for a single goal. In this
specific experiment, reaching multiple goals requires the re-fabrication of a
real robot since the link length and pulley radii cannot be changed after
assembly without fabricating new parts and reassembling the manipulator.

  2. 2.

Goal reaching via online hardware updates: In this experiment, we optimize the
hardware design in two stages: we first optimize design parameters that
require re-fabrication (i.e., pulley radii) to a specific goal. Then, we fix
all parameters except the pre-tension elongation and co-optimize for a
different goal.

  3. 3.

Goal reaching via multi-goal control policies: We optimize all parameters to
learn a shared design and control that can achieve multiple goals with a
single control policy. In this case, we do not need to adapt any of the
hardware parameters to reach multiple goals after the policy is trained.

The first experiment is our design’s most limited application due to the need
for re-fabrication. However, it demonstrates that our transmission can be
optimized to reach different areas of the fully actuated workspace. We also
train control policies with fixed initial design parameters to compare task
performance. Fig.5 shows the goal distribution for the first experiment set.
The second set of experiments leverages a key aspect of our design - the
ability to easily adapt some design parameters without having to go through
tedious re-design, fabrication, and assembly. For this second set, we keep the
same pulley radii and only optimize the pretension elongation for a goal
different from the first experiment. In the first and second experiments, we
use the state space described in Section.III for inputs of our control policy.

The last experiment demonstrates that our robot, while being underactuated,
can achieve different goals with the same design parameters and control
policy—if the desired goals are not too far away. The state space is extended
from the first two experiment sets with the addition of the goal locations. To
establish the specific set of goals we hope to achieve, We sample goals
randomly around a center location $(-0.16,0)$ within a radius of
$50\mathrm{mm}$. For evaluation, we randomly sample $20$ goals and execute the
policy $10$ times for each goal. This experiment requires re-assembly of the
entire tentacle as we optimize all the design parameters. Therefore, we only
run this experiment in simulation.

For all experiments, We use a dense reward function:
$r=||\mathrm{x}_{e.e.}-\mathrm{x}_{goal}||^{2}+b$, where $\mathrm{x}_{e.e}$
and $\mathrm{x}_{goal}$ are the position of the end-effector and the goal, $b$
is a task completion bonus that is provided when the distance of E.E. and the
goal is less than $2\mathrm{mm}$. We initialize the hardware parameters to be
$R^{F}_{1}=[0.005,0.005,0.005]$ and $R^{F}_{2}=[0.005,0.01,0.02]$,
$L=[0.8,0.8,0.8]$ and $l^{pre}=[0,0,0]$. The units of all these parameters is
meters.

Figure 4: Average episode returns for goal reaching via re-fabrication.

###  IV-B Hardware implementation and sim-to-real transfer

Shown in Fig.1 and 6, we physically build the manipulator with the optimized
pulley radii in the first two experiment sets. To simplify both the
fabrication and assembly, we set $(R_{e})=15\mathrm{mm}$. Our transmission in
this robot is driven by two Dynamixel XM-430 servos that sense and control the
angle $(\theta_{a})$ of the actuated winches. During our experiments, we fix
the manipulator to a experimental rig, which has a mounted camera looking down
on the robot. The camera is used in lieu of joint angle encoders, as we
calculate the joint angles of our robot $(\theta_{J})$by tracking the pose of
AR tags attached to our robot as shown in Fig. 6. We use the joint angles
collected from the AR tags to execute closed-loop policies for the first two
experiment sets. In this case, we close the loop by taking observations from
the real robot during policy execution and use our optimized control policy to
determine the next action. In addition to closed-loop policies, we also simply
train a control policy to reach the goal in simulation, and then directly
transfer this policy to the real-robot without adjusting any of the actions
during runtime.

We run both the first and second experiment set on the real robot hardware for
both closed-loop and open-loop policies. In the next section, we evaluate the
accuracy and precision over 20 samples on the real robot for each set of data.

##  V Results and Analysis

###  V-A Co-optimization results

Figure 5: Qualitative results for optimizing underactuated robots to reach
different goals in simulation (goal reaching via re-fabrication). (a) shows
original design and goal locations. (b) and (c) show two optimized
underactuated tentacles reaching different goals.  Figure 6: Example of policy
rollout on the real robot. Top: real robot rollout. Bottom: comparisons
between simulated robots and real robot tracking. The real robot is shown in
purple in the bottom row of figures.

Our results from the first set of experiments show that our model learns
different design parameters for different goal locations with an error
$1.80\mathrm{mm}$. As shown in Fig.5, our model can adapt hardware parameters
to different goals. The initial design can achieve only $1$ out of $9$ goals,
achieving much lower average returns than our method (Fig.4) when we only
learn a control policy with a fixed design. Although the robot’s workspace may
contain the goal, the initial hardware parameters make learning a stable
control policy difficult. After the design parameters are optimized, the
control policy can learn a stable action sequence to arrive at the goal
position.

TABLE I: Final Distance to Goal Position of Optimized Robots for real-world reaching tasks. Open and closed represent open-loop and closed-loop execution, respectively. All units of distance in this table are $\mathrm{mm}$, and each value is calculated for 20 samples. | Avg. distance to goal | Std. distance
---|---|---
|  |
Stage | Open | Closed | Open | Closed
Pulley radii | 13.27 | 6.30 | 7.28 | 2.87
Preloads | 21.27 | 12.97 | 7.26 | 4.39
TABLE II: Quantitative results for sim-to-real hardware transfer. Average error between real hardware states and simulated states for both open-loop and closed-loop policy rollouts. Each value is calculated over a total of 40 samples. | Open-loop | Closed-loop
---|---|---
Joint angle (rad) | 0.176 | 0.0932
End-effector position (m) | 17.523 | 10.877

In the second experiment set, our results show we can optimize only the pulley
radii to achieve one goal $(0.13,0.3)$ with $1.8\mathrm{mm}$ error in
simulation. Then, we fix the pulley radii and only optimize the preloaded
extension to another goal $(0.16,-0.08)$ with $3.2\mathrm{mm}$ error in
simulation. Using derived design parameters, we build a real tentacle robot
and directly transfer the control policy from simulation to the real world.

Our results show that the transfer policy achieves, on average,
$6.3\mathrm{mm}$ error to the goal location. During training, since the
control policy is co-optimized with the non-stationary hardware proxy, it is
robust to hardware parameters perturbation. We test both open- and closed-loop
performance on the real robot. In the open-loop test case, we record action
sequences executed from the simulation and directly execute them on the real
robot. The robot can touch the goal with $12.97\mathrm{mm}$ errors.

###  V-B Sim-to-real accuracy

A key focus in our paper is the ability to transfer the optimized design and
control policy to physical robotic hardware. We achieve accurate transfer in
terms of task performance, as shown in Table I. Our closed-loop results show
that taking observations directly from the real robot improves task
performance substantially. The closed-loop policies’ standard deviation of the
distance to the goal is much lower, which indicates that they yield more
consistent results in the real world. Additionally, the average final distance
to the goal for the closed-loop policies is much lower, indicating higher
accuracy. Our policy is robust to actuation errors by reasoning about the
current robot state and producing actions that correct its errors. In Table
II, we calculate the average difference between the real robot state and the
simulated robot state for both the closed-loop and open-loop sets of
experiments. The average error in joint space for closed-loop policy execution
is about half the error for open-loop execution. Similarly, the average
distance between the real and simulated end-effector positions for the entire
action sequence of each experiment is much lower for closed-loop policies.

###  V-C Energy manifold optimization

Figure 7: Manifold comparisons before (a) and after (b) hardware optimization.
Each grid represents the constraint satisfaction (yellow for constraints
satisfied, dark purple for not satisfied) of its center. We set the joint
angles of $\\#2$ to be 0 for visualization purposes.

Our forward model is an optimization-based model. Each forward step requires
many optimization steps to find a global energy minimum, which is crucial for
accurately simulating our robot. Each set of hardware parameters, along with a
control action, corresponds to a manifold in the energy landscape. Although
the global energy landscape does not have any local optima, a manifold can be
discontinuous and has multiple local minimums. This introduces difficulties
for optimization and constrains our choice of optimization algorithm to global
optimization (e.g., genetic algorithms). However, given the time budget,
global optimization can be inefficient and slow down our simulated robot since
it generally searches in a large state space.

As mentioned in Section III-B, when we optimize the control and hardware
parameters of this robot, we optimize the manifold for energy optimization.
Our experiments show that our co-optimization process results in an energy
manifold that is easy to learn by task. As shown in Fig.7, the original
hardware design parameters provide manifolds that are discontinuous, spread in
different regions in the energy landscape, and have several local minimums.
After MORPH, the resulting manifold (see Fig.7b) is a continuous and smooth
manifold. To further analyze the optimized manifold’s effect, we apply
Sequential Least Squares Programming optimizer (SLSQP) [31], a local
optimization algorithm with a time budget of $500$ optimization steps and
$100$ random actions, to both the unoptimized and optimized manifolds and
compare their results to those of a global search algorithm implemented in-
house. When the hardware is unoptimized, given the same time budget, the
simulated results have a much higher error ($0.42$) than the optimized
hardware design ($0.18$). This means that the resulting manifold of our
optimization is more suitable for fast simulation with local optimization.

##  VI Conclusion

In this work, we present a method for co-optimizing the design and control of
underactuated kinematic chains. The key to our approach is an optimization-
based forward actuation model that effectively captures the behavior of our
robot design, and a co-optimization pipeline is capable of learning with non-
differentiable physics. Our experimental results from simulation and the real
world demonstrate that our design, along with the flexibility provided by
hardware optimization, results in flexible robot capabilities while enjoying
the benefits of underactuation. A key limitation of our current work is the
task complexity. While being general, our forward actuation model assumes
quasi-static, making contact-rich tasks difficult. In future works, we aim to
extend our work to more complex tasks. Another future direction is to further
utilize our online adaptable design and explore novel mechanisms to make part
of the design adaptable.

## References

  * [1] Shigeo Hirose and Yoji Umetani  “The development of soft gripper for the versatile robot hand”  In _Mechanism and Machine Theory_ 13.3, 1978, pp. 351–359  DOI: https://doi.org/10.1016/0094-114X(78)90059-9
  * [2] Matei Ciocarlie and Peter Allen  “Data-driven optimization for underactuated robotic hands”  In _2010 IEEE International Conference on Robotics and Automation_ , 2010, pp. 1292–1299  DOI: 10.1109/ROBOT.2010.5509793
  * [3] Matei Ciocarlie, Fernando Mier Hicks and Scott Stanford  “Kinetic and dimensional optimization for a tendon-driven gripper”  In _2013 IEEE International Conference on Robotics and Automation_ , 2013, pp. 2751–2758  DOI: 10.1109/ICRA.2013.6630956
  * [4] Angus B. Clark, Lois Liow and Nicolas Rojas  “Force Evaluation of Tendon Routing for Underactuated Grasping”  In _Journal of Mechanical Design_ 143.10, 2021, pp. 104502  DOI: 10.1115/1.4050382
  * [5] Priyanka Rao, Quentin Peyron, Sven Lilge and Jessica Burgner-Kahrs  “How to Model Tendon-Driven Continuum Robots and Benchmark Modelling Performance”  In _Frontiers in Robotics and AI_ 7, 2021  DOI: 10.3389/frobt.2020.630245
  * [6] Matteo Russo, Seyed Mohammad Hadi Sadati, Xin Dong, Abdelkhalick Mohammad, Ian D. Walker, Christos Bergeles, Kai Xu and Dragos A. Axinte  “Continuum Robots: An Overview”  In _Advanced Intelligent Systems_ 5.5, 2023, pp. 2200367  DOI: https://doi.org/10.1002/aisy.202200367
  * [7] Bowen Yi, Yeman Fan, Dikai Liu and Jose Guadalupe Romero  “Simultaneous Position-and-Stiffness Control of Underactuated Antagonistic Tendon-Driven Continuum Robots”, 2023  arXiv:2306.03865 [cs.RO]
  * [8] Raymond R. Ma and Aaron M. Dollar  “Linkage-Based Analysis and Optimization of an Underactuated Planar Manipulator for In-Hand Manipulation”  In _Journal of Mechanisms and Robotics_ 6.1, 2013, pp. 011002  DOI: 10.1115/1.4025620
  * [9] Raphael Deimel, Patrick Irmisch, Vincent Wall and Oliver Brock  “Automated co-design of soft hand morphology and control strategy for grasping”  In _2017 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2017, pp. 1213–1218  URL: https://api.semanticscholar.org/CorpusID:3037063
  * [10] Luca Barbazza, Damiano Zanotto, Giulio Rosati and Sunil K. Agrawal  “Design and Optimal Control of an Underactuated Cable-Driven Micro–Macro Robot”  In _IEEE Robotics and Automation Letters_ 2.2, 2017, pp. 896–903  DOI: 10.1109/LRA.2017.2651941
  * [11] Jie Xu, Tao Chen, Lara Zlokapa, Michael Foshey, Wojciech Matusik, Shinjiro Sueda and Pulkit Agrawal  “An End-to-End Differentiable Framework for Contact-Aware Robot Design”  In _ArXiv_ abs/2107.07501, 2021  URL: https://api.semanticscholar.org/CorpusID:235651322
  * [12] Tianjian Chen, Zhanpeng He and Matei Ciocarlie  “Hardware as Policy: Mechanical and ComputationalCo-Optimization using Deep Reinforcement Learning”  In _Conference on Robotic Learning (CoRL)_ , 2020
  * [13] Zhanpeng He and Matei T. Ciocarlie  “MORPH: Design Co-optimization with Reinforcement Learning via a Differentiable Hardware Model Proxy”  In _ArXiv_ abs/2309.17227, 2023  URL: https://api.semanticscholar.org/CorpusID:263310528
  * [14] Zhong-Ping Jiang  “Controlling Underactuated Mechanical Systems: A Review and Open Problems”  In _Advances in the Theory of Control, Signals and Systems with Physical Modeling_ Berlin, Heidelberg: Springer Berlin Heidelberg, 2011, pp. 77–88  DOI: 10.1007/978-3-642-16135-3˙7
  * [15] Ryuta Ozawa, Kazunori Hashirii and Hiroaki Kobayashi  “Design and control of underactuated tendon-driven mechanisms”  In _2009 IEEE International Conference on Robotics and Automation_ , 2009, pp. 1522–1527  DOI: 10.1109/ROBOT.2009.5152222
  * [16] Aaron M. Dollar and Robert D. Howe  “The SDM Hand: A Highly Adaptive Compliant Grasper for Unstructured Environments”  In _Experimental Robotics_ Berlin, Heidelberg: Springer Berlin Heidelberg, 2009, pp. 3–11
  * [17] Yanyu Su, Yan Wu, Kyuhwa Lee, Zhijiang Du and Yiannis Demiris  “Robust grasping for an under-actuated anthropomorphic hand under object position uncertainty”  In _2012 12th IEEE-RAS International Conference on Humanoid Robots (Humanoids 2012)_ , 2012, pp. 719–725  DOI: 10.1109/HUMANOIDS.2012.6651599
  * [18] Giorgio Grioli, Manuel Catalano, Emanuele Silvestro, Simone Tono and Antonio Bicchi  “Adaptive synergies: An approach to the design of under-actuated robotic hands”  In _2012 IEEE/RSJ International Conference on Intelligent Robots and Systems_ , 2012, pp. 1251–1256  DOI: 10.1109/IROS.2012.6385881
  * [19] Cosimo Della Santina, Cristina Piazza, Giorgio Grioli, Manuel G. Catalano and Antonio Bicchi  “Toward Dexterous Manipulation With Augmented Adaptive Synergies: The Pisa/IIT SoftHand 2”  In _IEEE Transactions on Robotics_ 34.5, 2018, pp. 1141–1156  DOI: 10.1109/TRO.2018.2830407
  * [20] Tianjian Chen, Tianyi Zhang and Matei Ciocarlie  “Design Paradigms Based on Spring Agonists for Underactuated Robot Hands: Concepts and Application”  In _2021 IEEE International Conference on Robotics and Automation (ICRA)_ , 2021, pp. 7100–7106  DOI: 10.1109/ICRA48506.2021.9561832
  * [21] L. Udawatta, K. Watanabe, K. Izumi and K. Kiguchi  “Control of underactuated robot manipulators using switching computed torque method: GA based approach”  In _Soft Computing_ 8.1 Springer Verlag, 2003, pp. 51–60  DOI: 10.1007/s00500-002-0257-8
  * [22] Lael U. Odhner and Aaron M. Dollar  “Dexterous manipulation with underactuated elastic hands”  In _2011 IEEE International Conference on Robotics and Automation_ , 2011, pp. 5254–5260  DOI: 10.1109/ICRA.2011.5980263
  * [23] Andrew S. Morgan, Kaiyu Hang, Bowen Wen, Kostas Bekris and Aaron M. Dollar  “Complex In-Hand Manipulation Via Compliance-Enabled Finger Gaiting and Multi-Modal Planning”  In _IEEE Robotics and Automation Letters_ 7.2, 2022, pp. 4821–4828  DOI: 10.1109/LRA.2022.3145961
  * [24] Lucy Jackson, Celyn Walters, Steve Eckersley, Pete Senior and Simon Hadfield  “ORCHID: Optimisation of Robotic Control and Hardware In Design using Reinforcement Learning”  In _2021 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2021, pp. 4911–4917  DOI: 10.1109/IROS51168.2021.9635865
  * [25] David R Ha  “Reinforcement Learning for Improving Agent Design”  In _Artificial Life_ 25, 2018, pp. 352–365  URL: https://api.semanticscholar.org/CorpusID:52945447
  * [26] Charles Schaff, David Yunis, Ayan Chakrabarti and Matthew R Walter  “Jointly learning to construct and control agents using deep reinforcement learning”  In _IEEE Intl. Conf. on Robotics and Automation_ , 2019, pp. 9798–9805  IEEE
  * [27] Tingwu Wang, Yuhao Zhou, Sanja Fidler and Jimmy Ba  “Neural Graph Evolution: Automatic Robot Design”  In _International Conference on Learning Representations_ , 2019  URL: https://openreview.net/forum?id=BkgWHnR5tm
  * [28] Kevin Sebastian Luck, Heni Ben Amor and Roberto Calandra  “Data-efficient Co-Adaptation of Morphology and Behaviour with Deep Reinforcement Learning”  In _Conference on Robot Learning_ , 2019  URL: https://api.semanticscholar.org/CorpusID:208139268
  * [29] Ye Yuan, Yuda Song, Zhengyi Luo, Wen Sun and Kris M. Kitani  “Transform2Act: Learning a Transform-and-Control Policy for Efficient Agent Design”  In _ArXiv_ abs/2110.03659, 2021  URL: https://api.semanticscholar.org/CorpusID:238419578
  * [30] Nikolaus Hansen and Andreas Ostermeier  “Completely derandomized self-adaptation in evolution strategies”  In _Evolutionary computation_ 9.2, 2001, pp. 159–195
  * [31] D. Kraft  “A Software Package for Sequential Quadratic Programming”, Deutsche Forschungs- und Versuchsanstalt für Luft- und Raumfahrt Köln: Forschungsbericht  Wiss. Berichtswesen d. DFVLR, 1988  URL: https://books.google.com/books?id=4rKaGwAACAAJ