∥

<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>aff1]Department of Mechanical Engineering, Bangladesh
University of Engineering and Technology aff2]Department of Mechanical
Engineering, Auburn University

aff3]Department of Electrical & Electronic Engineering, Bangladesh University
of Engineering and Technology

# Comparison of Different Control Theories on a Two Wheeled Self Balancing
Robot

MD Muhaimin Rahman  SM Hasanur Rashid  KM Rafidh Hassan  M. M. Hossain  [  [
[

###### Abstract

This paper is aimed to discuss and compare three of the most famous Control
Theories on a Two wheeled Self Balancing Robot Simulation using Robot
Operating System (ROS) and Gazebo. Two Wheeled Self Balancing Robots are one
of the most fascinating applications of Inverted Pendulum System. In this
paper, PID, LQR and Fuzzy logic controllers are discussed . Also,0 the
modeling and algorithms of the robot simulation is discussed. The primary
objectives of this paper is to discuss about the building of a robot model in
ROS and Gazebo , experimenting different control theories on them, documenting
the whole process with the analysis of the robot and comparison of different
control theories on the system.

##  1 INTRODUCTION

Inverted Pendulum is one of the most discussed and fascinating topics of Non-
linear control System in the academia. Many researches have been done on this
problem. According to Boubaker et al, no less than ten types of Control
theories have been researched on , explored using this problem [1]. One of the
most interesting applications of this theory is the Two Wheeled Self balancing
Robot. Many researches have been done on this robot also. Like, [3] discussed
about developing simple Self balancing robot from popular micro controller
board, Arduino. In [4], a ”Hopfield Neural Network” based controller. In [5] ,
same robot was developed using Sensor Fusion Algorithm. Again , in [2] , LQR
based controller was used for this robot. In [6] , self balancing robot was
developed using Boltzmann machine. This project is different than others in
main purpose. The purpose of this paper is to describe the development of
ready to use control algorithms using ROS and Gazebo simulation. ROS and
Gazebo are state-of-the-art Open Source Robotics Simulation platform. But
there is hardly any good resource on the ROS and control theory , let alone
complex non-linear controllers. So , this research opened a door for ROS and
Control Systems future Robotics researchers.
The paper is structured as following , first section discusses ROS and Gazebo,
second section discusses Experimenting Control system packages with ROS and
Gazebo. They are PID, Fuzzy Logic and LQR controllers. The fourth section
sheds light on comparison to the available works and literature on same topic
and the final section is the conclusion.

Table 1: Nomenclature M | Mass of the cart, $0.0754kg$ | I | Mass moment of inertia 0.01094 $kgm^{2}$
---|---|---|---
m | Mass of the pendulum | b | Dynamic Friction Coefficient = 0.65
l | Length of the pendulum , 0.157 $m$ | g | Gavitational Acceleration= 9.8 $ms^{-}2$
$P(s)$ | Transfer Function of the system | $\phi(s)$ | Transfer function of the yaw angle
$U(s)$ | Transfer Function of the Force |  |

##  2 ROBOT OPERATING SYSTEM AND GAZEBO

Robot Operating System (ROS) [8] is arguably the most powerful tool for modern
Robotics Software development. It is a collection of Open Source frameworks,
tools ,libraries for research and development of modern robots . Like
Universal Robots, PR2 etc are using ROS at present. In this project we are
using ROS and Gazebo simulation. Because of following reasons,

  * •

It is Open Source, hence it has a vast working frameworks available.

  * •

We want Robotics Community use and explore various control theories.

  * •

But, there is hardly any good resources on ROS Control Theory publicly
available.

For these reasons, we have developed the self balancing robot Gazebo [9]
model. Gazebo is also arguably the best Robotics Simulator at present. Its
widely used with ROS interface. It is also open source, which makes it widely
accessible. Figure 1(b) shows the Gazebo Model of The Robot.

###  2.1 Building the Robot Model

To build the robot Model, the required steps are,

  * •

Defining Links

  * •

Defining Joints

  * •

Adding ROS Plugins

  * •

Writing the ROS code to test the model

The detailed description is given below

####  2.1.1 Link

Link is the building block of every Gazebo simulation . Every robot model will
have several types of links. Links can be defined in two ways.

  * •

Writing xml code from beginning

  * •

Adding a CAD model in mesh file format.

  * •

Combination of Both

For building the model by writing xml code , we need to define every aspect of
links. Like Mass, geometry, moment of Inertia around every axis etc. It
becomes almost impossible for complex links. In that case, the best way is to
integrate CAD files of the links. Another way is to integrate the CAD files as
well as defining the simpler links in the xml code.[17]

In our model of robot, we have three links, Two wheels and a chassis. A sample
xml code for the Right wheel of the robot is given below,

    Ψ<link name = "Right_Wheel">
    ΨΨ<inertial>
    ΨΨΨ<origin xyz ="0.00091 0.01535 0.07450" rpy = "0 0 0"/>
    ΨΨΨ<mass value ="0.03770"/>
    ΨΨΨ<inertia ixx ="0.000233" ixz = "0.000002560" iyy = "0.000236407" iyz = "0.000043120"
     ixy = "0.000000528" izz = "0.000024435"/>
    ΨΨ</inertial>
    ΨΨ<collision>
    ΨΨΨ<geometry>
    ΨΨΨ<cylinder radius="0.0316" length ="0.037"/>Ψ
    ΨΨΨ
    ΨΨΨ</geometry>
    ΨΨΨ<!--
    ΨΨΨ<origin rpy = "0 0 0" xyz = "0 0.1100 0.0316"/> -->
    ΨΨΨ<origin rpy = "0 0 0" xyz = ".005 0.160 0.010"/>

    ΨΨ</collision>
    ΨΨ<visual>
    ΨΨΨ<geometry>
    ΨΨΨ<mesh filename="package://self_balancing_robot/meshes/WheelZTiny.dae"/>
    ΨΨΨ</geometry>
    ΨΨΨ
    ΨΨΨ<origin rpy = "0 0 0" xyz = ".005 0.160 0.010"/>
    ΨΨΨ
    ΨΨ</visual>
    Ψ</link>

####  2.1.2 Joint

Joint tags define how the links are attached to the base link and other links.
[17] According to [18], There are 6 types of joints

  * •

_Revolute_ \- links can rotate around specified axis within limits.

  * •

_Continuous_ \- It is like revolute type of joint but it has no limits

  * •

_Prismatic_ \- Prismatic joint defines the joints where one link can slide on
other

  * •

_Fixed_ \- Fixed joints are defined for links which cannot move with respect
to other

  * •

_Floating_ \- Floating type of joints allow movement around six axis of
rotation.

  * •

_Planar_ \- Planar type of joints allow motion perpendicular to a plane .

In our robot, the wheels are attached to the main body using continuous joint.
The figure 1 Shows the Block Diagram and the Gazebo model.

---
(a) Block Diagram for Gazebo Model

---
(b) Gazebo Simulation

Figure 1: Gazebo Simulation Block Diagram and Model
Figure 2: ROS Controller Application Block Diagram

####  2.1.3 ROS Plugins

Plugins make the robot models function according to ROS commands.Without
plugins we cannot operate the actuators of the model and we cannot access the
sensor outputs. There are many built in plugins for ROS and Gazebo.[19] We
have used Differential Drive Plugin to provide differential drive movement for
the robot and Imu Plugin for Imu sensor of the robot.

##  3 EXPERIMENT WITH ROS AND OUR CONTROL SYSTEM PACKAGE

In this section the implementation and performance of three controllers ,
namely PID, Fuzzy Logic and LQR on the Robot simulation have been discussed
using ROS and third party Control System Packages like [23], [25] and [14].
The general algorithm for implementing Controllers is given below

  * •

Measure the Pitch Angle of the Robot and save it to value y

  * •

Set Previous Error to Zero

  * •

Set Sum of Errors to Zero

  * •

Current Error = Pitch value-Setpoint

  * •

Error Difference = Current Error - Previous Error

  * •

Sum of Errors = Current Error+Sum of Errors

  * •

Get Output from controllers for different Parameters

###  3.1 PID

Proportional Integral and Derivative (PID) is the most widely used Controller
in the world. Li Yun et al discussed about this controller extensively in
their paper [10]. In our previous paper [11] , we have discussed about
applying PID controllers in a real Self Balancing Robot using Simulink. In
this paper, we are showing the performance of PID packages [23] . It follows
the equation 1 . Figure 3 show the respose curves for various PID gains. The
curves are plotted the $rqt\\_plot$ , which also can be used for Real Robot’s
response plotting.

$U=K_{p}*e+K_{i}\int edt+K_{d}\frac{de}{dt}$ (1)

---
(a) Kp=25, Ki =0.8, Kd =0.1

---
(b) Kp=50 Ki=0.8 Kd =0.05

---
(c) Kp=100 Ki=0.8 Kd =0.1

---
(d) Kp=1000, Ki =0.8, Kd =0.05

Figure 3: PID Response Curves for Four different values

###  3.2 Fuzzy Logic

Fuzzy logic controllers are one of the most robust control systems in the
world. In the academia , it was first described by Chuen Chien Lee [12]. There
are many advantages of Fuzzy logic controllers. First of all, they can work on
Multiple Input Multiple Output (MIMO) systems. Secondly, one doesn’t need to
derive the plant models for tuning the controllers. In this paper we have
developed fuzzy PD and fuzzy PD+I controllers based on [14]. The rule base is
given in Table 2. The membership functions for the velocity and yaw angle of
the robot are given in Figure 4

---
(a) Pitch Angle of the robot

---
(b) Pitch velocity of the robot

Figure 4: Fuzzy Membership Function Table 2: Fuzzy Rule Base $e^{\prime}/e$ | Long Negative | Short Negative | Medium | Short Positive | Long Positive
---|---|---|---|---|---
Long Negative | High Negative | Negative | Zero | Positive | Positive
Short Negative | High Negative | Negative | Zero | High Positive | High Positive
Medium | High Negative | High Negative | Zero | High Positive | High Poitive
Short Positive | High Negative | Negative | Zero | Positive | High Positive
Long Positive | Negative | Negative | Zero | Positive | High Positive

---
(a) Fuzzy PD

---
(b) Fuzzy PD+I

Figure 5: Performance Curves for Fuzzy Logic Controller

---
(a) $Q_{1}$ and $R_{1}$

---
(b) $Q_{2}$ and $R_{2}$

Figure 6: Performance Curves for LQR Controller

###  3.3 LQR

Another Robust Controller is Linear Quadratic Regulator Controller. H. Wang et
al [15] had discussed about the Implementation of LQR controllers on Inverted
pendulum. Also [7] has an extensive section about how to implement using
Matlab. But , to the best of our knowledge , it has not been implemented in
ROS simulation yet. The state space equations are given in [7] as following,

$X=A+Bu$ (2)

Here,

$u=-KX$ (3)

, so the equation 2 becomes

$X=(A-BK)X$ (4)

where ,

$A=\left[\matrix{0&1&0&0\cr
0&\frac{-(I+ml^{2})b}{I(M+m)+Mm^{2}}&\frac{m^{2}gl^{2}}{I(M+m)+Mm^{2}}&0\cr
0&0&0&1\cr
0&\frac{-mlb}{I(M+m)+Mm^{2}}&\frac{mgl(M+m)}{I(M+m)+Mm^{2}}&0}\right]$ (5)
$B=\left[\matrix{0\cr\frac{I+ml^{2}}{I(M+m)+Mm^{2}}\cr
0\cr\frac{ml}{I(M+m)+Mm^{2}}}\right]$ (6)
$[X=\left[\matrix{\dot{x}\cr\ddot{x}\cr\dot{\phi}\cr\ddot{\phi}}\right]$ (7)

As we are more concerned with the yaw angles of the robot, not the linear
movement and/or velocity , in our case the values of A,B and X will change. In
our case , the values will become,

$A=\left[\matrix{0&1\cr\frac{mgl(M+m)}{I(M+m)+Mm^{2}}&0}\right]$ (8)
$B=\left[\matrix{0\cr\frac{ml}{I(M+m)+Mm^{2}}}\right]$ (9)
$X=\left[\matrix{\dot{\phi}\cr\ddot{\phi}}\right]$ (10)

From [16] we get the cost function

$J=\int(x^{T}Qx+u^{T}Ru)dt$ (11)

Where , Q is the penalty matrix for the system, which penalizes the cost
function for any error in the input. R penalizes the force input in the
system. Now, the goal in the LQR controller is to choose $K$ matrix such that,
we will have minimum cost function according to our need , i.e. according to
our matrices $Q$ and $R$. In this project, we have experimented on two values
of Q and R. For first one ,

$Q_{1}=\left[\matrix{10&0\cr 0&100}\right]$ (12)

and

$R_{1}=0.001$ (13)

The Second values are

$Q_{2}=\left[\matrix{100&0\cr 0&1000}\right]$ (14)

and

$R_{2}=0.0001$ (15)

Figure 6 shows the Response curves .

##  4 COMPARISON TO PREVIOUS WORKS

###  4.1 Available ROS Controller

In the available literature on how to use control systems with ROS , only PID
control systems are available.They focus on how to use for different purpose,
like Joint Controller, Position Controller, Effort Controller etc.[20] But to
the best of our knowledge no work exists on non PID controllers with ROS.
Moreover, the existing work on PID Control system is not user friendly.
Impossible to tune the controller gains while simulation is running. The
application of available ROS controller system is shown in Fig. 2 which is
retrieved from [20]

###  4.2 MATLAB and Simulink Simulation

Matlab and Simulink are two of the most famous simulation softwares in Control
systems and Robotics. But they have some problems

  * •

To work on Matlab/Simulink simulation, one has to derive the control system of
the robot before hand using physical laws and some assumptions. There is a
high chance of error in this mathematical modeling , as it depends largely on
the assumptions of those physical laws. [21] In gazebo, there’s no need to go
through the mathematical model of the robot. The properties needed from
systems like, mass, moment of inertia, center of gravity are used to build the
robot model. The Physics Engine of Gazebo itself adjusts and derives the
mathematical model from the robot model’s properties.[17]

  * •

However, in Simmechanics one can import CAD models for simulation.[22] But it
has it’s own limits.

    * –

It is suitable for simple mechanical models.

    * –

It can only simulate mechanical properties , not the simulation for other
sensors , like, camera, imu, laser sensors etc.

  * •

MATLAB and SIMULINK can only be used for control system simulation. But ROS
and GAZEBO can be used like real world problems where simulation of control
system, computer vision, Three dimensional positionin, path planning of the
robots are needed.

  * •

The codes used in ROS simulation can directly be used in Real world robots,
which is not the case for MATLAB and SIMULINK

##  5 CONCLUSION

The main purpose of this paper is develop explore and experiment different
control theory on ROS simulation of Self Balancing Robot, which has been
achieved successfully. The Gazebo model were built successfully. The real time
data of the pitch angles vs time were plotted. From the Response curves, it
can be noticed that , the most stable performance was achieved using PID
controller with the gains 50,0.8 and 0.05 . The faster response is the
response curve of LQR controllers with First Q matrix and R value. The Fuzzy
logic controller based performance is highly unstable. Which leads to our next
research topic. Tuning the fuzzy rule base. In this project, we have written
the rule base , based on intuition. Several algorithms exist in the literature
to tune the correct rule base. Nevertheless, we hope that this project will be
beneficial to the Robotics Researchers to learn more about the Control theory.

## References

  * [1] Olfa Boubaker, _The Inverted Pendulum Benchmark in Nonlinear Control Theory: A Survey_ ,International Journal of Advanced Robotic Systems, 2013
  * [2] Jian Fang,_The LQR Controller Design of Two-Wheeled Self-Balancing Robot Based on the Particle Swarm Optimization Algorithm_ , Mathematical Problems in Engineering, Volume 2014, Article ID 729095
  * [3] Hau-Shiue Juang , Kai-Yew Lurrr Design and control of a two-wheel self-balancing robot using the arduino microcontroller board,_10th IEEE International Conference on Control and Automation (ICCA)_ , 2013 .
  * [4] Xian-gang RuanJianwei Zhao, ”The Flexible Two-Wheeled Self-balancing Robot Based on Hopfield”,_International Conference on Intelligent Robotics and Applications ICIRA 2009_ : Intelligent Robotics and Applications pp 1196-1204
  * [5] Jianhai Han, Xiangpan Li, and Qi Qin , _Design of Two-Wheeled Self-Balancing Robot Based on Sensor Fusion Algorithm_ ,International Journal of Automation Technology, Vol8 No.2 ,pp. 216-221, 2014
  * [6] Jianwei Zhao, Xiaogang Ruan, The flexible two-wheeled self-balancing robot intelligence controlling based on Boltzmann, _ROBIO’09 Proceedings of the 2009 international conference on Robotics and biomimetics_
  * [7] _Control Tutorials using Simulink and Matlab_. Link http://ctms.engin.umich.edu/CTMS/index.php?example=
InvertedPendulum&section=SystemAnalysis , retrieved 25 August, 2017

  * [8] Morgan Quigley, Brian Gerkey, Ken Conley, Josh Faust, Tully Foote, Jeremy Leibs, Eric Berger, Rob Wheeler, Andrew Ng, ROS: an open-source Robot Operating System, _ICRA workshop on open source software. Vol. 3. No. 3.2. 2009._
  * [9] Nathan Koenig, Andrew Howard, Design and Use Paradigms for Gazebo, An Open-Source Multi-Robot Simulator, _IEEE/RSJ International Conference on Intelligent Robots and systems 2004_
  * [10] Li, Yun, Kiam Heong Ang, and Gregory CY Chong, _PID control system analysis and design._ IEEE Control Systems 26.1 (2006): 32-41.
  * [11] Muhaimin Rahman, Md & , Ashik-E-Rasul & Haq, Nowab Md Aminul & Hassan, Mehedi & Mohammad Al Hasib, Irfan & M. Rafidh Hassan, K. (2016). Development of a two wheeled self balancing robot with speech recognition and navigation algorithm. _Proceedings of International Conference on Mechanical Engineering, 2015, Dhaka, BUET_ , American Institute of Physics, Vol. 1754. 060005. 10.1063/1.4958446.
  * [12] Chuen Chien Lee, _Fuzzy Logic in Control Systems: Fuzzy Logic Controller, Part II_ , IEEE Transactions on Systems, Man and Cybernetics, volume 20, No 2, March/ April 1990
  * [13] Jan Jantzen, _Tuning Of Fuzzy PID Controllers_
  * [14] Joshua D. Warner; Jason Sexauer; twmeggs; Alexandre M. S.; Aishwarya Unnikrishnan; Guilherme Castelão; Fernando Batista; Himanshu Mishra,_JDWarner/scikit-fuzzy: Scikit-Fuzzy 0.2_ , website https://github.com/scikit-fuzzy/scikit-fuzzy , retrieved at 26th August, 2017
  * [15] Wang, H., Dong, H., He, L., Shi, Y., and Zhang, Y. (2010). Design and simulation of lqr controller with the linear inverted pendulum. _In International Conference on Electrical and Control Engineering, China._
  * [16] Richard M. Murray , _Optimization-Based Control_ , Control and Dynamic Systems, California Institute of Technology, Page 12
  * [17] Nathan Koenig, Andrew Howard, Design and Use Paradigms for Gazebo, An Open-Source Multi-Robot Simulator, _International Conference on Intelligent Robots and Systems, 2004_
  * [18] _ROS URDF Joint Tutorial_ , http://wiki.ros.org/urdf/XML/joint , Retrieved at 15th November 2017
  * [19] _ROS Gazebo Plugins_ , http://gazebosim.org/tutorials? tut=ros_gzplugins , retrieved on 15th November
  * [20] _ROS control official documentation_ , https://github.com/ros-controls/ros_control/wiki , retrieved at 15th November, 2017
  * [21] V.K. Singh _MATLAB Simulation for Control System_ , International Journal of Innovative Research in Electrical , Electronics, Instrumentation and Control Engineering, Vol.2 Issue1, January 2014
  * [22] Hroncova, D. , & Pastor, M. (2013). _Mechanical System and SimMechanics Simulation_ , American Journal of Mechanical Engineering, 1(7), 251-255.
  * [23] Bipeen Acharya, Fred Gisa, and Simon D. Levy, _Python Library for PID control_ , https://github.com/simondlevy/PyQuadSim/blob/master/pidcontrol.py , retrieved at 15th November, 2017
  * [24] ROS _RQT Plot_ http://wiki.ros.org/rqt_plot, retrieved at 15th November,2017
  * [25] _Python Control Package_ , https://github.com/python-control/python-control, retrieved at 15th November, 2017


# Lecture Notes on Stationary Gamma Processes

Robert L. Wolpert

###### Abstract

For each $\lambda>0$ and every square-integrable infinitely-divisible (ID)
distribution there exists at least one stationary stochastic process $t\mapsto
X_{t}$ with the specified distribution for $X_{1}$ and with first-order
autoregressive (AR(1)) structure in the sense that the autocorrelation of
$X_{s}$ and $X_{t}$ is $\exp(-\lambda|s-t|)$ for all indices $s,t$. For the
special case of the standard Normal distribution, the process $X_{t}$ is
unique— namely, the first-order autoregressive Ornstein-Uhlenbeck velocity
process. The process $X_{t}$ is also uniquely determined if $X_{1}$ is
accorded the unit rate Poisson distribution.

For the Gamma distribution, however, $X_{t}$ is _not_ determined uniquely. In
these lecture notes we describe six distinct processes with the same
univariate marginal distributions and AR(1) autocorrelation function. We
explore a few of their properties and describe methods of simulating their
sample paths.

##  1 Introduction

For fixed $\alpha,\beta>0$ these notes present six different stationary time
series, each with Gamma $X_{t}\sim\mathsf{Ga}(\alpha,\beta)$ univariate
marginal distributions and autocorrelation function $\rho^{|s-t|}$ for
$X_{s},~{}X_{t}$. Each will be defined on some time index set ${\mathcal{T}}$,
either ${\mathcal{T}}={\mathbb{Z}}$ or ${\mathcal{T}}={\mathbb{R}}$.

Five of the six constructions can be applied to other Infinitely Divisible
(ID) distributions as well, both continuous ones (normal, $\alpha$-stable,
_etc_.) and discrete (Poisson, negative binomial, _etc_.). For specifically
the Poisson and Gaussian distributions, all but one of them (the Markov
change-point construction) coincide— essentially, there is just one
“AR(1)-like” Gaussian process (namely, the AR(1) process in discrete time, or
the Ornstein-Uhlenbeck process in continuous time), and there is just one
AR(1)-like Poisson process. For other ID distributions, however, and in
particular for the Gamma, each of these constructions yields a process with
the same univariate marginal distributions and the same autocorrelation but
with different joint distributions at three or more times.

First, by “$\mathsf{Ga}(\alpha,\beta)$” we mean the Gamma distribution with
mean $\alpha/\beta$ and variance $\alpha/\beta^{2}$, _i.e._ , with shape
parameter $\alpha$ and _rate_ parameter $\beta$. The pdf and chf are given by:

$\displaystyle f(x\mid\alpha,\beta)$
$\displaystyle=\frac{\beta^{\alpha}}{\Gamma(\alpha)}~{}x^{\alpha-1}e^{-\beta
x}\mathbf{1}_{\\{x>0\\}}$ $\displaystyle\chi(\omega\mid\alpha,\beta)$
$\displaystyle=(1-i\omega/\beta)^{-\alpha}$
$\displaystyle=\exp\left\\{-\int_{{\mathbb{R}}_{+}}e^{i\omega u}~{}\alpha
e^{-\beta u}\,u^{-1}\,du\right\\}.$ (1)

The sum $\xi_{+}:=\sum\xi_{j}$ of independent random variables
$\xi_{j}\mathrel{\mathop{\sim}\limits^{\mathrm{ind}}}\mathsf{Ga}(\alpha_{j},\beta)$
with the same rate parameter also has a Gamma distribution,
$\xi_{+}\sim\mathsf{Ga}(\alpha_{+},\beta)$, if
$\left\\{\alpha_{j}\right\\}\subset{\mathbb{R}}_{+}$ are summable with sum
$\alpha_{+}:=\Sigma\alpha_{j}<\infty$. Eqn (1) shows that the “Lévy measure”
for this distribution is $\nu(du)=\alpha e^{-\beta
u}\,u^{-1}\mathbf{1}_{\\{u>0\\}}du$.

From the Poisson representation of ID distributions

$X=\int_{\mathbb{R}}u\mathcal{N}(du)$

for $\mathcal{N}(du)\sim\mathsf{Po}\big{(}\nu(du)\big{)}$ we can show that the
extremal properties of ID random variables depend only on the Lévy measure
$\nu(du)$: for large $u\in{\mathbb{R}}_{+}$,

$\displaystyle{\mathsf{P}}[X>u]$
$\displaystyle\approx{\mathsf{P}}\big{[}\mathcal{N}\big{(}(u,\infty)\big{)}\neq\emptyset]$
$\displaystyle=1-\exp\big{(}-\nu\big{(}(u,\infty)\big{)}$
$\displaystyle\approx\nu\big{(}(u,\infty)\big{)}\approx(\alpha/\beta
u)e^{-\beta u},$ (2)

from which one might mount a study of the multivariate extreme properties of
the six processes presented below.

##  2 Six Stationary Gamma Processes

###  2.1 The Gamma $\textsf{AR}(1)$ Process

Fix $\alpha,\beta>0$ and $0\leq\rho<1$. Let
$X_{0}\sim\mathsf{Ga}(\alpha,\beta)$ and for $t\in{\mathbb{N}}$ define $X_{t}$
recursively by

$\displaystyle X_{t}$ $\displaystyle:=\rho X_{t-1}+\zeta_{t}$ (3)

for iid $\left\\{\zeta_{t}\right\\}$ with chf
${\mathsf{E}}e^{i\omega\zeta_{t}}=(1-i\omega/\beta)^{-\alpha}(1-i\rho\omega/\beta)^{\alpha}=\big{[}\frac{\beta-i\omega}{\beta-i\rho\omega}\big{]}^{-\alpha}$
(easily seen to be positive-definite, with Lévy measure
$\nu_{\zeta}(du)=\alpha\left[e^{-\beta u}-e^{-\beta
u/\rho}\right]u^{-1}\mathbf{1}_{\\{u>0\\}}du$). A simple way of generating
$\left\\{\zeta_{t}\right\\}$ with this distribution from [Walker, 2000] is
presented in the Appendix. The process $\\{X_{t}\\}$ has Gamma univariate
marginal distribution $X_{t}\sim\mathsf{Ga}(\alpha,\beta)$ for every
$t\in{\mathbb{R}}_{+}$ and, at consecutive times $0,1$, joint chf

$\displaystyle\chi(s,t)$ $\displaystyle={\mathsf{E}}\exp(isX_{0}+itX_{1})$
$\displaystyle={\mathsf{E}}\exp(i(s+\rho t)X_{0}+it\zeta_{1})$
$\displaystyle=\Big{[}\frac{(1-i(s+\rho
t)/\beta)~{}(1-it/\beta)}{1-it\rho/\beta}\Big{]}^{-\alpha}.$ (4)

Since this is asymmetric in $s,t$, the process is not time-reversible; this is
also evident from the observation that

${\mathsf{P}}[X_{t}\geq\rho X_{t-1}]={\mathsf{P}}[\zeta_{t}\geq
0]=1>{\mathsf{P}}[\zeta_{t}\leq
X_{t-1}(1-\rho^{2})/\rho]={\mathsf{P}}[X_{t-1}\geq\rho X_{t}].$

This process has marginal distribution $X_{t}\sim\mathsf{Ga}(\alpha,\beta)$ at
all times $t$, and has autocorrelation
${\mathsf{Corr}}\big{(}X_{s},X_{t})=\rho^{|s-t|}$ (easily found from either
(3) or (4)). It is clearly Markov (from (3)), and can be shown to have
infinitely-divisible (ID) multivariate marginal distributions of all orders,
since the $\\{\zeta_{t}\\}$ are ID.

###  2.2 Thinned Gamma Process

If $X\sim\mathsf{Ga}(\alpha_{1},\beta)$ and
$Y\sim\mathsf{Ga}(\alpha_{2},\beta)$ are independent, then
$Z:=X+Y\sim\mathsf{Ga}(\alpha_{+},\beta)$ and
$U:=X/Z\sim\mathsf{Be}(\alpha_{1},\alpha_{2})$ are also independent (where
$\alpha_{+}:=\alpha_{1}+\alpha_{2}$), because under the change of variables
$x=uz$, $y=(1-u)z$ with Jacobian $J(u,z)=z$,

$\displaystyle f(u,z)$
$\displaystyle=\left\\{\frac{\beta^{\alpha_{1}}}{\Gamma(\alpha_{1})}x^{\alpha_{1}-1}e^{-\beta
x}\right\\}\left\\{\frac{\beta^{\alpha_{2}}}{\Gamma(\alpha_{2})}y^{\alpha_{2}-1}e^{-\beta
y}\right\\}\,J(u,z)$
$\displaystyle=\left\\{\frac{\beta^{\alpha_{+}}}{\Gamma(\alpha_{1})\Gamma(\alpha_{2})}\right\\}u^{\alpha_{1}-1}(1-u)^{\alpha_{2}-1}z^{\alpha_{+}-2}e^{-\beta
z}z$
$\displaystyle=\left\\{\frac{\Gamma(\alpha_{+})}{\Gamma(\alpha_{1})\Gamma(\alpha_{2})}u^{\alpha_{1}-1}(1-u)^{\alpha_{2}-1}\right\\}\quad\left\\{\frac{\beta^{\alpha_{+}}}{\Gamma(\alpha_{+})}z^{\alpha_{+}-1}e^{-\beta
z}\right\\}.$

Thus if $Z\sim\mathsf{Ga}(\alpha_{1}+\alpha_{2},\beta)$ and
$U\sim\mathsf{Be}(\alpha_{1},\alpha_{2})$ are independent then
$X=UZ\sim\mathsf{Ga}(\alpha_{1},\beta)$ and
$Y=(1-U)Z\sim\mathsf{Ga}(\alpha_{2},\beta)$ are independent too. Let

$\displaystyle X_{0}$ $\displaystyle\sim\mathsf{Ga}(\alpha,\beta)$ and, for
$t\in{\mathbb{N}}$ (or $t\in-{\mathbb{N}}$, resp.) set $\displaystyle X_{t}$
$\displaystyle:=\xi_{t}+\zeta_{t}$ where $\displaystyle\xi_{t}$
$\displaystyle:=B_{t}\cdot X_{t-1},\qquad
B_{t}\sim\mathsf{Be}(\alpha\rho,\alpha{\bar{\rho}}),$ $\displaystyle\zeta_{t}$
$\displaystyle\sim\mathsf{Ga}(\alpha{\bar{\rho}},\beta)$ (or
$\xi_{t}=B_{t}\cdot X_{t+1}$, resp.), where ${\bar{\rho}}:=(1{-}\rho)$ and all
the $\left\\{B_{t}\right\\}$ and $\left\\{\zeta_{t}\right\\}$ are independent.

Then, by induction, $\xi_{t}\sim\mathsf{Ga}(\alpha\rho,\beta)$ and
$\zeta_{t}\sim\mathsf{Ga}(\alpha{\bar{\rho}},\beta)$ are independent, with sum
$X_{t}\sim\mathsf{Ga}(\alpha,\beta)$. Thus $\\{X_{t}\\}$ is a Markov process
with Gamma univariate marginal distribution
$X_{t}\sim\mathsf{Ga}(\alpha,\beta)$, now with symmetric joint chf

$\displaystyle\chi(s,t)$ $\displaystyle={\mathsf{E}}\exp(isX_{0}+itX_{1})$
$\displaystyle={\mathsf{E}}\exp\left\\{is(X_{0}-\xi_{1})+i(s+t)\xi_{1}+it\zeta_{1}\right\\}$
$\displaystyle=(1-is/\beta)^{-\alpha{\bar{\rho}}}(1-i(s+t)/\beta)^{-\alpha\rho}(1-it/\beta)^{-\alpha{\bar{\rho}}}$
(5)

and autocorrelation ${\mathsf{Corr}}\big{(}X_{s},X_{t})=\rho^{|s-t|}$. The
process of passing from $X_{t-1}$ to $\xi_{t}=X_{t-1}B_{t}$ is called
_thinning_ , so $X_{t}$ is called the _thinned gamma process_. A similar
construction is available for any ID marginal distribution.

###  2.3 Random Measure Gamma Process

Let ${\mathcal{G}}(dx\,dy)$ be a countably additive random measure that
assigns independent random variables
${\mathcal{G}}(A_{i})\sim\mathsf{Ga}(\alpha|A_{i}|,\beta)$ to disjoint Borel
sets $A_{i}\in{\mathcal{B}}({\mathbb{R}}^{2})$ of finite area $|A_{i}|$ (this
is possible by the Kolmogorov consistency conditions, and is illustrated in
the Appendix) and, for $\lambda:=-\log\rho$, consider the collection of sets:

$G_{t}:=\left\\{(x,y):~{}x\in{\mathbb{R}},~{}0\leq y<\lambda
e^{-2\lambda|t-x|}\right\\}$

Figure 1: Random measure construction of process $X_{t}={\mathcal{G}}(G_{t})$

shown in Figure (1) whose intersections have area $|G_{s}\cap
G_{t}|=e^{-\lambda|s-t|}$. For $t\in{\mathcal{T}}={\mathbb{R}}$, set

$X_{t}:={\mathcal{G}}(G_{t}).$ (6)

For any $n$ times $t_{1}<t_{2}<\dots<t_{n}$ the sets
$\left\\{G_{t_{i}}\right\\}$ partition ${\mathbb{R}}^{2}$ into $n(n+1)/2$ sets
of finite area (and one with infinite area, $(\cup G_{t_{i}})^{c}$), so each
$X_{t_{i}}$ can be written as the sum of some subset of $n(n+1)/2$ independent
Gamma random variables. In particular, any $n=2$ variables $X_{s}$ and $X_{t}$
can be written as

$X_{s}={\mathcal{G}}(G_{s}\backslash G_{t})+{\mathcal{G}}(G_{s}\cap
G_{t}),\qquad X_{t}={\mathcal{G}}(G_{t}\backslash
G_{s})+{\mathcal{G}}(G_{s}\cap G_{t})$

just as in the thinning approach of Section (2.2), so both 1-dimensional and
2-dimensional marginal distributions for the random measure process coincide
with those for the thinning process. Again the joint chf is

$\displaystyle\chi(s,t)$ $\displaystyle={\mathsf{E}}\exp(isX_{0}+itX_{1})$
$\displaystyle=(1-is/\beta)^{-\alpha{\bar{\rho}}}(1-i(s+t)/\beta)^{-\alpha\rho}(1-it/\beta)^{-\alpha{\bar{\rho}}}$
(7)

and the autocorrelation is
${\mathsf{Corr}}\big{(}X_{s},X_{t})=\exp\big{(}-\lambda|s-t|\big{)}$ or, for
integer times, $\rho^{|s-t|}$ for $\rho:=\exp(-\lambda)$. The distribution for
consecutive _triplets_ differs from those of the Thinned Gamma Process,
however, an illustration that the thinning process is Markov but the random
measure is not. The Random Measure process does feature infinitely-divisible
(ID) marginal distributions of all orders, while the thinned process does not.

###  2.4 The Markov change-point Gamma Process

Let
$\\{\zeta_{n}:~{}n\in{\mathbb{Z}}\\}\mathrel{\mathop{\sim}\limits^{\mathrm{iid}}}\mathsf{Ga}(\alpha,\beta)$
be iid Gamma random variables and let $N_{t}$ be a standard Poisson process
indexed by $t\in{\mathbb{R}}$ (so $N_{0}=0$ and
$(N_{t}-N_{s})\sim\mathsf{Po}(t-s)$ for all $-\infty<s<t<\infty$, with
independent increments), and set

$X_{t}:=\zeta_{n},\qquad n={N_{\lambda t}}.$

Then each $X_{t}\sim\mathsf{Ga}(\alpha,\beta)$ and, for $s,t\in{\mathbb{R}}$,
$X_{s}$ and $X_{t}$ are either identical (with probability $\rho^{|s-t|}$) or
independent— reminiscent of a Metropolis MCMC chain. The chf is

$\displaystyle\chi(s,t)$ $\displaystyle={\mathsf{E}}\exp(isX_{0}+itX_{1})$
$\displaystyle=\rho\big{(}1-i(s+t)/\beta\big{)}^{-\alpha}+{\bar{\rho}}(1-is/\beta)^{-\alpha}(1-it/\beta)^{-\alpha}$
(8)

and once again the marginal distribution is
$X_{t}\sim\mathsf{Ga}(\alpha,\beta)$ and the autocorrelation function is
${\mathsf{Corr}}\big{(}X_{s},X_{t})=\rho^{|s-t|}$.

###  2.5 The Squared O-U Gamma Diffusion

Let
$\\{Z_{i}\\}\mathrel{\mathop{\sim}\limits^{\mathrm{iid}}}\textsf{OU}(\lambda/2,1)$
be independent Ornstein-Uhlenbeck velocity processes, mean-zero Gaussian
processes with covariance
${\mathsf{Cov}}\big{(}Z_{i}(s),Z_{j}(t)\big{)}=\exp\big{(}-{\frac{\lambda}{2}}|s-t|\big{)}\delta_{ij}$,
and set

$X_{t}:=\frac{1}{2\beta}\sum_{i=1}^{n}Z_{i}(t)^{2}$

for $n\in{\mathbb{N}}$ and $\beta\in{\mathbb{R}}_{+}$. Note
$Z_{i}(t)\sim\mathsf{No}(0,1)$, so ${\mathsf{E}}Z_{i}(t)^{2}=1$ and
${\mathsf{E}}Z_{i}(t)^{4}=3$; it follows that
${\mathsf{E}}Z_{i}(s)^{2}Z_{i}(t)^{2}=1+2\exp(-\lambda|s-t|)$. Then
$X_{t}\sim\mathsf{Ga}(\alpha,\beta)$ for $\alpha=n/2$, with
${\mathsf{E}}X_{s}=\alpha/\beta$ and

$\displaystyle{\mathsf{E}}X_{s}X_{t}$
$\displaystyle=\frac{1}{4\beta^{2}}\left\\{n\big{(}1+2\exp(-\lambda|s-t|)\big{)}+n(n-1)\right\\}$
$\displaystyle=\frac{\alpha^{2}+\alpha\exp\big{(}-\lambda|s-t|\big{)}}{\beta^{2}}$

so the autocovariance is
${\mathsf{Cov}}\big{(}X_{s},X_{t})=\frac{\alpha}{\beta^{2}}~{}e^{-\lambda|s-t|}$
and the autocorrelation at integer times is
${\mathsf{Corr}}\big{(}X_{s},X_{t})=\rho^{|s-t|}$ for $\rho:=\exp(-\lambda)$.
The chf at consecutive integer times is

$\displaystyle\chi(s,t)$ $\displaystyle={\mathsf{E}}\exp(isX_{0}+itX_{1})$
$\displaystyle=\big{(}1-i(s+t)/\beta-st(1-\rho)/\beta^{2}\big{)}^{-\alpha},$
(9)

distinct from (4), (5)=(7), and (8), so this process is new. Itô’s formula is
used in [Wolpert and Brown, 2021] to show that $X_{t}$ has stochastic
differential equation (SDE) representation

$\displaystyle X_{t}$
$\displaystyle=X_{0}-\int_{0}^{t}2\lambda\big{(}X_{s}-\alpha/\beta\big{)}\,ds+2\sqrt{\lambda/\beta}\int_{0}^{t}\sqrt{X_{s}}\,dW_{s}$
(10)

and hence has generator
$\mathfrak{A}\phi(x)=(\partial/\partial\epsilon){\mathsf{E}}[\phi(X_{t+\epsilon})\mid
X_{t}=x]\big{|}_{\epsilon=0}$ given by

$\displaystyle\mathfrak{A}\phi(x)$
$\displaystyle=-2\lambda(x-\alpha/\beta)\phi^{\prime}(x)+(2\lambda/\beta)x\phi^{\prime\prime}(x),$
(11)

which we will use to distinguish this process from that of Section (2.6).
While the construction above required half-integer values for $\alpha$, (9) is
positive-definite and the SDE (10) has a unique strong solution for all
$\alpha>0$, so a time-reversible stationary Markov diffusion process exists
with this distribution. Cox et al. [1985, Eqn (17)], building on [Feller,
1951], found the transition kernel for this process:

$\displaystyle f(y,t\mid x,s)$
$\displaystyle=ce^{-u-v}\Big{(}\frac{v}{u}\Big{)}^{(\alpha-1)/2}I_{(\alpha-1)}\big{(}2\sqrt{uv}\big{)},$

where

$c:=\beta/\big{[}1-e^{-2\lambda|t-s|}\big{]},\qquad
u:=cye^{-2\lambda|t-s|},\qquad v:=cx.$

where $I_{q}(x)$ denotes the modified Bessel function of the first kind of
order $q$. With this one can construct likelihood functions and generate
posterior samples, MLEs, _etc_. for the parameters of this process. Cox et al.
show that the process $X_{t}$ is strictly positive at all times if $\alpha\geq
1$, but occasionally reaches zero for $\alpha<1$.

###  2.6 Continuously Thinned Gamma Process

Pick a large integer $n$ and set $\epsilon:=1/n$, $q:=\exp(-\lambda\epsilon)$,
and $p:=1-q=\lambda\epsilon+o(\epsilon)$. Draw
$X_{0}\sim\mathsf{Ga}(\alpha,\beta)$ and, for integers $i,j\in{\mathbb{N}}$,
draw independently

$\zeta_{i}\sim\mathsf{Ga}(\alpha p,\beta)\qquad b_{j}\sim\mathsf{Be}(\alpha
p,\alpha q).$

Set:

$\displaystyle X_{0}$ $\displaystyle=X_{0}$ $\displaystyle X_{\epsilon}$
$\displaystyle=X_{0}(1-b_{1})$ $\displaystyle\,+\,$ $\displaystyle\zeta_{1}$
$\displaystyle X_{2\epsilon}$ $\displaystyle=X_{\epsilon}\,(1-b_{2})$
$\displaystyle\,+\,$ $\displaystyle\zeta_{2}$
$\displaystyle=X_{0}(1-b_{1})(1-b_{2})$ $\displaystyle\,+\,$
$\displaystyle\zeta_{1}(1-b_{2})$ $\displaystyle\,+\,$
$\displaystyle\zeta_{2}$ $\displaystyle X_{3\epsilon}$
$\displaystyle=X_{0}(1-b_{1})(1-b_{2})(1-b_{3})$ $\displaystyle\,+\,$
$\displaystyle\zeta_{1}(1-b_{2})(1-b_{3})$ $\displaystyle\,+\,$
$\displaystyle\zeta_{2}(1-b_{3})$ $\displaystyle\,+\,$
$\displaystyle\zeta_{3}$

and, in general,

$\displaystyle X_{k\epsilon}$
$\displaystyle=X_{0}\prod_{j=1}^{k}(1-b_{j})+\sum_{i=1}^{k}\left\\{\zeta_{i}\prod_{i<j\leq
k}(1-b_{j})\right\\}.$ (12a) In the limit as $n\to\infty$ and $k\epsilon\to t$
the products converge to the product integral of the beta process introduced
by Hjort [1980, §3] [and described lucidly by Thibaux and Jordan, 2007, §2]
and the sum to an ordinary gamma stochastic integral, $\displaystyle X_{t}$
$\displaystyle=X_{0}\prod_{s\in(0,t]}[1-dB(s)]+\int_{0}^{t}\left\\{\prod_{s\in(r,t]}[1-dB(s)]\right\\}\,\zeta(dr),$
(12b)

where $\zeta(dr)\sim\mathsf{Ga}\big{(}\alpha\lambda\,dr,\beta\big{)}$ is a
Gamma random measure and $B(s)\sim\textsf{BP}\big{(}\alpha,\lambda\,ds\big{)}$
is a Beta process, _i.e._ , an SII Lévy process with Lévy measure

$\nu_{B}(du)=\lambda\alpha\,u^{-1}(1-u)^{\alpha-1}\,\mathbf{1}_{\\{0<u<1\\}}\,du$

with constant “concentration function” $\alpha(s)\equiv\alpha$ and
translation-invariant “base measure” $\lambda(ds)=\lambda ds$. The product
integral can be written as a ratio

$\displaystyle\prod_{s\in(r,t]}[1-dB(s)]$
$\displaystyle=\frac{1-F(t)}{1-F(r)}$

where $F(t)=\prod_{s\in(t,\infty)}[1-dB(s)]$ satisfies

$\displaystyle\frac{dF(t)}{1-F(t)}$ $\displaystyle=dB(t),\qquad
B_{t}=\int_{(0,t]}\frac{dF(s)}{1-F(s)}.$ (13)

####  2.6.1 Generator

The Gamma process of Eqn (12b) is stationary and Markov, with generator

$\displaystyle\mathfrak{A}\phi(x)=$
$\displaystyle\int_{0}^{\infty}\big{[}\phi(x+u)-\phi(x)\big{]}\,\alpha\lambda
u^{-1}e^{-\beta u}\,du$ $\displaystyle+$
$\displaystyle\int_{0}^{x}\big{[}\phi(x-u)-\phi(x)\big{]}\,\alpha\lambda
u^{-1}(1-u/x)^{\alpha-1}\,du$ (14)

Because this differs from (11) (and in particular because it is a non-local
operator, showing $X_{t}$ has jumps), this process is new. Once again the one-
dimensional marginal distributions are $X_{t}\sim\mathsf{Ga}(\alpha,\beta)$
and the autocorrelation is
${\mathsf{Corr}}(X_{s},X_{t})=\exp\big{(}-\lambda|s-t|\big{)}$.

##  3 Discussion

We have now constructed six distinct processes that share the same univariate
marginal distribution and autocorrelation function, but which all differ in
their $n$-variate marginal distributions for $n\geq 3$. Some are Markov, some
not; some are time-reversible, some not; some have ID marginal distributions
of all orders, some don’t. Similar methods can be used to construct
$\textsf{AR}(1)$-like processes with _any_ ID marginal distribution, such as
those listed on _p._ ID Distributions; only in the two cases of Gaussian and
Poisson do all these constructions coincide. Many of these are useful for
modeling time-dependent phenomena whose dependence falls off over time, but
for which traditional Gaussian methods are unsuitable because of heavy tails,
or integer values, or positivity, or for other reasons. I know of very little
work (yet!) exploring inference for processes like these; a beginning appears
in [Wang, 2013, §3]. Code in R to generate samples from each of these six
processes is available on request.

The author would like to thank Lawrence D. Brown for many fruitful discussions
about ID processes, and Victor Peña for helping develop methods of
distinguishing these six processes from observations.

## Appendix

######  Proposition 1 (Walker [2000]).

The innovations $\zeta_{t}$ in Eqn (3) can be constructed successively as
follows:

$\displaystyle\lambda_{t}\sim\mathsf{Ga}(\alpha,1),\qquad
N_{t}\mid\lambda_{t}\sim\mathsf{Po}\big{(}{\textstyle\frac{1-\rho}{\rho}}\lambda_{t}\big{)},\qquad\zeta_{t}\mid
N_{t}$
$\displaystyle\sim\mathsf{Ga}\big{(}N_{t},{\textstyle\frac{\beta}{\rho}}\big{)}.$

###### Proof.

For $\omega\in{\mathbb{R}}$,

$\displaystyle{\mathsf{E}}e^{i\zeta_{t}\omega}$ $\displaystyle\mid$
$\displaystyle N_{t}$ $\displaystyle=(1-i\omega\rho/\beta)^{-N_{t}}$
$\displaystyle{\mathsf{E}}e^{i\zeta_{t}\omega}$ $\displaystyle\mid$
$\displaystyle\lambda_{t}$
$\displaystyle=\sum_{n=0}^{\infty}(1-i\omega\rho/\beta)^{-n}\big{(}{\textstyle\frac{1-\rho}{\rho}}\lambda_{t}\big{)}^{n}\exp\big{(}-{\textstyle\frac{1-\rho}{\rho}}\lambda_{t}\big{)}/n!$
$\displaystyle=\exp\Big{(}i\frac{(1-\rho)\omega}{\beta-i\rho\omega}\lambda_{t}\Big{)}$
$\displaystyle{\mathsf{E}}e^{i\zeta_{t}\omega}$
$\displaystyle=\Big{[}1-i\frac{(1-\rho)\omega}{\beta-i\rho\omega}\Big{]}^{-\alpha}=\Big{[}\frac{\beta-i\omega}{\beta-i\rho\omega}\Big{]}^{-\alpha}.$

∎

### Poisson and Gamma SII Processes

The chf for a Poisson random variable $X\sim\mathsf{Po}(\nu)$ is

$\displaystyle\chi_{X}(\theta)={\mathsf{E}}[e^{i\theta X}]$
$\displaystyle=\sum_{k=0}^{\infty}e^{i\theta
k}\left\\{\frac{\nu^{k}}{k!}e^{-\nu}\right\\}=e^{(e^{i\theta}-1)\nu}$ so for
any $u\in{\mathbb{R}}$ the re-scaled random variable $Y:=uX$ has chf
$\displaystyle\chi_{uX}(\theta)$ $\displaystyle={\mathsf{E}}[e^{i\theta
uX}]=\chi_{X}(u\theta)$ $\displaystyle=e^{(e^{i\theta u}-1)\nu}$ and a linear
combination $Y:=\sum u_{j}X_{j}$ of independent
$X_{j}\sim\mathsf{Po}(\nu_{j})$ has chf $\displaystyle\chi_{Y}(\theta)$
$\displaystyle=\prod_{j}\left\\{e^{(e^{i\theta u_{j}}-1)\nu_{j}}\right\\}$
$\displaystyle=\exp\left\\{\sum_{j}{(e^{i\theta u_{j}}-1)\nu_{j}}\right\\}$
(15a) $\displaystyle=\exp\left\\{\int_{\mathbb{R}}{(e^{i\theta
u}-1)\nu(du)}\right\\}$ (15b)

for the discrete measure

$\nu(du)=\sum\nu_{j}\delta_{u_{j}}(du)$

that assigns mass $\nu_{j}$ to each point $u_{j}$, provided the sum in (15a)
converges. Of course the sum converges if it has only finitely-many terms, or
even if there are infinitely-many with $\sum\nu_{j}<\infty$ (because
$|e^{i\theta u}-1|\leq 2$), but that condition isn’t actually necessary. Since
also $|e^{i\theta u}-1|\leq|\theta u|$, the random variable $Y$ will be well-
defined and finite provided that

$\displaystyle\sum_{j}\big{(}1\wedge|u_{j}|\big{)}\nu_{j}$
$\displaystyle<\infty$ (16a) or, in integral form,
$\displaystyle\int_{\mathbb{R}}\big{(}1\wedge|u|\big{)}\nu(du)$
$\displaystyle<\infty.$ (16b)

A random variable with chf of form (15a) for a sequence satisfying (16a) has a
“compound Poisson” distribution; one with the more general chf of form (15b)
for a measure satisfying (16b) is called Infinitely Divisible, or ID.

#### ID Distributions

For any $\sigma$-finite measure satisfying (16b) it’s easy to make a
stochastic process with stationary independent increments of the more general
form of (15b), beginning with a Poisson random measure $\mathcal{N}(du\,ds)$
on ${\mathbb{R}}\times{\mathbb{R}}_{+}$ with intensity measure
${\mathsf{E}}\mathcal{N}(du\,ds)=\nu(du)\,ds$:

$\displaystyle X_{t}$
$\displaystyle:=\iint_{{\mathbb{R}}\times(0,t]}u\mathcal{N}(du\,ds).$ (17)

This is a right-continuous independent-increment nondecreasing process that
begins at $X_{0}=0$ with jumps
$\Delta_{t}=[X_{t}-X_{t-}]\equiv[X_{t}-\lim_{s\nearrow t}X_{s}]$ of magnitudes
$\Delta_{t}\in E$ at rate $\nu(E)$ for any Borel $E\subset{\mathbb{R}}$. The
Poisson process itself is the special case where
$\nu(E)=\nu\mathbf{1}_{\\{1\in E\\}}$, with jumps of magnitude $\Delta_{t}=1$
at rate $\nu\in{\mathbb{R}}_{+}$.

Khinchine and Lévy [1936] showed that a random variable $Y$ has a chf of the
form (15b) if and only111For nonnegative random variables this is true as
stated, but a slightly more general form is necessary for real-valued ID
random variables, with a condition on $\nu(du)$ somewhat weaker than (16b)
(see (18c))— if you get interested, ask me about “compensation”. This is
needed for the $\mathsf{Ca}(\delta,\gamma)$ example and, for $\alpha\geq 1$,
the $\mathsf{St}_{0}(\alpha,\beta,\gamma,\delta)$ and
$\textsf{S$\alpha$S}(\alpha,\gamma)$ examples below. if, for every
$n\in{\mathbb{N}}$, one can write $Y=\zeta_{1}+\cdots\zeta_{n}$ as the sum of
$n$ iid random variables $\zeta_{j}$. This property is called “Infinite
Divisibility” (abbreviated ID), and the processes we have constructed whose
increments have this property are called “SII” processes for their stationary
independent increments. Examples of ID distributions (or SII processes) and
their Lévy measures include:

Poisson $\mathsf{Po}(\lambda)$ $\nu(du)=\lambda\delta_{1}(du)$ Negative
Binomial $\mathsf{NB}(\alpha,p)$
$\nu(du)=\sum_{k\in{\mathbb{N}}}\alpha\frac{q^{k}}{k}\delta_{k}(du),\qquad
q:=(1-p)$ Gamma $\mathsf{Ga}(\alpha,\beta)$ $\nu(du)=\alpha u^{-1}e^{-\beta
u}\mathbf{1}_{\\{u>0\\}}\,du$ $\alpha$-Stable
$\mathsf{St}_{0}(\alpha,\beta,\gamma,\delta)$
$\nu(du)=\frac{\alpha\gamma}{\pi}\Gamma(\alpha)\sin\frac{\pi\alpha}{2}|u|^{-\alpha-1}(1+\beta\mathop{\mathrm{sgn}}u)\,du$
Symmetric $\alpha$-Stable $\textsf{S$\alpha$S}(\alpha,\gamma)$
$\nu(du)=\frac{\alpha\gamma}{\pi}\Gamma(\alpha)\sin\frac{\pi\alpha}{2}|u|^{-\alpha-1}\,du$
Cauchy $\mathsf{Ca}(\delta,\gamma)$ $\nu(du)=\frac{\gamma}{\pi}|u|^{-2}\,du.$

The defining condition for a random variable $Y$ to be “Infinitely Divisible”
(ID) is that for each $n\in{\mathbb{N}}$ there must exist iid random variables
$\\{\zeta_{j}:~{}1\leq j\leq n\\}$ such that $Y$ and $\sum_{j=1}^{n}\zeta_{j}$
have the same distribution. This is clearly equivalent to the condition that
every power $\chi^{\alpha}(\omega)$ of the characteristic function
$\chi(\theta):={\mathsf{E}}\exp(i\theta Y)$ must _also_ be a characteristic
function (_i.e._ , must be positive-definite) for each inverse integer
$\alpha=1/n$, because we can just take $\\{\zeta_{j}\\}$ to be iid with chf
$\chi^{1/n}(\omega)$ and verify that their sum has chf $\chi(\omega)$. Less
obvious but also true is that $Y$ is ID if and only if $\chi^{\alpha}(\omega)$
is positive-definite for _all_ real $\alpha>0$, and even less obvious is the
Khinchine and Lévy theorem that $\chi$ must take the specific form

$\displaystyle\chi(\theta)$
$\displaystyle=\exp\left\\{i\theta\delta-\theta^{2}\sigma^{2}/2+\int_{\mathbb{R}}\left[e^{i\theta
u}-1\right]\,\nu(du)\right\\}$ (18a) for some $\delta\in{\mathbb{R}}$,
$\sigma^{2}\geq 0$, and Borel measure $\nu$ satisfying $\nu(\\{0\\})=0$ and
(16b) or, a little more generally, the form $\displaystyle\chi(\theta)$
$\displaystyle=\exp\left\\{i\theta\delta-\theta^{2}\sigma^{2}/2+\int_{\mathbb{R}}\left[e^{i\theta
u}-1-i\theta h(u)\right]\,\nu(du)\right\\}$ (18b) for any bounded function $h$
that satisfies $h(u)=u+O(u^{2})$ near $u\approx 0$ (like $\arctan u$ or
$u\mathbf{1}_{\\{|u|<1\\}}$ or $u/(1+u^{2})$) and a Borel measure $\nu$
satisfying the weaker restriction
$\displaystyle\int_{\mathbb{R}}\big{(}1\wedge u^{2}\big{)}\nu(du)$
$\displaystyle<\infty.$ (18c)

Some properties of ID distributions beyond the scope of these notes, but
covered in [Steutel and van Harn, 2004] (see also [Bose et al., 2002]),
include:

######  Theorem 1 (S&vH, Thm 2.13).

Let $\chi(\theta)$ be the chf of an Infinitely Divisible distribution. Then
$(\forall\theta\in{\mathbb{R}})\left\\{\chi(\theta)\neq 0\right\\}$. Also, if
$\chi(\theta)$ is analytic in some open domain $\Omega\subset{\mathbb{C}}$,
then $(\forall\theta\in\Omega)\left\\{\chi(\theta)\neq 0\right\\}$. Thus, ID
chfs do not vanish on ${\mathbb{R}}$ or anywhere in ${\mathbb{C}}$ that
$\chi(\theta)$ is analytic.

######  Theorem 2 (S&vH, Thm 9.8).

Let $X$ be an infinitely divisible random variable that is not normal or
degenerate. Then the two-sided tail of $X$ satisfies

$\lim_{x\to\infty}\frac{-\log{\mathsf{P}}[|X|>x]}{x\log x}=c$

for a number $0\leq c<\infty$ given by
$c^{-1}=\max[\nu({\mathbb{R}}_{+}),\nu({\mathbb{R}}_{-})]$. An ID random
variable that is not degenerate has a normal distribution if and only if the
same limit is $c=\infty$.

######  Theorem 3 (S&vH, Prop 2.3).

No non-degenerate bounded random variable is ID

This one’s easy enough to prove. If $\|X\|_{\infty}=B<\infty$ and $X$ has the
same distribution as $\sum_{j=1}^{n}\zeta_{j}$ for iid $\\{\zeta_{j}\\}$ then
$\|\zeta_{j}\|_{\infty}=B/n$ and
$\sigma^{2}:={\mathsf{V}}(X)=n{\mathsf{V}}(\zeta_{j})\leq
n{\mathsf{E}}\zeta_{j}^{2}\leq B^{2}/n$, so $\sigma^{2}=0$ and $X$ must be
degenerate.

#### Gamma Variables & Processes

The Gamma distribution $X\sim\mathsf{Ga}(\alpha,\beta)$ with mean
$\alpha/\beta$ and variance $\alpha/\beta^{2}$ has chf

$\displaystyle\chi_{X}(\theta)$ $\displaystyle={\mathsf{E}}[e^{i\theta X}]$
$\displaystyle=\int_{0}^{\infty}e^{i\theta
x}\left\\{\frac{\beta^{\alpha}}{\Gamma(\alpha)}x^{\alpha-1}e^{-\beta
x}\right\\}\,dx$ $\displaystyle=\big{(}1-i\theta/\beta\big{)}^{-\alpha}$
$\displaystyle=\exp\left\\{-\alpha\log(1-i\theta/\beta)\right\\}$
$\displaystyle=\exp\left\\{-\int_{0}^{\infty}\big{(}e^{i\theta
u}-1\big{)}\alpha u^{-1}e^{-\beta u}\,du\right\\},$ exactly of the form of
(15b) with Lévy measure $\displaystyle\nu(du)$ $\displaystyle=\alpha
u^{-1}e^{-\beta u}\mathbf{1}_{\\{u>0\\}}\,du.$

This measure, while infinite, does satisfy condition (16b):

$\displaystyle\int_{\mathbb{R}}\big{(}1\wedge|u|\big{)}\nu(du)$
$\displaystyle=\int_{0}^{1}\big{(}|u|\big{)}\alpha u^{-1}e^{-\beta
u}\,du+\int_{1}^{\infty}\big{(}1\big{)}\alpha u^{-1}e^{-\beta u}\,du$
$\displaystyle\leq\int_{0}^{\infty}\alpha e^{-\beta
u}\,du=\alpha/\beta<\infty.$

Thus gamma-distributed random variables are ID and the SII gamma _process_

$X_{t}=\iint_{{\mathbb{R}}\times(0,t]}u\mathcal{N}(du\,ds)$

has infinitely-many non-negative “jumps” $\Delta_{t}=X_{t}-X_{t-}$ in any time
interval $a<t\leq b$. Their sum $\sum\\{\Delta_{t}:~{}a<t\leq
b\\}=X_{b}-X_{a}$ is finite, however, with probability distribution

$X_{b}-X_{a}\sim\mathsf{Ga}\big{(}\alpha(b-a),~{}\beta\big{)}.$

The gamma _random measure_
${\mathcal{G}}(dx\,dy)\sim\mathsf{Ga}(\alpha\,dx\,dy,~{}\beta)$ of Section
(2.3) has a similar construction,

${\mathcal{G}}(A)=\iiint_{{\mathbb{R}}\times A}u\mathcal{N}(du\,dx\,dy)$

as the sum of the heights $u_{j}$ of a Poisson cloud of points
$(u_{j},x_{j},y_{j})$ for which $(x_{j},y_{j})\in A$. Wolpert and Ickstadt
[1998] show how to simulate such random measures very efficiently, drawing the
jumps $\\{u_{j}\\}$ in monotone decreasing order.

## References

  * Bose et al. [2002] Arup Bose, Anirban DasGupta, and Herman Rubin.  A contemporary review and bibliography of infinitely divisible distributions and processes.  _Sankhyā, Ser. A_, 64(3):763–819, 2002\.  doi: 10.2307/25051430.  Special issue in memory of D. Basu.
  * Cox et al. [1985] John C. Cox, Jonathan E. Ingersoll, Jr., and Stephen A. Ross.  A theory of the term structure of interest rates.  _Econometrica_ , 53(2):385–408, 1985.  doi: 10.2307/1911242.
  * Feller [1951] William Feller.  Two singular diffusion problems.  _Annals of Mathematics_ , 54(1):173–182, 1951\.  doi: 10.2307/1969318.
  * Hjort [1980] Nils Lid Hjort.  Nonparametric Bayes estimators based on beta processes in models for life history data.  _Ann. Stat._ , 18(3):1259–1294, 1980.  doi: 10.1214/aos/1176347749.
  * Khinchine and Lévy [1936] Alexander Ya. Khinchine and Paul Lévy.  Sur les lois stables.  _Comptes rendus hebdomadaires des seances de l’Académie des sciences. Académie des science (France). Serie A. Paris_ , 202:374–376, 1936.
  * Steutel and van Harn [2004] Frederik W. Steutel and Klaas van Harn.  _Infinite Divisibility of Probability Distributions on the Real Line_ , volume 259 of _Monographs and Textbooks in Pure and Applied Mathematics_.  Marcel Dekker, New York, NY, 2004.  ISBN 0-8247-0724-9.  doi: 10.1201/9780203014127.
  * Thibaux and Jordan [2007] Romain Thibaux and Michael I. Jordan.  Hierarchical beta processes and the Indian buffet process.  In Marina Meila and Xiaotong Shen, editors, _Proceedings of the Eleventh International Conference on Artificial Intelligence and Statistics (AISTATS 2007), March 21–24, San Juan, Puerto Rico_ , volume 2, pages 564–571. PMLR, 2007.  URL http://proceedings.mlr.press/v2/thibaux07a.html.
  * Walker [2000] Stephen G. Walker.  A note on the innovation distribution of a gamma distributed autoregressive process.  _Scand. J. Stat._ , 27(4):575–576, 2000.  doi: 10.1111/1467-9469.00208.
  * Wang [2013] Jianyu Wang.  _Bayesian Modeling and Adaptive Monte Carlo with Geophysics Applications_.  PhD thesis, Duke University Statistical Science Department, 2013.
  * Wolpert and Brown [2021] Robert L. Wolpert and Lawrence D. Brown.  Markov infinitely-divisible stationary time-reversible integer-valued processess.  On arXiv.
  * Wolpert and Ickstadt [1998] Robert L. Wolpert and Katja Ickstadt.  Poisson/gamma random field models for spatial statistics.  _Biometrika_ , 85(2):251–267, 1998.  doi: 10.1093/biomet/85.2.251.


# Fast counter-diabatic Thouless pumping in the Rice-Mele model

Joshua Chiel  Department of Physics, University of Maryland, College Park,
Maryland 20742-4111, US  Christopher Jarzynski  Department of Physics,
University of Maryland, College Park, Maryland 20742-4111, US  Department of
Chemistry and Biochemistry, University of Maryland, College Park, MD 20742 USA
Institute for Physical Science and Technology, University of Maryland, College
Park, MD 20742 USA  Jay Sau  Joint Quantum Institute and Condensed Matter
Theory Center, Department of Physics, University of Maryland, College Park,
Maryland 20742-4111, US

###### Abstract

Thouless pumping is a transport phenomenon where a periodically varying
Hamiltonian can transfer a quantized amount of charge when the time-dependence
of the Hamiltonian is quasi-adiabatic. Past proposals to speed up this process
involving Floquet techniques lead to a subtle problem of setting the initial
state of the system. In this work we apply counter-diabatic driving to the
Rice-Mele model, which is one of the simplest models for Thouless pumping, to
ensure that the system remains in the ground state for any driving speed. We
show that the pumped charge across each bond of the Rice-Mele model is given
by a topologically quantized Chern number in this case. However, the counter-
diabatic driving in a general case turns out to involve long-range hopping. We
show that this can be mitigated either by choosing a very specific example of
the Rice-Mele model or by numerical optimization of the Hamiltonian to create
experimentally realizable variants of fast pumping in the Rice-Mele model.

##  I Introduction

Topological Thouless pumps provide dissipationless transport that is robustly
and exactly quantized [1]. A Thouless pump is driven by a time-periodic
Hamiltonian, and the charge that is pumped per driving cycle is given by a
Chern number associated with the system’s ground state. In its original
formulation, a Thouless pump must be driven adiabatically, i.e. slowly
relative to the time scale set by the ground state gap. It is natural to ask
whether the same result – quantized, dissipationless transport – can be
achieved without the requirement of slow driving.

A growing field of study addresses this question. One approach is to take
advantage of the time-periodicity of the Hamiltonian in the Thouless scheme
and view it as a Floquet Hamiltonian. In this case, it was shown that states
in a gapped Floquet band can exhibit quantized charge pumping [2]. Similar
fast pumping can occur in the anomalous Floquet-Anderson insulator [3].
However, these cases require fine-tuning of initial conditions to occupy
quasi-energy bands, which generally differ from the ground state of the
Hamiltonian. Explicitly including damping can also lead to fast pumping in
open dissipative systems [4] governed by a non-Hermitian generalization of the
Rice-Mele Hamiltonian. Both the Floquet Hamiltonian approach and the open
system setting differ from the case considered originally by Thouless, in
which a closed quantum system begins in the ground state of a static
Hamiltonian, then undergoes quantized pumping during one cycle of driving, and
finally returns to the ground state.

Shortcuts to adiabaticity (STA) offer a natural, additional line of attack to
speed up transport and preserve a given state (in this case the ground state)
of a system throughout the system’s evolution [5]. In the counterdiabatic (CD)
approach, given a time-dependent Hamiltonian $H_{0}(t)$ and a system
initialized in an arbitrary instantaneous eigenstate of $H_{0}(t)$, adding an
appropriately designed auxiliary term $H_{CD}(t)$ causes the system to remain
in this eigenstate of $H_{0}(t)$ at all times. That is, the CD Hamiltonian
$H_{CD}(t)$ suppesses transitions to other eigenstates [6, 7, 8]. CD methods
have also been employed to speed up classical geometric pumps [9, 10] and have
been applied to a Cooper-pair pump, which is a bosonic system with an
effectively single phase quantum degree of freedom [11].

In this work, we apply CD methods to the standard Thouless pump i.e. a non-
interacting fermionic chain in one dimension where the pumping protocol starts
in the ground state of an initial Hamiltonian. The system then evolves under
the Hamiltonian $H_{0}(t)+H_{CD}(t)$ as $H_{0}(t)$ varies periodically in time
for a finite number of periods. By adding $H_{CD}$ to $H_{0}$ we preserve the
filled single particle Bloch states of $H_{0}$ throughout the system’s
evolution. We will show that for any gapped “original” Hamiltonian $H_{0}(t)$,
the charge pumped across each bond of the system (see Fig. 1) is given by a
Chern number of the ground state density matrix of $H_{0}(t)$ and is therefore
quantized in a topologically robust way. While CD driving can produce fast
quantized particle transport arbitrarily fast, for an appropriate choice of
$H_{CD}(t)$, the resulting $H_{CD}$ can potentially have arbitrarily long-
range hopping terms, that limits its experimental realizability. We show that
this problem can be avoided using the so-called “control-freak” protocol [12]
of the Rice-Mele (RM) model [13] shown in Fig. 1. This protocol provides a
specific, though fine-tuned, example of $H_{0}(t)$ where $H_{CD}(t)$ contains
only nearest neighbor hopping terms and is therefore realizable. We show that
the charge pumping in this protocol is given by a quantized Chern number. For
the case of more general Hamiltonians $H_{0}(t)$, we find an approximate
solution to $H_{CD}$ which involves only nearest neighbor hopping using a
numerical optimization procedure.

This paper is structured as follows. In section II we review the Rice-Mele
model and derive an exact CD protocol for it. In section III we construct a
current density operator and show that the exact CD protocol preserves the
pumped charge over a cycle. In section IV we develop a general approach to
construct nearest neighbor Hamiltonians that will produce finite-time Thouless
pumping, and we conclude in section V by discussing our results.

##  II CD Rice-Mele Model

The Rice-Mele model provides a useful lattice model of non-interacting
Thouless pumping. When driven through a quasi-adiabatic cycle, the system
produces quantized particle transport. The model consists of a $1$d chain of
atoms, two per unit cell, with two types of bonds – intracell bonds and
intercell bonds. See Figure 1.

Figure 1: The Rice-Mele $1$d lattice model consists of a $1$d chain of two
types of atoms, with positions $x_{j,A}=aj$ (blue) and $x_{j,B}=a(j+1)$ (red),
connected by intracell bonds (d, double green lines) at positions
$y_{j,d}=a(j+1/2)$ and intercell bonds (s, single green lines) at positions
$y_{j,s}=a(j-1/2)$. The dashed box shows a unit cell. The intracell/intercell
hopping amplitudes are given by $J_{1/2}$.

In real space we have:

$\hat{H}_{0}(t)=-\sum_{j}\left(J_{1}(t)\hat{a}_{j}^{\dagger}\hat{b}_{j}+J_{2}(t)\hat{a}_{j+1}\hat{b}_{j}+h.c.\right)+\Delta(t)\sum_{j}\left(\hat{a}_{j}^{\dagger}\hat{a}_{j}-\hat{b}_{j}^{\dagger}\hat{b}_{j}\right)$
(1)

where $\hat{a}_{j}^{\dagger}(\hat{a}_{j})$ and
$\hat{b}_{j}^{\dagger}(\hat{b}_{j})$ are the creation (annihilation) operators
acting on the left and right sites of the $j$th unit cell. In the position
representation, the two locations of physical atoms per unit cell are
$x_{j,A}=aj$ and $x_{j,B}=a(j+1)$. We also write the locations of the
intracell bonds as $y_{j,d}=a(j+1/2)$ and the intercell bonds as
$y_{j,s}=a(j-1/2)$. The time dependent coefficients $J_{1}(t)$, $J_{2}(t)$ and
$\Delta(t)$ tune the intracell hopping amplitude, intercell hopping amplitude,
and a staggered on-site potential respectively. In momentum space, the
Hamiltonian can be written as

$\hat{H}_{0}(k,t)=\bm{R}(k,t)\cdot\bm{\sigma},$ (2)

where $\bm{\sigma}$ are $2\times 2$ Pauli matrices acting on the space of
Bloch states $(\tilde{a}_{k}^{\dagger},\tilde{b}_{k}^{\dagger})^{T}$ with
$\tilde{a}_{k}^{\dagger},\tilde{b}_{k}^{\dagger}$ being Fourier transforms of
the orbital operators $a_{j}^{\dagger}$ and $b_{j}^{\dagger}$ respectively and

$\bm{R}(k,t)=\left\\{-J_{1}(t)-J_{2}(t)\cos(ka),-J_{2}(t)\sin(ka),\Delta(t)\right\\},$
(3)

where $J_{1/2}(t)$ and $\Delta(t)$ are periodic on the interval $[0,T]$. We
moreover assume that $R(k,t)$ is constant outside of $(0,T)$. We note that
this gives us a Floquet-Bloch Hamiltonian, which is gapped for all $k,t$ away
from $J_{1}(t)-J_{2}(t)=\Delta(t)=0$. We emphasize that our system will thus
begin and end in the ground state of $H_{0}$, rather than requiring an
initialization into a particular Floquet-Bloch state in the infinite past.

Eq. (2) has the same form as that considered by Berry[8] (see Eq. 3.1
therein). Since this is a two-band system, we can immediately evoke Berry’s
results and write down the corresponding $H_{CD}(k,t)$ in lattice-momentum
space:

$\hat{H}(k,t)=\hat{H}_{0}(k,t)+\hat{H}_{CD}(k,t)=\left(\bm{R}+\frac{1}{R^{2}}\bm{R}\crossproduct\partial_{t}\bm{R}\right)\cdot\bm{\sigma}\equiv\bm{u}\cdot\bm{\sigma}.$
(4)

where we define $\bm{u}(k,t)$, a vector which is periodic in $k$ and $t$. For
a system initialized in its single particle ground Bloch state of $H_{0}$, we
will show that in the presence of $H_{CD}$, the particle transport remains
quantized under arbitrarily fast driving. Note that $\hat{H}_{CD}$ is
manifestly periodic and smooth in $k$ and thus, by a well known theorem of
Fourier analysis, $\hat{H}_{CD}$ is exponentially localized in real space.

##  III Pumped Charge

Even granting the spatial, exponential-locality of $\hat{H}_{CD}$, it is not a
priori obvious how to define a current operator for an operator whose matrix
elements involve two distinct spatial points. We will employ a generalized
Peierls substitution method to define such a current and show it gives the
desired result.

We consider an arbitrary non-local Hamiltonian in one spatial dimension:

$i\partial_{t}\Psi(x)=\sum_{x^{\prime}}H(x,x^{\prime})\Psi(x^{\prime})$ (5)

where $x,x^{\prime}\in\\{aj:j\in\mathbb{Z}\\}$ represent positions of
atoms/sites. Motivated by the Peierls substitution method [14] framework, we
introduce a vector potential $A(y)$ on the bonds as follows:

$H(x,x^{\prime})[A]=H(x,x^{\prime})e^{i\textrm{sign}(x-x^{\prime})\sum_{x<y<x^{\prime}}A(y)}$
(6)

where $H(x,x^{\prime})$ is the non-local Hamiltonian introduced above, and $y$
live on bonds, i.e. $y_{s/d}\in\\{a(j/2\pm 1):j\in\mathbb{Z}\\}$ The current
operator in the tight-binding approximation can be written as:

$\bra{x}\hat{J}(y)\ket{x^{\prime}}=\left.\frac{\delta
H(x,x^{\prime})[A]}{\delta A(y)}\right|_{A\to
0}=i\Bigl{(}\Theta(y-x)-\Theta(y-x^{\prime})\Bigr{)}H(x,x^{\prime})$ (7)

where $\Theta$ is the Heaviside step function. Note that, consistent with Fig.
1, we label the position of the bonds with $y$, which is where the vector
potential connects sites $x,x^{\prime}$. In appendix A we show this current
density operator satisfies a continuity equation.

The total current operator of the system (with periodic boundary conditions)
$\hat{j}=a\sum_{y}\hat{J}(y)$ can be computed from the local current density
calculated in Eq. 7 as:

$\displaystyle\bra{x}\hat{j}\ket{x^{\prime}}=ia\sum_{y}\Bigl{(}\Theta(y-x)-\Theta(y-x^{\prime})\Bigr{)}H(x,x^{\prime})=i(x-x^{\prime})H(x,x^{\prime})$
(8)

for any pair of points $x,x^{\prime}$ away from the boundary $x=0,L$, which
matches the expression for average current in the crystal [15]. Using lattice
translational invariance, the total current over a unit-cell is written as

$\begin{split}j_{cell}(t)&=\sum_{R,x^{\prime}}\bra{x}\hat{j}\ket{x^{\prime}+R}\rho(x^{\prime}+R,x)\\\
&=\sum_{R,x^{\prime}}i(x-x^{\prime}+R)H(x,x^{\prime}+R)\rho(x^{\prime}+R,x)\end{split}$
(9)

where $\rho(x,x^{\prime}+R)$ is the density matrix of the occupied states.
Writing the real space Hamiltonian in terms of the Bloch Hamiltonian
$H(x,x^{\prime}+R)=\int\frac{dk}{2\pi}e^{ik(x-x^{\prime}+R)}H_{k}(x,x^{\prime})$
as well as the density matrix in the above equation and then using $H(k,t)$ in
Eq. (4) for the Bloch Hamiltonian together with
$\hat{\rho}_{k}(t)=\frac{1}{2}\left(\sigma_{0}-\hat{\bm{R}}(k,t)\cdot\bm{\sigma}\right)$
as the density matrix, the average current in the system simplifies to:

$\begin{split}j_{cell}(t)&=\int\frac{dk}{2\pi}Tr[\partial_{k}H(k,t)\rho_{k}^{\dagger}]=\int\frac{dk}{2\pi}\partial_{k}\bm{u}\cdot\hat{\bm{R}},\end{split}$
(10)

where $\hat{\bm{R}}=\bm{R}/R$.

The current $j_{cell}(t)$ is actually the total current over all bonds in the
unit-cell. The physical charge pumped across a unit-cell is the average over
the unit-cell in terms of bonds i.e.

$\displaystyle Q_{pump}(t)=\int_{0}^{t}d\tau j_{cell}(\tau)/N_{cell},$ (11)

where $N_{cell}$ is the number of sites in the unit-cell. For the Rice-Mele
model, $N_{cell}=2$. Integrating (11) by parts, and substituting $\bm{u}$ from
Eq.(4) we find

$Q_{pump}(t)=\frac{1}{4\pi}\int_{0}^{t}d\tau\int_{-\frac{\pi}{a}}^{\frac{\pi}{a}}dk\,\hat{\bm{R}}\cdot\left(\partial_{k}\hat{\bm{R}}\crossproduct\partial_{t}\hat{\bm{R}}\right).$
(12)

It is manifest from the final line of Eq. 12 that the pumped charge at the end
of the cycle $Q_{pump}(T)$ is the Chern number $C$ of the vector field
$\hat{\bm{R}}(k,t)$ over the torus and is an integer as long as
$\hat{\bm{R}}(k,t)$ is a continuous vector field on the torus due to
$\hat{\bm{R}}$’s periodic boundary conditions in $k$ and $t$[16].

Note that Eq. (8) matches the standard current operator form[15] which is for
the average current over a unit cell. In the quasi-adiabatic regime, charge
does not build up in the cell, so average current is the only relevant
quantity. Indeed, $Q_{pump}(t)$ is the charge pumped across the unit-cell and
does not represent the complete motion of the charge across bonds. This
information is provided by the time-dependent density of charge in the unit-
cell

$\begin{split}Q(x,t)=\int_{-\pi/a}^{\pi/a}\frac{dk}{2\pi}\rho_{k}(x,x)&=\int\frac{dk}{2\pi}(1-\hat{R}_{z}(k,t)s_{\sigma(x)})\end{split}$
(13)

where $\hat{R}_{z}$ is the $z$th component of $\hat{\bm{R}}$. For a given
cell, we take the atom sites to be at $x\in\\{0,a\\}$ and bond positions to be
$x_{b}\in\\{-a/2,a/2\\}$ and so $s_{\sigma}(0)=-1$ and $s_{\sigma}(a)=1$. The
charge pumped across any bond $b$, is $Q_{pump,b}(t)=\int
d\tau\bra{x}J(y_{b})\ket{x^{\prime}}$, where $J(y_{b})=Tr(\hat{\rho}\hat{J})$
is the expectation value of the current operator defined in Eq. (7). Using the
conservation law $\partial_{t}Q(x,t)=a^{-1}(J(x+a/2)-J(x-a/2))$ together with
Eq. (11), we find

$\displaystyle
Q_{pump,b}(t)=Q_{pump}(t)+\sum_{x\in\\{0,a\\}}\Phi(x,x_{b})(Q(x,t)-Q(x,0)),$
(14)

which rearranges charge in the unit-cell in addition to pumping charge across
the cell. Moreover, $\Phi(x,x_{b})=\Theta(x-x_{b})-x/a$ is determined from
$a^{-1}\int dx_{b}Q_{pump,b}=Q_{pump}(t)$. The latter contribution vanishes in
the adiabatic case since the density matrix retains the ground state form.
Note that $Q_{pump,b}(T)=Q_{pump}(T)$.

##  IV Finite-time, nearest-neighbor Thouless pumps

Equations (4) and (12) imply that a time-dependent, Bloch Hamiltonian
$H(k,t)=H_{0}(k,t)$ in the quasi-adiabatic regime gives rise to a Thouless
pump if and only if $H(k,t)=H_{0}(k,t)+H_{CD}(k,t)$ gives rise to a Thouless
pump for arbitrary fast driving. Hence, for the purpose of demonstrating a
“fast” Thouless pump, we can choose a form for $H(k,t)$ that is experimentally
realizable, and reverse-engineer a corresponding $H_{0}(k,t)$, so long as
$H_{0}$ is periodic in $t$ and gapped for all $k$ and $t$ ($H_{0}$ is
automatically periodic in $k$ as a Bloch Hamiltonian). For the purposes of
this work, “an experimentally realizable” protocol refers to a RM model with
complex hopping terms i.e. can be described by Eq. (4) with

$\displaystyle\bm{u}(k,t)=Re[\left\\{-u_{1}(t)-u_{2}(t)e^{ika},-iu_{1}(t)-iu_{2}(t)e^{i(ka)},u_{0}(t)\right\\}]$
(15)

where the hopping parameters $J_{1,2}(t)$ in Fig. 1 are replaced by
$u_{1,2}(t)$, which are now allowed to be complex. The total Hamiltonian $H$
generated using Eq. 4 starting with $H_{0}$ given by an RM model is typically
not a nearest neighbor model (i.e. does not satisfy the constraint above) and
therefore experimentally challenging to implement. Therefore, we choose
$\bm{u}(k,t)$ to satisfy the above constraint and solve for $\bm{R}$ using Eq.
4 to determine $H_{0}(k,t)$, which is now not necessarily a RM model.
Furthermore, the system must start in the ground state of the Hamiltonian
$H(k,t<0)$ and return to the ground state at the end of the protocol, which
leads to the constraint

$\bm{R}(k,0)=\bm{u}(k,t<0)=\bm{u}(k,t>T)=\bm{R}(k,T).$ (16)

This condition is critical to ensuring that $\bm{R}$ is periodic in $k$ and
$t$ during the drive, which in turn is required for the periodicity condition
that leads to quantization of the pumped charge. While the first half of this
condition is simply set as an initial condition for $\bm{R}$ in Eq. 4 viewed
as a time-evolution equation for $\bm{R}(k,t)$, the second half constrains the
final value $\bm{R}(k,t)$ of the solution. As discussed below this constraint
is satisfied by choosing $\bm{u}(k,t)$ according to Eq. (15).

To determine the appropriate $\bm{R}(k,t)$, we start by noting that the
equation for $H_{0}(k,t)$ (i.e. Eq. (4)) is explicitly written as a
differential equation
$\bm{u}=\left(\bm{R}+\frac{1}{R^{2}}\bm{R}\crossproduct\partial_{t}\bm{R}\right)$.
This equation can be simplified into a differential equation

$\displaystyle\partial_{t}\hat{\bm{R}}=\bm{u}\crossproduct\hat{\bm{R}}$ (17)

with the magnitude of $\bm{R}$ determined by

$R=\bm{u}\cdot\hat{\bm{R}}.$ (18)

The above problem can be solved numerically to find solutions near a RM model
with particular parameters. However, before that we discuss an analytically
soluble limit that is based on Thouless pumping using the so-called ”control
freak” protocol [12]. This protocol was shown to allow quantized pumping at
arbitrary speeds but is somewhat limited because it requires one of
$J_{1,2}(t)$ in Fig. 1 to vanish at any instant of time, so that the lattice
is never completely connected. At the same time, this limit provides a nice
analytic limit to study pumping [12] and therefore we use this as an analytic
test example. For such a limit, we use an ansatz for the vector

$\displaystyle\bm{R}(k,0<t<T/2)=(R_{\perp}(t),0,R_{||}(t))$
$\displaystyle\bm{R}(k,T/2<t<T)=(R_{\perp}(t)\cos{k},R_{\perp}(t)\sin{k},R_{||}(t)).$
(19)

This ansatz corresponding to the so-called “control freak” protocol [12] where
the first phase $0<t<T/2$ corresponds to a state where electrons are
delocalized on one set of bonds and the second phase corresponds to electrons
on the other set $T/2<t<T$. The latter set of bonds are entirely inside the
unit cell (i.e. double bonds in Fig. 1) and therefore lead to a Hamiltonian
with no $k-$dependence. The total Hamiltonian to implement this protocol is
given by a vector $\bm{u}(k,0<t<T/2)=\bm{R}(k,t)+\dot{\theta}(0,1,0)$ and
$\bm{u}(k,T/2<t<T)=\bm{R}(k,t)+\dot{\theta}(-\sin{k},\cos{k},0)$, which
clearly satsifies the constraint Eq. 15. Here we have reparameterized the
time-dependence of $R_{\perp,||}(t)$ as
$R_{||}(t)=e^{\lambda(t)}\cos{\theta(t)}$ and
$R_{\perp}(t)=e^{\lambda(t)}\sin{\theta(t)}$. Continuity and time-periodicity
of the vector $\bm{R}(k,t)$ require $\theta(t=0)=0$, $\theta(T/2)=\pi$ and
$\theta(t=T)=2\pi$. Substituting $\bm{R}$ from Eq. IV into the pumped charge
Eq. 12 we get $Q_{pump}(T/2<t<T)=\cos{\theta(t)}$ (and $Q_{pump}(0<t<T/2)=0$).
Substituting into Eq. 13, $Q(0,0,t)-Q(0,0,0)=(1-\cos{\theta(t)})/2$. These can
be combined, using Eq. 14, into equations for charge pumped across each bond:

$\displaystyle
Q_{pump,d}(0<t<T/2)=0,\quad\,Q_{pump,d}(0<t<T/2)=\frac{(1+\cos{\theta(t)})}{2}$
$\displaystyle
Q_{pump,s}(0<t<T/2)=\frac{(1-\cos{\theta(t)})}{2},\quad\,Q_{pump,s}(0<t<T/2)=1,$
(20)

where $\theta(t)$ varies from $0$ to $2\pi$ over the protocol. Based on these
boundary conditions for $\theta$ it is clear that both $Q_{pump,s}$ and
$Q_{pump,d}$ increase from $0$ to $1$ over a time-period. Furthermore, a
single electron is pumped across the $d-$bond in Fig. 1 from site A to B over
the first half of the cycle followed by the electron being pumped across the
$s-$bond from site B to A in the second half of the protocol.

In the case of a generic starting RM model, we numerically solve these
differential equations for $\bm{R}(k,t)$ with the initial condition
$\bm{R}(k,0)=\bm{u}(k,t<0)$ (i.e. Eq. 16). For these results we choose
$J_{1/2}=J_{0}\pm\delta_{0}\cos(\varphi(t)+\phi_{shift})$ and
$\Delta(t)=\Delta_{0}\sin(\varphi(t)+\phi_{shift})$, with lattice spacing $a$
and $\omega=2\pi/T$, for some period $T$, so that $H_{0}(t)$ comes back to
itself after a period $T$. Here $\phi_{shift}$ is a time shift to allow one to
choose general initial starting points. The function
$\varphi(t)=\pi(1-\cos{(\omega t/2)})$ ensures that the derivative of
$\partial_{t}\rm{R}(k,t)|_{t=0,T}=0$. We then vary $\bm{u}(k,t)$ by choosing
$u_{j}(t)=u_{j}^{(0)}(t)+\sum_{n\geq 0}\left(u_{j,n}e^{in\varphi(t)}\right)$
where the coefficients $u_{j,n\neq 0}$ are chosen by numerical optimization
(fminunc in MATLAB) to minimize the least squares error in the final condition
i.e. $\mathcal{E}=\sum_{k}|\bm{R}(k,T)-\bm{R}(k,0)|^{2}$ where the sum is
taken over a discrete grid of $k$. Here the functions $u_{j}^{(0)}(t)$ set the
starting point of the iteration to the Rice-Mele vector vector
$\bm{u}(k,t)=\bm{R}(k,t)$ in Eq. 3.

The components of the optimized solution $\hat{\bm{R}}(k,t)$ for the protocol
are plotted in Figure 2 along with the components of $\bm{u}$, the input to
the optimization method, for fixed $\omega=0.5$ and $k=1.1$. The parameters
for the simulation are chosen to be $a=1$, $J_{0}=1.1$, $\delta_{0}=0.9$,
$\Delta_{0}=1$, with a uniform mesh of one hundred $t$ points between $(0,T)$.
Notice the approximate periodicity of $\hat{R}_{i}$ and their being non-zero.
We see from Figure that $\min_{k}\left(R(k,t)\right)>0$ for all $t$, thus
ensuring that the reverse-engineered $H_{0}(k,t)$ will be gapped for all $k$
and $t$. The error metric $\mathcal{E}$ sums over a uniform mesh of one
hundred $k$ points between $(-\pi/a,\pi/a)$.

(a)

(b)

Figure 2: (a) The components of $\bm{u}(k,t)$ from Eq. (15) which are inputs
to the numerical optimization method plotted with the components of
$\hat{R}(k,t)$ which are the output of the method and solutions to (17)
subject the boundary conditions described in the main text. Plots produced for
$\omega=0.5$ and $k=1.1$. Notice the approximate periodicity of $\hat{R}_{i}$
and their being non-zero. (b) $\min_{k}(R(k,t))$ calculated from the output of
the optimization method as a function of $t$. Here again, $\omega=0.5$. We see
that the magnitude of $R(k,t)$ is greater than $0$ for all $t$ and $k$ thus
ensuring that the reverse-engineered $H_{0}(k,t)=\bm{R}(k,t)\cdot\bm{\sigma}$
will be gapped for all $k$ and $t$.

From the plot of $Q_{pump}$ and $Q_{pump,b}$ of Eq. (12) and Eq. (14) in
Figure 3 a simple physical picture emerges: the particle at site $A$ is
sequentially pumped to site $B$ in the same unit cell and then pumped to site
$A$ in the next unit cell. This follows from the pumped charge initially being
larger across the bond linking $A$ to $B$ within the same unit cell to then
flipping to be larger on the bond linking $B$ to $A$ on the next unit cell.
This is qualitatively similar with the result for the ideal fast pumping model
described by Eq. IV.

Figure 3: $Q_{pump}(t)$ of Eq. (12) and $Q_{pump,b}$ of Eq. (12) for
$\omega=0.5$. Note the quantization of pumped charged over a cycle $T$.
Physically, the particle at site $A$ is sequentially pumped to site $B$ in the
same unit cell and then pumped to site $A$ in the next unit cell.

##  V Discussion

Starting from a Rice-Mele model, we derive an exact $H_{CD}$ protocol and
analytically show that it preserves the charge pumped across each bond in the
model (see Fig. 1) as a Chern number. We do this by determining the pumped
charge $Q_{pump,b}(t)$ for the two types of bond $s$ and $d$ in the lattice
shown in the RM model shown in Fig. 1. However, this exact $H_{CD}$ involves
higher harmonics of $k$ which in real space correspond to long-range hopping
terms. Experimentally, such terms are challenging if not impossible to
actualize. The total Hamiltonian $H=H_{0}+H_{CD}$ is what is actually
implemented in an experiment. Based on this we considered the possibility of
restricting $H(t)$ to only include nearest neighbor hopping terms through Eq.
15.

With this constraint we found that the “control freak” protocol on the RM
model [12] allowed for an analytically soluble example of fast pumping in an
RM model as was clear from Eq. IV. We then relaxed the requirements by
numerically finding a more general Hamiltonian $\bm{u}(t)$, which led to
quantized charge pumping as was seen in Fig. 3. In both these case, the time-
dependence of the charge pumped across bonds provided a physical
interpretation of the charge being first pumped from site $A$ to $B$ and then
back to site $A$ in the next unit-cell. These results on fast quantized
pumping over a class of Hamiltonians provide a wealth of opportunities for
experimental realizations of finite-time Thouless pumps.

The Rice-Mele model we considered is a lattice model, that in $k$ space is a
two-band system, without disorder or interactions. Applying STA methods to
multiband lattice systems and systems involving interactions and disorder are
important future steps. Our work contributes to the promise of actualizing
finite-time, robust and dissipationless charge transport.

Note added: In the process of completing this manuscript, we became aware of
the preprint by Liu et al [17], who also apply counterdiabatic methods to the
Rice-Mele model. Their work focuses on Wannier centers, and they obtain a
counterdiabatic field $\bm{u}$ that generally includes non-nearest-neighbor
interactions.

Acknowledgements: JC is supported by the National Science Foundation Graduate
Research Fellowship Program under Grant No. DGE 2236417\. CJ acknowledges
support by the U.S. National Science Foundation under Grant No. 2127900. Any
opinions, findings, and conclusions or recommendations expressed in this
material are those of the authors and do not necessarily reflect the views of
the National Science Foundation. J.S. acknowledges support from the Joint
Quantum Institute and support by the Laboratory for Physical Sciences through
its continuous support of the Condensed Matter Theory Center at the University
of Maryland. JC gratefully acknowledges Michael Hinczewski and Sebastian
Deffner for helpful conversations.

##  Appendix A Current density operator

Starting from (7), we consider the right hand side of the continuity equation.
For this purpose we define a discrete derivative
$\delta_{x}f(x)=f(x+a/2)-f(x-a/2)$, where the function $f$ is defined over
bond when $x$ is on a site. The discrete analog of the continuity equation is
then given by:

$\begin{split}-\delta_{x}\langle\hat{J}\rangle&=-\delta_{x}\left(\sum_{x_{1},x^{\prime}}\Psi^{*}(x_{1})\bra{x_{1}}\hat{J}(x)\ket{x^{\prime}}\Psi(x^{\prime})\right)\\\
&=-i\sum_{x_{1},x^{\prime}}\Psi^{*}(x_{1})H(x_{1},x^{\prime})\Psi(x^{\prime})\Bigl{(}\delta_{x,x_{1}}-\delta_{x-x^{\prime}}\Bigr{)}\\\
&=-i\left(\sum_{x^{\prime}}\Psi^{*}(x)H(x,x^{\prime})\Psi(x^{\prime})-\sum_{x_{1}}\Psi^{*}(x_{1})H(x_{1},x)\Psi(x)\right)\\\
&=-i\left(\sum_{x_{1}}\Psi^{*}(x)H(x,x)\Psi(x)-\sum_{x_{1}}(\Psi(x)(H(x,x_{1}))^{*}\Psi(x_{1})\right)\\\
&=\Psi^{*}(x)\partial_{t}\Psi(x)+\Psi(x)\partial_{t}\Psi^{*}(x)\\\
&=\partial_{t}\rho(x)\end{split}$ (21)

Thus, the current density operator defined in (7) satisfies a continuity
equation.

## References

  * Thouless [1983] D. Thouless, Physical Review B 27, 6083 (1983).
  * Kitagawa  _et al._ [2010] T. Kitagawa, E. Berg, M. Rudner, and E. Demler, Physical Review B 82, 235114 (2010).
  * Titum  _et al._ [2016] P. Titum, E. Berg, M. S. Rudner, G. Refael, and N. H. Lindner, Physical Review X 6, 021013 (2016).
  * Fedorova  _et al._ [2020] Z. Fedorova, H. Qiu, S. Linden, and J. Kroha, Nature communications 11, 3758 (2020).
  * Guéry-Odelin  _et al._ [2019] D. Guéry-Odelin, A. Ruschhaupt, A. Kiely, E. Torrontegui, S. Martínez-Garaot, and J. G. Muga, Reviews of Modern Physics 91, 045001 (2019).
  * Demirplak and Rice [2003] M. Demirplak and S. A. Rice, The Journal of Physical Chemistry A 107, 9937 (2003).
  * Demirplak and Rice [2005] M. Demirplak and S. A. Rice, The Journal of Physical Chemistry B 109, 6838 (2005).
  * Berry [2009] M. V. Berry, Journal of Physics A: Mathematical and Theoretical 42, 365303 (2009).
  * Funo  _et al._ [2020] K. Funo, N. Lambert, F. Nori, and C. Flindt, Physical Review Letters 124, 150603 (2020).
  * Takahashi  _et al._ [2020] K. Takahashi, K. Fujii, Y. Hino, and H. Hayakawa, Physical Review Letters 124, 150602 (2020).
  * Erdman  _et al._ [2019] P. A. Erdman, F. Taddei, J. T. Peltonen, R. Fazio, and J. P. Pekola, Physical Review B 100, 235428 (2019).
  * Asbóth  _et al._ [2016] J. K. Asbóth, L. Oroszlány, and A. Pályi, Lecture notes in physics 919, 166 (2016).
  * Rice and Mele [1982] M. Rice and E. Mele, Physical Review Letters 49, 1455 (1982).
  * Qi and Zhang [2011] X.-L. Qi and S.-C. Zhang, Reviews of Modern Physics 83, 1057 (2011).
  * Vanderbilt [2018] D. Vanderbilt,  _Berry phases in electronic structure theory: electric polarization, orbital magnetization and topological insulators_ (Cambridge University Press, 2018).
  * Hasan and Kane [2010] M. Z. Hasan and C. L. Kane, Reviews of modern physics 82, 3045 (2010).
  * Liu  _et al._ [2024] W. Liu, Y. Ke, and C. Lee, arXiv preprint arXiv:2401.17081 (2024).