# Solitons with Self-induced Topological Nonreciprocity

Pedro Fittipaldi de Castro and Wladimir A. Benalcazar  Department of Physics,
Emory University, Atlanta, GA 30322, USA

###### Abstract

The nonlinear Schrödinger equation supports solitons – self-interacting,
localized states that behave as nearly independent objects. We show solitons
with self-induced nonreciprocal dynamics in a discrete nonlinear Schrödinger
equation. This nonreciprocal behavior, dependent on soliton power, arises from
the interplay between linear and nonlinear terms in the equations of motion.
Initially stable at high power, solitons exhibit nonreciprocal instabilities
as power decreases, leading to unidirectional acceleration and amplification.
This behavior is topologically protected by winding numbers on the solitons’
mean-field Hamiltonian and their stability matrix, linking nonlinear dynamics
and point gap topology in non-Hermitian Hamiltonians.

Nonreciprocity, or the unequal behavior of channels moving in opposite
directions, is a peculiar feature of physical systems that lack spatial
reflection and time-reversal symmetries. It has wide applications in optical
and electrical circuit components [1, 2], such as isolators [3] and
circulators [4], which are crucial for source protection and prevention of
spurious interferences. Some platforms rely on a magnetic field to break the
necessary symmetries [5, 6, 7, 8, 9], while others achieve nonreciprocity, for
example, by spatially and temporally modulating resonator arrays [10, 11, 12,
13, 14].

When defined in a crystalline structure, some systems display a robust,
topologically protected form of nonreciprocity that occurs, for instance, at
the 1D edges of systems hosting the quantum Hall effect [15, 16, 17, 18] and
in non-Hermitian lattices with a “skin effect” [19, 20, 21, 22, 23]. In the
former case, edge states exhibit unidirectional propagation due to a
topological bulk-boundary correspondence [24], while in the latter, the energy
spectra take complex values $E_{n}(k)\in\mathbb{C}$ that wind around a certain
point in the complex energy plane across the Brillouin zone [25, 26, 27]. The
resulting winding number is a topological invariant associated with fields
that exhibit unidirectional amplification and acceleration [28, 21, 29, 30].

In parallel with the development of non-Hermitian quantum mechanics, advances
in condensed matter, such as the recent observation of fractional Chern
insulators [31, 32, 33, 34], have highlighted the importance of interactions
in enabling novel topological phases. Still, the role of interactions and
nonlinearities as catalysts of topologically protected nonreciprocity remains
largely unexplored. One way of attacking the problem is to consider nonlinear
lattices supporting discrete solitons or breathers – collective self-
interacting localized states that behave as individual particles, with well-
defined quantized properties [35, 36], and analyze their dynamical properties.

In this work, we combine the Ablowitz-Ladik (AL) model [37, 38], a particular
discretization of the nonlinear Schödinger equation, with a parity-broken
Hermitian linear lattice – neither of which is nonreciprocal on its own – and
show that the resulting model supports solitons with self-induced topological
nonreciprocity. We show that while the asymmetry in the spatial distribution
of our solitons about their “center of mass” is necessary for their
nonreciprocal dynamics, it is not sufficient. Additionally, their
nonreciprocal behavior is topologically protected by point gaps with
nontrivial windings in the solitons’ (non-Hermitian) stability and mean-field
(MF) energy spectra under periodic boundary conditions (PBC), thus unveiling
an intimate connection between nonlinear, nonreciprocal dynamics and point gap
topology in non-Hermitian linear Hamiltonians. Indeed, as we shall see, the
nonreciprocal characteristics of the solitons described herein are akin to
those found in the wave packet dynamics of linear systems exhibiting a skin
effect [39, 28].

The type of nonlinearity we describe appears in the approximate descriptions
of the dynamics of small-amplitude excitons in atomic chains with exchange or
dipole-dipole interactions [40], energy transport in biological systems [41],
or even solitons in nonlinear waveguides [42].

Figure 1: Dynamics of nonreciprocal solitons. (a,b) Time evolution of the
square amplitude per site of solitons in the parity-broken AL model, Eq. (1).
In (a), the initial configuration is a static soliton with norm $|\phi|=1.05$,
energy $\omega=-2.32$, period of oscillation $T=2\pi/|\omega|=2.71$, and
random complex perturbations $|\delta\phi_{a}|<0.0001$ added across the
lattice at each time step under PBC. In (b), the initial configuration is an
unperturbed static soliton with norm $|\phi|=0.68$, energy $\omega=-2.09$, and
period of oscillation $T=2\pi/|\omega|=3.01$ under PBC. (c) Evolution of the
square amplitude of a soliton with $|\phi|=0.68$, energy $-2.09$, and period
of oscillation $T=2\pi/\omega=3.01$, but now located next to the boundary
under open boundary conditions. (d-f) Squared norm as a function of time for
the solitons displayed in (a-c). The insets in (d-f) show the field profiles
$\phi_{R}$ at $t=0$. All simulations use the 4th order Runge-Kutta algorithm
with a time step $\Delta t=0.01$ and lattice parameters $t_{1}=1.35$,
$t_{2}=0.65$, $\Delta=0.35$ and $L=40$ unit cells ($80$ sites).

_Model –_ As a minimal example of a nonlinear system supporting self-induced
nonreciprocal solitons, we consider a field $\Phi$ defined in a 1D lattice and
governed by the equations of motion

$\displaystyle{\mathrm{i}}{\frac{\mathrm{d}}{\mathrm{d}t}}\Phi_{a}=-\sum_{b}(h_{ab}+|\Phi_{a}|^{2}f_{ab})\Phi_{b},$
(1)

where the composite index $a=(R,\alpha)$ runs over unit cells $R=1,\dots,L$
and their internal degrees of freedom $\alpha=A,B$ (the lattice basis). The
linear part of (1) is determined by the tight-binding matrix

$\displaystyle h_{(R,\alpha)(R^{\prime},\beta)}$
$\displaystyle=\left(t_{1}\sigma^{x}_{\alpha\beta}+\Delta\sigma^{z}_{\alpha\beta}\right)\delta_{R,R^{\prime}}$
$\displaystyle\ \
+t_{2}\left(\sigma^{+}_{\alpha\beta}\delta_{R-1,R^{\prime}}+\sigma^{-}_{\alpha\beta}\delta_{R+1,R^{\prime}}\right),$
(2)

where $\sigma^{x},\sigma^{y},\sigma^{z}$ are the Pauli matrices and
$\sigma^{\pm}=(\sigma^{x}\pm i\sigma^{y})/2$. Here, $t_{1}$ and $t_{2}>0$ are
intra- and inter-unit cell NN hopping parameters, and $\Delta$ represents the
staggered on-site energies of the $A$ and $B$ sublattices.

The nonlinear contribution to (1) is proportional to

$f_{(R,\alpha)(R^{\prime},\beta)}=(\sigma^{x}_{\alpha\beta}\delta_{R,R^{\prime}}+\sigma^{+}_{\alpha\beta}\delta_{R-1,R^{\prime}}+\sigma^{-}_{\alpha\beta}\delta_{R+1,R^{\prime}}),$
(3)

and also connects NNs. The probability density function (PDF)
$|\Phi_{R,\alpha}|^{2}$ controls the strength of the nonlinear hopping from
the site $(R,\alpha)$ to its NN destination. The system breaks parity symmetry
by setting $t_{1}\neq t_{2}$ and $\Delta\neq 0$, which is a necessary
condition for nonreciprocal solitons. Still, their nonreciprocal dynamics
additionally depend on the power of the soliton’s field. At high enough
powers, parity-broken solitons are always static (and hence reciprocal). At
lower powers, nonreciprocal linear instabilities appear on static solitons,
Fig. 1(a), while at even lower powers, the static soliton solutions disappear,
replaced instead by solitons with sustained unidirectional acceleration and
amplification, Fig. 1(b), except at the one end of the lattice to which the
nonreciprocal solitons accelerate, where a single static soliton solution
persists, Fig. 1(c).

_Static solitons and their linear stability –_ To expose the nonreciprocal
behavior in our model, we investigate the fate of static solitons in the bulk
of a lattice obeying (1) as we vary the parity-breaking parameters
$t_{2}/t_{1}$ and $\Delta$. We find these solutions by changing variables to
the rotating frame $\Phi_{a}(t)=e^{-{\mathrm{i}}\omega t}\phi_{a}(t)$, where
$\omega\in\mathbb{R}$ is some real frequency to be determined [43]. In the
rotating frame, the equations of motion (1) take the form

${\mathrm{i}}{\frac{\mathrm{d}}{\mathrm{d}t}}\phi_{a}=-\sum_{b}\left(h_{ab}+\omega\delta_{ab}+|\phi_{a}|^{2}f_{ab}\right)\phi_{b}.$
(4)

The static solitons we are after must be stationary in the rotating frame,
${\frac{\mathrm{d}}{\mathrm{d}t}}\phi_{a}=0$ [44], yielding the time-
independent expression

$\displaystyle\omega\phi_{a}$
$\displaystyle=-\sum_{b}(h_{ab}+|\phi_{a}|^{2}f_{ab})\phi_{b}\equiv\sum_{b}H_{ab}(\phi)\phi_{b},$
(5)

where $H(\phi)$ is a mean-field (MF) Hamiltonian built using the state $\phi$
with frequency $\omega$. We self-consistently solve for $\phi$, $\omega$, and
$H(\phi)$ as detailed in the Supplemental Material.

In addition to $\phi$, $H(\phi)$ has a multitude of extended eigenvectors
$\psi$ with energies $\varepsilon$. In the $|\phi|^{2}\to 0$ limit, the mean-
field Hamiltonian reduces to the tight-binding matrix $h$ in (1), so the
eigenstates $\psi$ are simply Bloch waves forming the two energy bands

$\varepsilon_{\pm}(k)=\pm\sqrt{(t_{1}+t_{2}\cos k)^{2}+(t_{2}\sin
k)^{2}+\Delta^{2}}$

of a linear Rice-Mele chain [45]. The extreme levels of the linear dispersion
are $E_{\pm}=\pm\sqrt{(t_{1}+t_{2})^{2}+\Delta^{2}}$. As $|\phi|^{2}$
increases, the Bloch waves at the edges $E_{\pm}$ of the two energy bands
continuously deform to give rise to two localized eigenvectors of $H(\phi)$:
the soliton $\phi$ itself and a state $\phi^{\prime}$ with frequency
$-\omega$. Thus, we say that a given soliton “bifurcates” from the edge of a
particular energy band. From now on, we will focus on the properties of the
soliton $\phi$ bifurcating from the lower band, corresponding to the ground
state of the MF Hamiltonian. Although the energy $\omega$ of the soliton is
real, $H(\phi)$ is non-Hermitian; thus, the energies $\varepsilon$ of the
extended states may be complex. From (5), we see that $H^{*}(\phi)=H(\phi)\neq
H^{T}(\phi)$. Therefore, the MF Hamiltonian has time-reversal symmetry (TRS)
and breaks the ramified version of TRS for non-Hermitian Hamiltonians, called
$\mathrm{TRS}^{\dagger}$, lying in the symmetry class AI of the 38-fold
classification of non-Hermitian Hamiltonians [46], which has a $\mathbb{Z}$
topological invariant in one spatial dimension in the presence of a point gap.

Along with its mean-field spectrum, we probe the soliton’s linear stability by
adding arbitrary time-dependent complex perturbations $\delta\phi(t)$ to them,
so the field becomes $\tilde{\phi}_{a}(t)=\phi_{a}+\delta\phi_{a}(t)$. Without
loss of generality, we write the complex perturbations as
$\delta\phi_{a}=v_{a}e^{{\mathrm{i}}\lambda
t}+w_{a}^{*}e^{-{\mathrm{i}}\lambda^{*}t}$ [47]. Substituting
$\tilde{\phi}_{a}(t)$ into (4) and expanding up to linear order in
$\delta\phi_{a}$, we obtain the linear equation for the fluctuations around
the soliton $\phi$

$\lambda\begin{pmatrix}v\\\ w\end{pmatrix}=M(\phi)\begin{pmatrix}v\\\
w\end{pmatrix},\quad M(\phi)\equiv\begin{pmatrix}\mathcal{A}&\mathcal{B}\\\
-\mathcal{B}^{*}&-\mathcal{A}^{*}\end{pmatrix},$ (6)

where $v=[v_{(1,A)},v_{(1,B)},\dots,v_{(L,A)},v_{(L,B)}]^{T}$ is the vector
containing the component $\\{v_{a}\\}$ for all choices of $a$ in the lattice
basis, and similarly for $w$. The $\mathcal{A}$ and $\mathcal{B}$ matrices
forming the non-Hermitian stability matrix $M(\phi)$ are given in Eq. 18 of
the Supplemental Material. The linear stability of the soliton $\phi$ is
determined by the spectrum of $M(\phi)$: if at least one $\lambda$ has a
positive imaginary part, the corresponding perturbation $\delta\phi(t)$ grows
exponentially as a function of $t$, indicating that the original solution is
unstable. By construction, the stability matrix enjoys the ramified version of
particle-hole symmetry for non-Hermitian systems (called
$\mathrm{PHS}^{\dagger}$), $\Gamma M(\phi)^{*}\Gamma=-M(\phi)$, with
$\Gamma=[0\;1;1\;0]$ and as a consequence, the eigenvalues of $M(\phi)$ must
come in pairs $(\lambda,-\lambda^{*})$ or lie in the imaginary axis. Also,
since the tight-binding matrices $h$ and $f$ are real and we can choose
$\phi_{a}\in\mathbb{R}$ without loss of generality for static solitons, we get
$M(\phi)^{*}=M(\phi)$, which means that the eigenvalues must come in complex
conjugate pairs $(\lambda,\lambda^{*})$ or lie on the real axis. Together with
$\mathrm{PHS}^{\dagger}$, the reality of $M(\phi)$ implies that complex
stability eigenvalues must generally come in quadruplets
$(\lambda,\lambda^{*},-\lambda,-\lambda^{*})$. Therefore, if there is one
$\lambda$ with a negative imaginary part, there will necessarily be other
eigenvalues with positive imaginary parts. Therefore, we arrive at our
criterion: A soliton is linearly stable if and only if the spectrum of the
stability matrix $M(\phi)$ is purely real. Lastly, $\mathrm{PHS}^{\dagger}$
makes $M(\phi)$ fall into symmetry class $D^{\dagger}$ of the 38-fold
classification, which is also characterized by a $\mathbb{Z}$ topological
invariant in the presence of a point gap.

_Nonreciprocal effects –_ The spectra of the MF and stability Hamiltonians,
Eqs. 5 and 6, fall into three different regimes under PBC, illustrated in Fig.
2 (for the sake of simplicity, in the present discussion, we set the lattice
parameters to $t_{1}=1.35$, $t_{2}=0.65$, and $\Delta=0.35$ for all
simulations, and vary the squared norm $|\phi|^{2}$ of the self-consistent
solitons, which encode the nonlinearity strength). First, there is a strongly
nonlinear regime (I) at $|\phi|^{2}>1.25$ supporting static bulk solitons that
are linearly stable (all stability eigenvalues $\lambda$ are real) and whose
mean-field Hamiltonians have purely real spectra consisting of two bands
stemming from the extended eigenvectors $\psi$ and two localized states
outside the gap, the soliton $\phi$ and its chiral partner, bifurcating from
the lower and upper band, respectively. Then, at $|\phi|^{2}=1.25$, some
stability eigenvalues $\lambda$ become complex, opening _point gaps_ in the
interiors of two out of the four stability bands. As we continue to reduce
power, point gaps eventually open in intervals within both MF energy bands,
but remain closed near the outermost edges $E_{\pm}$ as long as
$|\phi|^{2}>0.94$. Thus, $0.94<|\phi|^{2}<1.25$ defines an intermediate regime
(II). When $|\phi|^{2}=0.94$, the point gaps are open at points in the real
axis arbitrarily close to $E_{\pm}$ (see Section VI of the Supplemental
Material for detailed plots of the MF spectrum close to this transition).
Below this value, we obtain a weakly nonlinear domain (III) where the
iterative sequence $\\{\phi^{(i)}\\}$ in the self-consistent mean-field
solution stops converging to machine precision and instead reaches a lower
bound $\delta\sim 10^{-3}$, suggesting the absence of _stationary_ solitons in
that regime.

Figure 2: Complex mean-field and stability spectra of static solitons as a
function of their norm squared $|\phi|^{2}$. (a) Eigenvalues $E$ (small black
dots) of the mean-field Hamiltonian $H(\phi)$ showing two bands associated
with extended states and two energies outside the gap, corresponding to the
soliton $\phi$ (ground state) and its chiral partner. (b) Eigenvalues (small
black dots) $\lambda$ of the stability matrix $M(\phi)$ stemming from the same
solitons as in figure (a). Now, there are four bands because $M$ is twice the
size of the mean-field Hamiltonian. In both (a) and (b), the gray, red, and
blue regions correspond to regimes (I), (II), and (III), and the thick dots
are the stability spectra of solitons at representative points (one in regime
(I) and two in regime (II)). Parameters used in the numerical calculations:
$t_{1}=1.35$, $t_{2}=0.65$, $\Delta=0.35$, and $L=200$ unit cells.

The complex stability spectrum indicates that the solitons in regime (II) are
linearly unstable and that two of the stability bands exhibit point gaps in
their interior. Our Supplemental Material demonstrates that the stability
spectrum has a winding number $\nu=\pm 1$ relative to any point in these gaps.
We evaluate $\nu$ employing the method from Refs. [46, 30] for determining the
winding number of U(1)-symmetric systems without translation invariance, which
involves the deformation of the stability matrix $M(\phi)$ as
$M_{ab}(\phi)\mapsto e^{{\mathrm{i}}\theta(a-b)/2L}M_{ab}(\phi)$. Then, the
quantity $\det[M_{\theta}(\phi)-\lambda_{0}]$ is periodic in $\theta$ and the
winding of the stability matrix relative to a point $\lambda_{0}\in\mathbb{C}$
is given by

$\nu(\lambda_{0})\equiv\frac{1}{2\pi{\mathrm{i}}}\int_{0}^{2\pi}{\mathrm{d}}\theta\
\partial_{\theta}\ln\det\left[M_{\theta}(\phi)-\lambda_{0}\right]\;\;\in\mathbb{Z},$
(7)

provided that $\det[M_{\theta}(\phi)-\lambda_{0}]\neq 0$ for all
$\theta\in(0,2\pi)$. We can analogously define the winding of the MF
Hamiltonian $H(\phi)$, as shown in Eq. 27 of the Supplemental Material.
Whenever the point gaps in the energy spectrum are open, we also observe that
$H(\phi)$ has a winding number $\nu=\pm 1$.

Since the winding number is odd under parity, the overall effect of growth and
decay of perturbations must favor motion in a particular direction, as in non-
Hermitian systems with a skin effect. Figure 1(a) shows the time evolution of
an initially stationary soliton in regime (II) in the presence of random
perturbations $\delta\phi_{a}$, confirming our prediction that the perturbed
soliton exhibits an overall displacement in the direction determined by the
winding number of the stability bands and the mean field spectrum below
$|\phi|^{2}<1.25$. Moreover, the squared norm $|\phi|^{2}$ is no longer
constant over time, in contrast with solitons in the usual AL model [48, 49].
Figure 1(d) shows that the norm of the unstable soliton grows when it drifts
in the preferred direction 111Solitons in regime (I) can also display the same
kind of nonreciprocal instabilities for strong enough perturbations that
become relevant beyond our linear stability analysis, i.e., to higher orders
in perturbation theory..

The most extreme behavior occurs in the weakly nonlinear regime (III).
Although not static, the resulting states have essentially the same form as
the solitons in the previous regimes. Under time evolution, they remain
localized, but now smoothly accelerate and amplify as they move in the
direction determined by the winding number [Figs. 1(b) and 1(e)]. Such a
phenomenon is consistent with the fact that, at low powers, the soliton $\phi$
inherits the properties of the extended modes $\psi$ with energies near the
band edge $E_{-}$ from which it bifurcates [50, 43].

A regime having no static bulk solutions and, at the same time, supporting
unidirectional acceleration and amplification of states is reminiscent of the
skin effect [25, 51, 52]. By further pushing the analogy with non-Hermitian
linear systems, we expect that, in regime (III), stationary solitons will be
present exclusively at one of the edges of the lattice, just as a linear
system with a skin effect has eigenstates only at one of its edges. Indeed,
when we adopt open boundary conditions (OBC) and look for MF solutions living
next to the right boundary, we obtain states that converge to machine
precision and are linearly stable. The dynamical simulation in Fig. 1(c)
confirms that such configurations remain static over long periods, and the
plot in Fig. 1(f) shows that their norm is constant – not by necessity out of
the equations of motion (1), but as a property of the solution. At the same
time, there are no converging static solitons near the left boundary.

_Discussion –_

Figure 3: Parity breaking and nonreciprocity in our model, Eq. 1. (a-c)
Circles and horizontal lines represent sites and linear couplings in the
lattice. Blue vertical bars indicate the field amplitudes $|\Phi_{R,\alpha}|$
of representative solitons. (a) $t_{1}\neq t_{2}$ (black and orange lines) and
$\Delta=0$ (white circles); (b) $t_{1}=t_{2}$ (black lines) and $\Delta\neq 0$
(white and purple circles); (c) $t_{1}\neq t_{2}$ and $\Delta\neq 0$. Only (c)
breaks all parity symmetries, generating solitons susceptible to nonreciprocal
dynamics. (d) Displacement of the center of mass (measured in unit cells) of
initially static solitons corresponding to a time evolution from $t=0$ to
$t=100$ as a function of the dimerization $t_{2}/t_{1}$ and the staggered
energy $\Delta$. Colors represent the displacement for an initial
configuration with $|\phi|^{2}=0.64$. Simulations are done with the RK4 method
and a time step $\delta t=0.02$.

Setting $t_{1}=t_{2}$ and $\Delta=0$ in our system recovers the Ablowitz-Ladik
(AL) equation,
${\mathrm{i}}{\frac{\mathrm{d}}{\mathrm{d}t}}\Phi_{R}=-(\epsilon+|\Phi_{R}|^{2})(\Phi_{R+1}+\Phi_{R-1})$,
where $R\in\mathbb{Z}$ and $\epsilon\in\mathbb{R}$, which is an integrable
discretization of the nonlinear Schödinger equation known to support stable
solitons [37, 38]. Studies have produced edge solitons in 1D lattices by
introducing unequal parameters to the left- and right-nonlinear hopping terms
of the AL equation [53, 54, 55]. Another recent work [56] generated
nonreciprocal solitons in a mechanical metamaterial using on-site
nonlinearities and a nonreciprocal linear driving mechanism, similar to the
Hatano-Nelson model [51, 52]. In both cases, reciprocity is explicitly broken.

In contrast, our work shows the emergence of topological nonreciprocity from
the interplay of linear and nonlinear terms, none of which are non-reciprocal
a priori. Since the nonlinear contribution of the AL equation for hopping from
a site $R$ to $R+1$ is proportional to $|\Phi_{R}|^{2}$ whereas that for going
in the opposite direction is proportional to $|\Phi_{R+1}|^{2}$ (Fig. 3), only
the MF Hamiltonian of a non-uniformly distributed configuration of the field
$\Phi$ – as in a soliton – becomes non-Hermitian. Even then, parity symmetry
hinders the manifestation of nonreciprocity. Indeed, all solitons in the
original (parity-symmetric) AL model are inversion-symmetric about their
center of mass [47], leading to nonlinear hopping amplitudes that cannot
imbalance their motion in a particular direction; hence, these solitons are
always reciprocal. Figure 3 illustrates that setting $t_{1}\neq t_{2}$ breaks
the parity symmetry of the lattice about sites $A$ or $B$, and $\Delta\neq 0$
does the same about the center of the unit cell. The simultaneous dimerization
and staggered on-site energies ensure that no solution with well-defined
energy will be inversion-symmetric. Then, the resulting unbalanced
distribution of the field may accelerate the state as a whole. Figure 3(d)
shows how choosing $t_{1}/t_{2}$ and $\Delta$ determines the self-induced
unidirectional propagation of initially static solitons. The lower white
regions in Fig. 3(d) represent the points of parameter space where there is no
self-acceleration as the solitons belong to regimes (I) and (II).

Regardless of the particular details of our model, there are two critical
elements for generating nonreciprocal dynamics in nonlinear periodic lattice
models: (i) non-Hermiticity, which automatically breaks TRS or
$\mathrm{TRS}^{\dagger}$, and (ii) parity symmetry breaking. Thus, more
generally, any lattice that contains the same sort of nonlinear coupling and
breaks parity symmetry could also exhibit nonreciprocal solitons. In the
Supplemental Material, we show this to be the case for the Salerno model [57],
which contains AL and Gross-Pitaevskii contributions to the nonlinear
Schödinger equation. In addition, since the nonreciprocal solitons appear in
regime (III) where $|\phi|^{2}$ is small, it is natural to ask if such states
endure as the field amplitude decreases enough to enter the quantum realm. The
quantized version of the AL model [58] is an example of a q-boson model [59,
60, 61] and the effects of parity breaking on this kind of strongly
interacting bosonic system are still unknown.

###### Acknowledgements.

_Acknowledgements –_ We thank Vir Bulchandani, Nicolas Regnault, Shinsei Ryu,
Frank Schindler, and Ajit Srivastava for discussions on related topics. W.A.B.
thanks the support of the startup funds from Emory University and by grant NSF
PHY-1748958 to the Kavli Institute for Theoretical Physics (KITP) during the
program “A Quantum Universe in a Crystal: Symmetry and Topology across the
Correlation Spectrum”.

## References

  * Sounas and Alù [2017] D. L. Sounas and A. Alù, Nature Photonics 11, 774 (2017).
  * Nagulu  _et al._ [2020] A. Nagulu, N. Reiskarimian, and H. Krishnaswamy, Nature Electronics 3, 241 (2020).
  * Jalas  _et al._ [2013] D. Jalas, A. Petrov, M. Eich, W. Freude, S. Fan, Z. Yu, R. Baets, M. Popović, A. Melloni, J. D. Joannopoulos, _et al._ , Nature Photonics 7, 579 (2013).
  * Tanaka  _et al._ [1965] S. Tanaka, N. Shimomura, and K. Ohtake, Proceedings of the IEEE 53, 260 (1965).
  * doi [1991] Polarization and crystal optics, in  _Fundamentals of Photonics_ (John Wiley & Sons, Ltd, 1991) Chap. 6, pp. 193–237.
  * Dötsch  _et al._ [2005] H. Dötsch, N. Bahlmann, O. Zhuromskyy, M. Hammer, L. Wilkens, R. Gerhardt, P. Hertel, and A. F. Popkov, JOSA B 22, 240 (2005).
  * Shoji  _et al._ [2012] Y. Shoji, M. Ito, Y. Shirato, and T. Mizumoto, Optics express 20, 18440 (2012).
  * Wang and Fan [2005] Z. Wang and S. Fan, Optics letters 30, 1989 (2005).
  * Bi  _et al._ [2011] L. Bi, J. Hu, P. Jiang, D. H. Kim, G. F. Dionne, L. C. Kimerling, and C. Ross, Nature Photonics 5, 758 (2011).
  * Kamal  _et al._ [2011] A. Kamal, J. Clarke, and M. Devoret, Nature Physics 7, 311 (2011).
  * Yu and Fan [2009] Z. Yu and S. Fan, Nature photonics 3, 91 (2009).
  * Hafezi and Rabl [2012] M. Hafezi and P. Rabl, Optics express 20, 7672 (2012).
  * Peterson  _et al._ [2019] C. W. Peterson, W. A. Benalcazar, M. Lin, T. L. Hughes, and G. Bahl, Phys. Rev. Lett. 123, 063901 (2019).
  * Kim  _et al._ [2015] J. Kim, M. C. Kuzyk, K. Han, H. Wang, and G. Bahl, Nature Physics 11, 275 (2015).
  * Wang and Goldman [1991] J. K. Wang and V. J. Goldman, Phys. Rev. Lett. 67, 749 (1991).
  * MacDonald [1990] A. H. MacDonald, Phys. Rev. Lett. 64, 220 (1990).
  * Hatsugai [1993] Y. Hatsugai, Phys. Rev. Lett. 71, 3697 (1993).
  * Haldane [1988] F. D. M. Haldane, Phys. Rev. Lett. 61, 2015 (1988).
  * Martinez Alvarez  _et al._ [2018] V. M. Martinez Alvarez, J. E. Barrios Vargas, and L. E. F. Foa Torres, Phys. Rev. B 97, 121401(R) (2018).
  * Rotter [2007] I. Rotter, Non-hermitian hamilton operator in open quantum systems (2007), arXiv:0711.2926 [quant-ph] .
  * Song  _et al._ [2019] F. Song, S. Yao, and Z. Wang, Phys. Rev. Lett. 123, 170401 (2019).
  * Silveirinha [2019] M. G. Silveirinha, Phys. Rev. B 99, 125155 (2019).
  * Kunst  _et al._ [2018] F. K. Kunst, E. Edvardsson, J. C. Budich, and E. J. Bergholtz, Phys. Rev. Lett. 121, 026808 (2018).
  * Chang  _et al._ [2023] C.-Z. Chang, C.-X. Liu, and A. H. MacDonald, Rev. Mod. Phys. 95, 011002 (2023).
  * Yao and Wang [2018] S. Yao and Z. Wang, Phys. Rev. Lett. 121, 086803 (2018).
  * Gong  _et al._ [2018] Z. Gong, Y. Ashida, K. Kawabata, K. Takasan, S. Higashikawa, and M. Ueda, Phys. Rev. X 8, 031079 (2018).
  * Hu and Zhao [2021] H. Hu and E. Zhao, Phys. Rev. Lett. 126, 010401 (2021).
  * Longhi [2022] S. Longhi, Phys. Rev. B 105, 245143 (2022).
  * Zhang  _et al._ [2023] Z.-X. Zhang, J. Cao, J.-Q. Li, Y. Zhang, S. Liu, S. Zhang, and H.-F. Wang, Phys. Rev. B 108, 125402 (2023).
  * Kawabata  _et al._ [2022] K. Kawabata, K. Shiozaki, and S. Ryu, Phys. Rev. B 105, 165137 (2022).
  * Regnault and Bernevig [2011] N. Regnault and B. A. Bernevig, Physical Review X 1, 021014 (2011).
  * Sheng  _et al._ [2011] D. Sheng, Z.-C. Gu, K. Sun, and L. Sheng, Nature Communications 2, 10.1038/ncomms1380 (2011).
  * Neupert  _et al._ [2011] T. Neupert, L. Santos, C. Chamon, and C. Mudry, Phys. Rev. Lett. 106, 236804 (2011).
  * Wang  _et al._ [2011] Y.-F. Wang, Z.-C. Gu, C.-D. Gong, and D. N. Sheng, Phys. Rev. Lett. 107, 146803 (2011).
  * Kartashov  _et al._ [2011] Y. V. Kartashov, B. A. Malomed, and L. Torner, Rev. Mod. Phys. 83, 247 (2011).
  * Lederer  _et al._ [2008] F. Lederer, G. I. Stegeman, D. N. Christodoulides, G. Assanto, M. Segev, and Y. Silberberg, Physics Reports 463, 1 (2008).
  * Ablowitz and Ladik [2008a] M. J. Ablowitz and J. F. Ladik, Journal of Mathematical Physics 16, 598 (2008a).
  * Ablowitz and Ladik [2008b] M. J. Ablowitz and J. F. Ladik, Journal of Mathematical Physics 17, 1011 (2008b).
  * Li and Wan [2022] H. Li and S. Wan, Phys. Rev. B 106, L241112 (2022).
  * Konotop and Takeno [1997] V. V. Konotop and S. Takeno, Phys. Rev. B 55, 11342 (1997).
  * Kundu [2000] K. Kundu, Phys. Rev. E 61, 5839 (2000).
  * Aceves  _et al._ [1996] A. B. Aceves, C. De Angelis, T. Peschel, R. Muschall, F. Lederer, S. Trillo, and S. Wabnitz, Phys. Rev. E 53, 1172 (1996).
  * Schindler  _et al._ [2023] F. Schindler, V. B. Bulchandani, and W. A. Benalcazar, Nonlinear breathers with crystalline symmetries (2023), arXiv:2309.07244 [cond-mat.mes-hall] .
  * Flach and Willis [1998] S. Flach and C. Willis, Physics Reports 295, 181 (1998).
  * Rice and Mele [1982] M. J. Rice and E. J. Mele, Phys. Rev. Lett. 49, 1455 (1982).
  * Kawabata  _et al._ [2019] K. Kawabata, K. Shiozaki, M. Ueda, and M. Sato, Phys. Rev. X 9, 041015 (2019).
  * Kevrekidis [2009] P. G. Kevrekidis,  _The Discrete Nonlinear Schrodinger Equation_ (Springer, 2009).
  * Wang and Wen [2020] H.-T. Wang and X.-Y. Wen, Applied Mathematics Letters 100, 106013 (2020).
  * Xie  _et al._ [2017] X.-Y. Xie, B. Tian, J. Chai, X.-Y. Wu, and Y. Jiang, Optical and Quantum Electronics 49, 1 (2017).
  * Jürgensen  _et al._ [2023] M. Jürgensen, S. Mukherjee, C. Jörg, and M. C. Rechtsman, Nature Physics 19, 420 (2023).
  * Hatano and Nelson [1996] N. Hatano and D. R. Nelson, Phys. Rev. Lett. 77, 570 (1996).
  * Hatano and Nelson [1997] N. Hatano and D. R. Nelson, Phys. Rev. B 56, 8651 (1997).
  * Lang  _et al._ [2021] L.-J. Lang, S.-L. Zhu, and Y. D. Chong, Phys. Rev. B 104, L020303 (2021).
  * Ezawa [2022] M. Ezawa, Phys. Rev. B 105, 125421 (2022).
  * Yuce [2021] C. Yuce, Physics Letters A 408, 127484 (2021).
  * Veenstra  _et al._ [2024] J. Veenstra, O. Gamayun, X. Guo, A. Sarvi, C. V. Meinersen, and C. Coulais, Nature 627, 528 (2024).
  * Salerno [1992] M. Salerno, Phys. Rev. A 46, 6856 (1992).
  * Chowdhury and Dasgupta [1993] A. R. Chowdhury and N. Dasgupta, Journal of Physics A: Mathematical and General 26, L583 (1993).
  * Arik and Coon [2008] M. Arik and D. D. Coon, Journal of Mathematical Physics 17, 524 (2008).
  * Biedenharn [1989] L. C. Biedenharn, Journal of Physics A: Mathematical and General 22, L873 (1989).
  * Macfarlane [1989] A. J. Macfarlane, Journal of Physics A: Mathematical and General 22, 4581 (1989).


# Noise-resilient phase estimation with randomized compiling

Yanwu Gu<EMAIL_ADDRESS>Beijing Academy of Quantum Information Sciences,
Beijing 100193, China  State Key Laboratory of Low Dimensional Quantum
Physics, Department of Physics, Tsinghua University, Beijing, 100084, China
Yunheng Ma  Beijing Academy of Quantum Information Sciences, Beijing 100193,
China  State Key Laboratory of Low Dimensional Quantum Physics, Department of
Physics, Tsinghua University, Beijing, 100084, China  Nicolò Forcellini
Beijing Academy of Quantum Information Sciences, Beijing 100193, China  Dong
E. Liu<EMAIL_ADDRESS>State Key Laboratory of Low Dimensional
Quantum Physics, Department of Physics, Tsinghua University, Beijing, 100084,
China  Beijing Academy of Quantum Information Sciences, Beijing 100193, China
Frontier Science Center for Quantum Information, Beijing 100184, China  Hefei
National Laboratory, Hefei 230088, China

###### Abstract

We develop an error mitigation method for the control-free phase estimation.
We prove a theorem that under the first-order correction, the phases of a
unitary operator are immune to the noise channels with only Hermitian Kraus
operators, and therefore, certain benign types of noise for phase estimation
are identified. By further incorporating the randomized compiling protocol, we
can convert the generic noise in the phase estimation circuits into stochastic
Pauli noise, which satisfies the condition of our theorem. Thus we achieve a
noise-resilient phase estimation without any quantum resource overhead. The
simulated experiments show that our method can significantly reduce the
estimation error of the phases by up to two orders of magnitude. Our method
paves the way for the utilization of quantum phase estimation before the
advent of fault-tolerant quantum computers.

Introduction.— Quantum phase estimation (QPE) Nielsen and Chuang (2004); Cleve
_et al._ (1998) is a crucial component in quantum algorithms which are
believed to achieve an exponential speedup over their classical counterparts
for solving problems, such as integer factorization Shor (1994), linear
systems of equations Harrow  _et al._ (2009), the Hamiltonian spectrum Abrams
and Lloyd (1999). However, the conventional phase estimation algorithm, which
is based on the quantum Fourier transform (QFT), requires ancillary qubits to
perform the controlled unitary operators and quantum error correction (QEC)
Shor (1996); Aharonov and Ben-Or (2008); Preskill (1998); Knill  _et al._
(1998); Kitaev (2003); Fowler  _et al._ (2012) to combat noise. These
requirements surpass the capability of the current noisy intermediate-scale
quantum (NISQ) devices Preskill (2018); Arute  _et al._ (2019); Wu and et al.
(2021); Pino  _et al._ (2021).

Regarding the controlled unitary operations in the conventional schemes, their
implementation typically necessitates the employment of many more native gates
than the unitary itself. Some variants of the phase estimation algorithm have
been devised to reduce the number of control qubits to one Kitaev (1995);
Aspuru-Guzik  _et al._ (2005); Knill  _et al._ (2007); Higgins  _et al._
(2007); Svore  _et al._ (2014); Wiebe and Granade (2016); O’Brien  _et al._
(2019a); Somma (2019). Furthermore, the control-free phase estimation Kimmel
_et al._ (2015); Roushan  _et al._ (2017); Russo  _et al._ (2021); Neill  _et
al._ (2021); Lu  _et al._ (2021); Zintchenko and Wiebe (2016); O’Brien  _et
al._ (2021); Lin and Tong (2022); Dutkiewicz  _et al._ (2022) is a more
appealing variant in the NISQ era, where the unitary of interest is repeatedly
applied instead of its control version. Control-free phase estimation was
first applied to calibrate the unitary errors in single-qubit gates as “robust
phase estimation” Kimmel  _et al._ (2015), and later adapted to two-qubit gate
calibration as “Floquet calibration” Neill  _et al._ (2021).

Noise affects all quantum operations in NISQ devices. Understanding and
correcting the effect of noise is necessary for meaningful quantum
computations. In the NISQ era, quantum error mitigation (QEM) O’Brien  _et
al._ (2021); Kandala  _et al._ (2019); Temme  _et al._ (2017); Li and Benjamin
(2017); Temme  _et al._ (2017); Endo  _et al._ (2018); McClean  _et al._
(2017); Bonet-Monroig  _et al._ (2018); Sagastizabal  _et al._ (2019); McArdle
_et al._ (2019); McClean  _et al._ (2020); Arute  _et al._ (2020); Mi  _et
al._ (2021, 2022); Montanaro and Stanisic (2021) is a more feasible technique
than the resource-expensive QEC. Among the QEM schemes, verified phase
estimation (VPE) O’Brien  _et al._ (2021) is a technique that is directly
linked to the QPE. This method achieves error mitigation by postselecting the
shots in which the system register is recovered to its initial state. Due to
the postselection, each circuit requires many extra repetitions to collect
enough statistics. The number of repetitions is inversely proportional to the
circuit fidelity, which decreases exponentially as the number of gate
operations increases. In addition, VPE cannot mitigate the unitary noise,
which is one of the most common types of noise in multi-qubit circuits, due to
ineluctable miscalibrations and crosstalks Sarovar  _et al._ (2020). As a
result, it is urgently necessary to develop a more practical scheme for
general noise kinds that uses less resources in NISQ devices.

In this work, we develop an error mitigation technique for control-free phase
estimation using randomized compiling, with no resource overhead other than
some efficient classical computations. We theoretically analyze the noise
types that the phase estimation is insensitive to. We prove a theorem that if
all Kraus operators of the associated noise channel are Hermitian, the phases
extracted from noisy QPE circuits do not change under the first-order
correction in comparison to those from ideal circuits. In addition, the
randomized compiling method is applied to convert general noise types
including coherent noise in the circuit to stochastic Pauli noise, which fits
the criterion in our theorem. As a result, we achieve a practical error
reduction in phase estimation. The validity of our error mitigation method is
tested using two simulated experiments. The results show that our method can
reduce errors in both unitary and stochastic noise, particularly in unitary
noise by up to two orders of magnitude.

Figure 1: The circuit structure of control-free phase estimation. The initial
state is prepared by applying a unitary operator $U_{s}$. The unitary $U$ is
then repeated for $L$ times. Finally, the expectation value of an operator $O$
is measured by performing a measurement unitary $U_{O}$ prior to the
computational basis measurement.

Control-free phase estimation.— Let us briefly review the control-free phase
estimation Kimmel  _et al._ (2015); Roushan  _et al._ (2017); Russo  _et al._
(2021); Neill  _et al._ (2021) where its circuit structure is shown in Fig. 1.
The spectral decomposition of the unitary operator $U$ is

$U=\sum_{a}e^{i\lambda_{a}}|\phi_{a}\rangle\langle\phi_{a}|$ (1)

where $\lambda_{a}$ is a phase and $|\phi_{a}\rangle$ is the corresponding
eigenstate. The system is prepared to a state
$|\psi\rangle=\sum_{a}c_{a}|\phi_{a}\rangle$ by conducting a unitary operator
$U_{s}$. To obtain the phases of interest, the initial state $|\psi\rangle$
should be carefully chosen so that it only has a few relevant eigenstates of
$U$. The unitary operator $U$ is then repeatedly applied to $|\psi\rangle$ for
$L$ times, where $L$ is an integer ranging from $1$ to $L_{\textrm{max}}$.
Finally, an operator $O$ is measured by performing a unitary $U_{O}$ before
the computational basis measurement. The expectation value of $O$ is

$\displaystyle{\left\langle{O}\right\rangle}_{L}$ $\displaystyle=$
$\displaystyle\langle\psi|(U^{\dagger})^{L}OU^{L}|\psi\rangle$ (2)
$\displaystyle=$
$\displaystyle\sum_{a,b}c_{a}c^{*}_{b}\langle\phi_{b}|O|\phi_{a}\rangle
e^{i(\lambda_{a}-\lambda_{b})L}\,.$

The difference $\lambda_{a}-\lambda_{b}$ can be retrieved by conducting a
discrete Fourier transform on the vector of
${\left\langle{O}\right\rangle}_{L}$ or by performing a function fitting.
Control-free phase estimation, like randomized benchmarking Knill  _et al._
(2008); Magesan  _et al._ (2011), is robust to state preparation and
measurement (SPAM) errors, because SPAM errors only affect the coefficients
$c_{a}c_{b}^{*}\langle\phi_{b}|O|\phi_{a}\rangle$ but not phases (this is also
true for other phase estimation algorithms with repeated control version of
target unitary).

To get the individual phases, we assume a reference eigenstate
$|\phi_{0}\rangle$ exists with known phase $\lambda_{0}$. Then, we prepare an
initial state
$|\psi\rangle=\frac{1}{\sqrt{2}}\left(|\phi_{0}\rangle+|t\rangle\right)$,
where $|t\rangle=\sum_{n=1}^{N_{p}}c_{n}|\phi_{n}\rangle$ is a state
containing $N_{p}$ eigenstates whose phases are to be estimated. The
measurement operator is chosen as non-Hermitian $O=2|\phi_{0}\rangle\langle
t|$, which can be always decomposed as the sum of Hermitian operators Nielsen
and Chuang (2004). In this case, the expectation value reads

$\displaystyle{\left\langle{O}\right\rangle}_{L}$ $\displaystyle=$
$\displaystyle 2\langle\psi|(U^{\dagger})^{L}|\phi_{0}\rangle\langle
t|U^{L}|\psi\rangle$ (3) $\displaystyle=$
$\displaystyle\sum_{n}c_{n}c^{*}_{n}e^{i(\lambda_{n}-\lambda_{0})L}\,.$

Finally, the measurements of the decomposed Hermitian operators yield the
phases $\\{\lambda_{n}\\}$.

Benign type of noise for phase estimation.— Any quantum algorithms running on
the current devices are affected by noise. If we can identify the noise types
that are mild to phase estimation, we may tailor the noise in quantum devices
to the desired type and achieve error mitigation.

We use the language of quantum channels (represented with Calligraphical
symbols) to describe the noise effect. The corresponding unitary channel
$\mathcal{U}$ of the unitary operator $U$ in Eq. (1) has the effect

$\mathcal{U}\left(|\phi_{a}\rangle\langle\phi_{b}|\right)=U|\phi_{a}\rangle\langle\phi_{b}|U^{\dagger}=e^{i(\lambda_{a}-\lambda_{b})}|\phi_{a}\rangle\langle\phi_{b}|\,.$
(4)

Thus $|\phi_{a}\rangle\langle\phi_{b}|$ is an eigen-operator of the unitary
channel $\mathcal{U}$ with eigenvalue $e^{i(\lambda_{a}-\lambda_{b})}$. If a
noisy channel $\mathcal{E}(\rho)=\sum_{k}E_{k}\rho E_{k}^{\dagger}$ is
appended to the unitary channel $\mathcal{U}$, then the resultant noisy gate
is $\widetilde{\mathcal{U}}=\mathcal{E}\,\mathcal{U}$. In Sec.II of
supplementary information (SI) sup , we prove that the expectation value in
Eq. (3) is transformed into

${\left\langle{O}\right\rangle}_{L}\approx\sum_{n}p_{n}\,(g_{n0})^{L}\,e^{i\lambda_{n0}L}$
(5)

under weak noise ${\mathcal{E}}$, where $p_{n}\approx c_{n}c_{n}^{*}$ is a
real number close to the proportion of $|\phi_{n}\rangle$ in $|t\rangle$. The
eigenvalue of the noisy operation $\widetilde{\mathcal{U}}$, that is deviated
from the ideal unitary operation ${\mathcal{U}}$, is modified to
$g_{n0}e^{i\lambda_{n0}}$ with the noisy amplitude $g_{n0}$ and phase
$\lambda_{n0}$, respectively. The amplitude has the constraint $g_{n0}\leq 1$
for completely positive trace preserving (CPTP) maps describing physical
channels Wolf (2012); Watrous (2018). Thus Eq. (5) is a damping oscillating
model. One can still conduct a Fourier transform or fit Eq. (5) to obtain an
estimation of the noisy phases $\lambda_{n0}$. To reduce the inaccuracy in the
phase estimation, we can identify noise types benign for the phase estimation,
and then design a scheme to convert the general noise to those desired types.
Following this strategy, we first provide a theorem below.

###### Theorem 1

If every Kraus operator $E_{k}$ of a noise $\mathcal{E}$ is Hermitian, then
the noisy version $\widetilde{\mathcal{U}}=\mathcal{E}\,\mathcal{U}$ of a
unitary channel $\mathcal{U}$ keeps the phases unchanged up to the leading
order of noise strength.

This theorem is based on the first-order perturbation theory (see Sec.I of SI
sup  for the detailed proof, which includes Refs. Sakurai and Napolitano
(2010); Kato (1995); Sarkar and Pereira (1995); Potts and Tasche (2013);
Helsen  _et al._ (2019); O’Brien  _et al._ (2019b)). The proof makes no
assumptions regarding the form of the unitary operator $U$, and therefore, the
theorem applies to any phase estimation problem. Our theorem provides a
sufficient condition regarding the benign types of noise, which already
encompasses a variety of noise types, including stochastic Pauli noise.
Particularly, the phase damping channel ($T_{2}$ error) is a stochastic Pauli
channel.

The noise in a real device is exceedingly complex and rarely meets the
criterion in Theorem 1. For example, the amplitude damping channel ($T_{1}$
error) and the unitary noise cannot satisfy the condition of our theorem.
Fortunately, there is a technique known as randomized compiling (RC) Wallman
and Emerson (2016); Hashim  _et al._ (2021) that can turn the noise in the
corresponding circuit of $\mathcal{U}$ into stochastic Pauli noise. Note RC
can only convert the noise in each cycle of $U$ to stochastic Pauli noise, but
we suppose the noise channel ${\mathcal{E}}$ occurs at the end of $U$ in
Theorem 1. As a result, the condition in Theorem 1 is stronger than what RC
can achieve; see the detailed discussion below.

The effect of randomized compiling.— Randomized compiling (RC) Wallman and
Emerson (2016); Hashim  _et al._ (2021) is a technique to reduce the general
and complex noise in circuits to a specific simple noise type, namely
stochastic Pauli noise. Initially, the bare circuit $U$ is partitioned into
many cycles and each cycle contains a layer of single-qubit gates and a layer
of two-qubit gates, that is

$U=G_{K}C_{K}\cdots G_{k}C_{k}\cdots G_{1}C_{1}$ (6)

where $C_{k}$ and $G_{k}$ represent single-qubit and two-qubit layers,
respectively. Each layer of single-qubit gates $C_{k}$ is replaced with a
round of randomized dressed gates $\widetilde{C}_{k}=T_{k}C_{k}T_{k-1}^{c}$
where $T_{k}$ are chosen uniformly at random from the Pauli group and the
correction operators are set to
$T_{k}^{c}=G_{k}T_{k}^{\dagger}G_{k}^{\dagger}$ Wallman and Emerson (2016);
Hashim  _et al._ (2021). In the RC protocol, each dressed gate layer
$\widetilde{C}_{k}$ should be implemented as a single layer of elementary
gates. Thus, randomized compiling can create a set of randomized circuits that
are logically equivalent to the original circuit without increasing circuit
depth. The noise in circuits is customized to stochastic Pauli noise after
averaging the outputs of these randomized circuits.

For phase estimation in Fig. 1, we perform RC on the circuit segments $U^{L}$.
Running these randomized circuits with the same initial state and final
measurement, and taking the average of results tailor the noise in each cycle
to stochastic Pauli noise. Now the actual operations by the circuit $U^{L}$
can be written as

$\left(\mathcal{G}_{K}\mathcal{E}_{p}^{(K)}\mathcal{C}_{K}\cdots\mathcal{G}_{1}\mathcal{E}_{p}^{(1)}\mathcal{C}_{1}\right)^{L}=\left(\mathcal{E^{\prime}}\,\mathcal{U}\right)^{L}$
(7)

where $\mathcal{E}_{p}$ is a stochastic Pauli channel. The noise
$\mathcal{E}^{\prime}$ is defined by moving all the $\mathcal{E}_{p}$s to the
end of each segment $U$. We ask whether the composite noise channel
${\mathcal{E}}^{\prime}$ after RC fits the condition in Theorem 1?

If $U$ is a Clifford circuit, the final error is still a Pauli error, which
meets our criterion. However, if the circuit $U$ includes some non-Clifford
gates, the condition in Theorem 1 is generally not satisfied. For the
stochastic Pauli noise after the RC protocol, we can perform a probabilistic
analysis. First, if we assume the probability of error in each cycle is
relatively small, there is a high possibility that only one Pauli error
(denoted as $P$) occurs in $U$. We assume the unitary operator of the circuit
segments before and after the error $P$ to be $U_{1}$ and $U_{2}$, and thus
the ideal unitary is $U=U_{2}U_{1}$. With error $P$, the whole operation
becomes $U_{2}PU_{1}$. Let us move $P$ to the end of $U$, and obtain
$U_{2}PU_{2}^{\dagger}U_{2}U_{1}=(U_{2}PU_{2}^{\dagger})U$. Now, the new error
operator $U_{2}PU_{2}^{\dagger}$ is Hermitian, which clearly satisfies the
condition of Theorem 1. For the cases of two errors and beyond, many cases
still satisfy the criteria in Theorem 1, for example, if the circuit segments
between the first $P_{1}$ and the last $P_{n}$ error are all Clifford.

Now, we provide more rigorous results for the effect of RC (see detailed
analysis in Sec.III of SI sup ). We show that the phase error in bare circuits
is proportional to the $\|\mathcal{E}-\mathcal{I}\|_{\diamond}$, the diamond
norm distance between the noise ${\mathcal{E}}$ in bare circuit and identity
channel ${\mathcal{I}}$. After RC, the phase error is proportional to the
$\|\mathcal{E}^{\prime}-\mathcal{I}\|_{\diamond}^{\alpha}$ where $1<\alpha\leq
2$. For Clifford circuit $U$, $\alpha=2$; for the non-Clifford circuit,
$\alpha<2$ because the criterion in Theorem 1 can only partially be met. If
the gate noise in bare circuit $U$ is stochastic with a characteristic noise
probability $p$, the phase error is $\sim p$ in the bare circuits and $\sim
p^{\alpha}$ in the RC circuits. For the case of unitary noise with the
characteristic rotation angle $\theta$, the phase error in the bare circuits
is $\sim\theta$. After RC, the unitary error is converted into stochastic
Pauli noise with some noise probability $p\sim\theta^{2}$ Kliesch and Roth
(2021). The phase error by RC circuits should be $\sim
p^{\alpha}\sim\theta^{2\alpha}$. These results show that our method can reduce
the estimation error of the phases for both stochastic and unitary noise,
particularly has a stronger effect on unitary noise. Note the theoretical
scaling may not be obtained in actual experiments due to the finite number of
shots, random circuits, and repetitions of $U$. Nevertheless, our method shows
strong error mitigation power under practical experimental settings as shown
by the simulated experiments.

Simulated experiments.—To demonstrate the performance of our method, we
present two simulated experiments: 1) The estimation of the quasi-energies of
a Floquet system Neill  _et al._ (2021), 2) an order finding problem Nielsen
and Chuang (2004). The measured signals
${\left\langle{X}\right\rangle}_{L}+i{\left\langle{Y}\right\rangle}_{L}$ are
first Fourier transformed to the frequency domain, and then the peak locations
in the frequency spectrum provide a rough estimate of phases. To obtain a more
accurate estimate, one can fit the data to the expression in Eq. (5). The
estimation error is defined as the average distance between the estimates
$\widehat{\lambda}_{n}$ using our method and the actual phases $\lambda_{n}$,
i.e.,
$\frac{1}{N}\sum_{n=1}^{N}\left|\widehat{\lambda}_{n}-\lambda_{n}\right|$
where $N$ is the number of non-degenerate actual phases. See Sec. IV of SI sup
for the circuits and more simulation details.

Firstly, we use a 6-qubit version of the Floquet system in Ref. Neill  _et
al._ (2021) to validate the phase estimation error scaling with the noise
strength after RC. But this scaling can only be obtained in a very ideal case,
i.e., an infinite number of shots $N_{s}$ (to eliminate sampling error), an
infinite number of random circuits $N_{r}$ for each bare circuit (to perfectly
transform noise into stochastic Pauli noise) and enough number of repetitions
$L$ of target unitary (to remove the possible effect of not enough data).
Thus, we directly compute the expectation values from the final density matrix
of each circuit. Furthermore, we set the noise of each gate to be stochastic
Pauli noise and use a very large circuit length $L_{\textrm{max}}$. As shown
in Fig. 2, when the noise probability of each gate is very low, the
probability of two or more Pauli errors occurring in the circuit $U$ is
negligible, and thus the scaling power $\alpha$ approaches 2. However, when
noise probability is relatively large, the cases of two and more Pauli errors
cause $\alpha$ to be less than 2, which coincides with our theoretical
analysis.

Figure 2: Phase estimation error under stochastic Pauli noise. Each qubit at
each gate location is subjected to a stochastic Pauli channel with Pauli noise
probabilities $p_{x},p_{y},p_{z}$. Under very weak noise, the scaling power
approaches 2. But for relatively large noise, the scaling power is less than 2
because the non-negligible cases of two or more Pauli errors in circuit $U$
may violate the condition in Theorem 1.

Then, we show the performance of our method under practical experimental
settings, i.e., a finite number of $N_{s},N_{r}$ and a small
$L_{\textrm{max}}$. We add unitary noise or stochastic noise (with non-
Hermitian Kraus operators) to each gate. The shot for each bare circuit is set
to $N_{s}$. We use RC to generate $N_{r}=20$ randomized circuits for each bare
circuit. To maintain the same resource cost, we set the shot for each
randomized circuit to $N_{s}/N_{r}$.

Figure 3: Simulated results on estimating Floquet quasi-energies Neill  _et
al._ (2021) ($10$-qubit system), under unitary noise (a) and stochastic noise
(b). At small $\theta$, the estimation error decreases slowly due to a small
number of $N_{r}$. One can increase $N_{r}$ to get more error mitigation power
because the phase estimation error scales with $N_{r}$ as $1/\sqrt{N_{r}}$,
see Sec. IV of SI sup .

In Fig. 3, we show the results of simulated experiments on estimating Floquet
quasi-energies with $N_{s}=10^{7},L_{\textrm{max}}=50$. The unitary operation
$U$ is a 10-qubit non-Clifford unitary operator. As shown in Fig. 3(a), our
technique has strong error mitigation power for unitary noise. The estimation
error is reduced by up to two orders of magnitude. If there is no error
mitigation, it is linearly increased with the rotation angle of unitary noise.
For the error mitigated case in Fig. 3(a), we fit the results in the strong
noise regime with a power law function and obtain its power of almost $2.73$.
Both the scaling behaviors of estimation error with or without error
mitigation are consistent with our analytical analysis. The absence of
theoretical scaling behaviour in the weak noise regime is a result of the
small number of random circuits $N_{r}$. In Sec. IV of SI sup , we show the
phase estimation error scales with $N_{r}$ as $1/\sqrt{N_{r}}$. In actual
experiments, $N_{r}$ can be increased up to $N_{s}$ in order to increase error
mitigation capability without incurring experimental overhead. Though not as
significant as in the case of unitary noise, our method provides a nearly
$60\%$ reduction in estimation error for stochastic noise as shown in Fig.
3(b).

Figure 4: Simulated results of order finding $x=4$ and $N=255$, for two types
of noise: unitary noise (a)(b) and stochastic noise (c)(d). The correct peaks
in the frequency domain are located at 0, 25, 50, 75. Some extra peaks appear
in (a)(c) due to the splitting of degenerated eigenvalues under noise. Our
method completely eliminates these extra peaks and the remaining peaks are
located at the correct positions in (b)(d).

Finally, we consider the order finding problem, that is to find the least
positive integer $r$ such that for two specified coprime numbers $x$ and $N$
$(x<N)$ we have $x^{r}=1(\mathrm{mod}\,N)$. Here we simulate an order finding
problem with $x=4$ and $N=255$. The corresponding unitary is an 8-qubit
Clifford unitary with only SWAP gates sup . The main numerical results are
shown in Fig. 4 (with $N_{s}=10^{5}$), where we present the Fourier transform
of the measured signals under unitary noise $(a,b)$ and stochastic noise
$(c,d)$. The maximum number of repetitions of the unitary of interest is
$L_{\textrm{max}}=100$. In this particular order finding problem, the order
should be $r=4$, which induces four peaks at 0, 25, 50, 75 in the frequency
domain. Under unitary noise, degenerated eigenvalues split, resulting in extra
peaks in Fig. 4(a). With our error mitigation method, these extra peaks are
completely removed and the surviving peaks are located at the correct
positions in Fig. 4(b). Under stochastic noise, some peaks are almost buried
by noise in Fig. 4(c). However, our method recovers the right peaks from noise
in Fig. 4(d).

Conclusion— We have developed an efficient method to mitigate the errors in
control-free phase estimation. The method is based on Theorem 1 that the noise
with Hermitian Kraus operators causes only mild changes in the phases of a
unitary operator. To achieve the desired type of noise in actual devices,
randomized compiling can be used, which has no resource overhead. The two
simulated experiments on Hamiltonian eigenvalue estimation and order finding
demonstrate the strong error mitigation power of our method. We emphasize that
our scheme requires an assumption that the noise is time independent, and this
assumption is reasonable, at least in the “Good devices” Arute  _et al._
(2019); Wu and et al. (2021); Neill  _et al._ (2021). In addition, the effect
of the noise drift can be controlled and calibrated, for example, using the
method Proctor  _et al._ (2020). Finally, our method can be readily applied to
the phase estimation scheme with a single control qubit.

###### Acknowledgements.

Acknowledgments.— This work was supported by the Innovation Program for
Quantum Science and Technology (Grant No. 2021ZD0302400), the National Natural
Science Foundation of China (Grant No. 12147123 and 11974198), and the Beijing
Natural Science Foundation (No. Z220002). The source code for the simulated
experiments is available at this site https://github.com/yanwu-gu/noise-
resilient-phase-estimation.

## References

  * Nielsen and Chuang (2004) M. A. Nielsen and I. L. Chuang,  _Quantum computation and quantum information_ , 1st ed. (Cambridge University Press, 2004).
  * Cleve  _et al._ (1998) R. Cleve, A. Ekert, C. Macchiavello, and M. Mosca, Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences 454, 339 (1998).
  * Shor (1994) P. Shor, in  _Proceedings 35th Annual Symposium on Foundations of Computer Science_ (1994) pp. 124–134.
  * Harrow  _et al._ (2009) A. W. Harrow, A. Hassidim, and S. Lloyd, Phys. Rev. Lett. 103, 150502 (2009).
  * Abrams and Lloyd (1999) D. S. Abrams and S. Lloyd, Physical Review Letters 83, 5162 (1999).
  * Shor (1996) P. Shor, in  _Proceedings of 37th Conference on Foundations of Computer Science_ (1996) pp. 56–65.
  * Aharonov and Ben-Or (2008) D. Aharonov and M. Ben-Or, SIAM Journal on Computing 38, 1207 (2008).
  * Preskill (1998) J. Preskill, Proc. R. Soc. London, Ser. A 454, 385 (1998).
  * Knill  _et al._ (1998) E. Knill, R. Laflamme, and W. H. Zurek, Proc. R. Soc. London, Ser. A 454, 365 (1998).
  * Kitaev (2003) A. Kitaev, Annals of Physics 303, 2 (2003).
  * Fowler  _et al._ (2012) A. G. Fowler, M. Mariantoni, J. M. Martinis, and A. N. Cleland, Phys. Rev. A 86, 032324 (2012).
  * Preskill (2018) J. Preskill, Quantum 2, 79 (2018).
  * Arute  _et al._ (2019) F. Arute, K. Arya, R. Babbush, and et al., Nature 574 (2019), 10.1038/s41586-019-1666-5.
  * Wu and et al. (2021) Y. Wu and et al., Phys. Rev. Lett. 127, 180501 (2021).
  * Pino  _et al._ (2021) J. M. Pino, J. M. Dreiling, C. Figgatt, and et al., Nature 592, 209 (2021).
  * Kitaev (1995) A. Y. Kitaev, arXiv:quant-ph/9511026 (1995).
  * Aspuru-Guzik  _et al._ (2005) A. Aspuru-Guzik, A. D. Dutoi, P. J. Love, and M. Head-Gordon, Science 309, 1704 (2005).
  * Knill  _et al._ (2007) E. Knill, G. Ortiz, and R. D. Somma, Phys. Rev. A 75, 012328 (2007).
  * Higgins  _et al._ (2007) B. L. Higgins, D. W. Berry, S. D. Bartlett, H. M. Wiseman, and G. J. Pryde, Nature 450, 393 (2007).
  * Svore  _et al._ (2014) K. M. Svore, M. B. Hastings, and M. Freedman, Quantum Info. Comput. 14, 306–328 (2014).
  * Wiebe and Granade (2016) N. Wiebe and C. Granade, Phys. Rev. Lett. 117, 010503 (2016).
  * O’Brien  _et al._ (2019a) T. E. O’Brien, B. Tarasinski, and B. M. Terhal, New Journal of Physics 21, 023022 (2019a).
  * Somma (2019) R. D. Somma, New Journal of Physics 21, 123025 (2019).
  * Kimmel  _et al._ (2015) S. Kimmel, G. H. Low, and T. J. Yoder, Phys. Rev. A 92, 062315 (2015).
  * Roushan  _et al._ (2017) P. Roushan, C. Neill, J. Tangpanitanon, V. M. Bastidas, A. Megrant, R. Barends, Y. Chen, Z. Chen, B. Chiaro, A. Dunsworth, A. Fowler, B. Foxen, M. Giustina, E. Jeffrey, J. Kelly, E. Lucero, J. Mutus, M. Neeley, C. Quintana, D. Sank, A. Vainsencher, J. Wenner, T. White, H. Neven, D. G. Angelakis, and J. Martinis, Science 358, 1175 (2017).
  * Russo  _et al._ (2021) A. E. Russo, K. M. Rudinger, B. C. A. Morrison, and A. D. Baczewski, Phys. Rev. Lett. 126, 210501 (2021).
  * Neill  _et al._ (2021) C. Neill, T. McCourt, X. Mi, Z. Jiang, M. Niu, W. Mruczkiewicz, I. Aleiner, F. Arute, K. Arya, J. Atalaya, _et al._ , Nature 594, 508 (2021).
  * Lu  _et al._ (2021) S. Lu, M. C. Bañuls, and J. I. Cirac, PRX Quantum 2, 020321 (2021).
  * Zintchenko and Wiebe (2016) I. Zintchenko and N. Wiebe, Phys. Rev. A 93, 062306 (2016).
  * O’Brien  _et al._ (2021) T. E. O’Brien, S. Polla, N. C. Rubin, W. J. Huggins, S. McArdle, S. Boixo, J. R. McClean, and R. Babbush, PRX Quantum 2, 020317 (2021).
  * Lin and Tong (2022) L. Lin and Y. Tong, PRX Quantum 3, 010318 (2022).
  * Dutkiewicz  _et al._ (2022) A. Dutkiewicz, B. M. Terhal, and T. E. O’Brien, Quantum 6, 830 (2022).
  * Kandala  _et al._ (2019) A. Kandala, K. Temme, A. D. Córcoles, A. Mezzacapo, J. M. Chow, and J. M. Gambetta, Nature 567, 491 (2019).
  * Temme  _et al._ (2017) K. Temme, S. Bravyi, and J. M. Gambetta, Phys. Rev. Lett. 119, 180509 (2017).
  * Li and Benjamin (2017) Y. Li and S. C. Benjamin, Phys. Rev. X 7, 021050 (2017).
  * Endo  _et al._ (2018) S. Endo, S. C. Benjamin, and Y. Li, Phys. Rev. X 8, 031027 (2018).
  * McClean  _et al._ (2017) J. R. McClean, M. E. Kimchi-Schwartz, J. Carter, and W. A. de Jong, Phys. Rev. A 95, 042308 (2017).
  * Bonet-Monroig  _et al._ (2018) X. Bonet-Monroig, R. Sagastizabal, M. Singh, and T. E. O’Brien, Phys. Rev. A 98, 062339 (2018).
  * Sagastizabal  _et al._ (2019) R. Sagastizabal, X. Bonet-Monroig, M. Singh, M. A. Rol, C. C. Bultink, X. Fu, C. H. Price, V. P. Ostroukh, N. Muthusubramanian, A. Bruno, M. Beekman, N. Haider, T. E. O’Brien, and L. DiCarlo, Phys. Rev. A 100, 010302(R) (2019).
  * McArdle  _et al._ (2019) S. McArdle, X. Yuan, and S. Benjamin, Phys. Rev. Lett. 122, 180501 (2019).
  * McClean  _et al._ (2020) J. R. McClean, Z. Jiang, N. C. Rubin, R. Babbush, and H. Neven, Nature communications 11, 1 (2020).
  * Arute  _et al._ (2020) F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, A. Bengtsson, S. Boixo, M. Broughton, B. B. Buckley, _et al._ , arXiv:2010.07965 (2020).
  * Mi  _et al._ (2021) X. Mi, P. Roushan, C. Quintana, S. Mandrà, J. Marshall, C. Neill, F. Arute, K. Arya, J. Atalaya, R. Babbush, _et al._ , Science 374, 1479 (2021).
  * Mi  _et al._ (2022) X. Mi, M. Ippoliti, C. Quintana, A. Greene, Z. Chen, J. Gross, F. Arute, K. Arya, J. Atalaya, R. Babbush, _et al._ , Nature 601, 531 (2022).
  * Montanaro and Stanisic (2021) A. Montanaro and S. Stanisic, arXiv:2102.02120 [quant-ph] (2021).
  * Sarovar  _et al._ (2020) M. Sarovar, T. Proctor, K. Rudinger, K. Young, E. Nielsen, and R. Blume-Kohout, Quantum 4, 321 (2020).
  * Knill  _et al._ (2008) E. Knill, D. Leibfried, R. Reichle, J. Britton, R. B. Blakestad, J. D. Jost, C. Langer, R. Ozeri, S. Seidelin, and D. J. Wineland, Phys. Rev. A 77, 012307 (2008).
  * Magesan  _et al._ (2011) E. Magesan, J. M. Gambetta, and J. Emerson, Phys. Rev. Lett. 106, 180504 (2011).
  * (49) See Supplemental Information for more details.
  * Wolf (2012) M. M. Wolf, “Quantum channels and operations - guided tour,”  (2012), graue Literatur.
  * Watrous (2018) J. Watrous,  _The Theory of Quantum Information_ (Cambridge University Press, 2018).
  * Sakurai and Napolitano (2010) J. J. Sakurai and J. J. Napolitano,  _Modern quantum mechanics_ (Pearson, 2010).
  * Kato (1995) T. Kato,  _Perturbation theory for linear operators_ , Vol. 132 (Springer Berlin, Heidelberg, 1995).
  * Sarkar and Pereira (1995) T. Sarkar and O. Pereira, IEEE Antennas and Propagation Magazine 37, 48 (1995).
  * Potts and Tasche (2013) D. Potts and M. Tasche, Linear Algebra and its Applications 439, 1024 (2013), 17th Conference of the International Linear Algebra Society, Braunschweig, Germany, August 2011.
  * Helsen  _et al._ (2019) J. Helsen, F. Battistel, and B. M. Terhal, npj Quantum Information 5, 1 (2019).
  * O’Brien  _et al._ (2019b) T. E. O’Brien, B. Tarasinski, and B. M. Terhal, New Journal of Physics 21, 023022 (2019b).
  * Wallman and Emerson (2016) J. J. Wallman and J. Emerson, Phys. Rev. A 94, 052325 (2016).
  * Hashim  _et al._ (2021) A. Hashim, R. K. Naik, A. Morvan, J.-L. Ville, B. Mitchell, J. M. Kreikebaum, M. Davis, E. Smith, C. Iancu, K. P. O’Brien, I. Hincks, J. J. Wallman, J. Emerson, and I. Siddiqi, Phys. Rev. X 11, 041039 (2021).
  * Kliesch and Roth (2021) M. Kliesch and I. Roth, PRX Quantum 2, 010201 (2021).
  * Proctor  _et al._ (2020) T. Proctor, M. Revelle, E. Nielsen, K. Rudinger, D. Lobser, P. Maunz, R. Blume-Kohout, and K. Young, Nature communications 11, 1 (2020).

See pages 1 of NRPE_SI.pdf See pages 2 of NRPE_SI.pdf See pages 3 of
NRPE_SI.pdf See pages 4 of NRPE_SI.pdf See pages 5 of NRPE_SI.pdf See pages 6
of NRPE_SI.pdf See pages 7 of NRPE_SI.pdf See pages 8 of NRPE_SI.pdf See pages
9 of NRPE_SI.pdf See pages 10 of NRPE_SI.pdf See pages 11 of NRPE_SI.pdf


# Valveless pumping at low Reynolds numbers

Gabriel Amselem<EMAIL_ADDRESS>Christophe Clanet  Michael
Benzaquen  Laboratoire d’Hydrodynamique, CNRS, École polytechnique, Institut
Polytechnique de Paris, 91120 Palaiseau, France

###### Abstract

Pumping at low Reynolds number is a ubiquitously encountered feature, both in
biological organisms and engineered devices. Generating net flow requires the
presence of an asymmetry in the system, which traditionally comes from
geometric flow rectifiers. Here, we study a valveless system of $N$
oscillating pumps in series, where the asymmetry comes not from the geometry
but from time, that is the phase shifts between the pumps. Experimental and
theoretical results are in very good agreement. We provide the optimal phase
shifts leading to the maximal net flow in the continuous $N\rightarrow\infty$
limit. The maximal flow rate is larger by 25% than that of a traditional
peristaltic sinusoidal wave, for the same amplitude of actuation. Our results
pave the way for the design of more efficient microfluidic pumps.

##  I Introduction

Pumping at low Reynolds number is a broadly encountered feature, both in
nature and engineered systems. On the biological side, insects such as
butterflies and mosquitoes feed essentially on liquids such as plant nectar or
blood, and many rely on the action of one or more muscular pumps to intake
their food [1]. On the man-made side, the design of micropumps for MEMS
applications started with the seminal work of Smits 30 years ago [2], and
continued with the Quake-valves now ubiquitous in microfluidic designs [3],
triggering a large variety of novel designs [4, 5, 6, 7].

In one way or another, these natural and engineered devices rely on the
presence of an asymmetry in the system, which is essential to pump at low
Reynolds number. Some biological and man-made pump designs rely on valves, a
role played by the pharynx in mosquitoes [8, 9, 10]. Yet, valves are prone to
mechanical failure and clogging, and valveless systems of micropumps were
conceived to make them easier to fabricate and more resilient [5, 11]. Their
mode of operation relies on geometrically asymmetric elements having flow-
rectifying properties. Valveless pumping at low Reynolds number is also
prevalent in biology, throughout the animal kingdom: in the early stages of
embryonic development of a number of species, the heart consists in two
chambers that pump fluid without a valve [12, 13, 14].

In this communication we study experimentally and theoretically a valveless
pumping system consisting of $N$ identical contracting pumps in series,
operating at low Reynolds number. The system does not contain any flow
rectifier (valveless), and the asymmetry comes not from the built-in geometry
but from the phase shifts between the oscillating pumps, i.e. from the
dynamics of actuation. The average flow rate through the system is computed
from the system geometry and the phases of the oscillating pumps. Theoretical
results are in very good agreement with experiments. We find numerically the
phase shifts between the $N$ pumps leading to maximal flow rate, for arbitrary
values of $N$. When the number of pumps is increased to infinity (continuous
limit), and for a given amplitude of deformation, the maximum flow rate
obtained is larger by 25% than that of a traditional peristaltic, sinusoidal
wave. Our results thereby cut the first turf for more efficient designs of
microfluidic pumps.

Figure 1: Sketch of the experimental setup: two cylindrical wooden rods of
diameter $1\;\rm cm$ placed on computer-controlled syringe pumps are used to
press periodically on an elastic silicone tubing. The tubing is partially
filled with V1000 silicon oil (gray). The parameterization of the system is as
follows: two pumps of length $\ell$ delimit sections of a tubing of total
length $L_{\rm tot}=2\ell+L_{1}+L_{2}+L_{3}$. The tubing radius is
$R=0.75\;\rm mm$. The varying tubing radius under the pumps is denoted
$r_{1,2}(t)$.  Figure 2: (a) Oscillating position of the oil/air meniscus when
one pump presses periodically on the silicone tubing, with a period $T=5\;\rm
s$. (b) When two pumps press periodically on the silicone tubing, the position
of the meniscus shows a net drift superposed on oscillations at the forcing
period. The drift direction depends on the sign of the phase shift
$\Delta\varphi=\varphi_{1}-\varphi_{2}$ between the two pumps. (c) Evolution
of the mean drift speed $v_{\rm exp}$ as a function of the phase shift
$\Delta\varphi$. Points: experimental data. Line: fit to a sine function.

##  II Pumping with two pumps

###  II.1 Experimental setup

The experimental setup consists in a silicon tubing of internal diameter
$2R=1.5$ mm (Thermo Scientific TM Sterilin, Thermo Fischer, France), partially
filled with silicon oil (Rhodorsil ® 47V1000, kinematic viscosity
$\nu=10^{-3}\;\rm{m^{2}\cdot s^{-1}}$). Cylindrical wooden rods of diameter 1
cm are placed on computer-controlled syringe pumps (Nemesys, Cetoni GMBH,
Germany), and press periodically on the silicon tubing, see Fig. 1. These
wooden rods constitute the pumps. The evolution of the tubing diameter at the
location of rod $i$ can be approximated as
$r_{i}(t)=R[1-\frac{\varepsilon}{2}+\frac{\varepsilon}{2}\sin(\omega
t+\varphi_{i})]$, where $\varepsilon$ is the normalized pumping amplitude,
$\omega$ the pumping angular frequency and $\varphi_{i}$ the phase shift of
pump $i$, see Fig. 1. The column of oil of length $L_{\rm tot}$ is imaged at
10 Hz with a CMOS camera (Lt225c, Teledyne Lumenera, Canada), and the position
of the oil/air interface is followed over time.

When only one pump is actuated, the oil/air meniscus moves periodically and
symmetrically left and right at the forcing frequency, leading to on average
zero net flow through the system, see Fig. 2a. This is typical of the
reversible behavior of flows at low-Reynolds number: in the absence of
inertia, a symmmetric forcing on the fluid does not lead to any net flow [15].
To estimate the value of $\varepsilon$, we measure the length of the column of
liquid $L_{\rm tot,0}$ in the absence forcing. We then measure its length
$L_{\rm tot,max}$ when a single pump is actuated, and maximally pressing on
the tube. The resulting effective radius of the deformed cross-section is then
$r_{\rm min}=r(1-\varepsilon)$. By volume conservation, we find that
$\varepsilon=1-\sqrt{1-(L_{\rm tot,max}-L_{\rm tot,0})/\ell}$, where
$\ell=1\;\rm cm$ is the diameter of the wooden rod, and so the size of the
pumping region.

When two pumps are actuated, the meniscus still moves periodically left and
right at the forcing frequency, but in the presence of an additional, slower
drift. The amplitude and direction of the drift depend on the phase shift
$\Delta\varphi=\varphi_{1}-\varphi_{2}$ between the two pumps, see Fig. 2b.
The drift velocities of the meniscus are extracted from kymographs of the
experiment, for all values of the phase shifts tested. The drift velocity $v$
shows a sinusoidal dependency on the phase shift $\Delta\varphi$, see Fig. 2c,
with no drift for $\Delta\varphi=0,\pm\pi$, and a maximum drift for
$\Delta\varphi=\pm\pi/2$. This can be intuitively understood by seeing that
when the two pumps are in phase opposition ($\Delta\varphi=\pm\pi$), the fluid
pushed by one pump fills up the space of the tube under the second pump,
resulting in a reversible situation and so an absence of net flow. When the
phase shift is $\Delta\varphi=\pm\pi/2$, the pump in phase retardation pushes
the fluid preferentially away from the other pump, in a discretized version of
peristaltic pumping.

###  II.2 Theory

The flow generated by the pumps can be calculated under the assumption of
small deformations $\varepsilon\ll 1$. Consider the valveless system of two
pumps pushing on a tubing of radius $R$ and total length $L_{\rm tot}$, as
shown in Fig. 1. The pumps have a length $\ell$, and the radius of the tubing
at the location of the pumps is $r_{i}(t)$, $i=1,2$. The tubing is delimited
by the pumps into three sections of length $L_{1}$, $L_{2}$ and $L_{3}$. The
total length of the tubing is then $L_{\rm tot}=L_{1}+L_{2}+L_{3}+2\ell$. The
pressures at the inlet and outlet of the $i$-th section are $p_{i}^{\rm
in}(t)$ and $p_{i}^{\rm out}(t)$, respectively (see bottom panel in Fig. 1).
At low Reynolds numbers and in the cylindrical geometry presented here, flows
obey Poiseuille’s law, so that the flow rate $Q_{i}(t)$ in channel $i$ is
linked to the pressures $p_{i}^{\rm in}(t)$ and $p_{i}^{\rm out}(t)$ through:

$\forall i\in[1,3],\qquad Q_{i}(t)=k_{i}[p_{i}^{\rm in}(t)-p_{i}^{\rm
out}(t)],$ (1)

where $k_{i}=\pi R^{4}/8\eta L_{i}$ is the resistivity of the channel. The
flow rates between adjacent channels are related by mass conservation:

$\forall i\in[1,2],\qquad Q_{i}(t)-Q_{i+1}(t)=2\ell\pi
r_{i}(t)\partial_{t}r_{i}.$ (2)

Last, the flow in channel $i$ can be expressed as a function of the
resistivity $\kappa_{i}(t)=\pi r_{i}(t)^{4}/8\eta\ell$ of pump $i\in[1,2]$:

$Q_{i}(t)={\kappa_{i}}(t)[p_{i}^{\rm out}(t)-p_{i+1}^{\rm in}(t)]+\ell\pi
r_{i}(t)\partial_{t}r_{i}.$ (3)

Equations (1)-(3) constitute a linear system of 7 scalar equations. The
pressures are imposed to be identical at the inlet and outlet of the system:
$p_{1}^{\rm in}=p_{3}^{\rm out}=p_{0}$. We solve for the $7$ unknown variables
$(Q_{1},Q_{2},Q_{3},p_{1}^{\rm out},p_{2}^{\rm in},p_{2}^{\rm out},p_{3}^{\rm
in})$, looking for a solution to the equations when pumps are forced
sinusoidically:
$r_{i}(t)=R\left[1-\frac{\varepsilon}{2}+\frac{\varepsilon}{2}\sin(\omega
t+\varphi_{i})\right]$.

The time-averaged flow $\bar{Q}_{1}=\bar{Q}_{2}=\bar{Q}_{3}:=\bar{Q}$ is found
to depend on the phase difference $\Delta\varphi=\varphi_{1}-\varphi_{2}$
between the two pumps through:

$\bar{Q}(\Delta\varphi)=\frac{\pi R^{2}\ell^{2}(L_{\rm tot}-\ell-
L_{2})}{L^{2}_{\rm
tot}}\omega\varepsilon^{2}\sin\Delta\varphi+\mathcal{O}(\varepsilon^{3}).$ (4)

The average flow velocity $v_{\rm theor,\;2\;pumps}=\bar{Q}/\pi R^{2}$ in the
tubing is then given by:

$v_{\rm theor,\;2\;pumps}=(L_{\rm tot}-\ell-L_{2})\frac{\ell^{2}}{L^{2}_{\rm
tot}}\omega\varepsilon^{2}\sin\Delta\varphi+\mathcal{O}(\varepsilon^{3}).$ (5)

We recover the experimentally observed dependency of the flow speed on the
sine of the phase difference $\sin\Delta\varphi$. The flow speed is
proportional to the forcing frequency $\omega$, as can be expected for low-
Reynolds flows, and depends on a geometric prefactor that involves the ratio
of the pump size $\ell$ to the total length of the fluid column $L_{\rm tot}$:
the larger the fraction of the fluid being pressed by the pumps, the more
efficient the pumping. Note that $\ell=L_{\rm tot}/2$ at most, for which the
flow speed saturates to its maximum value. Last, note that pumping is a
second-order effect in $\varepsilon$.

###  II.3 Comparison between experiments and theory

Experiments are repeated for different values of $\varepsilon$, $L_{\rm tot}$
and $\Delta\varphi$, and the experimental drift velocities of the oil are
plotted as a function of the theoretical velocities in Fig. 3a. We find good
agreement in the scaling between the second-order theory and experiments for
small values of $\varepsilon\lesssim 0.2$, see dashed line in Fig. 3a.
Deviation of the data from the dashed line at larger values of $\varepsilon$
indicate that higher-order terms are needed to fully describe the experimental
data. Solving the equations of fluid motion up to order
$\mathcal{O}(\varepsilon^{3})$ and comparing the experiments to the third-
order theory gives very good agreement for all parameters tested, up to a
numerical prefactor of $\approx 1.5$, see Fig. 3b. The expression for the
theoretical speed up to third order in $\varepsilon$ is given in 111see
Supplemental Material at [URL will be inserted by publisher], which includes
Ref. [33], for detailed expressions of the pumping speed, estimates of the
tubing deformation and pumping energy, and a movie of the optimal waveform.

The fact that the experimental speed is $\approx 1.5$ times larger than
predicted by theory can most likely be attributed to the complex deformation
of the tube when the pumps are actuated. Indeed, under an actuator, the cross-
section of the tubing does not stay circular but becomes approximately
elliptical, which leads to an increased hydrodynamic resistance at the
location of the actuator [17]. Such an increased resistance itself leads to a
larger pumping efficiency. Moreover, due to elasticity, the deformation of the
tubing extends beyond the region of the actuator, which also affects the
hydrodynamic resistance of the tube. A discussion of these two effects can be
found in [16]. Last, note that the resistance to flow is inversely
proportional to the square of the cross-sectional area [17], and so a small
error in the estimation of $\varepsilon$ can lead to large differences in the
hydrodynamic resistance.

Figure 3: Comparison between the experimental and theoretical flow speeds
(points), for a system of two pumps. The normalized pumping amplitude
$\varepsilon$ is color-coded, different symbols are for different lengths of
liquid column $L_{\rm tot}$. Other parameters: $2R=1.5\;\rm mm$, $l=1\;\rm
cm$, $L_{2}=3\;\rm cm$. (a) Experimental vs. theoretical speed calculated up
to $\mathcal{O}(\varepsilon^{2})$. The dashed line is the best fit to data
with $\varepsilon\leq 0.2$. (b) Experimental vs. theoretical speed calculated
up to $\mathcal{O}(\varepsilon^{3})$. The dashed line is the best fit to all
data, and shows the line $y=1.46x$. For a discussion on the existence of a
numerical prefactor in the linear fit of the data, see main text. Figure 4:
Comparison between the more than 200 experimental flow speed values measured
and the theoretical speed values (points), for a system of three pumps. The
normalized pumping amplitude $\varepsilon$ is color-coded. The liquid length
$L_{\rm tot}$ was varied between 13 and 18 cm. Other parameters: $2R=1.5\;\rm
mm$, $l=1\;\rm cm$. Theoretical speed calculated up to
$\mathcal{O}(\varepsilon^{2})$, according to Eq. (6). The dashed line is the
best fit to all data. The existence of a numerical prefactor between theory
and experiment is discussed in the main text.

##  III Pumping with $N\geq 3$ pumps

###  III.1 3 pumps

The theory can be extended to a system of $N$ pumps in series (see below). In
particular, for a system of 3 identical pumps, and using the same notations as
above, the average flow velocity is found to be:

$\displaystyle v_{\rm theor,\;3\;pumps}$ $\displaystyle=$
$\displaystyle\frac{\ell^{2}}{L^{2}_{\rm
tot}}\omega\varepsilon^{2}\big{[}a_{1}\sin(\varphi_{1}-\varphi_{2})+a_{2}\sin(\varphi_{1}-\varphi_{3})$
(6) $\displaystyle\qquad\qquad+\
a_{3}\sin(\varphi_{2}-\varphi_{3})\big{]}+\mathcal{O}(\varepsilon^{3}),$

where the $a_{i}$ are geometrical prefactors that have the dimension of a
length, and $\varphi_{i}$ is the phase shift of pump $i$ (see [16] for the
expressions of the $a_{i}$). The validity of this prediction is tested by
mounting a system with 3 pumps, where the forcing amplitude $\varepsilon$, the
length of the column of oil $L_{\rm tot}$, and the phase shifts $\varphi_{1}$
and $\varphi_{2}$ are systematically varied. Again, experimental results are
in good agreement with theoretical predictions, see Fig. 4. As for the system
with two pumps, the experimental velocities for the three-pump device are
larger than theoretical predictions, by a factor $\approx 2$ (see black dashed
line in Fig. 4), likely for the same reasons as cited above.

###  III.2 $N$ pumps

Is pumping more efficient when more pumps are added to the system? To answer
this question, consider a system of $N$ identical pumps of length $\ell$ in
series, actuated sinusoidally. For simplicity, we now consider that all
subsections of the tubing have the same length $L$, and so the same
resistivity $k=\pi R^{4}/8\eta L$. The total length of the tubing is then
$L_{\rm tot}=N\ell+(N+1)L$. Similarly to Eq. (1)-(3), there are now $3N+1$
equations to solve in order to compute the flow rate:

$&\forall i\in[1,N+1],\quad Q_{i}(t)=k[p_{i}^{\rm in}(t)-p_{i}^{\rm
out}(t)].\\\ \forall i\in[1,N],\quad Q_{i}(t)-Q_{i+1}(t)=2\ell\pi
r_{i}(t)\partial_{t}r_{i},\\\ \quad\quad\quad
Q_{i}(t)={\kappa_{i}}(t)[p_{i}^{\rm out}(t)-p_{i+1}^{\rm in}(t)]+\ell\pi
r_{i}(t)\partial_{t}r_{i}.\quad\quad$ (7)

The same pressure is imposed at the inlet and the outlet ($p_{1}^{\rm
in}=p_{N+1}^{\rm out}=p_{0}$), and the system of equations is solved for the
$3N+1$ unknown variables $(Q_{1},\cdots,Q_{N},p_{1}^{\rm out},p_{2}^{\rm
in},p_{2}^{\rm out},\cdots,p_{N}^{\rm out},p_{N+1}^{\rm in})$, when pumps are
forced sinusoidically
$r_{i}(t)=R\left[1-\frac{\varepsilon}{2}+\frac{\varepsilon}{2}\sin(\omega
t+\varphi_{i})\right]$.

The average flow rate over one period $T$ through a system of $N$ pumps is
found to be:

$\displaystyle\bar{Q}(\\{\varphi_{i}\\})$ $\displaystyle=$ $\displaystyle 4\pi
R^{2}\ell\omega
f_{N}(\varepsilon,b,\\{\varphi_{i}\\})+\mathcal{O}(\varepsilon^{3}),$ (8)

where we have introduced $b=\ell/L$, the ratio of the pump dimension to that
of the tubing sections, and

$f_{N}(\varepsilon,b,\\{\varphi_{i}\\})=\frac{b}{(Nb+N+1)^{2}}\times\\\
\left(\sum_{M=1}^{N-1}(M(b+1)+1)\sum_{i=1}^{M}\sin(\varphi_{i}-\varphi_{i+N-M})\right)\varepsilon^{2}+O(\varepsilon^{3})\
.$ (9)

The flow rate therefore depends on the phase shifts between the pumps, and on
the aspect ratio $b=\ell/L$; it increases linearly with $b$ for small $b$ and
saturates for $b\gg 1$, when the pumps essentially press on the entire length
of the tubing.

The phase shifts giving the maximal flow rate for a given aspect ratio $b$
were obtained by solving Eq. (9) numerically, using Mathematica (Wolfram
Research, Inc., IL, USA). For the particular case of $N=2$ pumps, the maximum
flow rate is obtained for $\Delta\varphi=\pm\pi/2$, independent of the aspect
ratio $b$. For a system of $N\geq 3$ pumps, the optimal phase shifts depend on
the aspect ratio $b$, albeit only weakly (for $N=3$ pumps, varying $b$ over 10
orders of magnitude changes the optimal values of the phase shifts by $\approx
5\%$, see Supp. Fig. 2 in [16]). Actually, for all tested values of $N\geq 3$,
the optimal phases shifts are essentially constant for $b\ll 1$ and $b\gg 1$,
and undergo a small variation when the aspect ratio $b$ goes from $0.1$ to
$10$.

We now fix the aspect ratio $b=1$ as well as the tube size $L_{\rm tot}$, and
consider that pumps of size $\ell=L_{\rm tot}/(2N+1)$ are actuated. The flow
rate in the system is then simply given by
$\tilde{Q}(\\{\varphi_{i}\\})=\bar{Q}(\\{\varphi_{i}\\})/(2N+1)$. How does the
optimal flow rate depend on the number of pumps $N$? We find that the maximal
flow rate in the system $\tilde{Q}_{\rm max}(N)$ increases with the number of
pumps $N$ until saturation in the continuous limit ($N\rightarrow\infty$), see
black markers in Fig. 5a. The maximal flow rate in a system of 2 pumps is
$\approx 50\%$ as large as in the continuous limit. The increase in flow rate
is very steep when $N$ initially increases: the flow rate in a system of 4
pumps is already $\approx 80\%$ that of the continuous limit.

Regardless of the number of pumps $N$ in the system, the optimal flow rate
$\tilde{Q}_{\rm max}$ is always larger than the flow rate $\tilde{Q}_{\rm
sine}$ obtained by naively discretizing a sine wave, for which the phase
shifts are equally spaced by $2\pi/N$ (see the orange markers in Fig. 5a, and
the dashed orange line in Fig. 5b). This is most obvious for $N=2$, where
discretizing a sine wave leads to the two pumps being actuated in phase
opposition, which does not induce any net flow. For $N=3$ pumps, the optimal
flow rate is $\approx 65\%$ larger than $\tilde{Q}_{\rm sine}$ (see inset in
Fig. 5a). The ratio $\tilde{Q}_{\rm max}/\tilde{Q}_{\rm sine}$ decreases when
$N$ increases, and in the continuous limit, the optimal flow rate is 25%
larger than the flow rate obtained with a sinusoidal peristaltic wave, see the
difference between black and orange points in Fig. 5a and inset in Fig. 5a.

The waveform leading to the maximal flow rate in the continuous limit is shown
in black in Fig. 5b for large $N$, at time $t=0$. This optimal wave travels
along the tube over one period of oscillation. The optimal waveform is
different from the traditional sinusoidal peristaltic wave; in particular the
phase values of the first and last pumps are different, see the animated phase
profiles in the Supplementary Movie [16].

##  IV Discussion

The literature on miniaturized pumps, which operate at low Reynolds number, is
extremely vast, and covers aspects from fluid dispensing in lab-on-a-chip
devices for biological analysis [18] to biomedical applications such as
microdosing for drug delivery [7]. Yet, the broad majority of the existing
work is concerned with the flow generated by at least 3 pumps in series, with
very few exceptions on pumping with two pumps [19, 20, 21]. A physical model
of the flow due to pumping with 3 pumps was described in Goulpeau et al. [22].
In our work, we studied in detail the influence of the geometric parameters of
the fluidic circuit, and of the phase shifts between the pumps, on the flow
rate, for $N=2$ and 3 pumps, and extended the theory of pumping to an
arbitrary number $N$ of pumps.

Figure 5: (a) Increasing the number of pumps in the system while keeping the
tube length constant increases the optimal flow rate $\tilde{Q}_{\rm max}(N)$
(black dots), up to a limit corresponding to the continuous case
$N\rightarrow\infty$. The optimal flow rate is always larger than the flow
rate obtained by naively discretizing a sine wave, where the phase shifts are
equally spaced by $2\pi/N$ (orange dots). Inset: ratio between the optimal
flow rate $\tilde{Q}_{\rm max}(N)$, and the flow rate $\tilde{Q}_{\rm
sine}(N)$ obtained by discretizing a sine wave, as a function of $N$. (b)
Waveform leading to the maximum flow in a system of $N$ pumps with $N=30$.
Increasing $N$ further does not change the waveform.

To achieve the highest flow rate possible, existing work usually sets the
phase shifts to arbitrary, intuitive values, and then measures the flow rate
as a function of the frequency of actuation of the pumps, see e.g. [23, 24,
25, 26, 21]. At low Reynolds number however, both quantities are proportional;
deviation from proportionality indicates either a change of regime to high
Reynolds number, or reveals the presence of a compliance in the microsystem,
which induces a phase shift between the applied actuation signal and the real
mechanical actuation. This difference between the applied actuation and the
real actuation is probably at the origin of surprising results, such as
obtaining a maximum flow rate in a system of two pumps when they are
supposedly actuated in phase ($\Delta\varphi=0$) or in phase opposition
($\Delta\varphi=\pi$) [21], in direct contradiction with the reversibility of
flows at low Reynolds number. Our results show that, in the low Reynolds
regime, valveless pumping with two pumps is optimal for phase shifts between
the pumps $\Delta\varphi=\pm\pi/2$, and we provide the exact analytical
expression of pumping efficiency as a function of the geometric parameters of
the problem. Theoretical and experimental results are in good agreement, see
Figs. 3 and 4.

The flow due to a simple peristaltic wave travelling along a tube can be
analytically calculated, see [27, 28], with a number of applications. Original
ideas making use of peristaltic waves were put forth, e.g. for the swimming of
microorganisms [29, 30]. The optimization of peristaltic wave shapes at low
Reynolds number has been studied numerically using variational methods [31]
and, more recently, boundary integral methods [32]. Here, we obtained
numerical results on the optimal peristaltic wave shape at low Reynolds number
by studying a system of $N$ discrete pumps, and let $N\rightarrow\infty$, see
Fig. 5. The non-trivial optimal peristaltic shape we find leads to a flow rate
larger by 25% than that due to a sinusoidal peristaltic wave. Moreover, in the
discrete case of a small number of pumps, the difference between the optimal
flow rate and the flow rate obtained by naively discretizing a sine wave is
even larger: for $N=3$ pumps, the optimal flow rate is $\approx 65\%$ larger
than that obtained with a discretized sine wave.

Note that we have found the wave shape leading to maximal flow rate,
considering $N$ pumps actuated at a fixed amplitude $\varepsilon$, and with a
pump length scaling as $1/N$. Numerical results indicate that, under these
constraints, the energy needed to actuate $N$ pumps with the optimal waveform
is larger than for a sinusoidal wave with the same amplitude of pumping
$\varepsilon$, see Supp. Fig. 4a in [16]. For a given energy, and relaxing the
constraint on the amplitude of pumping, a sinusoidal wave leads to a larger
flow than our optimal wave form, see Supp. Fig. 4b in [16]. We plan to
investigate in the future what the optimal phases and amplitudes of actuation
are in a system of $N$ pumps, under the constraint of a fixed energy.

Finally, the results presented here are valid for a system of $N$ pumps
actuated independently of each other, in a regime where the imposed actuation
is instantaneously transmitted to the fluid. They do not take into account
effects due to the elasticity of the tubing enclosing the fluid. Elasticity
introduces a compliance in the system as mentioned above, and can also lead to
two neighboring pumps not being independent of each other: if the pumps are
too close to each other, closing one pump will open the neighboring one. The
analysis of such intricate effects is left for future work.

###### Acknowledgements.

We thank Roméo Antier for his contribution to the early stages of this work,
and Caroline Frot for her valuable help with fabrication.

## References

  * Borrell and Krenn [2006] B. J. Borrell and H. W. Krenn, Nectar Feeding in Long-Proboscid Insects, in  _Ecology and Biomechanics: A Mechanical Approach to the Ecology of Animals and Plants_ , edited by A. Herrel, T. Speck, and N. P. Rowe (CRC Press, 2006) Chap. 9, pp. 185 – 212\.
  * Smits [1990] J. G. Smits, Piezoelectric micropump with three valves working peristaltically, Sensors and Actuators: A. Physical 21, 203 (1990).
  * Unger  _et al._ [2000] M. a. Unger, H. P. Chou, T. Thorsen, a. Scherer, and S. R. Quake, Monolithic microfabricated valves and pumps by multilayer soft lithography., Science (New York, N.Y.) 288, 113 (2000).
  * Laser and Santiago [2004] D. J. Laser and J. G. Santiago, A review of micropumps, Journal of Micromechanics and Microengineering 14, R35 (2004).
  * Woias [2005] P. Woias, Micropumps—past, progress and future prospects, Sensors and Actuators B: Chemical 105, 28 (2005).
  * Iverson and Garimella [2008] B. D. Iverson and S. V. Garimella, Recent advances in microscale pumping technologies: A review and evaluation, Microfluidics and Nanofluidics 5, 145 (2008).
  * Bußmann  _et al._ [2021] A. B. Bußmann, L. M. Grünerbel, C. P. Durasiewicz, T. A. Thalhofer, A. Wille, and M. Richter, Microdosing for drug delivery application—A review, Sensors and Actuators, A: Physical 330, (2021).
  * Uchida [1979] K. Uchida, Cibarial sensilla and pharyngeal valves in Aedes albopictus (skuse) and Culex pipiens pallens coquillett (Diptera : Culicidae), International Journal of Insect Morphology and Embryology 8, 159 (1979).
  * Kim  _et al._ [2012] B. H. Kim, E. S. Seo, J. H. Lim, and S. J. Lee, Synchrotron X-ray microscopic computed tomography of the pump system of a female mosquito, Microscopy Research and Technique 75, 1051 (2012).
  * Kikuchi  _et al._ [2018] K. Kikuchi, M. A. Stremler, S. Chatterjee, W. K. Lee, O. Mochizuki, and J. J. Socha, Burst mode pumping: A new mechanism of drinking in mosquitoes, Scientific Reports 8, 1 (2018).
  * Stemme and Stemme [1993] E. Stemme and G. Stemme, A valveless diffuser/nozzle-based fluid pump, Sensors and Actuators: A. Physical 39, 159 (1993).
  * Glickman and Yelon [2002] N. S. Glickman and D. Yelon, Cardiac development in zebrafish: coordination of form and function, Seminars in Cell & Developmental Biology 13, 507 (2002).
  * Olson [2006] E. N. Olson, Gene regulatory networks in the evolution and development of the heart, Science 313, 1922 (2006).
  * Scherz  _et al._ [2008] P. J. Scherz, J. Huisken, P. Sahai-Hernandez, and D. Y. Stainier, High-speed imaging of developing heart valves reveals interplay of morphogenesis and function, Development 135, 1179 (2008).
  * Purcell [1977] E. M. Purcell, Life at low reynolds number, American journal of physics 45, 3 (1977).
  * Note [1] See Supplemental Material at [URL will be inserted by publisher], which includes Ref. [33], for detailed expressions of the pumping speed, estimates of the tubing deformation and pumping energy, and a movie of the optimal waveform.
  * Mortensen  _et al._ [2005] N. A. Mortensen, F. Okkels, and H. Bruus, Reexamination of hagen-poiseuille flow: Shape dependence of the hydraulic resistance in microchannels, Physical Review E 71, 057301 (2005).
  * Wang and Fu [2018] Y. N. Wang and L. M. Fu, Micropumps and biomedical applications – A review, Microelectronic Engineering 195, 121 (2018).
  * Berg  _et al._ [2003] J. M. Berg, R. Anderson, M. Anaya, B. Lahlouh, M. Holtz, and T. Dallas, A two-stage discrete peristaltic micropump, Sensors and Actuators, A: Physical 104, 6 (2003).
  * Geipel  _et al._ [2007] A. Geipel, A. Doll, P. Jantscheff, N. Esser, U. Massing, P. Woias, and F. Goldschmidtboeing, A novel two-stage backpressure-independent micropump: Modeling and characterization, Journal of Micromechanics and Microengineering 17, 949 (2007).
  * Kim  _et al._ [2014] B. H. Kim, I. C. Kim, Y. J. Kang, J. Ryu, and S. J. Lee, Effect of phase shift on optimal operation of serial-connected valveless micropumps, Sensors and Actuators A: Physical 209, 133 (2014).
  * Goulpeau  _et al._ [2005] J. Goulpeau, D. Trouchet, A. Ajdari, and P. Tabeling, Experimental study and modeling of polydimethylsiloxane peristaltic micropumps, Journal of Applied Physics 98, 10.1063/1.1947893 (2005).
  * Xie  _et al._ [2004] J. Xie, J. Shih, Q. Lin, B. Yang, and Y. C. Tai, Surface micromachined electrostatically actuated micro peristaltic pump, Lab on a Chip 4, 495 (2004).
  * Jang and Kan [2007] L. S. Jang and W. H. Kan, Peristaltic piezoelectric micropump system for biomedical applications, Biomedical Microdevices 9, 619 (2007).
  * Liu  _et al._ [2010] G. Liu, C. Shen, Z. Yang, X. Cai, and H. Zhang, A disposable piezoelectric micropump with high performance for closed-loop insulin therapy system, Sensors and Actuators, A: Physical 163, 291 (2010).
  * Fuchs  _et al._ [2012] O. Fuchs, Y. Fouillet, S. Maubert, M. Cochet, C. Chabrol, N. David, X. Médal, and R. Campagnolo, A novel volumetric silicon micropump with integrated sensors, Microelectronic Engineering 97, 375 (2012).
  * Shapiro  _et al._ [1969] A. Shapiro, M. Jaffrin, and S. Weinberg, Peristaltic pumping with long wavelengths at low Reynolds number, J. Fluid Mech 37, 799 (1969).
  * Jaffrin and Shapiro [1971] M. Y. Jaffrin and a. H. Shapiro, Peristaltic Pumping, Annual Review of Fluid Mechanics 3, 13 (1971).
  * Ajdari  _et al._ [1994] A. Ajdari, P.-G. De Gennes, J.-P. Hulin, L. Leibler, and J. Prost, Propulsion péristaltique, Comptes rendus de l’Académie des sciences. Série II, Mécanique, physique, chimie, astronomie 319, 861 (1994).
  * Ajdari and Stone [1999] A. Ajdari and H. A. Stone, A note on swimming using internally generated traveling waves, Physics of Fluids 11, 1275 (1999).
  * Walker and Shelley [2010] S. W. Walker and M. J. Shelley, Shape optimization of peristaltic pumping, Journal of Computational Physics 229, 1260 (2010).
  * Bonnet  _et al._ [2020] M. Bonnet, R. Liu, and S. Veerapaneni, Shape optimization of Stokesian peristaltic pumps using boundary integral methods, Advances in Computational Mathematics 46, 1 (2020).
  * Ramanujan [1914] S. Ramanujan, Modular equations and approximations to $\pi$, Quart. J. Math 45, 350 (1914).