# Nonadiabatic Correction and Adiabatic Criteria
of Noninteracting Quantum Dot Systems

Masahiro Hasegawa1 and Takeo Kato1<EMAIL_ADDRESS>for Solid
State Physics1Institute for Solid State Physics The University of Tokyo The
University of Tokyo Kashiwa Kashiwa Chiba 277-8581 Chiba 277-8581 Japan Japan

###### Abstract

We theoretically study nonadiabatic corrections for charge pumping in a
noninteracting electron model of a single-level quantum dot. We derive a
formula for the velocity limit of parameter driving to realize adiabatic
pumping and illustrate its features in the wide-band limit. We discuss the
effect of the nonadiabatic corrections in terms of a typical velocity limit,
which is defined by averaging the velocity limits on the parameter contour. We
also show that the typical velocity limit vanishes when the adiabaticity
breakdown occurs in a quantum dot coupled to reservoirs with a band edge.

##  1 Introduction

Adiabatic pumping, a transport phenomenon induced by quasi-static periodic
driving of model parameters, has been investigated for a long time in various
physical systems, since its concept was first proposed by Thouless [1]. In the
field of mesoscopic physics, electron transport induced by adiabatic pumping
has been studied in a number of theoretical and experimental studies [2, 3, 4,
5, 6]. Charge transferred per cycle of adiabatic pumping can be expressed in
terms of the Berry curvature for non-interacting quantum systems [7], and a
similar geometrical framework has been generalized to interacting electron
systems [8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24].
Adiabatic pumping has also been examined in the context of quasi-static
operations to characterize non-equilibrium states in nanoscale quantum heat
engines [25, 26, 27, 28, 29, 30, 31, 32, 33].

The geometric structure of adiabatic pumping is helpful for understanding non-
equilibrium properties of nanoscale quantum engines. However, in order to
justify a geometric description of adiabatic pumping, its pumping protocol has
to satisfy an adiabatic (or quasi-static) condition in the target systems. In
previous theoretical work, the adiabatic condition was not studied seriously
and was roughly estimated intuitively or empirically. However, in order to
estimate the non-equilibrium properties of nanoscale engines, for example, the
upper bound of efficiency of pumping, one has to evaluate the maximum velocity
of parameter driving to satisfy the adiabatic condition because the parameter
driving should be as fast as possible to increase the pumping power. In
particular, the upper limit of the driving velocity generally depends on the
position in the parameter space, which is crucial to geometrical optimization
of the pumping power [34]. This indicates that estimation of the adiabatic
condition is not a simple problem. Therefore, to be able to understand the
non-equilibrium nature of nanoscale engines, the adiabatic condition should be
defined in a systematic way.

To this end, we theoretically derive a detailed condition for adiabaticity by
taking the dependence of the parameter driving contour into account in this
work. As a minimum setup for adiabatic pumping, we choose a non-interacting
electron system composed of a quantum dot and two reservoirs as shown in Fig.
1 and consider adiabatic pumping induced by periodic modulation of the dot-
reservoir coupling. Our formula constitutes a general procedure to calculate
the upper limit of the driving velocity as a function of parameters and can be
used to detect bottleneck regions in the parameter space at which the limit
velocity is suppressed. We examine the adiabatic condition by analyzing the
leading nonadiabatic corrections for a quantum dot system coupled to
reservoirs with a constant density of states by assuming an infinite
bandwidth. Then we examine a quantum dot coupled to reservoirs with band edges
and discuss the problem of adiabaticity breakdown [35] in terms of our
formulation of the adiabatic condition.

This paper is organized as follows. In Sect. 2, we introduce a model of a
quantum dot and formulate the transferred charge per cycle of parameter
driving in terms of Keldysh Green’s functions. In Sect. 3, we derive a
nonadiabatic correction for the transferred charge and derive an adiabatic
condition for the velocity of parameter driving. In Sect. 4, we evaluate the
adiabatic condition for a quantum dot coupled to reservoirs with a constant
density of states in the wide-band limit. In Sect. 5, we examine adiabaticity
breakdown in a quantum dot coupled to the reservoirs with band edges. We
summarize our results in Sect. 6. The appendices present details on certain
relevant equations.

Throughout this paper, we employ the unit of $\hbar=1$ and denote the electron
charge as $e$ ($<0$).

Figure 1: Schematics of our model. A quantum dot is coupled to two electron
reservoirs, L and R, whose chemical potentials are denoted by $\mu_{\rm L}$
and $\mu_{\rm R}$, respectively. The dot-reservoir couplings denoted by
$\Gamma_{\rm L}(t)$ and $\Gamma_{\rm R}(t)$ are varied periodically by
applying ac modulation to the potential barriers.

##  2 Model

###  2.1 Model Hamiltonian

We consider a noninteracting electron model of a single-level quantum dot
coupled to two electric reservoirs (see Fig. 1) described by the Hamiltonian:

$\displaystyle
H(t)=\epsilon_{d}d^{\dagger}d+\sum_{r,k}\left[\epsilon_{k}c_{rk}^{\dagger}c_{rk}+\gamma_{r}(t)(d^{\dagger}c_{rk}+c_{rk}^{\dagger}d)\right].$
(1)

Here, $d$ and $d^{\dagger}$ are the creation and annihilation operators of an
electron in the quantum dot with an energy level $\epsilon_{d}$, while
$c_{rk}$ and $c_{rk}^{\dagger}$ are those of electrons with a wavenumber $k$
and an energy $\epsilon_{k}$ in the reservoir $r=\mathrm{L},\mathrm{R}$,
respectively. The time-dependent coupling constant between the dot and the
reservoir $r$ is denoted as $\gamma_{r}(t)$. We assume that the reservoir
temperature is zero and set the reservoir electrochemical potentials as
$\mu_{r}$. Then, the distribution function of the reservoir $r$ can be
described as

$\displaystyle f_{r}(\omega)=\Theta(\mu_{r}-\omega),$ (2)

where $\Theta(\omega)$ is the Heaviside step function.

Throughout this paper, we assume that both reservoirs have the same continuous
density of states $\rho(\omega)$. In the wide-band limit, the strength of the
dot-reservoir coupling is described by a time-dependent linewidth, as

$\displaystyle\Gamma_{r}(t)=2\pi\rho|\gamma_{r}(t)|^{2},$ (3)

where $\rho$ is the density of states at the Fermi energy. In Sect. 5, where
we discuss the breakdown of adiabaticity, we need to consider the dependence
of the density of states carefully and use a time-dependent linewidth, defined
as

$\displaystyle\Gamma_{r}(\omega;t)=2\pi\rho(\omega)|\gamma_{r}(t)|^{2}.$ (4)

###  2.2 Transferred charge in one cycle

We consider charge pumping induced by a time-dependent dot-reservoir coupling.
The charge transferred from the reservoir $r$ into the quantum dot in one
pumping cycle, denoted by $\delta Q_{r}$, is calculated as

$\displaystyle\delta Q_{r}$ $\displaystyle=\int_{0}^{T}dt\,\Braket{I_{r}(t)},$
(5) $\displaystyle I_{r}(t)$
$\displaystyle=-e\frac{d}{dt}\left[\sum_{k}c_{rk}^{\dagger}(t)c_{rk}(t)\right]$
$\displaystyle=(-i)e\gamma_{r}(t)\left(d^{\dagger}(t)c_{rk}(t)-c_{rk}^{\dagger}(t)d(t)\right),$
(6)

where $T$ is the pumping period. The average of the charge current operator
can be expressed in terms of Keldysh Green’s functions as

$\displaystyle\Braket{I_{r}(t)}$ $\displaystyle=2e\int dt_{1}\
\mathrm{Re}\biggl{[}G^{R}(t,t_{1})\Sigma_{r}^{<}(t_{1},t)$
$\displaystyle\biggl{.}\hskip
28.45274pt+G^{<}(t,t_{1})\Sigma_{r}^{A}(t_{1},t)\biggr{]}.$ (7)

Here, $G$ and $\Sigma$ denote the Green’s function and the one-particle-
irreducible self-energy, respectively, and the superscripts, $R$, $A$, and $<$
indicate the retarded, advanced, and lesser components (see Appendix A).

##  3 Adiabatic pumping and nonadiabatic correction

###  3.1 Average-time approximation

To calculate the charge pumped in one cycle in the slow driving region, we
employ the adiabatic approximation, that is called the average-time
approximation in Ref. Splettstoesser2005Dec. In this approximation, the time-
dependent Hamiltonian is expanded up to the first order at frozen time $t$:

$\displaystyle H(t^{\prime})$ $\displaystyle\simeq
H(t)+(t^{\prime}-t)\frac{dH(t)}{dt}$
$\displaystyle=H(t)+(t^{\prime}-t)\sum_{r}F_{r}\frac{d\gamma_{r}(t)}{dt},$ (8)

where $F_{r}=\partial H/\partial\gamma_{r}$ is the response operator. Pumped
charge in one cycle is approximated into two terms,

$\displaystyle\delta Q_{r}$ $\displaystyle\simeq\delta Q_{r}^{(0)}+\delta
Q_{r}^{(1)},$ (9) $\displaystyle\delta Q_{r}^{(a)}$
$\displaystyle=\int_{0}^{T}dt\ \Braket{I_{r}^{(a)}(t)}_{t},\quad(a=0,1),$ (10)

where $\Braket{\cdot}_{t}$ denotes the steady-state ensemble average described
by the Hamiltonian at frozen time $t$. The current operators, $I_{r}^{(0)}$
and $I_{r}^{(1)}$, are calculated as

$\displaystyle I_{r}^{(0)}(t)$ $\displaystyle=I_{r}(t),$ (11) $\displaystyle
I_{r}^{(1)}(t)$ $\displaystyle=\sum_{r_{1}}\int dt_{1}\
\mathcal{T}\left[I_{r}(t)F_{r_{1}}(t_{1})\right](t_{1}-t)\frac{d\gamma_{r_{1}}(t)}{dt}.$
(12)

Here, $\mathcal{T}$ is the time-ordering operator. The former term, $\delta
Q_{r}^{(0)}$, is referred to as the steady-state term, because it is just the
integral of the steady-state charge current in one cycle. The latter term,
$\delta Q_{r}^{(1)}$, is referred to as the adiabatic pumping term, because it
is independent of the pumping period $T$, while the steady state contribution
is of order $O(T)$. This term can be interpreted in a geometrical manner by
substituting Eq. (12) into Eq. (10) as

$\displaystyle\delta Q^{(1)}_{r}$ $\displaystyle=\int_{0}^{T}dt\
\sum_{r_{1}}A_{r,r_{1}}(t)\frac{d\gamma_{r_{1}}(t)}{dt},$ (13) $\displaystyle
A_{r,r_{1}}(t)$ $\displaystyle=\int dt_{1}\,\langle{\cal
T}[I_{r}(t)F_{r_{1}}(t_{1})]\rangle(t_{1}-t).$ (14)

###  3.2 Adiabatic expansion

In obtaining the pumped charge for adiabatic pumping in the previous
subsection, the higher-order corrections, such as, $\frac{d^{2}H(t)}{dt^{2}}$
and $\left(\frac{dH(t)}{dt}\right)^{2}$, were neglected. In other words,
adiabatic charge pumping is realized under the condition that these higher-
order corrections are small enough compared to the adiabatic pumping term.
Here, in order to see the details of the condition for adiabatic pumping, we
derive the higher-order corrections explicitly.

The pumped charge can be generally written as a series expansion with respect
to the inverse of the pumping period, $1/T$:

$\displaystyle\delta Q_{r}$ $\displaystyle=\sum_{n=0}^{\infty}\delta
Q_{r}^{(n)}$ (15) $\displaystyle\delta Q_{r}^{(n)}$
$\displaystyle=\int_{0}^{T}dt\ \Braket{I_{r}^{(n)}(t)}_{t}$ (16)

Here, $\delta Q^{(n)}$ is of order $O((1/T)^{n-1})$ and $I^{(n)}_{r}(t)$ is
referred to as the $n$th-order current operator. The zeroth- and first-order
current operators correspond to Eqs. (11) and (12), respectively.

The second-order current operator is calculated as

$\displaystyle I^{(2)}_{r}(t)=I^{(2,1)}_{r}(t)+I^{(2,2)}_{r}(t),$ (17)
$\displaystyle I^{(2,1)}_{r}(t)$ $\displaystyle=\frac{1}{2}\sum_{r_{1}}\int
dt_{1}\ \mathcal{T}\left[I_{r}(t)F_{r_{1}}(t_{1})\right]$ $\displaystyle\hskip
28.45274pt\times(t_{1}-t)^{2}\frac{d^{2}\gamma_{r_{1}}(t)}{dt^{2}},$ (18)
$\displaystyle I^{(2,2)}_{r}(t)$
$\displaystyle=\frac{1}{2}\sum_{r_{1},r_{2}}\int dt_{1}dt_{2}\
\mathcal{T}\left[I_{r}(t)F_{r_{1}}(t_{1})F_{r_{2}}(t_{2})\right]$
$\displaystyle\hskip
28.45274pt\times(t_{1}-t)(t_{2}-t)\frac{d\gamma_{r_{1}}(t)}{dt}\frac{d\gamma_{r_{2}}(t)}{dt}.$
(19)

Substituting Eqs. (17)-(19) into Eq. (16), the second-order correction can be
expressed in quadratic form,

$\displaystyle\delta Q_{r}^{(2)}$ $\displaystyle=\int_{0}^{T}dt\
\Biggl{(}\sum_{r_{1}}B^{(1)}_{r,r_{1}}(t)\frac{d^{2}\gamma_{r_{1}}(t)}{dt^{2}}$
$\displaystyle\hskip
28.45274pt+\sum_{r_{1},r_{2}}B^{(2)}_{r,r_{1}r_{2}}(t)\frac{d\gamma_{r_{1}}(t)}{dt}\frac{d\gamma_{r_{2}}(t)}{dt}\Biggr{)}$
(20) $\displaystyle B_{r,r_{1}}^{(1)}(t)$ $\displaystyle=\int
dt_{1}\,\langle{\cal T}[I_{r}(t)F_{r_{1}}(t_{1})]\rangle(t_{1}-t)^{2},$ (21)
$\displaystyle B_{r,r_{1}r_{2}}^{(2)}(t)$ $\displaystyle=\int
dt_{1}dt_{2}\,\langle{\cal
T}[I_{r}(t)F_{r_{1}}(t_{1})F_{r_{2}}(t_{2})]\rangle$ $\displaystyle\hskip
28.45274pt\times(t_{1}-t)(t_{2}-t).$ (22)

Using integration by parts, the second-order correction becomes

$\displaystyle\braket{I_{r}^{(2)}(t)}=\sum_{r_{1},r_{2}}B_{r,r_{1}r_{2}}(t)\frac{d\gamma_{r_{1}}(t)}{dt}\frac{d\gamma_{r_{2}}(t)}{dt}.,$
(23) $\displaystyle
B_{r,r_{1}r_{2}}(t)=B^{(2)}_{r,r_{1}r_{2}}(t)-\frac{\partial
B^{(1)}_{r,r_{1}}(t)}{\partial\gamma_{r_{2}}}.$ (24)

The third-order current operator is calculated in the same way:

$\displaystyle I^{(3)}_{r}(t)$
$\displaystyle=I^{(3,1)}_{r}(t)+I^{(3,2)}_{r}(t)+I^{(3,3)}_{r}(t),$ (25)
$\displaystyle I^{(3,1)}_{r}(t)$
$\displaystyle=\frac{1}{3!}\sum_{r_{1}}\int\\!dt_{1}\,\mathcal{T}\left[I_{r}(t)F_{r_{1}}(t_{1})\right]$
$\displaystyle\times(t_{1}-t)^{3}\frac{d^{3}\gamma_{r_{1}}(t)}{dt^{3}},$ (26)
$\displaystyle I^{(3,2)}_{r}(t)$
$\displaystyle=\frac{1}{2!}\sum_{r_{1},r_{2}}\int\\!dt_{1}\,\mathcal{T}\left[I_{r}(t)F_{r_{1}}(t_{1})F_{r_{2}}(t_{2})\right]$
$\displaystyle\times(t_{1}-t)^{2}(t_{2}-t)\frac{d^{2}\gamma_{r_{1}}(t)}{dt^{2}}\frac{d\gamma_{r_{2}}(t)}{dt},$
(27) $\displaystyle I^{(3,3)}_{r}(t)$
$\displaystyle=\frac{1}{3!}\sum_{r_{1},r_{2},r_{3}}\int\\!dt_{1}\,\mathcal{T}\left[I_{r}(t)F_{r_{1}}(t_{1})F_{r_{2}}(t_{2})F_{r_{3}}(t_{3})\right]$
$\displaystyle\times(t_{1}-t)(t_{2}-t)(t_{3}-t)\frac{d\gamma_{r_{1}}(t)}{dt}\frac{d\gamma_{r_{2}}(t)}{dt}\frac{d\gamma_{r_{3}}(t)}{dt}.$
(28)

Using integration by parts, the third-order correction can be summarized into
two terms:

$\displaystyle\braket{I_{r}^{(3)}(t)}=\braket{I_{r}^{(3v)}(t)}+\braket{I_{r,a}^{(3a)}(t)},$
(29)
$\displaystyle\braket{I_{r}^{(3v)}(t)}=\\!\\!\sum_{r_{1},r_{2},r_{3}}\\!\\!C_{r,r_{1}r_{2}r_{3}}(t)\frac{d\gamma_{r_{1}}(t)}{dt}\frac{d\gamma_{r_{2}}(t)}{dt}\frac{d\gamma_{r_{3}}(t)}{dt},$
(30)
$\displaystyle\braket{I_{r}^{(3a)}(t)}=\sum_{r_{1},r_{2}}D_{r,r_{1}r_{2}}(t)\frac{d^{2}\gamma_{r_{1}}(t)}{dt^{2}}\frac{d\gamma_{r_{2}}(t)}{dt}.$
(31)

The first term, $\langle I_{r}^{(3v)}(t)\rangle$, consists only of the first
derivative of the driving parameter, whereas the second term, $\langle
I_{r}^{(3a)}(t)\rangle$, includes the second derivative of the driving
parameter. In general, we can derive higher-order corrections by using a
similar procedure for the $n$th-order operator, defined as

$\displaystyle I^{(n)}_{r}(t)=\int dt_{1}\cdots dt_{n}\
\mathcal{T}\left[I_{r}(t)F_{n}(t_{1},\cdots,t_{n})\right],$ (32)
$\displaystyle
F_{n}(t_{1},\cdots,t_{n})=\prod_{m=1}^{n}\sum_{l_{m}}\sum_{r^{\prime}}\frac{\partial
H}{\partial\gamma_{r^{\prime}}}(t_{m})\left(\frac{d}{dt}\right)^{l_{m}}\gamma_{r^{\prime}}(t),$
(33)

though their expressions become much more complex.

###  3.3 Adiabatic condition

We define the driving velocity of the parameters as

$\displaystyle{\bm{v}}(t)=(v_{\rm L}(t),v_{\rm
R}(t))=\biggl{(}\frac{d\gamma_{\rm L}}{dt},\frac{d\gamma_{\rm
R}}{dt}\biggr{)}.$ (34)

To write down the adiabatic condition, it is convenient to introduce the
normalized velocity, defined as

$\displaystyle\tilde{\bm{v}}(t)=(\tilde{v}_{\rm L}(t),\tilde{v}_{\rm
R}(t))=\frac{{\bm{v}}(t)}{|{\bm{v}}(t)|}.$ (35)

Let us discuss how the driving velocity is restricted in order to realize
adiabatic pumping. Charge pumping can be regarded as ‘adiabatic’ when the
higher-order (nonadiabatic) corrections of the charge current are much smaller
than the adiabatic current:

$\displaystyle|\langle{I_{r}^{(2)}(t)}\rangle|,|\langle{I_{r}^{(3v)}(t)}\rangle|,|\langle{I_{r}^{(3a)}(t)}\rangle|,\cdots\ll|\langle
I_{r}^{(1)}(t)\rangle|.$ (36)

Note that this condition is always satisfied in the slow-driving limit of
${\cal T}\rightarrow\infty$ since $\langle I_{r}^{(n)}(t)\rangle$ is
proportional to $(1/T)^{n-1}$. For a finite speed of driving, we introduce a
small tolerance factor $\eta$ ($\ll 1$) and describe the extent to which the
ratio of the higher-order current corrections to the adiabatic current is
allowed as follows:

$\displaystyle|\langle I_{r}^{(2)}(t)\rangle|\leq\eta|\langle
I_{r}^{(1)}(t)\rangle|,$ (37) $\displaystyle|\langle
I_{r}^{(3v)}(t)\rangle|\leq\eta|\langle I_{r}^{(1)}(t)\rangle|,$ (38)
$\displaystyle|\langle I_{r}^{(3a)}(t)\rangle|\leq\eta|\langle
I_{r}^{(1)}(t)\rangle|.$ (39)

We call these the adiabatic conditions.

The adiabatic condition for the second-order correction, Eq. (37), can be
rewritten using Eqs. (13) and (23) as

$\displaystyle|v_{r}(t)|\leq v^{(2)}_{{\rm lim},r}(t),$ (40)

where the ‘velocity limit’ $v^{(2)}_{{\rm lim},r}(t)$ is given as

$\displaystyle v_{{\rm lim},r}^{(2)}(t)=\eta\frac{|A_{r}(t)|}{|B_{r}(t)|},$
(41) $\displaystyle A_{r}(t)=\sum_{r_{1}}A_{r,r_{1}}(t)\tilde{v}_{r_{1}}(t),$
(42) $\displaystyle
B_{r}(t)=\sum_{r_{1},r_{2}}B_{r,r_{1},r_{2}}(t)\tilde{v}_{r_{1}}(t)\tilde{v}_{r_{2}}(t).$
(43)

In a similar way, the adiabatic condition for the third-order correction, Eq.
(38), can be rewritten as

$\displaystyle|v_{r}(t)|\leq v_{\mathrm{lim},r}^{(3)}(t),$ (44) $\displaystyle
v_{\mathrm{lim},r}^{(3)}(t)=\sqrt{\eta\frac{|A_{r}(t)|}{|C_{r}(t)|}},$ (45)
$\displaystyle
C_{r}(t)=\sum_{r_{1},r_{2},r_{3}}C_{r,r_{1}r_{2}r_{3}}(t)\tilde{v}_{r_{1}}(t)\tilde{v}_{r_{2}}(t)\tilde{v}_{r_{3}}(t).$
(46)

The last condition, Eq. (39), imposes a restriction on the acceleration of the
pumping parameter (see Appendix B). For simplicity, we will not discuss this
condition and instead assume that the driving velocity is nearly constant.

Combining these conditions, we obtain

$\displaystyle
v(t)\leq\mathrm{min}\left[v_{\mathrm{lim},r}^{(2)}(t),v_{\mathrm{lim},r}^{(3)}(t),\cdots\right].$
(47)

One might guess that, as long as one considers slow driving, the strictest
condition is imposed by a second-order correction, since the third- and
higher-order corrections are much smaller than the second-order one in the
adiabatic expansion. Although this is true in most cases, there is an
exception: When the system is in equilibrium, that is, when all the reservoirs
have the same temperature and chemical potential, the second-order correction
is always zero for any driving of $\gamma_{r}(t)$ (this exception will be
discussed in the next subsection). In this case, the velocity limit
$v^{(2)}_{{\rm lim},r}(t)$ becomes infinite, and there is no restriction on
the pumping velocity. Therefore, the third correction should be included in
the estimate of the velocity limit.

###  3.4 Disappearance of the second-order correction

Comparing the charge transfer induced by a driving operation and its time-
reversed operation, one can see that the $2n$-th-order correction disappears
when the system is in equilibrium. Let us consider a driving protocol $C$
defined in the parameter space, such as,

$\displaystyle C=\\{(\gamma_{L}(t),\gamma_{R}(t))|t\in[0,T]\\},$ (48)

and its time-reversed protocol $\bar{C}$, such as,

$\displaystyle\bar{C}=\\{(\gamma_{L}(T-t),\gamma_{R}(T-t))|t\in[0,T]\\}.$ (49)

In general, the amount of charge transfer induced by $C$ and $\bar{C}$ can
vary, because steady charge current flows between the reservoirs. This
indicates that when the system is in equilibrium, that is, when the
temperature and chemical potential are the same among all the reservoirs,
there is no steady charge current and the amount of charge transfer induced by
$C$ and $\bar{C}$ should satisfy the following relation:

$\displaystyle\delta Q_{r}[C]=-\delta Q_{r}[\bar{C}].$ (50)

Here, $\delta Q[C]$ denotes the amount of charge transfer under the driving
protocol $C$. Let us consider the adiabatic expansion of $\delta Q[C]$ and
$\delta Q[\bar{C}]$.

$\displaystyle\delta Q_{r}[C]=\sum_{n=0}\delta Q^{(n)}_{r}[C],$ (51)
$\displaystyle\delta Q_{r}[\bar{C}]=\sum_{n=0}\delta Q^{(n)}_{r}[\bar{C}]$
(52)

Focusing on the first-order corrections (adiabatic pumping terms), they
satisfy an anti-symmetric relation,

$\displaystyle\delta Q^{(1)}_{r}[\bar{C}]$ $\displaystyle=\int_{0}^{T}dt\
\sum_{r_{1}}A_{r,r_{1}}(T-t)\frac{d\gamma_{r_{1}}(T-t)}{dt}$
$\displaystyle=-\int_{0}^{T}dt\
\sum_{r_{1}}A_{r,r_{1}}(t)\frac{d\gamma_{r_{1}}(t)}{dt}$
$\displaystyle=-\delta Q^{(1)}[C],$ (53)

while the second-order corrections satisfy a symmetric relation,

$\displaystyle\delta Q^{(2)}_{r}[\bar{C}]$
$\displaystyle=\int_{0}^{T}dt\sum_{r_{1},r_{2}}B_{r,r_{1}r_{2}}(T-t)\frac{d\gamma_{r_{1}}(T-t)}{dt}\frac{d\gamma_{r_{2}}(T-t)}{dt}$
$\displaystyle=\int_{0}^{T}dt\sum_{r_{1},r_{2}}B_{r,r_{1}r_{2}}(t)\frac{d\gamma_{r_{1}}(t)}{dt}\frac{d\gamma_{r_{2}}(t)}{dt}=\delta
Q^{(2)}_{r}[C].$ (54)

In general, the $(2n)$-th-order corrections and $(2n+1)$-th-order corrections
satisfy symmetric and anti-symmetric relations, respectively,

$\displaystyle\delta Q^{(2n)}[C]$ $\displaystyle=\delta Q^{(2n)}[\bar{C}],$
(55) $\displaystyle\delta Q^{(2n+1)}[C]$ $\displaystyle=-\delta
Q^{(2n+1)}[\bar{C}].$ (56)

Combining Eq. (50) and Eq. (56), we can conclude that the $(2n)$-th-order
correction should be zero for any driving protocols when the system is in
equilibrium. In this case, the nonadiabatic correction starts with the third-
order correction $\delta Q^{(3)}$, so the velocity limit is imposed by the
third-order correction, not the second-order correction.

###  3.5 Typical velocity limit

As shown in Eqs. (41) and (45), the velocity limit may drop to zero under
certain conditions, which is known as adiabaticity breakdown. There are two
cases in which the velocity limit is zero: (i) the adiabatic pumping term
(numerator of Eqs. (41) and (45)) is zero; (ii) the nonadiabatic corrections
(denominator of Eqs. (41) and (45)) are infinitely large. The former case is
not serious for adiabaticity, because the nonadiabatic corrections are finite
and one can pass through such adiabaticity breaking points in finite velocity
as long as the nonadiabatic corrections are negligible compared to the
adiabatic pumping term in the whole cycle. We refer to this case as ‘non-
serious adiabaticity breakdown’. On the other hand, the latter case is
serious, because the nonadiabatic corrections diverge, so one should not pass
through such points at any finite velocity. We refer to this case as ‘serious
adiabaticity breakdown’.

Unfortunately, it is difficult to distinguish serious and non-serious
adiabaticity breakdown only by the velocity limit itself, and it is necessary
to define an indicator that can distinguish them. Here, we define a ‘typical’
velocity limit as such an indicator,

$\displaystyle\bar{v}_{\mathrm{lim},r}^{(2/3)}=\int_{0}^{T}dt\
v_{\mathrm{lim},r}^{(2/3)}(t)p^{(2/3)}(t),$ (57)

where $p^{(2)}(t)$ and $p^{(3)}(t)$ are weight functions describing how large
the nonadiabatic corrections are on the driving contour:

$\displaystyle
p^{(2)}(t)=\frac{\left|\sum_{r_{1},r_{2}}B_{r,r_{1}r_{2}}(t)\tilde{v}_{r_{1}}(t)\tilde{v}_{r_{2}}(t)\right|}{\int_{0}^{T}dt\left|\sum_{r_{1},r_{2}}B_{r,r_{1}r_{2}}(t)\tilde{v}_{r_{1}}(t)\tilde{v}_{r_{2}}(t)\right|},$
(58) $\displaystyle
p^{(3)}(t)=\frac{\left|\sum_{r_{1},r_{2},r_{3}}C_{r,r_{1}r_{2}r_{3}}(t)\tilde{v}_{r_{1}}(t)\tilde{v}_{r_{2}}(t)\tilde{v}_{r_{3}}(t)\right|}{\int_{0}^{T}dt\left|\sum_{r_{1},r_{2},r_{3}}C_{r,r_{1}r_{2}r_{3}}(t)\tilde{v}_{r_{1}}(t)\tilde{v}_{r_{2}}(t)\tilde{v}_{r_{3}}(t)\right|}.$
(59)

The typical velocity limit becomes zero for the serious adiabaticity
breakdown, while it remains finite for the non-serious adiabaticity breakdown.

##  4 Velocity limit in wide-band limit

Here, we present the velocity limit calculated by the second-order correction
and the third-order correction in the wide-band limit. First, we consider the
quantum dot system with a finite electrochemical potential bias and a square
driving contour, as illustrated in Fig. 2. Next, we estimate the
electrochemical potential bias dependence for a typical velocity limit.

Figure 2: Square contour.

###  4.1 Velocity limit with finite electrochemical potential bias

Figure 3: Velocity limit on the line (a) $A\to B$ and (b) $B\to C$. A finite
electrochemical potential bias is considered, $\mu_{L}=-\mu_{R}=\mu$, and
$\mu$ is taken to be the unit of energy. The rest of the parameters are set as
follows: $\epsilon_{d}=0$, $\Gamma_{\mathrm{min}}=0$,
$\Gamma_{\mathrm{max}}/\mu=1.0$, $\eta=0.01$. The vertical axis denotes the
velocity limit in units of $\mu^{2}$. $v_{\mathrm{lim,L}}^{(2)}$ and
$v_{\mathrm{lim,R}}^{(2)}$ are plotted ten times larger in order to draw the
lines on the same graph.

Figure 3 shows four velocity limits, $v_{\mathrm{lim,L}}^{(2)}$,
$v_{\mathrm{lim,R}}^{(2)}$, $v_{\mathrm{lim,L}}^{(3)}$, and
$v_{\mathrm{lim,R}}^{(3)}$ on the contour shown in Fig. 2. The limits are for
the symmetric bias case, with $\mu_{L}=-\mu_{R}=\mu>0$, and $\mu$ is taken to
be the unit of energy. The parameters of the contour, $\Gamma_{\mathrm{min}}$
and $\Gamma_{\mathrm{max}}$, are taken as 0 and $\mu$, respectively (see Fig.
2). The velocity limits are plotted for (a) line $A\to B$ and (b) line $B\to
C$. The velocity limits on lines $C\to D$ and $D\to A$ are not presented
because they can be obtained by replacing the reservoir indexes
($\mathrm{L}\leftrightarrow\mathrm{R}$) in Fig. 3. The horizontal axis denotes
a dimensionless parameter $s\in[0,1]$,

$\displaystyle(\Gamma_{L},\Gamma_{R})$
$\displaystyle=\left\\{\begin{array}[]{cc}s(\Gamma_{\mathrm{max}},\Gamma_{\mathrm{min}})+(1-s)(\Gamma_{\mathrm{min}},\Gamma_{\mathrm{min}}),&(A\to
B),\\\
s(\Gamma_{\mathrm{max}},\Gamma_{\mathrm{max}})+(1-s)(\Gamma_{\mathrm{max}},\Gamma_{\mathrm{min}}),&(B\to
C),\end{array}\right.$ (62)

Only $v_{\mathrm{lim,L}}^{(3)}$ is plotted in Fig. 3 (a). On this line, we
obtain $v_{\mathrm{lim,R}}^{(2)},v_{\mathrm{lim,R}}^{(3)}=\infty$ because
reservoir R is always disconnected. Also, when only one reservoir is connected
to the dot, $v_{\mathrm{lim,L}}^{(2)}$ is infinite due to the disappearance of
the second-order correction discussed in Sect. 3.4. At the origin ($s=0$), the
velocity limit drops to zero. This is because the lifetime of the electron
state in the quantum dot diverges and violates adiabaticity as the quantum dot
becomes isolated from the reservoirs. The peak at $s\simeq 0.4$ reflects the
local maximum of the adiabatic pumping term $A_{\mathrm{L,L}}$.

In Fig. 3 (b), reservoir R is connected to the dot, so all
$v_{\mathrm{lim,L}}^{(2)}$, $v_{\mathrm{lim,R}}^{(2)}$,
$v_{\mathrm{lim,L}}^{(3)}$, and $v_{\mathrm{lim,R}}^{(3)}$ are plotted. As
shown, $v_{\mathrm{lim,L}}^{(2)}$ and $v_{\mathrm{lim,R}}^{(2)}$ are
respectively smaller than $v_{\mathrm{lim,L}}^{(3)}$ and
$v_{\mathrm{lim,R}}^{(3)}$. At the origin ($s=0$), all the velocity limits
diverge for the following reason: around the origin, reservoir R is almost
disconnected from the dot while reservoir L is strongly connected so that a
small change in $\gamma_{R}$ does not affect charge transport at all. This
indicates that both the adiabatic pumping term and the other nonadiabatic
corrections are almost zero. The divergence occurs because the nonadiabatic
corrections more rapidly decrease to zero than the adiabatic pumping term.

As shown in Fig. 3 (b), $v_{\mathrm{lim,R}}^{(2)}$ and
$v_{\mathrm{lim,R}}^{(3)}$ drop to zero around $s\simeq 0.2$. This drop
happens just because the adiabatic pumping term accidentally becomes zero at
this point. This is different from the decrease in the velocity limit brought
about by the isolation of the quantum dot because the nonadiabatic corrections
do not diverge and remain finite. In terms of the velocity limit, the decrease
of velocity limit toward $(\Gamma_{\mathrm{L}},\Gamma_{\mathrm{R}})=(0,0)$
shown in Fig. 3 (a) should be more serious than this drop in Fig. 3 (b),
because the nonadiabatic corrections diverge in the former case.

Figure 4: Typical velocity limit calculated by (a) the second-order correction
and (b) the third-order correction. The electrochemical potential bias is
symmetric, $\mu_{\mathrm{L}}=-\mu_{\mathrm{R}}=\mu>0$. The other parameters
are set as, $\epsilon_{d}=0$ and $\eta=0.01$. All quantities are plotted in
energy units of $\mu$.

Now let us examine the typical velocity limit defined in Sect. 4.1. It is
plotted in Figure 4 (a) $\bar{v}_{\mathrm{lim,L}}^{(2)}$ and (b)
$\bar{v}_{\mathrm{lim,L}}^{(3)}$ on different square contours. We assume a
symmetric electrochemical potential bias, i.e.,
$\mu_{\mathrm{L}}=-\mu_{\mathrm{R}}=\mu>0$. Six cases are plotted:
$\Gamma_{\mathrm{min}}/\mu=0,0.1,0.2,0.3,0.4,0.5$. As shown in Figs. 4 (a) and
(b), the velocity limit decreases to zero and becomes insensitive to
$\Gamma_{\mathrm{max}}$ as the driving contour reaches the origin
$(\Gamma_{\mathrm{L}},\Gamma_{\mathrm{R}})=(0,0)$. This reflects the isolation
of the quantum dot. Around the point
$(\Gamma_{\mathrm{L}},\Gamma_{\mathrm{R}})=(0,0)$, due to the isolation
effect, the nonadiabatic corrections diverge, so that the typical velocity
limit strongly depends on the speed limit around the origin. On the other
hand, once the driving contour is far from the origin, nonadiabatic
corrections remain finite so that the typical velocity limit becomes sensitive
to $\Gamma_{\mathrm{max}}$. Also, comparing Fig. 4 (a) with Fig. 4 (b), we can
see that $\bar{v}_{\mathrm{lim,L}}^{(2)}$ is almost ten times smaller than
$\bar{v}_{\mathrm{lim,L}}^{(3)}$. This supports our conclusion that the
velocity limit is imposed by the second-order correction in most cases.

###  4.2 Electrochemical potential bias dependence of velocity limit

Figure 5: Electrochemical potential bias dependence of typical velocity
limits, $\bar{v}_{\mathrm{lim,L}}^{(2)}$ (purple) and
$\bar{v}_{\mathrm{lim,L}}^{(3)}$ (green). The red line shows the minimum value
between $\bar{v}_{\mathrm{lim,L}}^{(2)}$ and $\bar{v}_{\mathrm{lim,L}}^{(3)}$.
The electrochemical potential bias is
$\mu_{\mathrm{L}}=-\mu_{\mathrm{R}}=\mu>0$. The other parameters are set as
$\Gamma_{\mathrm{max}}/\epsilon_{d}=2.0$,
$\Gamma_{\mathrm{max}}/\epsilon_{d}=0.1$, and $\eta=0.01$. All quantities are
plotted in energy units of $\epsilon_{d}$.

The previous section showed that the velocity limit is mostly determined by
the second-order correction, but there is an exception as discussed in Sect.
3.4; i.e., when the electrochemical potential bias is zero, the second-order
correction disappears and the typical velocity limit
$\bar{v}_{\mathrm{lim,L}}^{(2)}$ diverges. However, for a finite
electrochemical potential bias, the second-order correction remains finite and
$\bar{v}_{\mathrm{lim,L}}^{(2)}$ becomes smaller than
$\bar{v}_{\mathrm{lim,L}}^{(3)}$. Thus, it is expected that the velocity limit
should be determined by $\bar{v}_{\mathrm{lim,L}}^{(3)}$ for small
electrochemical potential biases, while it should be determined by
$\bar{v}_{\mathrm{lim,L}}^{(2)}$ for large electrochemical potential biases.

To see this transition, let us examine the electrochemical potential bias
dependence of the typical velocity limit. Figure 5 shows the velocity limits,
$\bar{v}_{\mathrm{lim,L}}^{(2)}$ (purple line) and
$\bar{v}_{\mathrm{lim,L}}^{(3)}$ (green line). The red line plots
$\mathrm{min}[\bar{v}_{\mathrm{lim,L}}^{(2)},\bar{v}_{\mathrm{lim,L}}^{(3)}]$.
Let us consider a fixed square contour,
$\Gamma_{\mathrm{max}}/\epsilon_{d}=1.0$ and
$\Gamma_{\mathrm{min}}/\epsilon_{d}=0.1$, and estimate typical velocity
limits, $\bar{v}_{\mathrm{lim,L}}^{(2)}$ and $\bar{v}_{\mathrm{lim,L}}^{(3)}$,
for different electrochemical potential biases,
$\mu_{\mathrm{L}}=-\mu_{\mathrm{R}}=\mu>0$. Different from Fig. 4, we will set
the dot level $\epsilon_{d}\neq 0$ and employ it as the energy unit. This is
because we want to avoid the special point, $\epsilon_{d}=\mu=0$, where no
charge can be pumped so that all the higher-order corrections are zero.

As shown in the figure, for small biases, $\bar{v}_{\mathrm{lim,L}}^{(2)}$
diverges, while $\bar{v}_{\mathrm{lim,L}}^{(3)}$ remains finite. For large
biases, $\bar{v}_{\mathrm{lim,L}}^{(2)}$ converges to a finite value which is
sufficiently smaller than $\bar{v}_{\mathrm{lim,L}}^{(3)}$. This result
indicates that the velocity limit should be estimated by
$\bar{v}_{\mathrm{lim,L}}^{(2)}$ when the electrochemical potential bias is
finite, while it should be estimated by $\bar{v}_{\mathrm{lim,L}}^{(3)}$ when
the bias is small.

##  5 Adiabaticity breakdown in band-edge system

When the system shows critical phenomena, for example, phase transition,
transport changes drastically at the critical point and the adiabaticity
easily breaks down in this regime. This adiabaticity breakdown can be observed
as a divergence of nonadiabatic corrections. In this section, we examine this
adiabaticity breakdown in a quantum dot system with the electron reservoirs
that have a band edge. There are two cases in which adiabaticity breaks down:
(i) the quantum dot is isolated from all the reservoirs,
$\gamma_{L}=\gamma_{R}=0$, so adiabaticity trivially breaks down because the
quantum dot never relaxes; (ii) when the dot level is located near the band
edge of the electron reservoirs, a bound state with an infinite lifetime
appears in the quantum dot system and it violates adiabaticity.

To see the details, let us assume the reservois have the following density of
states:

$\displaystyle\rho(\omega)=\Theta(\omega)\rho\left(\frac{\omega}{\omega_{c}}\right)^{\alpha}\exp\left[-\frac{\omega}{\omega_{c}}\right]$
(63)

Here, we assume that the band edge is at $\omega=0$ and $\omega_{c}$ and
$\rho$ are the cutoff frequency and amplitude of the density of states,
respectively. The bound state appears when the coupling constants between the
dot and the reservoirs are larger than a critical value:

$\displaystyle\pi\rho\sum_{r}|\gamma_{r}(t)|^{2}=\frac{\epsilon_{d}}{\Gamma(\alpha)},$
(64)

where $\Gamma(\alpha)$ is the Gamma function. As a quantum state with an
infinite lifetime emerges, the system cannot relax and any finite-velocity
driving violates adiabaticity. Such a violation would appear as a divergence
of nonadiabatic corrections.

Figure 6: Coupling strength dependence of the second-order correction,
$B_{r,\mathrm{LL}}(\gamma_{\mathrm{L}},\gamma_{\mathrm{R}})$. The parameters
are set as follows: $\mu_{\mathrm{L}}/\epsilon_{d}=0$,
$\mu_{\mathrm{R}}/\epsilon_{d}=2.0$, $\omega_{c}/\epsilon_{d}=10.0$,
$\alpha=1.0$. The vertical axis is a logarithmic scale.

To see this, let us consider small $\gamma_{L}$ driving,

$\displaystyle(\gamma_{\mathrm{L}},\gamma_{\mathrm{R}})\to(\gamma_{\mathrm{L}}+\Delta\gamma,\gamma_{\mathrm{R}})$
(65)

over an elapsed time $\Delta t$. The second-order correction for this driving
is calculated as

$\displaystyle\delta
Q_{r}^{(2)}=B_{r,\mathrm{LL}}(\gamma_{\mathrm{L}},\gamma_{\mathrm{R}})\frac{(\Delta\gamma)^{2}}{\Delta
t}.$ (66)

Figure 6 presents coupling constant dependence of
$B_{r,\mathrm{LL}}(\gamma_{\mathrm{L}},\gamma_{\mathrm{R}})$. A symmetric
coupling, $\gamma_{\mathrm{L}}=\gamma_{\mathrm{R}}=\gamma$, is considered. As
shown in the figure, there are two divergences: the one at
$2\pi\rho\gamma^{2}/\epsilon_{d}=0$ corresponds to the isolation of the
quantum dot and the one at $2\pi\rho\gamma^{2}/\epsilon_{d}=1.0$ corresponds
to the emergence of the bound state. The singularity at
$2\pi\rho\gamma^{2}/\epsilon_{d}=1.0$ is much stronger than that at
$2\pi\rho\gamma^{2}/\epsilon_{d}=0$. This indicates that the emergence of the
bound state with an infinite lifetime is fatal to adiabaticity and adiabatic
operations cannot be performed around this critical point.

##  6 Summary

We studied nonadiabatic corrections in charge pumping and derived a formula
for adiabatic criteria by employing a noninteracting electron model of a
single-level quantum dot as a minimum model. The upper velocity limit for the
parameter for adiabatic pumping was derived from the $n$th-order term of the
pumped charge, which is proportional to $T^{n-1}$ where $T$ is the period of
the pumping cycle. The leading nonadiabatic correction derived from the second
term ($n=2$) imposes the strongest constraint in most cases. However, the next
leading correction derived from the third term ($n=3$) becomes important when
the two reservoirs have the same electrochemical potential and the same
temperature because the leading nonadiabatic correction vanishes in such a
situation. As a demonstration, we considered charge pumping induced by the
time-dependent dot-reservoir coupling and calculated the velocity limit in the
wide-band limit. We found that the velocity limit indeed depends on the point
in the parameter space; this feature is important for optimization to maximize
the pumping power. We also showed how the velocity limits imposed from the
second- and third-order terms behave near thermal equilibrium.

To detect adiabatic breakdown, we defined a typical velocity limit by
averaging the velocity limit on the parameter contour. We showed that the
typical velocity limit approaches zero in the wide-band limit when the dot-
reservoir coupling is weakened. We also calculated the typical velocity limit
for a quantum dot coupled to reservoirs with a band edge and showed that it
approaches zero, corresponding to adiabatic breakdown caused by the appearance
of a bound state in this model.

We expect that our findings will be helpful not only for accurate control of
adiabatic pumping in quantum devices but also for the foundation of non-
equilibrium statistics such as in the study of quantum heat engines. While our
work illustrates a general procedure for obtaining adiabatic criteria just in
a simple model, it is a straightforward exercise to extend our formulation to
more complex models. An important future problem would be to extend our work
to take the effect of Coulomb interactions between electrons into account.

M.H. acknowledges the financial support provided by the Advanced Leading
Graduate Course for Photon Science. T.K. was supported by JSPS Grants-in-Aid
for Scientific Research (No. JP24540316, JP26220711, and JP20K03831).

##  Appendix A Keldysh Green’s functions

In this appendix, we define the Keldysh Green’s functions (GFs) and the one-
particle-irreducible self-energies. The retarded and advanced GFs of an
electron in the isolated quantum dot are defined as

$\displaystyle g^{R}(t_{1},t_{2})$
$\displaystyle=-i\Theta(t_{1}-t_{2})\Braket{[d(t_{1}),d^{\dagger}(t_{2})]_{+}}_{\gamma_{r}=0},$
(67) $\displaystyle g^{A}(t_{1},t_{2})$
$\displaystyle=i\Theta(t_{2}-t_{1})\Braket{[d(t_{1}),d^{\dagger}(t_{2})]_{+}}_{\gamma_{r}=0},$
(68)

and those of electrons in the reservoir $r$ with wavenumber $k$ are defined as

$\displaystyle g_{rk}^{R}(t_{1},t_{2})$
$\displaystyle=-i\Theta(t_{1}-t_{2})\Braket{[c_{rk}(t_{1}),c_{rk}^{\dagger}(t_{2})]_{+}}_{\gamma_{r}=0},$
(69) $\displaystyle g_{rk}^{A}(t_{1},t_{2})$
$\displaystyle=i\Theta(t_{2}-t_{1})\Braket{[c_{rk}(t_{1}),c_{rk}^{\dagger}(t_{2})]_{+}}_{\gamma_{r}=0}.$
(70)

Here $[\cdot,\cdot]_{+}$ is the anti-commutator, and $\braket{\cdot}$ denotes
the ensemble average. The retarded GF of a dot electron for a finite dot-
reservoir coupling is calculated using the Dyson equation:

$\displaystyle G^{R}(t_{1},t_{2})=g^{R}(t_{1},t_{2})$ $\displaystyle\hskip
14.22636pt+\\!\int\\!dt_{3}dt_{4}\,g^{R}(t_{1},t_{3})\Sigma^{R}(t_{3},t_{4})G^{R}(t_{4},t_{2}),$
(71)

where $\Sigma^{R}(t_{1},t_{2})$ is the one-particle-irreducible self-energy
for the retarded GF, defined as

$\displaystyle\Sigma^{R}(t_{1},t_{2})$
$\displaystyle=\sum_{r}\Sigma_{r}^{R}(t_{1},t_{2}),$ (72)
$\displaystyle\Sigma_{r}^{R}(t_{1},t_{2})$
$\displaystyle=\sum_{k}\gamma_{r}(t_{1})g_{rk}^{R}(t_{1},t_{2})\gamma_{r}(t_{2}).$
(73)

The advanced GF and its self-energy are obtained from the retarded ones as
$G^{A}(t_{1},t_{2})=-G^{R}(t_{2},t_{1})$ and
$\Sigma^{A}(t_{1},t_{2})=-\Sigma^{R}(t_{2},t_{1})$, respectively. The lesser
GF is calculated from the retarded and advanced GFs:

$\displaystyle
G^{<}(t_{1},t_{2})=\\!\int\\!dt_{3}dt_{4}\,G^{R}(t_{1},t_{3})\Sigma^{<}(t_{3},t_{4})G^{A}(t_{4},t_{2}).$
(74)

where $\Sigma^{<}(t_{1},t_{2})$ is the lesser component of the self-energy
defined as

$\displaystyle\Sigma^{<}(t_{1},t_{2})$
$\displaystyle=\sum_{r}\Sigma_{r}^{<}(t_{1},t_{2}),$ (75)
$\displaystyle\Sigma_{r}^{<}(t_{1},t_{2})$
$\displaystyle=\sum_{k}\gamma_{r}(t_{1})g_{rk}^{<}(t_{1},t_{2})\gamma_{r}(t_{2}).$
(76)

##  Appendix B Adiabatic condition for the acceleration

We define the acceleration of the driving parameters and its normalized value
as

$\displaystyle{\bm{a}}(t)=(a_{\rm L}(t),a_{\rm
R}(t))=\biggl{(}\frac{d^{2}\gamma_{\rm L}(t)}{dt^{2}},\frac{d^{2}\gamma_{\rm
R}(t)}{dt^{2}}\biggr{)},$ (77) $\displaystyle\tilde{\bm{a}}(t)=(\tilde{a}_{\rm
L}(t),\tilde{a}_{\rm L}(t))=\frac{{\bm{a}}(t)}{|{\bm{a}}|},$ (78)

respectively. The adiabatic condition, Eq. (39), can be rewritten as

$\displaystyle|a_{r}(t)|\leq a_{\mathrm{lim},r}^{(3)}(t),$ (79) $\displaystyle
a_{\mathrm{lim},r}^{(3)}=\eta\frac{|A_{r}(t)|}{|D_{r}(t)|},$ (80)
$\displaystyle
D_{r}(t)=\sum_{r_{1},r_{2}}D_{r,r_{1}r_{2}}(t)\tilde{a}_{r_{1}}(t)\tilde{v}_{r_{2}}(t).$
(81)

Here $a_{\mathrm{lim},r}^{(3)}(t)$ is the limit of the acceleration of the
driving parameters. We do not consider this condition and instead assume that
the acceleration of the parameters is sufficiently small.

## References

  * [1] D. J. Thouless: Phys. Rev. B 27 (1983) 6083.
  * [2] L. P. Kouwenhoven, A. T. Johnson, N. C. van der Vaart, C. J. P. M. Harmans, and C. T. Foxon: Phys. Rev. Lett. 67 (1991) 1626.
  * [3] H. Pothier, P. Lafarge, C. Urbina, D. Esteve, and M. H. Devoret: EPL 17 (1992) 249.
  * [4] M. Büttiker, A. Prêtre, and H. Thomas: Phys. Rev. Lett. 70 (1993) 4114.
  * [5] M. Büttiker, H. Thomas, and A. Prêtre: Z. Phys. B 94 (1994) 133\.
  * [6] A. Prêtre, H. Thomas, and M. Büttiker: Phys. Rev. B 54 (1996) 8130\.
  * [7] P. W. Brouwer: Phys. Rev. B 58 (1998) R10135.
  * [8] I. L. Aleiner and A. V. Andreev: Phys. Rev. Lett. 81 (1998) 1286.
  * [9] O. Entin-Wohlman, A. Aharony, and Y. Levinson: Phys. Rev. B 65 (2002) 195411.
  * [10] P. W. Brouwer, A. Lamacraft, and K. Flensberg: Phys. Rev. B 72 (2005) 075316.
  * [11] J. Splettstoesser, M. Governale, J. König, and R. Fazio: Phys. Rev. Lett. 95 (2005) 246803.
  * [12] J. Splettstoesser, M. Governale, J. König, and R. Fazio: Phys. Rev. B 74 (2006) 085305.
  * [13] A. R. Hernández, F. A. Pinheiro, C. H. Lewenkopf, and E. R. Mucciolo: Phys. Rev. B 80 (2009) 115311.
  * [14] F. Reckermann, J. Splettstoesser, and M. R. Wegewijs: Phys. Rev. Lett. 104 (2010) 226803.
  * [15] S. Andergassen, V. Meden, H. Schoeller, J. Splettstoesser, and M. R. Wegewijs: Nanotechnology 21 (2010) 272001.
  * [16] H. L. Calvo, L. Classen, J. Splettstoesser, and M. R. Wegewijs: Phys. Rev. B 86 (2012) 245308.
  * [17] J. P. Pekola, O.-P. Saira, V. F. Maisi, A. Kemppinen, M. Möttönen, Y. A. Pashkin, and D. V. Averin: Rev. Mod. Phys. 85 (2013) 1421.
  * [18] M. Hasegawa and T. Kato: J. Phys. Soc. Jpn. 86 (2017) 024710.
  * [19] T. Aono: Phys. Rev. Lett. 93 (2004) 116601.
  * [20] J. Splettstoesser, M. Governale, J. König, and R. Fazio: Phys. Rev. Lett. 95 (2005) 246803.
  * [21] E. Sela and Y. Oreg: Phys. Rev. Lett. 96 (2006) 166802.
  * [22] A. K. Eissing, V. Meden, and D. M. Kennes: Phys. Rev. Lett. 116 (2016) 026801.
  * [23] J. I. Romero, P. Roura-Bas, A. A. Aligia, and L. Arrachea: Phys. Rev. B 95 (2017) 235117.
  * [24] M. Hasegawa and T. Kato: J. Phys. Soc. Jpn. 87 (2018) 044709.
  * [25] G. Benenti, G. Casati, K. Saito, and R. Whitney: Phys. Rep. 694 (2017) 1 .
  * [26] T. Yuge, T. Sagawa, A. Sugita, and H. Hayakawa: J. Stat. Phys. 153 (2013) 412.
  * [27] S. Nakajima and Y. Tokura: J. Stat. Phys. 169 (2017) 902.
  * [28] M. Esposito, M. A. Ochoa, and M. Galperin: Phys. Rev. Lett. 114 (2015) 080602.
  * [29] M. A. Ochoa, N. Zimbovskaya, and A. Nitzan: Phys. Rev. B 97 (2018) 085434\.
  * [30] C. Chamon, E. R. Mucciolo, L. Arrachea, and R. B. Capaz: Phys. Rev. Lett. 106 (2011) 135504.
  * [31] S. Juergens, F. Haupt, M. Moskalets, and J. Splettstoesser: Phys. Rev. B 87 (2013) 245423.
  * [32] R. S. Whitney: Phys. Rev. Lett. 112 (2014) 130601.
  * [33] N. Shiraishi, K. Saito, and H. Tasaki: Phys. Rev. Lett. 117 (2016) 190601\.
  * [34] M. Hasegawa and T. Kato: J. Phys. Soc. Jpn. 89 (2020) 064706.
  * [35] M. Hasegawa, E. Jussiau, and R. S. Whitney: Phys. Rev. B 100 (2019) 125420\.


# On the Dynamic Similarity between Bipeds and Quadrupeds:
a Case Study on Bounding

Zhenyu Gan, Ziyuan Jiao, and C. David Remy This material is based upon work
supported by the National Science Foundation under Grant No. 1453346.The
authors are with the Robotics and Motion Laboratory (RAMlab), Department of
Mechanical Engineering, University of Michigan, Ann Arbor, MI
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

In this paper, we identify passive dynamic bounding gaits of quadrupedal
robots and investigate the dynamic similarities between these gaits and those
of a bipedal model. For this purpose, we develop a passive dynamic bounding
model of a quadruped with an extended main body. This model is based on the
established SLIP model, yet has dynamic swing leg motions instead of a
predefined angle of attack. We present a gait-identification framework that is
based on a numerical continuation approach and that allows us to find all
possible passive periodic motions that emerge from pronking and bounding in
place. These gaits include forward pronking as well as forward bounding with
gathered suspension, extended suspension, and with two suspensions. By
conducting a parameter study on the body inertia, we can draw parallels
between bipedal and quadrupedal gaits, and can also show how the structure of
possible solutions changes with the main body inertia. In particular, we see
sudden changes in the gait structure at certain critical values of the main
body inertia. The work shows the connections between different types of
quadrupedal gaits and may serve as a guideline in the design and control of
gaits for energetically economical legged robots. It also sheds some light at
the question why animals with different morphology use different gaits at the
same Froude numbers.

##  I Introduction

Legged animals employ different _gaits_ when moving at different speeds.
Different gaits can be characterized by a specific footfall sequence
[alexander1984gaits, hildebrand1989quadrupedal], by the way kinetic and
potential energy exchange over time [cavagna1976sources], or by a typical
ground reaction force (GRF) profile [alexander1980optimum]. Broadly speaking,
gaits can be categorized into two classes: symmetrical and asymmetrical gaits
[hildebrand1977analysis, alexander1984gaits]. In symmetrical gaits, the left
and right half of the animal performs the exact same motion, half a stride out
of phase. This includes, for example, walking and running for bipeds and
trotting for quadrupeds. In contrast, in asymmetrical gaits, the motions of
the left and right legs show a different phase shift.

For bipeds, these asymmetrical gaits include _hopping_ and _bipedal galloping_
or _leaping_ [franz2005gait]. Red kangaroos, for example, use hopping as their
main gait over a large range of speeds [dawson1973energetic]. Galloping is
used by some bipedal primates such as lemurs [demes1999takeoff] and can also
be observed in humans, who use it in special situations such as fast turning,
stairs descending, or in low gravity environments [minetti1998biomechanics,
fiers2013biomechanics, pavei2015skipping].

In contrast to bipeds, quadrupeds can utilize a much larger variety of
asymmetrical gaits. Among other gaits, quadrupeds can _pronk_ , _bound_ , and
_gallop_. Within these broad definitions, sub-classes of gait exist.
Hildebrand [hildebrand1977analysis], for example, further classified the bound
into four categories according to the durations of the flight phases
(suspensions) between the front and back contact: 1) gaits with _no
suspension_ , 2) gaits with _gathered suspension_ , 3) gaits with _extended
suspension_ , and 4) gaits with _two suspensions_. Animals with varying body
weights and agility use different categories of bounding: heavier animals such
as the bison and giraffe tend to locomote in gaits with gathered suspension or
no suspension. While some fast moving animals like the cheetah and greyhound
prefer to use the gait with two suspensions at their highest speeds
[hudson2012high].

It is notable, that despite the vast differences in morphology, the gaits of
bipedal and quadrupedal animals share some important similarities. Heglund et
al. [heglund1982energetics] investigated the _dynamic similarity_ between
_walking_ in bipeds and quadrupeds and hypothesized that they utilize the same
mechanism similar to an inverted pendulum in which kinetic energy is exchanged
for potential (gravitational) energy and vice versa. This implies that
fluctuations in kinetic and potential energy happen out of phase. These
energy-based observations can be extended to other gaits: in bipedal _running_
or _hopping_ and in quadrupedal _trotting_ , fluctuations of potential and
kinetic energy happen in phase and both are exchanged for elastic energy
[farley1993running]. However, this analysis breaks down for asymmetrical gaits
of quadrupeds. Due to the lack of the additional pair of legs, a biped cannot
move in a fashion that is dynamically similar to a galloping quadruped
[alexander1983dynamic].

In this paper, we explore the dynamic similarity between bipedal gaits and
asymmetrical quadrupedal gaits by using _simplistic passive models_ as
illustrated in Fig. 1. These models are built on an extensive body of previous
work that investigates the passive dynamics of legged locomotion
[mochon1980ballistic, blickhan1989spring]. Most relevant in this context is
the work by Geyer et al. [Geyer2006] who first proposed a simple unified model
to explain the underlying dynamics of both bipedal walking and running. In the
present work, we employ our two models to reveal potential dynamic
relationships between bipedal gaits on the one side and quadrupedal
asymmetrical gaits on the other. For our quadrupedal model, we only focus on
_bounding_ gaits and neglect the phase difference within the fore/hind leg
pairs which is characteristic of different forms of quadrupedal _galloping_.
By letting the inertia of the torso in the quadrupedal model vary from zero to
infinitely large, we explicitly connect the two models and link all bounding
gaits of the quadrupedal model to the two-legged gaits of the bipedal model.

Figure 1: In this paper, we extend the well-known SLIP model to include
passive swing leg dynamics with constant swing leg frequency
$\omega_{\text{swing}}$ and apply it to a quadrupedal bounding model with main
body mass M and inertia J (a). Our results are contrasted against a similar
bipedal model with a point mass as the main body (b).

We believe these results to be of importance to the robotics community, as
roboticists have already developed a large number of legged machines that rely
on asymmetrical gaits. Talebi et al. [talebi2001quadruped], for example,
designed the quadruped robot Scout II and used a simple open-loop controller
to reset contact angles to demonstrate stable bounding at a range of speeds.
Wildcat developed by Boston Dynamics can bound and gallop at
$28\,\mathrm{mph}$. Similarly, MIT Cheetah 2 [park2017high] can bound at a
large range of speeds using impulse-based controller. Eckert et al. also found
that with different spine designs, Lynx-robot demonstrates bounding gaits with
different footfall patterns and duty factors [eckert2015comparing]. Most
closely related to the aims of this paper, De found different bounding and
pronking gaits by tuning the main body inertia of the Minitaur robot
[de2017modular].

Our approach to finding gaits is based on the idea of the continuation of
periodic motions. Once a gait (i.e., a periodic orbit in the state space of
the model) has been identified, we find adjacent gaits through numerical
continuation of the existing solutions. This approach is similar to the one
used by Rosa et al. to develop optimal gaits for actuated robots
[rosa2014extending]. In addition, we carefully monitor for bifurcation points
and track multiple solutions (each corresponding to a different gait) if they
arise (as done in [gan2017SLIP]). This method and the underlying model are
described in Section II. The results are presented in Section III, followed by
a brief discussion in Section IV.

##  II Methods

###  II-A Model Description

The quadrupedal model used in this study consists of an extended main body
with mass $M$ and inertia $J$ (Fig. 1a). The vertical and horizontal positions
of the COM of the main body, and the main body’s pitch angle were given by the
variables $x$, $y$, and $\phi$, respectively. Each pair of legs (with index
$i\in\left[H,F\right]$) was modeled as a massless linear spring with length
$l_{i}$, rest leg length $l_{o}$, total spring stiffness $k_{\text{leg}}$, no
damping, and a point mass $m$ at the foot. The leg pairs were connected to the
main body through frictionless rotational joints with joint angle
$\alpha_{i}$. The hip joint (back leg pair) was located at a distance of
$l_{b,H}$ behind the COM of the main body and the shoulder joint at a distance
of $l_{b,F}$ in front of the COM. We defined the value of $l_{b,H}$ to have a
negative sign. Torsional springs with stiffness $k_{\text{swing}}$ and no
damping were added to these joints. The springs are uncompressed at zero leg
angles when the main body and the respective leg are perpendicular to each
other. Rather than defining a value for $k_{\text{swing}}$ directly, we
prescribed a leg swing frequency $\omega_{\text{swing}}$, and derive:

$k_{\text{swing}}=\omega^{2}_{\text{swing}}{ml_{o}^{2}}.$ (1)

While this model is energetically conservative during flight, it would lose
energy when a foot’s velocity is brought to zero in a touch-down collision.
With such losses, it would be impossible to identify passive periodic gaits.
To resolve this issue, we took the limit of $m$ going towards zero ($m\to 0$),
similar to the method used in [garcia1998simplest] and [o2009relative]. This
implies that $k_{\text{swing}}\to 0$ as the swing frequency
$\omega_{\text{swing}}$ remains unchanged. During swing, dynamics along the
legs are ignored and the leg length remains as a constant value of their
uncompressed length $l_{o}$. To focus on the dynamics of the proposed model,
we assume the ground to be flat and rigid, and that there is no foot slipping
during the stance phase.

When the main body inertia of this quadrupedal model approaches infinity
($J\to\infty$), the torso would stop to perform any pitching motion. At this
point, and when started with a pitch of zero, the model can behave exactly as
a bipedal model in which the front and back legs of the bounding model move
like the left and right legs of the bipedal model. Such a bipedal model was
examined in detail in prior work [gan2017SLIP] (Fig. 1b).

###  II-B Equations of Motion

The equations of motion (EOM) are given for the main body as:

$\displaystyle\ddot{x}$
$\displaystyle={F_{x}}/M,\quad\ddot{y}={F_{y}}/M-g,\quad\ddot{\phi}=\tau/J,$
(2)

where $F_{x}$, $F_{y}$, and $\tau$ are the net forces and torques created by
the leg pairs. Depending on whether the leg pairs are in the air or in contact
with the ground, their dynamics differ. During the swing phase, the leg length
is set to $l_{i}=l_{o}$ and the leg angle accelerations are given by:

$\ddot{\alpha}_{i}=-\ddot{\phi}-\nicefrac{{1}}{{l_{o}}}\left(\cos\left({{\alpha_{i}+\phi}}\right)\ddot{x}+\sin\left({{\alpha_{i}+\phi}}\right)\left(g+\ddot{y}\right)\\\
-l_{b}\sin\left({{\alpha_{i}}}\right)\ddot{\phi}+l_{b}\cos\left({{\alpha_{i}}}\right)\dot{\phi}^{2}+\nicefrac{{1}}{{l_{o}}}k_{swing}\alpha_{i}\right).$
(3)

During stance, the leg angle and leg length are computed from kinematic
constraints, as we assume no sliding of the foot. The leg length is given by
$l_{i}=\frac{y+l_{b,i}\sin(\phi)}{\cos(\alpha_{i}+\phi)},$ and the leg angle
follows from

$x+l_{b}\cos(\phi)+\tan(\alpha_{i}+\phi)\left(y+l_{b}\sin(\phi)\right){\,=\,}x_{foot,i},$
(4)

where $x_{foot,i}$ is the horizontal position of the foot on the ground. We
solved this equation explicitly for leg velocities and accelerations, allowing
us to express the leg angle in both stance and swing phase as a second order
differential equation.

###  II-C Footfall Sequence

In order to determine the footfall pattern of a stride, we introduced four
timing variables $t_{j}$, ($j\in\left[Htd,Hlo,Ftd,Flo\right]$) for the touch-
down and lift-off events. To enable solutions with different footfall
patterns, the order of these events can be arbitrary, but their values are
confined within the time interval of one stride $\left[0,t_{stride}\right)$.
These timing variables are used to determine whether a leg pair is in the air
or on the ground at a given time $t$. The timing vector
$\mathbf{\boldsymbol{e}}^{T}=\left[t_{Htd},t_{Hlo},t_{Ftd},t_{Flo}\right]$
thus breaks the whole stride into 5 distinct intervals, each next interval
with a different contact configuration. Using the position and velocity vector
$\mathbf{\boldsymbol{q}}^{T}\mathrel{\mathop{:}}=[x,y,\phi,\alpha_{F},\alpha_{H}]$,
$\mathbf{\boldsymbol{\dot{q}}}^{T}\mathrel{\mathop{:}}=[\dot{x},\dot{y},\dot{\phi},\dot{\alpha}_{F},\dot{\alpha}_{H}]$
to describe the state of the system, we express the dynamics as a time variant
differential equation
$\ddot{\mathbf{\boldsymbol{q}}}=f\left(\mathbf{\boldsymbol{q}},\mathbf{\boldsymbol{\dot{q}}},t,\mathbf{\boldsymbol{e}}\right)$
which is parameterized by $\mathbf{\boldsymbol{e}}$. In addition, at the
moments of touch-down $t_{k}$, ($k\in\left[Htd,Ftd\right]$), we need to reset
the leg velocities according to the derivative of equation (4) resulting in
the additional discrete dynamics
$\mathbf{\boldsymbol{\dot{q}}}\left(t_{k}^{+}\right)=h\left(\mathbf{\boldsymbol{q}}\left(t_{k}^{-}\right),\mathbf{\boldsymbol{\dot{q}}}\left(t_{k}^{-}\right)\right)$.

###  II-D Gait Creation

Similar to [gan2017SLIP], we stated the gait creation as a _boundary value
problem_ (BVP) [hermannnonlinear]. We defined a gait as a periodic motion in
which all states except for the horizontal position $x$ return to their
original values after one full _stride_. Without loss of generality, the apex
transition ($\dot{y}=0$) was selected as the Poincaré section for this limit
cycle analysis. Finding a gait in this model was thus equivalent to solving
the following problem:

$\operatorname{T}(\mathbf{\boldsymbol{Z}}^{*}):=\begin{bmatrix}\ddot{\mathbf{\boldsymbol{q}}}-\operatorname{f}(\mathbf{\boldsymbol{q}},\mathbf{\boldsymbol{\dot{q}}},t,\mathbf{\boldsymbol{e}})\\\
\mathbf{\boldsymbol{\dot{q}}}\left(t_{k}^{+}\right)-h\left(\mathbf{\boldsymbol{q}}\left(t_{k}^{-}\right),\mathbf{\boldsymbol{\dot{q}}}\left(t_{k}^{-}\right)\right)\\\
\operatorname{R_{1-14}}(\mathbf{\boldsymbol{q}},\mathbf{\boldsymbol{\dot{q}}},\mathbf{\boldsymbol{e}},t_{stride})\end{bmatrix}=\mathbf{\boldsymbol{0}}.$
(5)

The solution vector of this equation,
$\mathbf{\boldsymbol{Z}}^{*}\mathrel{\mathop{:}}=\left(\mathbf{\boldsymbol{q}},\mathbf{\boldsymbol{\dot{q}}},\,\mathbf{\boldsymbol{e}},t_{stride}\right)$,
combines the generalized positions and velocities of the model
($\mathbf{\boldsymbol{q}}$ and $\mathbf{\boldsymbol{\dot{q}}}$) with the
vector of footfall timings $\mathbf{\boldsymbol{e}}$ and the overall stride
time $t_{stride}$. The boundary and interior-point conditions in
$\operatorname{R_{1-14}}$ were used to enforce (I) that all positions except
$x$ are periodic ($R_{1-4}$), (II) that all velocities are periodic
($R_{5-9}$), (III) that the stride ends at the Poincaré section ($R_{10}$),
and (IV) that the legs are fully extended at touchdown and liftoff
($R_{11-14}$):

$\displaystyle R_{1-4}{\,=\,}$
$\displaystyle\mathbf{\boldsymbol{\mathbf{\boldsymbol{\bar{q}}}}}\left(t_{stride}\right)-\mathbf{\boldsymbol{\mathbf{\boldsymbol{\bar{q}}}}}\left(t{\,=\,}0\right),$
(6) $\displaystyle R_{5-9}{\,=\,}$
$\displaystyle\mathbf{\boldsymbol{\mathbf{\boldsymbol{\dot{q}}}}}\left(t_{stride}\right)-\mathbf{\boldsymbol{\mathbf{\boldsymbol{\dot{q}}}}}\left(t{\,=\,}0\right),$
$\displaystyle R_{10}{\,=\,}$ $\displaystyle\dot{y}(t_{stride}),$
$\displaystyle R_{11-14}{\,=\,}$ $\displaystyle
y(t_{j})+l_{b,i}\sin(\phi(t_{j}))$ $\displaystyle-
l_{o}\cos(\phi(t_{j})+\alpha_{i}(t_{j})).$

###  II-E Continuation and Bifurcations

We found all solutions of eq. (5) through numerical continuation and
bifurcation analysis, similar to the continuation approach presented in
[rosa2014extending]. That is, starting from an initial solution
$\mathbf{\boldsymbol{Z}}_{1}$, we computed branches of solutions iteratively.
Assuming that $\mathbf{\boldsymbol{Z}}^{*}_{n}$ is the $n^{th}$ solution that
we found, we numerically search for an adjacent solution
$\mathbf{\boldsymbol{Z}}^{*}_{n+1}$ in a two-step predictor-corrector process
(Fig. 2): first we determine an initial guess
$\mathbf{\boldsymbol{Z}}^{o}_{n+1}{\,=\,}\mathbf{\boldsymbol{Z}}^{*}_{n}+d{\,\cdot\,}\left(\mathbf{\boldsymbol{Z}}^{*}_{n}-\mathbf{\boldsymbol{Z}}^{*}_{n-1}\right)$,
by continuing to move with a fixed step length of $d$ tangentially along the
existing solution branch. We then identify the exact solution
$\mathbf{\boldsymbol{Z}}^{*}_{n+1}$ by numerically solving the following
constrained problem:

$\displaystyle\operatorname{T}\left(\mathbf{\boldsymbol{Z}}^{*}_{n+1}\right)$
$\displaystyle=\mathbf{\boldsymbol{0}},$ (7a)
$\displaystyle\quad\left\lVert\mathbf{\boldsymbol{Z}}^{*}_{n+1}-\mathbf{\boldsymbol{Z}}^{*}_{n}\right\rVert$
$\displaystyle=d,$ (7b)
$\displaystyle\left(\mathbf{\boldsymbol{Z}}^{*}_{n+1}-\mathbf{\boldsymbol{Z}}^{*}_{n}\right)^{T}\left(\mathbf{\boldsymbol{Z}}^{*}_{n}-\mathbf{\boldsymbol{Z}}^{*}_{n-1}\right)$
$\displaystyle>0.$ (7c)

Here, eq. (7a) ensures that the new solution fulfills (5), eq. (7b) ensures
that we move with a constant step length along the branch, and eq. (7c)
ensures that we are not stepping backwards. By iteratively repeating these two
steps, we can identify complete solution branches.

Figure 2: The continuation algorithm used in this work employs a prediction-
correction process. We first determine an initial guess
$\mathbf{\boldsymbol{Z}}^{o}_{n+1}$ of a new solution by moving tangentially
along the branch with a step length of $d$. We then numerically solve eq. (7)
to find an exact new solution $\mathbf{\boldsymbol{Z}}^{*}_{n+1}$.

In some cases, eq. (7) has more then one solution. At such _bifurcation
points_ , new solution branches emerge, which we track in a similar fashion.
We explicitly detect these instances by tracking a change in the null space
dimension of the fréchet derivative of $\operatorname{T}$. We categorized
bifurcations by simply considering the changes in stability and the number of
solutions in the neighborhood of the bifurcation. Please see [gan2017SLIP,
merker2015numerical, hermannnonlinear, rosa2014extending] for more details of
this process.

Solutions obtained through this process can be characterized by the total
energy $E_{tot}$ stored in the system. Since the model is energetically
conservative, this value is constant throughout a stride and a good metric to
distinguish different solutions.

###  II-F Parameter Study

In order to keep our solutions general, we normalized all state and parameter
values with respect to $M$, $g$, and $l_{o}$ [Hof:1996]. We set the leg pair
stiffness to $k_{\text{leg}}{\,=\,}20\,\mathrm{\nicefrac{{Mg}}{{l_{o}}}}$ and
the swing frequency to
$\omega_{\text{swing}}{\,=\,}\sqrt{5}\,\mathrm{\sqrt{\nicefrac{{g}}{{l_{o}}}}}$.
The front leg pair and the hind leg pair were attached symmetrically to the
center of the body with $l_{b,F}=0.5\,\mathrm{l_{o}}$ and
$l_{b,H}=-0.5\,\mathrm{l_{o}}$. These dimensionless values are comparable to
our existing hardware [smith2017design]. The parameter of interest in this
paper was the main body inertia $J$. We started our investigation with the
extremal case of an infinite main body inertia $J=\infty$. While this case is
highly contrived, it allows us to draw the dynamic connections between the
gaits of a biped and the asymmetrical gaits of a quadruped, as the two models
are dynamically equal for $J=\infty$. When then gradually reduced the value of
$J$ and iteratively re-computed all primary solution branches, ignoring
motions with multiple leg-swing oscillations during swing in the process.

##  III Results

In this section, we visualize different solution branches by showing how
selected states $\mathbf{\boldsymbol{Z}}^{*}$ on the Poincaré section evolve
as a function of energy $E_{tot}$. Furthermore, we visualize and numerically
report bifurcation points (Table I). We also illustrate typical solutions
along branches by showing a series of key-frames of an animation of the
motion.

In order to highlight the dynamic similarities between the bipedal model and
quadrupedal model and show how do they diverge from each other, we also
conducted an detailed parameter study on the body inertia. In the bipedal
model (Fig. 1b) the legs are massless and the COM motion and the swing leg
rotations are fully decoupled. That is, for any finite value of the body
inertia $J$, the main body orientation will remain the same during the whole
stride. In other words, the solution branches of the bipedal gaits are
invariant under different values of $J$. As shown in the following sections,
this will not hold for the quadrupedal model. We first report on the extremal
case of $J=\infty$, then establish the correspondence between asymmetrical
quadrupedal gaits and the gaits of a biped, and finally report how these
solutions change as a function of $J$. For some selected periodic motions, we
also present a visualization of the motion111Please see the multimedia
attachment for a video showing these motions..

###  III-A Infinite Main Body Inertia, $J=\infty$

####  III-A1 Pronking and Bounding in Place

The most basic periodic motions of this model and the starting point for our
exploration are simple _pronking in place_ (_PP_) and _bounding in place_
(_BP_) gaits. In both gaits, the model moves purely vertically and the legs
are pointing straight down during the entire gait cycle. The primary
difference between pronking and bounding is the contact sequence. In pronking
(Fig. 3a), the front and hind leg pairs strike and leave the ground
simultaneously with a single flight phase in between, whereas in bounding
(Fig. 3b), there exist two flight phases and the front and hind leg pairs are
in contact at separate times. Strictly speaking, there are two versions of
this bounding gait, one in which the front contact happens first and one in
which it happens second in the stride that we consider. Since we assume an
infinite main body inertia, the off-center loading during bounding does not
induce any pitching motion and the main body remains horizontal. Increasing
$E_{tot}$ increases the apex height $y_{o}$ during flight and leads to more
compression of the legs during stance. The theoretical limit for the maximal
pronking height is $20\,\mathrm{l_{o}}$ and for the maximal bounding height it
is $10\,\mathrm{l_{o}}$. For higher apex heights, the legs would become fully
compressed during stance. Apart from the contact sequence, the solution
branches overlap. They are illustrated for $y_{o}<2\,\mathrm{l_{o}}$ in Fig.
4.

Figure 3:  Shown are the motions of pronking in place (_PP_ , a), bounding in
place (_BP_ , b), pronking with extended suspension (_PE_ , c), and pronking
with gathered suspension (_PG_ , d) of the quadrupedal model. The motions are
independent of the main body inertia $J$. Note that for _BP_ , only the first
half of a full stride is shown.  Figure 4: This visualization of the solution
branches shows the initial main body height $y_{o}$ and front leg angle
$\alpha_{F,o}$ at the Poincaré section for the _PP_ , _BP_ , _PE_ , and _PG_
gaits. The _PE_ and _PG_ branches bifurcate from the _PP_ branch at point $A$.
The solutions are independent of the main body inertia $J$. The labels (a)-(d)
refer to specific solutions that are shown in Fig. 3.

####  III-A2 Pronking in Place with Leg Swing (PE and PG)

Along the _PP_ branch, a pitchfork (PF) bifurcation occurs at point $A$ at an
apex height of $y_{o}=1.71\,\mathrm{l_{o}}$. The two new branches break the
front/hind symmetry of the _PP_ motions and correspond to two in-place gaits
in which the leg pairs swing against each other: _in-place pronking with
extended suspension_ (_PE_ , Fig. 3c) and _in-place pronking with gathered
suspension_ (_PG_ , Fig. 3d). As the _PE_ branch moves away from the
bifurcation point with increasing $E_{tot}$, the leg pairs in the
corresponding solutions rotate outward, reach their maximal excursion at apex,
and then move back inward again before hitting the ground pointing towards
each other. The opposite happens on the _PG_ branch, as the leg pairs gather
towards the center during the flight phase and are extended during stance. At
the same apex height, the main body motions of these two branches are exactly
the same. Yet, their leg motions are inverted. That is, at any instant of
motion, the angles and angle rates of the hind legs on _PE_ are identical to
those of the front legs on _PG_ , and vice versa. The two branches terminate
when the main body COM reaches an unstable equilibrium point during stance,
where the sum of the leg forces is equal to the weight of the main body,
$Mg{\,=\,}2k_{\text{leg}}\cos(\alpha_{i}+\phi)(l_{o}{\,-\,}l_{i})$. Increasing
the energy even further would result in the COM ‘dropping’ all the way into
the ground.

TABLE I: Initial states $\mathbf{\boldsymbol{X}}_{o}$ and event timings at the bifurcation points $A$ \- $F$. BP |  unit |  $A$ |  $B$ |  $C$ |  $D$ |  $E$ |  $F$
---|---|---|---|---|---|---|---
$\dot{x}_{o}$ |  $\mathrm{\sqrt{gl_{o}}}$ |  0 |  0 |  1.41 |  9.71 |  0 |  4.34
$y_{o}$ |  $\mathrm{l_{o}}$ |  1.71 |  1.61 |  1.58 |  0.45 |  1.02 |  1.09
$\dot{y}$ |  $\mathrm{\sqrt{gl_{o}}}$ |  0 |  0 |  0 |  0 |  0 |  0
$\phi_{o}$ |  $\left[\mathrm{\cdot}\right]$ |  0 |  0 |  0 |  0 |  0 |  0
$\dot{\phi}_{o}$ |  $\mathrm{\sqrt{g/l_{o}}}$ |  0 |  0 |  0 |  0 |  0 |  0.20
$\alpha_{H,o}$ |  $\left[\mathrm{\cdot}\right]$ |  0 |  0 |  0 |  0 |  0 |  1.30
$\dot{\alpha}_{H,o}$ |  $\mathrm{\sqrt{g/l_{o}}}$ |  0 |  0 |  1.22 |  3.25 |  0 |  1.38
$\alpha_{F,o}$ |  $\left[\mathrm{\cdot}\right]$ |  0 |  0 |  0 |  0 |  0 |  -1.30
$\dot{\alpha}_{F,o}$ |  $\mathrm{\sqrt{g/l_{o}}}$ |  0 |  0 |  1.22 |  3.25 |  0 |  1.38
$t_{htd}$ |  $\mathrm{\sqrt{l_{o}/g}}$ |  1.19 |  1.11 |  1.13 |  0.79 |  0.21 |  0.72
$t_{hlo}$ |  $\mathrm{\sqrt{l_{o}/g}}$ |  1.73 |  1.65 |  1.60 |  1.01 |  1.28 |  1.04
$t_{ftd}$ |  $\mathrm{\sqrt{l_{o}/g}}$ |  1.19 |  1.11 |  1.13 |  0.79 |  1.69 |  1.56
$t_{flo}$ |  $\mathrm{\sqrt{l_{o}/g}}$ |  1.73 |  1.65 |  1.60 |  1.01 |  2.76 |  1.88
$t_{apx}$ |  $\mathrm{\sqrt{l_{o}/g}}$ |  2.91 |  2.76 |  2.73 |  1.80 |  2.97 |  2.60

####  III-A3 Pronking Forward (PF)

Breaking the frontal symmetry, we found a _pronking forward_ (_PF_ , Fig. 5e)
gait that branches off the in-place pronking _PP_ solutions with a PF
bifurcation at point $B$ (Fig. 7). This bifurcation happens at an apex height
of $y_{o}{\,=\,}1.61\,\mathrm{l_{o}}$. In this gait, the two leg pairs
simultaneously strike and leave the ground and exhibit the same leg angles
throughout the motion. Along this branch, an increase in $E_{tot}$ leads to an
increased forward speed with larger contact angles, until a turning point with
the maximal energy $E_{tot}{\,=\,}61.43\,\mathrm{Mgl_{o}}$
($\dot{x}_{o}{\,=\,}11.02\,\mathrm{\sqrt{gl_{o}}}$) is reached. Solutions
further along the branch have lower energy and lower forward speeds. The gait
vanishes when the leg springs are fully compressed during stance, which
happens at a forward speed of
$\dot{x}_{o}{\,=\,}6.43\,\mathrm{\sqrt{gl_{o}}}$. One should note that there
also exists a corresponding pronking backward gait that branches off in the
opposite $\dot{x}_{o}$ direction.

Figure 5: Shown are the motions of pronking forward (_PF_ , e), bounding with
gathered suspension (_BG_ , f), bounding with extended suspension (_BE_ , g),
and bounding with two suspensions (_B2_ , h) of the quadrupedal model with
infinite inertia ($J{\,=\,}\infty$). Figure 6: This figure highlights the
correspondence between bipedal gaits and quadrupedal bounding with infinite
inertia (as shown in Fig. 5). In particular, bipedal hopping (e) corresponds
to quadrupedal pronking (_PF_), bipedal leaping (f, g) to bounding with
gathered/extended suspension (_BE_ , _BG_), and bipedal running (h) to
bounding with two suspensions (_B2_). These solution branches of the bipedal
gaits are invariant under different values of $J$ in the bipedal model.
Figure 7: This visualization of the solution branches shows the initial main
body height $y_{o}$, forward speed $\dot{x}_{o}$, and front leg angle
$\alpha_{F,o}$ at the Poincaré section for all gaits of the quadrupedal model
with infinite inertia. The bipedal branches are invariant with the main body
inertia $J$ and are identical to the quadrupedal branches with infinite
inertia. The PF branch (green) bifurcates from the _PP_ branch (grey) at point
$B$. The _BG_ and _BE_ branches (red) form a closed loop that bifurcates from
the _PF_ branch a points $C$ and $D$. Along the _BG_ and _BE_ branches, the
contact sequence changes and a second flight phase is introduced (dashed
section). Finally, the _B2_ branch (blue) bifurcates from the _BP_ branch at
point $E$. A second version of this branch (dotted line) can be obtained by
redefining the Poincaré section to the second apex transit. The labels e-h
refer to specific solutions that are shown in Fig. 5 and Fig. 6.

####  III-A4 Bounding with Gathered/Extended Suspension (BG and BE)

Along the _PF_ branch, we found a PF bifurcation at point $C$ (Fig. 7. On the
new branches, the front/hind symmetry of the _PF_ motions breaks, and the two
leg pairs loose their synchronization as they either gather inwards or extend
outwards during the flight phase. The corresponding gaits are called _bounding
with gathered suspension_ (_BG_ , Fig. 5f) and _bounding with extended
suspension_ (_BE_ , Fig. 5g). Moving along these branches, first with
increasing energy $E_{tot}$ and after a turning point with decreasing energy,
we can observe that they eventually merge back to the _PF_ branch at another
PF bifurcation at point $D$. These two bounding branches thus form a closed
loop, which is symmetrical with regard to the motion of the front and hind leg
pairs. That is, the motion of hind legs in _BG_ resembles the motion of front
legs in _BE_. Despite the differences in leg motions, the COM moves along the
same trajectories in both branches. One should note that along these branches,
the footfall sequence does change. In the central part of the branches (dashed
lines in Fig. 7), a second flight phases arises, as the order of touchdown and
lift-off of the leg pairs change in the middle of the stride.

####  III-A5 Bounding with Two Suspensions B2

The _BP_ branch has no front/hind symmetry due to the staggered contact
sequence. It thus bifurcates directly into branches that include
forward/backward motion (Fig. 7). Such a PF bifurcation was found at point $E$
at an apex height of $y_{o}{\,=\,}1.02\,\mathrm{l_{o}}$. The resulting motion
is _bounding with two suspensions_ (_B2_ , Fig. 5h). In this bounding gait,
after lift-off, the swing leg pair first continues to extend backward until it
reaches zero speed and the hip spring starts pulling it forward.. Then the leg
pair swings forward, reaches its maximal forward extension, and strikes the
ground during the backward retraction. Because of the infinite body inertia,
the main body remains horizontal within the whole stride cycle, the two flight
phases have equal durations, and the motion of the front and hind leg pairs is
exactly the same, just half a stride apart. Because of the two air-phases, the
_B2_ branch is represented by two curves in Fig. 7: the solid line corresponds
to the apex transit after the front leg stance (legs pointing inward at apex),
and the dotted line to the apex after the back leg stance (legs pointing
outward). Similar to the _PF_ branch, increasing $E_{tot}$ along the _B2_
branch increases forward speed until a turning point at
$E_{tot}{\,=\,}55.58\,\mathrm{Mgl_{o}}$
($\dot{x}_{o}{\,=\,}10.52\,\mathrm{\sqrt{gl_{o}}}$). The branch terminates
when full leg compression is reached at a speed of
$\dot{x}_{o}{\,=\,}5.35\,\mathrm{\sqrt{gl_{o}}}$. Similar to pronking, there
exist backward gaits that branch off in the negative $\dot{x}_{o}$ direction.

###  III-B Correspondence to Bipedal Gaits

As the main body does not pitch in the case of the infinite body inertia, we
can identify direct correspondences between the asymmetrical bounding gaits of
the quadrupedal model and the gaits of the bipedal model (Fig. 6). To this
end, we let the front leg pair of the quadrupedal model correspond to the left
leg of the bipedal model, and the hind leg pair correspond to the right leg.
In particular, _pronking_ corresponds to _bipedal hopping_ (Fig. 6e).
_Bounding with gathered suspension_ corresponds to _bipedal leaping_ with the
left leg leading (Fig. 6f) and _bounding with extended suspension_ to _bipedal
leaping_ with the right leg leading (Fig. 6g). Finally, _bounding with two
suspensions_ corresponds to _bipedal running_ (Fig. 6h).

###  III-C Finite Main Body Inertia, $J<\infty$

Since in the pronking gaits (_PG_ , _PE_ , and _PF_), all legs strike and
leave the ground simultaneously and the front and hind leg pairs have either
the same or opposite angles, the pitch torques on the main body are always
balanced. In these motions, the main body will remain horizontal, independent
of the value of $J$. In the following, we thus focus on the bounding branches
(_BG_ , _BE_ , and _B2_) which are indeed influenced by the choice of the main
body inertia $J$.

####  III-C1 Large Inertia,
$J\in\left[1.047\,\mathrm{Ml^{2}_{o}},\infty\right)$

As the value of $J$ gradually decreases from infinity, the structure of the
bounding gaits remains roughly the same until $J=1.047\,\mathrm{Ml^{2}_{o}}$.
The net torque on the torso creates an angular acceleration and the main body
pitches accordingly. This pitching motion is combined with the bounding
motion, yet the general motion and footfall sequence remain the same (Fig. 8).
What is lost, is the symmetry between the _BG_ and _BE_ branches. That is, the
motion of the hind legs in _BG_ does not longer resemble the motion of the
front legs in _BE_. This can also be observed in the branches (Fig. 9), where
the _BG_ and _BE_ are no longer symmetrical with respect to the _PF_ branch.
Notably, a new turning point arises along the _BG_ branch near Point $C$, as
the energy level along the branch now initially decreases. Similarly, we can
observe an asymmetry in the two different apex transits in _bounding with two
suspensions_ (_B2_).

Figure 8: Shown are the motions of bounding with gathered suspension (_BG_ ,
i), bounding with extended suspension (_BE_ , j), bounding with two
suspensions (_B2_ , k) of the quadrupedal model with a main body inertia of
$J{\,=\,}1.047\,\mathrm{Ml^{2}_{o}}$. Figure 9: This figure shows how the
quadrupedal gait branches evolve as the inertia goes from $J{\,=\,}\infty$
(faded lines) to $J{\,=\,}1.047\,\mathrm{Ml^{2}_{o}}$ (solid lines). At this
critical inertia value, the _B2_ and _BG_ branches merge at bifurcation point
$F$. The labels i-k refer to specific solutions that are shown in Fig. 8.

In the _bounding with two suspensions_ gait, we can also observe that, as the
inertia value decreases, the gathered suspension becomes longer than the
extended suspension. The main body motions along the _B2_ branch therefore
become similar to the solutions with two flight phases along the _BG_ branch;
in particular with regard to the timing of these flight phases and in the
sense that they both have a positive pitching velocity at the start of each
stride and a negative rotational velocity in the middle of that stride (Figs.
8h & 8j). As a result, the _BG_ and _B2_ branches approach each other for
decreasing values of $J$ and eventually merge at a critical inertia value of
$1.047\,\mathrm{Ml^{2}_{o}}$ (Point $F$).

####  III-C2 Intermediate Inertia,
$J\in\left[0.501\,\mathrm{Ml^{2}_{o}},1.047\,\mathrm{Ml^{2}_{o}}\right)$

As the main body inertia is decreased below $1.047\,\mathrm{Ml^{2}_{o}}$, the
_BG_ and _B2_ break apart at the bifurcation point $F$ (Fig. 10) but stay
connected with each other. As a result, we now get a seamless continuation
from the _B2_ branches onto the _BG_ branches at two new turning points
($F_{1}$ and $F_{2}$, Fig. 10). These two turning points move away from each
other as the inertia $J$ keeps decreasing. $F_{1}$ moves towards lower forward
speeds and $F_{2}$ towards higher forward speeds. As a result, bounding with
gathered suspension (_BG_) and bounding with two suspensions (_B2_) cease to
exist at the intermediate speed range and bounding with an extended suspension
(_BE_ branch) becomes the only passive bounding gait for certain speeds. The
pronking branch _PF_ and the bifurcation points ($C$ and $D$) remain in the
same locations on Poincaré section.

When the inertia value drops below $0.501\,\mathrm{Ml^{2}_{o}}$, the turning
points $F_{1}$ and $F_{2}$ disappear (Fig. 11). They merge with the existing
turning points on the _BG_ and _B2_ branches. Energy now monotonically
decreases when moving along _BG_ starting from point $C$, and only a single
turning point remains on the _BG_ branch starting from point $D$. With this,
it is not possible anymore to do a formal distinction between the _BG_ and
_B2_ branches. In a sense, the original bounding gait with double suspension
(_B2_) ceases to exist. In contrast, the _BE_ branch does not undergo any
fundamental changes and remains strikingly similar to bipedal leaping. The
branch leaves pronking at point $C$ and joins back to the pronking branch at
point $D$ after going through a turning point.

Figure 10: This figure shows how the quadrupedal gait branches evolve as the
inertia goes from $J{\,=\,}1.047\,\mathrm{Ml^{2}_{o}}$ (faded lines) to
$J{\,=\,}0.9\,\mathrm{Ml^{2}_{o}}$ (solid lines). The _BG_ and _B2_ branches
join and bifurcation point $F$ forms two turning points $F_{1}$ and $F_{2}$
that move away from each other. Figure 11: This figure shows how the
quadrupedal gait branches evolve as the inertia goes from
$0.9\,\mathrm{Ml^{2}_{o}}$ (faded lines) to $0.501\,\mathrm{Ml^{2}_{o}}$
(solid lines). The turning points $F_{1}$ and $F_{2}$ merge with other turning
points and disappear. at this point, the _B2_ branches cease to exist.

####  III-C3 Small Inertia, $J\in\left(0,0.501\,\mathrm{Ml^{2}_{o}}\right)$

As the inertia value keeps decreasing from $0.501\,\mathrm{Ml^{2}_{o}}$, the
rotations of the main body are getting smaller within one stride. In the
meantime, the branches _BG_ and _BE_ are approaching the pronking branch _PF_.
When the main body in the model has zero inertia, any time lag in the touch-
downs will cause infinite acceleration of the main body. Therefore all
bounding branches vanish, and the pronking branch _PF_ becomes the only gait
with nonzero speeds.

##  IV Conclusion and Discussion

In this paper, we have identified passive dynamic bounding gaits, and have
investigated the dynamic similarity between bipedal and asymmetrical
quadrupedal gaits. We did so by studying the passive motions of a simplistic
energetically-conservative quadrupedal model. Since there is no control in
this model, its motion is fully determined by its mechanical dynamics. Still,
by using numerical continuation techniques, we were able to identify the most
prevalent bounding gaits found in nature: pronking, bounding with extended
support, bounding with gathered support, and bounding with two flight phases.
These motions emerge completely passively, and may serve as a template to
develop energetically economical motions for legged robotic systems. While an
actual robot would need to compensate for energetic losses due to friction,
damping, and contact collisions, staying close to these passive motions would
ensure that the robot avoids negative actuator work. This is because the
energetic fluctuations that happen during locomotion can be fully realized by
the dynamics of the mechanical system.

As a starting point of our analysis, we assumed an infinite inertia of the
robot’s main body. Since this prevents any main body pitch, we were able to
draw direct connections between the gaits of a biped and those of a quadruped.
In fact, comparing the results presented in this work with our prior work on
bipeds [gan2017SLIP], showed that for the same set of parameters and the same
choice of the Poincaré section the two models behave exactly identical:
quadrupedal pronking is equal to bipedal hopping, quadrupedal bounding with a
gathered/extended suspension is equal to bipedal leaping, and bounding with
two suspensions is equal to bipedal running. Such similarities between bipedal
and quadrupedal gaits have been long postulated and have been studied
previously using the established SLIP model. However, due to the large
dimensions of solutions in the SLIP model, previous research only focused on
symmetrical gaits, such as the correspondence of running and trotting
[alexander1984gaits, farley1993running]. Knowing about the connections between
the bipedal and quadrupedal model, and also about the connections between the
different quadrupedal bounding gaits may help the design of robots and
controllers that can switch between different gaits.

Naturally, the similarities between bipeds and quadrupeds become less evident
as the main body inertia in the quadruped decreases and the solutions change.
In the bipedal model, because the legs are all connected to the COM, the main
body stays horizontal in all gaits and for all values of the main body
inertia. On the other hand, the two leg pairs in the quadrupedal model are
mounted at two ends of the torso, which results in a nonzero torque when the
footfalls of these two leg pairs are not simultaneous. Therefore, the
quadrupedal gaits start to differ as the main body begins to oscillate. Most
notably are the changes in bounding with a gathered suspension. For low values
of the main body inertia, these solutions disappear almost completely, leaving
bounding with an extended suspension and bounding with two suspensions as the
only passive modes of locomotion for a wide range of speeds.

These theoretical results are supported by recent experiments on hardware. De
[de2017modular] has shown that the stable phase relationship of two leg pairs
in the Minitaur robot are highly dependent on the body inertia. In the results
presented in this paper, when the body inertia is smaller than
$1\,\mathrm{Ml^{2}_{o}}$, the leg pairs are weakly coupled and bounding is
stable. As the body inertia further increases, the two leg pairs are strongly
coupled and the robot exhibits a stable pronking gait. Similar results were
also observed on the MIT Cheetah 2, when the gait pattern stabilizer was not
enforced, the phase difference in the resulting gait become smaller
[park2017high]. These results suggest that different morphology and varying
body inertia of robots may affect the stability and availability of their
gaits. As for roboticists, it is important to tune this value on the hardware
in order to achieve a specific robot gait.

In this work, we limited ourselves to bounding; that is, to gaits in which leg
pairs move together. While this greatly reduced the complexity of the
analysis, it should be noted that bounding is not a particularly common or
efficient gait. In fact, in our prior work, we have shown that bounding is
energetically inferior in quadrupedal robotic locomotion [xi2016selecting].
This holds both for quadrupeds with and without an articulated spine
[yesilevskiy2018spine]. A natural extension of this work is thus to expand our
approach to discover half-bound (three-beat) and galloping (four-beat) gaits
in which the footfalls of the front and hind legs are not necessarily
synchronized [hildebrand1977analysis]. As these gaits include bounding as a
special case, our preliminary results suggest that half-bound and galloping
gaits emerge through additional bifurcations from the solution branches that
have been presented in this paper.

Our work shows that all asymmetrical bounding can be generated with a single
passive model and a single set of parameters. In a sense, the different gaits
are just different oscillation modes of a single mechanical framework. Our
model therefore has the potential to serve as a general template
[full1999templates] for the study of quadrupedal locomotion. It might also
explain, why our past optimization studies [xi2016selecting] have converged to
similar motions with varying flight phases. Furthermore, the approach could be
extended to be used in the gait design of legged systems. In order to
reproduce desired gaits efficiently, one can tune the parameter values of
their specific hardware, such as the body mass, leg length, or the body
inertia, to ensure the gaits are passively available.