# Robust Feedback Motion Policy Design Using Reinforcement Learning on a 3D
Digit Bipedal Robot

Guillermo A. Castillo1, Bowen Weng1, Wei Zhang2, and Ayonga Hereid3 *This work
was supported in part by the National Science Foundation under grant
CNS-1552838, the OSU M&MS Discovery Theme Initiative, and the startup fund of
SUSTech. 1Electrical and Computer Engineering, Ohio State University,
Columbus, OH, USA; {castillomartinez.2<EMAIL_ADDRESS>Institute
of Robotics, Southern University of Science and Technology (SUSTech), China;
<EMAIL_ADDRESS>and Aerospace Engineering, Ohio State
University, Columbus, OH, USA<EMAIL_ADDRESS>
###### Abstract

In this paper, a hierarchical and robust framework for learning bipedal
locomotion is presented and successfully implemented on the 3D biped robot
Digit built by Agility Robotics. We propose a cascade-structure controller
that combines the learning process with intuitive feedback regulations. This
design allows the framework to realize robust and stable walking with a
reduced-dimension state and action spaces of the policy, significantly
simplifying the design and reducing the sampling efficiency of the learning
method. The inclusion of feedback regulation into the framework improves the
robustness of the learned walking gait and ensures the success of the sim-to-
real transfer of the proposed controller with minimal tuning. We specifically
present a learning pipeline that considers hardware-feasible initial poses of
the robot within the learning process to ensure the initial state of the
learning is replicated as close as possible to the initial state of the robot
in hardware experiments. Finally, we demonstrate the feasibility of our method
by successfully transferring the learned policy in simulation to the Digit
robot hardware, realizing sustained walking gaits under external force
disturbances and challenging terrains not included during the training
process. To the best of our knowledge, this is the first time a learning-based
policy is transferred successfully to the Digit robot in hardware experiments
without using dynamic randomization or curriculum learning.

##  I Introduction

For the bipedal locomotion problem, the policy _robustness_ is a critical
characteristic and remains one of the biggest challenges in the field. In
practice, the robustness of the control policy can be presented as (i) the
capability of handling various external disturbances (e.g., push recovery),
(ii) maintaining stable gaits while operating under various terrain
conditions, and (iii) accomplishing the sim-to-real transfer with as little
effort as possible. In this paper, we present a successful application of
designing a feedback motion policy on Digit, a challenging 3D bipedal robot,
and demonstrate robust performance among all of the aspects mentioned above.

Figure 1: Digit recovering from an external force disturbance applied in the
lateral direction.

In general, a robust bipedal locomotion policy can be obtained through (i) the
model-based approach, (ii) the learning-based approach, or (iii) a combination
of both. In the model-based regime, existing research typically relies on a
simplified template model [1, 2]. Some has also explored robust control
formulations of Control Lyapunov Function Control Barrier Functions (CLF-CBF)
[3] and Hybrid Zero Dynamics (HZD) inspired solutions [4]. In general, the
derived policy from the simplified model requires further tuning of low-level
control heuristics. On the other hand, a typical learning-based solution
involves supervised learning, reinforcement learning, and imitation learning.
The policy is typically constructed in an end-to-end manner [5] directly
working with the robot’ s full-body dynamics. Some existing work has shown
remarkable progress in push-recovery [6] and operating in various terrain
conditions [7]. While model-based methods often rely on a simplified model and
require extensive control gain tuning, most learning-based solutions require a
large amount of data.

A locomotion policy combining the learning-based and the model-based solution
is typically formulated in a cascade structure. The high-level controller
computes a reference trajectory of a selected anchor point (e.g., the center
of mass position). The lower-level controller then seeks to track such a
learned reference through basic model information such as kinematics. Morimoto
et.al. [8, 9] learned the Poincaré map of the periodic walking pattern and
applied the method to two 2D bipedal robots. Some recent work has proposed to
learn the joint-level trajectory for each joint as the reference motion
through supervised learning [10] or using reinforcement learning [11, 12, 13].
This approach simplifies the design of the lower-level tracking, which can be
as simple as a PD controller.

With the aforementioned approaches, the empirical successes are mostly
evaluated in simulations. Some recent work has seek to tackle the sim-to-real
problem through dynamics randomization [14], system identification [15], and
periodic reward composition [16]. However, many of these methods still suffer
from poor sampling efficiency. Moreover, the robustness is often partially
validated emphasizing one of the three mentioned phenomenons. To the best of
our knowledge, it remains a challenge to efficiently obtain a locomotion
policy that handles external disturbances, uneven terrains, with a little-
effort sim-to-real transfer. The challenge further escalates as one considers
a 3D, underactuated, and highly nonlinear bipedal system such as the Digit
(Figure 2).

Motivated by the existing challenges, this paper enhances the hybrid zero
dynamics inspired reinforcement learning framework presented in our previous
work [11, 13] emphasizing the policy robustness and sim-to-real transfer. The
proposed cascade control structure couples the high-level learning-based
controller with the low-level model-based controller. The obtained policy is
able to achieve robust bipedal locomotion on a challenging 3D Digit robot in
experiments to accomplish the following tasks.

  * •

Learn disturbance resistance without explicitly experiencing adversarial
attacks in training.

  * •

Adapt to various terrain conditions with the policy learning process only
limited to flat ground.

  * •

Transfer the learned policy between different simulation environments and from
simulation to the real robot without randomized dynamics in training or
extensive hyper-parameter tuning.

It is worth emphasizing that the policy is learned from scratch and does not
rely on demonstrations. To the best of our knowledge, this is the first
successful sim-to-real transfer of a learning-based policy working with Digit.
We further emphasize the following features that contribute to the
aforementioned robust performance of the proposed policy:

The high-level learning-based policy propagates reference motion trajectories
in real-time at low-frequency compatible with real hardware implementation.
This encourages the stability of the walking gait by modifying the reference
trajectories accordingly to the reduced state of the robot.

The low-level model-based regulation applies compensation to the reference
trajectories at a higher frequency based on instantaneous state feedback. This
feedback-based trajectory regulation compensates for uncertainty in the
environment, and based on our experiment observations, it is one of the key
factors leading to the effortless sim-to-real transfer (as illustrated in Fig.
6).

A model-based balancing controller is implemented that induces a set of
feasible initialization states for operations in the simulators and with the
real robot. While a feasible initialization is typically not of interest in
simulations given the state can be set arbitrarily, it is crucial to the
observed success of implementing the proposed locomotion policy to work with
Digit in real-world experiments.

The remainder of the paper is organized as follows. Section II presents
details of the Digit robot. Section III introduces the proposed framework
including the high-level learning module, the low-level model-based
controllers, the training, and the execution of the overall framework. A
series of experiments are included in Section V emphasizing the locomotion
robustness in terms of disturbance rejection, stable walking gaits on uneven
terrains, and effortless sim-to-real transfer. Section VI concludes the paper
and discusses future work.

##  II Kinematic Model of Digit Robot

In this section, we will briefly introduce the Digit robot and its kinematic
model and notations that will be used in the following sections.

###  II-A Robot Hardware

Digit is a versatile bipedal robot designed and built by Agility Robotics
[17]. The robot has an integrated perception system, 20 actuated joints, 1
IMU, and 30 degrees of freedom (DoF) that allows it to achieve robust and
dynamic locomotion. The robot design is based on its predecessor Cassie, whose
leg’s morphology is derived from an Orchid-like bird. Hence, the location of
the knee and ankle is not immediately obvious from the robot’s appearance.
Fig. 2. shows the Digit robot with a description of the body links and its
kinematic structure. The total weight of the robot is 48 kg. The torso of the
robot weighs 15 kg and contains Digit’s on-board computer, power source, and
vision sensors. Vision sensors include a LiDAR, an RGB camera, three
monochrome depth cameras and a RGB-depth camera.

Figure 2: Digit robot (left) and its kinematic tree structure (right). The
floating base coordinate is located at the pelvis link. Only left leg and left
arm joints are shown. Green represents actuated leg joints, blue represents
passive leg joints, and orange represents arm joints.

###  II-B Kinematic Model and Notations

Each arm of Digit has 4 joints for basic manipulation tasks. Each leg has
eight joints, from which four are directly actuated by electrical motors (hip
yaw, hip roll, hip pitch, knee), two are indirectly actuated by electrical
motors (toeA, toeB) and correspond to the actuated joints of the foot (foot
pitch, foot roll). The addition of the foot pitch and foot roll joints
improves the ability to balance stably on a wide variety of terrains. The
remaining two joints (ankle and shin) are passive and they are connected via
specially designed leaf-spring four-bar linkages for additional compliance.

As shown in Fig. 2, the variables $q_{x},q_{y},q_{z}$ denote the Cartesian
position of the robot’s pelvis, while the torso orientation is represented by
the x-y-z Euler angles (roll, pitch, yaw) as $q_{\psi},q_{\theta},q_{\phi}$.
Then, the coordinate of the floating base coordinate system at the pelvis is
denoted by

$\displaystyle\mathbf{q^{b}}=[q_{x},q_{y},q_{z},q_{\psi},q_{\theta},q_{\phi}]^{T}.$
(1)

The actuated joints are denoted by the vector given as,

$\displaystyle\mathbf{q^{a}}=[q_{1L},q_{2L},q_{3L},q_{4L},q_{7L},q_{8L},q_{9L},q_{10L},q_{11L},q_{12L},$
$\displaystyle
q_{1L},q_{2L},q_{3L},q_{4L},q_{7L},q_{8L},q_{9L},q_{10L},q_{11L},q_{12L}]^{T},$
(2)

where the joints $q_{1}$-$q_{4}$ correspond to the leg’s actuated joints (i.e.
hip roll, hip yaw, hip pitch, knee), joints $q_{7},q_{8}$ correspond to the
foot pitch and foot roll, and joints $q_{9}$-$q_{12}$ correspond to the arm’s
actuated joints (i.e. shoulder roll, shoulder pitch, shoulder yaw, elbow). The
passive joints are denoted by the vector:

$\displaystyle\mathbf{q^{p}}=[q_{5L},q_{6L},q_{5R},q_{6R}]^{T},$ (3)

where $q_{5},q_{6}$ correspond to the shin and tarsus joints. Hence, the
generalized coordinates of Digit robot are denoted by:

$\displaystyle\mathbf{q}=[\mathbf{q^{b}},\mathbf{q^{a}},\mathbf{q^{p}}]^{T}.$
(4)

In addition, the robot has 4 end-effectors, which are the left/right feet and
fists. We define the Cartesian position of any of these end-effectors as:

$\displaystyle\mathbf{x}_{ee}\mathbf{(q)}=[x_{ee},y_{ee},z_{ee}]^{T}$ (5)

The CoM Cartesian position is denoted by:

$\displaystyle\mathbf{x}_{CoM}\mathbf{(q)}=[x_{CoM},y_{CoM},z_{CoM}]^{T}$ (6)

The current position of the robot’s feet and CoM can be determined by applying
Forward Kinematics (FK). In particular, we use the the open-source package
FROST [18] to obtain symbolic expressions for
$\mathbf{x}_{footL}\mathbf{(q)}$, $\mathbf{x}_{footR}\mathbf{(q)}$, and
$\mathbf{x}_{CoM}\mathbf{(q)}$. We created a URDF model of Digit based on the
XML model provided by Agility Robotics.

##  III Learning Approach

In this section, we build upon our previous work proposed in [11, 13] to
implement a cascade-structure learning framework that realizes stable and
robust walking gaits for the 3D bipedal robots. The specific design ensures
successful transferring learned policies in simulation to robot hardware with
minimal turning.

Figure 3: Overall structure of the proposed Trajectory-based RL framework. The
trajectory planning phase is done by the neural network policy while the
feedback regulation block uses the robot’s sensor information to improve the
stability of the walking gait and the velocity tracking performance.

###  III-A Overview of the Learning Framework

As shown in Fig. 3, the proposed framework presents a cascade decoupled
structure to tackle the 3D walking problem through its two main components,
neural-network trajectory planning and feedback regulation.

The neural network trajectory planner uses a reduced-dimensional
representation of the robot states to compute a set of coefficients
$\mathbf{\alpha}$ that parameterize joint reference trajectories using
$5^{th}$-order Bézier Polynomials. The structure of the neural network is
provided in Section III-C. Since the forward and backward propagation of the
neural network could take significant time when compared with the running
frequency of the low-level controller (1 kHz), the trajectory planning phase
runs at a lower frequency (250 Hz). This consideration is especially relevant
for the implementation of the learned policy in the real robot. More
importantly, it opens the door to the real-time implementation of new learning
structures that may be computationally expensive in both simulation and
hardware.

The feedback regulation uses the robot’s pelvis velocities and torso
orientation to modify the joint reference trajectories. It compensates for the
uncertainty in the robot’s model and the environment, improving the robustness
of the walking gait. The joint-level PD controllers are used to track the
compensated reference joint trajectories. The compensations applied during the
feedback regulation are discussed in detail in Section IV-B. The low-level
feedback regulation is a key component of our control structure to close the
sim-to-real gap, enabling successful transfer of the learned policy to
hardware without the need for exhaustive tuning, dynamic randomization, or
curriculum learning.

###  III-B Dimension Reduction of the State and Action Space

By applying physical insights from the walking motion, we can significantly
reduce the number of inputs and outputs of the NN. Unlike other end-to-end RL
frameworks, we do not use all the available states of the robot to feed the
NN. Instead, we select states that provide insightful information such as the
torso orientation ($q_{\psi},q_{\theta},q_{\phi}$), pelvis linear velocity
($\dot{q}_{x},\dot{q}_{y},\dot{q}_{z}$), torso angular velocity
($\dot{q}_{\psi},\dot{q}_{\theta},\dot{q}_{\phi}$), and desired walking
direction, which is encoded by the desired forward and lateral walking
velocity ($\dot{q}_{x}^{d},\dot{q}_{y}^{d}$).

To reduce the dimension of the action space, we will keep the arms joints
fixed during the walking motion and keep the stance foot passive and the swing
foot parallel to the ground during the walking gait. Without the need for
determining reference trajectories for arm and foot joints, we reduce the
number of reference trajectories needed to be computed by the NN to eight,
with each leg having three hip joints and one knee joint. In addition, we
impose a symmetry condition between the right and left stance during the
walking gait. Therefore, given the set of coefficients for the right stance
$\mathbf{\alpha}^{R}\in\mathbb{R}^{M+1\times N}$, where $M=5$ is the degree of
the Bézier polynomials and $N=8$ is the number of reference trajectories, we
can obtain the set of coefficients for the left stance
$\mathbf{\alpha}_{L}\in\mathbb{R}^{M+1\times N}$ by the symmetry condition:

$\displaystyle\mathbf{\alpha}^{L}=\mathbf{T}\mathbf{\alpha}^{R},$ (7)

where $\mathbf{T}\in\mathbb{R}^{N\times N}$ is a sparse transformation matrix
that represents the symmetry between the joints of the right and left legs of
the robot.

To further reduce the action space and encourage the smoothness of the control
actions after the ground impact, we enforce that at the beginning of every
step, the initial point of the Bézier polynomial coincides with the current
position of the robot’s joints. That is, for each joint $i$ with Bézier
coefficients $\mathbf{\alpha}^{R}_{i}\in\mathbb{R}^{M+1}$, we have

$\displaystyle\alpha_{i}[0]=q_{i}(\tau(0)),\quad\tau\in[0,1],$ (8)

where $\tau$ is the time-based phase variable used to parameterize the Bézier
Polynomials. Finally, we enforce the position of the hip joints and knee
joints to be the same at the end of the current step and the beginning of the
next step ($\tau(t)=0$ and $\tau(t)=1$ respectively). This enforces continuity
in the joint position trajectories after switching the stance foot.

###  III-C Neural Network Structure

Given the considerations presented in Section III-B, the number of inputs for
the NN is $10$ and the number of outputs is $32$. We choose the number of
hidden layers of the NN to be $4$, each one with $32$ neurons. The activation
function for the hidden layers is ReLU, and the activation function for the
output layer is sigmoid. Finally, we scale the output of the NN within a range
of admissible motion for the robot’s joints. In particular, we use the convex
hull property of Bézier polynomials to translate the joint limits to the
corresponding coefficients limits. Given the number of inputs and outputs of
the NN, the number of trainable parameters is 4576. To the best of our
knowledge, this is the smallest NN implemented in hardware to realize 3D
walking locomotion.

###  III-D Walking Policy Learning Pipeline

Since the main purpose of this work is to implement a learning framework that
realizes a policy that is transferable to the real hardware, we need to create
a pipeline for the training process in simulation that renders working
conditions as close as possible to the real hardware. A diagram of the
training pipeline is shown if Fig. 4.

Figure 4: Pipeline of the proposed learning framework. Standing controller is
included within the pipeline of the learning process to obtain a feasible set
of initial states for the policy training using a customized Mujoco
environment. Then, the trained policy is tested in a more realistic real-time
simulator, and finally it is tested on hardware.

The initial state for each training episode is chosen randomly from a pool of
initial states whose kinematics and dynamics are feasible to be implemented on
the real robot. To achieve this, we implement a standing controller that
allows the robot to start up from an arbitrary resting position when the robot
is hanging up from the lab’s crane. We tested the controller successfully in
hardware and use it to replicate the standing process in simulation. More
details about the standing controller are provided in Section IV-A.

After the standing controller is activated and the robot is able to stand and
balance steadily in simulation, we capture the state of the robot and save it
in a pool of initial states. Then, we repeated the process 40 times to create
a diverse enough set of initial states that we can use for training our
controller. It is important to denote that this process does not provide a
fixed initial state for the walking gait, but a whole set of different initial
configurations that encourage the randomness of the initial state during the
training while providing feasible and safe starting conditions for the walking
gait. Moreover, this process does not provide any reference trajectories for
the walking gait that could bias the learning process to specific or
predefined walking gaits. That is, different from other RL approaches to
bipedal locomotion, our method learns walking gaits from scratch without the
need for predefined reference trajectories or imitation learning.

Given the pool of initial states, we use a customized environment using MuJoCo
[19] to start the training process using the Evolution Strategies (ES)
algorithm. When the training is finished, we proceed to test the trained
policy in the Agility Robotics proprietary simulation software. This
simulation software provides a more realistic representation of the robot’s
dynamic behavior since it runs in real-time and shares the same features with
the hardware such as low-level API, communication delay, and model dynamic
parameters.

Finally, once we verify the trained policy is deployed safely in the AR
simulation, we proceed to test it in the hardware and verify its performance
to do additional tuning of the low-level controller gains if needed. A
sequence of the whole learning process and transference to real hardware can
be seen in the accompanying video submission.

###  III-E Reward Function Design

Following the work in [13], we use a reward function of the form:

$\mathbf{r}=\mathbf{w}^{T}[r_{v_{x}},r_{v_{y}},r_{h},r_{a},r_{CoM},r_{ang},r_{angvel},r_{u},r_{fd}]^{T}.$
(9)

where $\mathbf{w}$ is a vector of weights corresponding to each of customized
rewards $r_{i}$. This reward function encourages forward and lateral velocity
tracking (through $r_{v_{x}},r_{v_{y}}$), height maintenance (with $r_{h}$),
energy efficiency (with $r_{u}$) and natural walking gaits. The episode length
is 10000 simulation steps, which are equivalent to 5 seconds, and it has an
early termination if any of the following conditions is violated:

$\displaystyle|q_{\psi}|<0.5,\quad|q_{\theta}|<0.5,\quad|q_{\phi}|<0.5,$ (10)
$\displaystyle|\dot{q}_{\psi}|<2,\quad|\dot{q}_{\theta}|<2,\quad|\dot{q}_{\phi}|<2,$
$\displaystyle 0.8<q_{z}<1.2,\quad\Delta_{f}<0.05,$

where $q_{z}$ is the height of the robot’s pelvis and $\Delta_{f}$ is the
distance between the feet.

##  IV Control

In this section, we provide details about the implementation of the high-level
and low-level controllers of the control-learning framework introduced in
Section III. Moreover, we provide details of the architecture of our
controller and its integration with Digit’s low-level API and communication
system. This architecture is presented in Fig. 5.

Figure 5: Our controller architecture allows asynchronous operation of the
high-level and low-level controllers. Moreover, it enables fast and reliable
communication between different control layers, external user inputs, and the
Digit simulator/hardware. Note that the same controller architecture can be
used for different tasks e.g. walking, standing, crouching.

The Agility Robotics low-level API streams data using the UDP protocol and
enables the user to access the low-level sensor data and to give commands
directly to the motor drives. To keep the connection of the low-level API
active to receive sensor data, commands must be sent periodically. Once the
client is connected and sending commands, torque control must be activated by
requesting transition to Low-level API operation using JSON messages sent
through Websocket protocol.

In addition, we use the Robot Operating System (ROS) to manage the
communication between different components in our system. This also adds
significant flexibility to our controller structure to include nodes for
additional tasks such as logging information, capturing external command
inputs, and accessing Digit’s perception system.

Since the low-level API needs to run at a very high frequency (2 kHz), we use
shared memory to enable fast communication between the low-level API and our
custom main control code (1 kHz). The low-level API reads the sensor
measurements and writes them into the shared memory. It also reads from the
shared memory the torque information and velocity commands written by the main
control code and sends the commands to the motors.

On the other hand, the main control code manages the integration and
synchronization of all the components in the overall control structure. It
reads the sensor information from the shared memory and publishes it in ROS
topics to make it available for the other components of the systems such as
the high level planner. In this work, the high-level planner is the trained NN
policy, which takes the information from the corresponding ROS topics to
compute the output (coefficients of the Bézier Polynomials) at a frequency of
250Hz. Then, the high-level planner publishes this information in its
respective ROS topics to make it available for the main control code.

The main control code reads the Bézier coefficients published by the high-
level planner and uses them to compute the corresponding Bézier Polynomials
that become the reference trajectories for the robot’s joints. In addition,
the main control code uses the sensor feedback information shared by the Low-
level API to compute the regulations needed to compensate for the model
mismatch between the simulation and hardware. The integration of these Low-
level feedback regulations into the learning framework is a key part that
makes the controller structure very robust to uncertainties in the model and
makes possible an almost zero effort transfer from the policy learned in
simulation into the real hardware. The detailed structure of the regulations
used in the low-level controller is discussed in Section IV-B.

###  IV-A Standing Controller

The standing controller implemented for Digit is based on the standing
controller implemented for Cassie in [20]. To ensure the robot keeps the
balance, the standing controller uses simple and intuitive regulations based
on the torso orientation and the position of the CoM with respect to the feet.
That is, keeping the robot CoM within the support polygon while both feet are
flat on the ground. The CoM position in the y-axis is adjusted by varying the
legs’ length, and the CoM in the x-axis is regulated by controlling the pitch
angle of the feet.

In our learning framework, the standing controller is used to bring the robot
to a stable configuration after initialization in arbitrary standing
positions. This process is repeated several times in simulation until we
obtain a diverse enough pool of stable configurations that are realistic and
feasible to be implemented in hardware. This set of feasible initial
configurations is then used to randomly choose the initial states during the
training process.

###  IV-B Feedback Regulations

The feedback regulation module in our controller structure is a key component
of the framework as it allows the controller to compensate the trajectories
obtained from the high level planner to adapt it to unknown disturbances such
as the mismatch between the simulation model and real robot, and environmental
factors like external disturbances or challenging terrains that the policy has
not experimented in simulation. These regulations are separated in three main
groups: foot placement, torso, and foot regulations.

Foot placement regulation controller has been widely used in 3D bipedal
walking robots with the objective of improving the speed tracking and the
stability and robustness of the walking gait [21, 22, 20]. Longitudinal speed
regulation, defined by (13), sets a target offset in the swing hip pitch
joint, whereas lateral speed regulation (11) do the same for the swing hip
roll angle. Direction regulation (12) add an offset to the yaw hip angle to
keep the torso yaw orientation at the desired angle.

$\displaystyle\delta_{q_{1i}}$
$\displaystyle=S_{y}(\tau(K_{p_{y}}(\dot{q}_{y}-\dot{q}_{y}^{d})+K_{d_{y}}(\dot{q}_{y}-\dot{q}_{y}^{ls}))+\beta_{y}),$
(11) $\displaystyle\delta_{q_{2i}}$
$\displaystyle=\tau(q_{\phi}-q_{\phi}^{d}),$ (12)
$\displaystyle\delta_{q_{3i}}$
$\displaystyle=-(\tau(K_{p_{x}}(\dot{q}_{x}-\dot{q}_{x}^{d})+K_{d_{x}}(\dot{q}_{x}-\dot{q}_{x}^{ls}))+\beta_{x}),$
(13)

where $i\in\\{L,R\\}$ depends on which foot is the swing foot, $S_{y}=1$ if
$i=L,S_{y}=-1$ if $i=R$, $\dot{q}_{x}$, $\dot{q}_{y}$ are the longitudinal and
lateral speeds of the robot, $\dot{q}_{x}^{ls}$, $\dot{q}_{y}^{ls}$ are the
speeds at the end of the previous step, $\dot{q}_{x}^{d}$, $\dot{q}_{y}^{d}$
are the reference speeds, and $K_{p_{x}},K_{d_{x}},K_{p_{y}},K_{d_{y}}$ are
the proportional and derivative gains.

(a) Limit walking cycle of trained policy in simulation.

(b) Limit walking cycle of trained policy in real robot.

Figure 6: Comparing limit walking cycles between the simulation trails and
real-robot tests: (i) For the same joint, simulated motion and real robot
motion share similar torque limits. (ii) For all of the hip joints, the
regulation term modifies the actual motion leading to the slight difference of
the limit walking cycles. (iii) As the regulation is not involved in any of
the knee joints, the limit walking cycles are almost identical.

The phase variable $\tau$ is used to smooth the regulation at the beginning of
each walking step and reduce torque overshoots. The term $\beta$ is the output
of an additional PI controller used to compensate for the accumulated error in
the velocity and prevent the robot to drift towards a non-desired direction.
Based in our experiments, the inclusion of $\beta$ is key in the successful
sim-to-real transfer of our controller.

Torso regulation is applied to keep the torso in an upright position, which is
desired for a stable walking gait. Assuming that the robot has a rigid body
torso, simple PD controllers defined by (14) and (15) can be applied
respectively to the hip roll $q_{1i}$ and hip pitch $q_{3i}$ angle of the
stance leg:

$\displaystyle u_{q_{1i}}$
$\displaystyle=-(K_{p_{troll}}(q_{\psi}-q_{\psi}^{d})+K_{d_{troll}}(\dot{q}_{\psi}-\dot{q}_{\psi}^{d})),$
(14) $\displaystyle u_{q_{3i}}$
$\displaystyle=S_{\theta}(K_{p_{tpitch}}(q_{\theta}-q_{\theta}^{d})+K_{d_{tpitch}}(\dot{q}_{\theta}-\dot{q}_{\theta}^{d})),$
(15)

where $i\in\\{L,R\\}$ depends on which foot is the stance foot, $S_{\theta}=1$
if $i=L,S_{\theta}=-1$ if $i=R$ where $q_{\phi}$ and $q_{\theta}$ are the
torso roll and pitch angles, and
$K_{p_{troll}},K_{d_{troll}},K_{p_{tpitch}},K_{d_{tpitch}}$ are manually tuned
gains.

Foot orientation regulation is applied to keep the swing foot flat during the
swinging phase to ensure a proper landing of the foot on the ground. Since
Digit has 2 DoF actuated foot, one regulation is needed for each of the roll
and pitch angles. By using forward kinematics, the swing foot roll and pitch
regulation are determined by

$\displaystyle q_{7i}^{d}$ $\displaystyle=q_{\psi}+q_{1i}+S_{f}(21\deg)$ (16)
$\displaystyle q_{8i}^{d}$ $\displaystyle=q_{\theta}+q_{3i}+S_{f}(6\deg),$
(17)

where $i\in\\{L,R\\}$ depends on which foot is the swing foot, $S_{f}=1$ if
$i=L,S_{f}=-1$ if $i=R$, and $q_{7i}^{d},q_{8i}^{d}$ are the desired pitch and
roll foot angle. Moreover, the stance foot is kept passive during the stance
phase, which helps to the stability of the walking gait, especially for soft
or irregular surfaces.

##  V Experimental Results

In this section, we validate the learning framework and controller structure
presented in this paper through both simulation and real hardware experiments
on the Digit robot. We prove that our controller structure allows the
transference of walking policies learned in simulation to the real hardware
with minimal tuning and enhanced robustness. To the best of our knowledge this
is the first time that a policy learned in simulation is successfully
transferred to hardware to realize stable and robust dynamic walking gaits for
Digit robot. In addition, we show that the improved controller structure
presented in this work is robust enough to mitigate the uncertainty in the
robot’s dynamics caused by the mismatch between the simulation and real
hardware. Finally, by thorough experiments on hardware, we show the robustness
of the learned controller against adversarial disturbances applied to the
robot as well as its capability to adapt to various terrains without need of
training for such challenging scenarios.

###  V-A Sim-to-real Transfer and Stability of the Walking Gait

We transferred the learned policy trained in simulation into the real robot
and we tested empirically the stability of the walking gait by analysing the
walking limit cycle described by the robot’s joint during the walking motion.
Fig. 6(a) shows the limit cycle described by some of the robot’s joints while
the robot is walking in place during about 1 minute in simulation, while Fig.
6(b) shows the limit cycles from the real experiment. The convergence of the
walking limit cycle to a stable periodic orbit demonstrates that the walking
gait is stable and symmetric, meaning that the policy is transferred
successfully from simulation to hardware.

###  V-B Robustness to Disturbances

To evaluate the robustness of the policy against disturbances, adversarial
forces are applied to the robot in different directions. Fig. 1 and Fig. 7
show respectively the transition of the robot recovering from a push in the
lateral direction, and the velocity profile of the robot during the
disturbance. In addition, Fig. 8 shows the limit cycles of the robot’s joints
when the disturbance is applied, where it can be seen the robot is able to
recover effectively from the push as the joint limit cycles return to a stable
periodic orbit after the disturbance.

Figure 7: Digit recovering from an external disturbance Figure 8: Disturbance
rejection of learned policy. The walking limit cycles described before, during
and after the lateral push show that the controller can handle effectively
large external disturbance forces.

###  V-C Robustness on Different Terrains

To further evaluate the robustness of the learned control policy, we make
Digit walk in various terrains with different levels of difficulty. These
terrains include flat vinyl ground, mulch, flat rubber ground, and irregular
rubber terrain. The controller is able to adapt to any of these terrains while
keeping an stable and robust walking gait. Snapshots of the walking gait over
the different terrains are shown in Fig. 9, and the complete motion can be
found in the accompanying video [23].

Figure 9: Digit walking on different terrains. From right to left we have
vinyl, mulch, flat rubber, and irregular rubber terrain. The stance foot is
kept passive during the stance phase to allow the foot to adapt easily to
different terrains.

##  VI Conclusions

In this paper, we present a framework for learning robust bipedal locomotion
policies that can be transferred to real hardware with minimal tuning. By
combining a sample efficient learning structure with intuitive but powerful
feedback regulations in a cascade structure, we decouple the learning problem
into 2 stages that work at a different frequency to facilitate the
implementation of the controller in the real hardware. While the trajectory
planning stage is handled by the neural network to produce reference
trajectories for the actuated joints of the robot at a lower frequency (250
Hz), the feedback regulation stage runs at a higher frequency (1 kHz) using
the sensor feedback to realize compensations of the reference trajectories
that guarantee the stability of the walking limit cycle. The end result are
policies learned from scratch that are transferred successfully to hardware
with minimal tuning. The controller is exhaustively tested in hardware
demonstrating stable walking gaits that are robust to external disturbances
and challenging terrain without being trained under those conditions.

## References

  * [1] M. Shafiee-Ashtiani, A. Yousefi-Koma, and M. Shariat-Panahi, “Robust bipedal locomotion control based on model predictive control and divergent component of motion,” in _2017 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2017, pp. 3505–3510.
  * [2] H. Chen, B. Wang, Z. Hong, C. Shen, P. M. Wensing, and W. Zhang, “Underactuated motion planning and control for jumping with wheeled-bipedal robots,” _IEEE Robotics and Automation Letters_ , 2020.
  * [3] Q. Nguyen and K. Sreenath, “Optimal robust safety-critical control for dynamic robotics,” _arXiv preprint arXiv:2005.07284_ , 2020.
  * [4] A. Hereid, E. A. Cousineau, C. M. Hubicki, and A. D. Ames, “3D dynamic walking with underactuated humanoid robots: a direct collocation framework for optimizing hybrid zero dynamics,” in _Proc. IEEE Int. Conf. Robotics and Automation (ICRA)_. Stockholm, Sweden: IEEE, May 2016, pp. 1447–1454.
  * [5] T. P. Lillicrap, J. J. Hunt, A. Pritzel, N. Heess, T. Erez, Y. Tassa, D. Silver, and D. Wierstra, “Continuous control with deep reinforcement learning,” _arXiv preprint arXiv:1509.02971_ , 2015.
  * [6] C. Yang, K. Yuan, W. Merkt, T. Komura, S. Vijayakumar, and Z. Li, “Learning whole-body motor skills for humanoids,” in _2018 IEEE-RAS 18th International Conference on Humanoid Robots (Humanoids)_. IEEE, 2018, pp. 270–276.
  * [7] Z. Xie, H. Y. Ling, N. H. Kim, and M. van de Panne, “Allsteps: Curriculum-driven learning of stepping stone skills,” _arXiv preprint arXiv:2005.04323_ , 2020.
  * [8] J. Morimoto, G. Cheng, C. G. Atkeson, and G. Zeglin, “A simple reinforcement learning algorithm for biped walking,” in _IEEE International Conference on Robotics and Automation_. IEEE, 2004, pp. 3030–3035.
  * [9] J. Morimoto, J. Nakanishi, G. Endo, G. Cheng, C. G. Atkeson, and G. Zeglin, “Poincare-map-based reinforcement learning for biped walking,” in _Proceedings of the 2005 IEEE International Conference on Robotics and Automation_. IEEE, 2005, pp. 2381–2386.
  * [10] X. Da, R. Hartley, and J. W. Grizzle, “Supervised learning for stabilizing underactuated bipedal robot locomotion, with outdoor experiments on the wave field,” in _2017 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2017, pp. 3476–3483.
  * [11] G. A. Castillo, B. Weng, A. Hereid, Z. Wang, and W. Zhang, “Reinforcement learning meets hybrid zero dynamics: A case study for rabbit,” in _2019 International Conference on Robotics and Automation (ICRA)_ , May 2019, pp. 284–290.
  * [12] T. Li, H. Geyer, C. G. Atkeson, and A. Rai, “Using deep reinforcement learning to learn high-level policies on the atrias biped,” in _2019 International Conference on Robotics and Automation (ICRA)_. IEEE, 2019, pp. 263–269.
  * [13] G. A. Castillo, B. Weng, W. Zhang, and A. Hereid, “Hybrid zero dynamics inspired feedback control policy design for 3d bipedal locomotion using reinforcement learning,” in _2020 IEEE International Conference on Robotics and Automation (ICRA)_ , 2020, pp. 8746–8752.
  * [14] J. Lee, J. Hwangbo, L. Wellhausen, V. Koltun, and M. Hutter, “Learning quadrupedal locomotion over challenging terrain,” _Science robotics_ , vol. 5, no. 47, 2020.
  * [15] W. Yu, V. C. Kumar, G. Turk, and C. K. Liu, “Sim-to-real transfer for biped locomotion,” _arXiv preprint arXiv:1903.01390_ , 2019.
  * [16] J. Siekmann, Y. Godse, A. Fern, and J. Hurst, “Sim-to-real learning of all common bipedal gaits via periodic reward composition,” _arXiv preprint arXiv:2011.01387_ , 2020.
  * [17] J. Hurst, “Building Robots That Can Go Where We Go,” _IEEE Spectrum: Technology, Engineering, and Science News_ , Feb. 2019. [Online]. Available: https://spectrum.ieee.org/robotics/humanoids/building-robots-that-can-go-where-we-go
  * [18] A. Hereid and A. D. Ames, “FROST: fast robot optimization and simulation toolkit,” in _Proc. IEEE/RSJ Int. Conf. Intelligent Robots and Systems (IROS)_. Vancouver, BC, Canada: IEEE/RSJ, Sept. 2017, pp. 719–726.
  * [19] E. Todorov, T. Erez, and Y. Tassa, “MuJoCo: a physics engine for model-based control,” in _2012 IEEE/RSJ International Conference on Intelligent Robots and Systems_ , Oct. 2012, pp. 5026–5033.
  * [20] Y. Gong, R. Hartley, X. Da, A. Hereid, O. Harib, J.-K. Huang, and J. Grizzle, “Feedback control of a Cassie bipedal robot: walking, standing, and riding a segway,” _American Control Conference (ACC)_ , 2019.
  * [21] X. Da, O. Harib, R. Hartley, B. Griffin, and J. W. Grizzle, “From 2D design of underactuated bipedal gaits to 3D implementation: Walking with speed tracking,” _IEEE Access_ , vol. 4, pp. 3469–3478, 2016.
  * [22] S. Rezazadeh, C. Hubicki, M. Jones, A. Peekema, J. Van Why, A. Abate, and J. Hurst, “Spring-mass walking with Atrias in 3D: robust gait control spanning zero to 4.3 kph on a heavily underactuated bipedal robot,” in _ASME 2015 Dynamic Systems and Control Conference_. American Society of Mechanical Engineers, 2015.
  * [23] “Robust Feedback Motion Policy Design Using Reinforcement Learning on a 3D Digit Bipedal Robot,” https://youtu.be/j8KbW-a9dbw, accessed: 2021-03-29.


# Phase Diagram and Crossover Phases of Topologically Ordered Graphene Zigzag
Nanoribbons: Role of Localization Effects

Hoang Anh Le  Department of Physics, Korea University, Seoul 02841, Korea
these authors contributed equally to this work  0000-0002-1668-8984 In Hwan
Lee  Department of Physics, Korea University, Seoul 02841, Korea  these
authors contributed equally to this work  0000-0002-7619-3780 Young Heon Kim
Department of Physics, Korea University, Seoul 02841, Korea
0000-0003-0000-539X S.-R. Eric Yang  Department of Physics, Korea University,
Seoul 02841, Korea  corresponding author<EMAIL_ADDRESS>0000-0003-3377-1859

###### Abstract

We computed the phase diagram of the zigzag graphene nanoribbons as a function
of on-site repulsion, doping, and disorder strength. The topologically ordered
phase undergoes topological phase transitions into crossover phases, which are
new disordered phases with a nonuniversal topological entanglement entropy
with significant variance. The topological order is destroyed by competition
between localization effects and on-site repulsion. We found that strong on-
site repulsion and/or doping weakens the nonlocal correlations between the
opposite zigzag edges. In one of the crossover phases, both $\frac{e^{-}}{2}$
fractional charges and spin-charge separation were absent; however, charge-
transfer correlations between the zigzag edges were possible. Another
crossover phase contains $\frac{e^{-}}{2}$ fractional charges, but no charge
transfer correlations. In low-doped zigzag ribbons the interplay between
electron localization and on-site repulsion contributes to the spatial
separation of quasi-degenerate gap-edge states and protects the charge
fractionalization against quantum fluctuations. In all these effects, mixed
chiral gap-edge states play an important role. The properties of
nontopological strongly disordered and strongly repulsive phases are also
observed. Each phase of the phase diagram has a different zigzag-edge
structure.

Is the formation of a fractional charge [1, 2, 3, 4, 5, 6] a necessary and
sufficient condition[7] for topologically ordered insulators [8] such as
fractional quantum Hall systems [9] and interacting disordered zigzag graphene
nanoribbons [10] (ZGNRs) with anyonic fractional charges? This issue is
related to whether the electron localization effects of doped systems destroy
or enhance the topological order [11, 12, 13] and quantization of fractional
charges. In a Laughlin state on a disordered sphere (no edges are present),
the added electrons fractionalize and form a quasi-degenerate peak in the gap
of the tunneling density of states (DOS). Electron localization [14] is
expected to suppress the quantum fluctuations of these fractional charges of
the quasi-degenerate gap states because these localized quasi-degenerate
energy states are spatially separated from each other, as explained in Ref.
[15] (if fractional charges are delocalized they overlap and become ill-
defined). However, excessive disorder is considered detrimental to topological
order.

In this study, we investigate similar issues with ZGNRs [16, 17, 18, 19, 20,
21, 22]. A recent study showed that weak randomness (disorder) in ZGNRs can
generate $e^{-}/2$ fractional charges [23, 24], which is a disorder effect
closely related to the change in the disorder-free symmetry-protected
topological insulator of ZGNRs to a topologically ordered [25] Mott-Anderson
insulator [26, 27, 28]. These systems have a universal value for topological
entanglement entropy (TEE) [11, 12] in the weak-disorder regime [29]. The
shape of entanglement spectrum is also found [30] to be similar to the DOS of
the edge states, as expected of topologically ordered systems [13]. In
interacting disordered ZGNRs, the gap is filled further by edge states [31,
23] with an increasing strength of the disorder potential. (We call these
states gap-edge states.)

The ground states have the opposite edge site spins in the absence of disorder
[16, 17, 18, 19, 20, 21, 22]. In the presence of disorder a spin
reconstruction of the zigzag edges can take place [25]. Nonetheless, a
topologically ordered ZGNR has two degenerate ground states, see Fig. 1(a).
Mixed chiral edge states play an important role in this effect. A short-range
disorder potential couples two nearly chiral gap-edge states residing on
opposite zigzag edges [23, 32], and mixed chiral gap-edge states with split
probability densities may form to display $e^{-}/2$ fractional semion charges
[33] (see Fig.1(b))(these states with midgap energies are solitonic with the
half of the spectral weight originating from the conduction band and the other
half from the valence band [34]). Note that a mixed chiral gap-edge state has
a nonzero fractional probability at the A- and B-carbon sites. In other words,
it is split into two nonlocal parts, each residing on the edges of the A or B
sublattice. The formation of mixed chiral gap-edge states is a nonperturbative
instanton effect [23]. (They are similar to the bonding and antibonding states
of a double quantum well.) It should be noted that well-defined $e^{-}/2$
fractional charges in the weak-disorder regime are emergent particles, i.e.,
they have new qualitative features and appear only in sufficiently long
ribbons. In a weak-disorder regime, the number of fractional charges is
proportional to the length of zigzag edges.

Figure 1:  (a) Two degenerate disordered ground states are illustrated on the
surface of a cylinder. The site spins of these ground states are reversed. The
expectation values of site spins that point outward (inward) normal to the
cylinder surface are solid red (dashed blue) vectors. Nearly zero-energy gap-
edge states exist on the zigzag edges. (b) Chiral edge states exist near the
Brillouin zone boundaries. $\textrm{u}(\textrm{l})$ stands for a chiral state
localized on the upper (lower) zigzag edge and $\sigma=\uparrow,\downarrow$
for the spin components along the $z$-axis. Short-range disorder potential
couples two nearly chiral edge states: a $(\textrm{u},\sigma)$ chiral edge
state is coupled to a $(\textrm{l},\sigma)$ chiral edge state This process
generates a mixed chiral state with two fractional charges.

Although weak disorder leads to formation of fractional charges strong
disorder may destroy them. Similar to fractional quantum Hall systems, the
topological order of a ZGNR is not immediately destroyed upon doping because
electron localization partially suppresses quantum fluctuations between quasi-
degenerate mid-gap states. The system may still be an insulator with a
fractional charge. However, in the presence of strong disorder or doping,
zigzag edge antiferromagnetism is expected to diminish, and thereby, the
topological order. (Away from the low doping region, a disordered anyon phase
with a distorted edge spin density wave was found [30].) These results suggest
that there may be several topological phase transitions in the zigzag ribbons.
What is the nature of these topological phase transitions and the physical
properties of the ground states? Does the presence of a fractional charge
imply a universal value of the TEE? Does the TEE become nonuniversal and vary
[35] with an increase in the disorder strength or doping level?

We explored the phase diagram of ZGNRs in the parameter space comprising on-
site repulsion ($U$), disorder strength ($\Gamma$), and doping concentration
($\delta N/N_{s}$) ($\delta N$ and $N_{s}$ are, respectively, the number of
doped electrons and the total number of sites in the ribbon). The competition
between localization and electron interactions can have detrimental effects on
the topological order and lead to several different phases, which includes
crossover phases. We found a number of different phases with a topological
order, quasi-topological order, and no order. Each of these phases is defined
by the value of TEE $\beta$ and its variance. These properties of $\beta$ are
related to the presence or absence of charge fractionalization and charge
transfer correlations between zigzag edges. When both of these properties are
present, in addition to correlations leading to spin-charge separation,
$\beta$ is universal, with small variances. In low-doped ZGNRs the interplay
between electron localization and on-site repulsion contributes to the spatial
separation of quasi-degenerate gap-edge states and protects the charge
fractionalization against quantum fluctuations. There are two other types of
phases with a quasi-topological order. We refer to these phases as crossover
phases, in which the variance of $\beta$ is significant. In one of these
phases, both $e^{-}/2$ fractional charges and spin-charge separation are
absent; however, the charge transfer ($\pm e^{-}/2$) correlations exist
between the zigzag edges. Another phase may contain stable $e^{-}/2$
fractional charges but no charge transfer correlations between the zigzag
edges. The ground state and zigzag edge properties of the various crossover
and nontopological phases are explored.

## Model Hamiltonian

The following mechanisms can all lead to fractional charges: the coupling
between the valleys mediated by short-range scatterers [23] and the sublattice
mixing facilitated by alternation of the nearest neighbor hopping parameters
[25]. Here we will consider only the effect of short-range scatterers. The
self-consistent Hartree-Fock (HF) approximation works well for graphene
systems [18, 19, 36]. The HF Hamiltonian of a ZGNR with length $L$ and width
$W$ is

$\displaystyle
H_{\text{MF}}=-t\sum_{n.n.,\sigma}c^{{\dagger}}_{i,\sigma}c_{j,\sigma}+\sum_{i,\sigma}V_{i}c_{i,\sigma}^{{\dagger}}c_{i,\sigma}+U\sum_{i}[n_{i,\uparrow}\langle
n_{i,\downarrow}\rangle+n_{i,\downarrow}\langle n_{i,\uparrow}\rangle-\langle
n_{i,\downarrow}\rangle\langle n_{i,\uparrow}\rangle]+\sum_{i}[s_{ix}\langle
h_{ix}\rangle+s_{iy}\langle h_{iy}\rangle],$ (1)

where the site index is given by $i=(k,l)$ ($k$ labels sites along the ribbon
direction and $l$ along the width), $c^{\dagger}_{i,\sigma}$ and
$n_{i_{,}\sigma}$ represent creation and occupation operators at site $i$ with
spin $\sigma=\left\\{\uparrow,\downarrow\right\\}$, respectively (periodic
boundary conditions are used along the ribbon direction). The site spin
operators are given by
$s_{ix(y)}=\frac{1}{2}(c^{\dagger}_{i,\uparrow},c^{\dagger}_{i,\downarrow})\sigma^{x(y)}(c_{i,\uparrow},c_{i,\downarrow})^{T}$,
where $\sigma^{x(y)}$ is the conventional Pauli matrix. The first term
represents the kinetic energy with hopping parameter $t$, $n.n$ implies the
summation over the nearest-neighbor sites. The second term represents the
short-range impurities parameterized by $V_{i}$, which is randomly chosen from
the energy interval $\left[-\Gamma,\Gamma\right]$. Throughout this study, the
density of the impure sites is fixed at $10\%$. $U$ denotes the on-site
repulsive strength. The last term in Eq. (1) represents self-consistent
“magnetic fields", where $\langle h_{ix}\rangle=-2U\langle s_{ix}\rangle$ and
$\langle h_{iy}\rangle=-2U\langle s_{iy}\rangle$. (These fields are present
only in doped ZGNRs. In the initial stage of the HF iteration, the values of
$\langle h_{ix}\rangle$ and $\langle h_{iy}\rangle$ can be selected from small
random numbers). In the presence of these fields, the HF eigenstates are mixed
spin states. The HF single-particle states $\left|k\right>$
($k=1,2,\ldots,2N_{s}$) can be written as a linear combination of site states
$\left|i,\sigma\right>$. In the language of second quantization this is
equivalent to

$a_{k}=\sum_{i,\sigma}A_{k,i,\sigma}c_{i,\sigma}.$ (2)

These magnetic fields are rather small for the disorder strength and doping
level considered in this study.

There may be several nearly degenerate HF ground states. We select the HF
initial ground state such that $\langle n_{i,\sigma}\rangle$ represents a
paramagnetic state with a small spin splitting. In addition, we choose small
random numbers of $\langle h_{ix}\rangle$ and $\langle h_{iy}\rangle$ (they do
not significantly affect the final results). The HF matrix dimension scales
with the number of carbon atoms, which is typically $<50000$. The HF
eigenstates and eigenenergies are self-consistently computed (this requires
approximately $20$ iterations). The TEE is computed using the disorder-
averaging results of numerous disorder realizations. Here, we used gpu to
speed up the solution of the HF matrix. The gpu calculations were intensive
and performed on a supercomputer. In the presence of disorder and in the low-
doping region, the obtained HF ground-state properties with solitons are in
qualitative agreement with those of the density matrix renormalization group
(DMRG) in the matrix product representation [30]. (In this work we do not
investigate the high doping region. The DMRG result is difficult to obtain in
this region because the computation is rather time consuming, and therefore,
it is not possible to determine which nearly degenerate HF ground state is the
true ground state.)

Note that the Mott gap $\Delta$ is well-developed only when $L\gg W$ (the
excitation spectrum of a ribbon with $L\sim W$ is similar to that of a gapless
two-dimensional graphene sheet [37]). The localization properties of ZGNRs are
unusual because both localized and delocalized states can exist [23, 32]. Gap-
edge states with energy $|E|\sim\Delta/2$ ($\Delta$ represents the Mott gap in
the absence of disorder) can have localization lengths $\sim W$ and overlap
significantly with each other.

## Effects of Anderson localization

Anderson localization plays a crucial role in the quantization of fractional
charges [15]. The effects of Anderson localization can be described using
self-consistent Hartree–Fock approximation (HFA) [38, 39]. The first important
effect is as follows: Anderson localization reduces the correlation length.
Therefore, in comparison to that of the nondisordered case, we can use a
smaller Wilson loop [12] to calculate the TEE of disordered interacting ZGNRs
. The correlation length may be computed from the entanglement entropy of an
area $A$. It is computed from the correlation function, which is also known as
the reduced density matrix of region $A$. The HF correlation function [29, 40]
between $i\in A$ and $j\in A$ decays exponentially as a function of distance
$x$ between $i$ and $j$

$C(x)=C_{i\uparrow,j\uparrow}=\langle\Psi|c^{\dagger}_{i\uparrow}c_{j\uparrow}|\Psi\rangle\sim\exp\left(-\frac{|x|}{\xi}\right),$
(3)

where $\Psi$ represents the HF ground state of the ZGNR and $\xi$ represents
the correlation length. By inverting the relation given in Eq.(2), we can
write $c_{i\sigma}$ as a linear combination of $a_{k}$. This makes it
straightforward to compute $C(x)$. To compute it accurately, the area must be
larger than the correlation length [7, 41].

We compute the correlation function and determine the correlation length
$\xi$, as shown in Fig. 2. For the disordered case $\Gamma\neq 0$, the
correlation length is obtained by averaging over several disorder
realizations. The Anderson localization reduces the correlation length
compared to that of disorder-free ribbons, as shown in Fig. 2(b). (Disorder-
free ZGNRs have a large correlation length for small $U$.) In contrast, doping
increases the correlation length, as shown in Fig. 2(c).

Figure 2: (a) ZGNRs consisting of two triangular sublattices A (black sites)
and B (grey sites). Each site of a ZGNR is labeled by two indices $(k,l)$ that
label blue and red carbon lines. The line $l=1$ comprises zigzag edge sites.
These sites are not connected to each other through the hopping term of the
Hamiltonian. (b) Correlation length $\xi$ is estimated by the fitting
correlation function with an exponential function along the zigzag edges.
Here, the ribbon length and width are $(L,W)=(100,64)$. At each data point,
$N_{D}=5$ disorder realizations were made. (c) Each point is calculated in
$(U,\Gamma)=(t,t)$ and $\delta N=0,2,5,10,20,30$. At each data point,
$N_{D}=6$ disorder realizations were created. (d) Fractional charge
probability densities of two quasi-degenerate mixed chiral states with
$E/(\Delta/2)=-0.005$ (left) and $E/(\Delta/2)=0.004$ (right) taken from a
single-disorder realization in the weak disorder regime at low doping. They
have the same spin. The Mott gap in this case is $\Delta=0.2t$. The green
(blue) surface refers to the probability density of sites A (B). On-site
repulsion, disorder strength, and number of doped electrons are
$(U,\Gamma,\delta N)=(t,0.01t,3)$. The geometry of the ribbon is
$(L,W)=(101,8)$.

Another important effect of Anderson localization in the presence of on-site
repulsion is that quasi-degenerate localized states are spatially separated
[14, 31], leading to well-defined fractional charges [15]. (In low-doped ZGNRs
added electrons fractionalize and form a narrow peak in the DOS near $E=0$
consisting of quasi-degenerate localized states [30].) The probability
densities of such two mid-gap states carrying fractional charges are shown in
Fig. 2(d). These gap-edge states are mixed chiral states [23], whose
probability densities peak at the two edges and rapidly decays inside the
ribbon. Note that these states do not overlap with each other. Non-interacting
electrons of disordered ZGNRs also display mixed chiral states near $E=0$.
However, although the overlap between nearly degenerate states in the weak-
disorder regime is small, it is not negligible. Thus well-defined fractional
charges do not readily form in non-interacting disordered ZGNRs.

## Phase Diagram

Topological order can be detected by investigating the TEE ($\beta$) [11, 12]
within the HFA [29]. We first select a set of values for
$(L,W,w,l_{\text{zig}},l_{\text{arm}})$, as defined in Fig. 3(a), to compute
$\beta$. Next, these quantities are increased by the same ratio and a new
$\beta$ is computed. This process is repeated several times (see Ref. [29] for
details). We apply finite-size scaling analysis to extract the value of the
TEE in the limit $L\rightarrow\infty$ (see Fig. 3(b)). We divide the parameter
space $(\Gamma,U,\delta N)$ into three-dimensional grid points, and at each
grid point, we compute $\beta$ (see Fig. 3(c)). The three-dimensional phase
diagram obtained is shown in Fig. 3(d). We find that $\beta$ can have three
types of values: (i) A universal value in the topologically ordered phase,
(ii) nonuniversal values of $\beta$ with large variances in the crossover
phases, and (iii) a zero value of $\beta$ in the normal-disordered phase.
Projections of the phase diagram, namely $U$-$\Gamma$, $\Gamma$-$\delta N$,
and $U$–$\delta N$ planes, are shown in Figs. 3(e)-(g).

In undoped ZGNRs a TO phase is found in regions $\Gamma/U\lesssim 1$ and
$U\lesssim t$, see Fig. 3(c). The topological phase transition into the
symmetric protected phase at $\Gamma=0$ is abrupt, consistent with the result
of Ref. [29] (The TEE of the symmetric phase is zero). There are also other
topological phase transitions but they are smooth transitions with crossover
regions [35]. Figs. 3(e)-(g) display the presence of crossover regimes lying
beyond the TO phase with an increase in the disorder, doping, and interaction
strength. The phase boundaries between topologically ordered and normal phases
are “blurred”, which indicate the presence of crossover phases (there are two
types of crossover (CO) phases, labeled COI and COII). The numerical results
of the TEE is shown in Fig. 3(c). The error bars in this figure include,
besides random fluctuations caused by disorder, the uncertainties that occur
in the extrapolation process of the finite scaling analysis. As $\Gamma/U$
increases, $\beta$ decreases (see the red line in Fig. 3(c)). The value of the
TEE thus changes across a crossover phase. In such a phase, $\beta$ has a
large variance, but the average values are not zero, which implies that the
topological order is not completely destroyed. In this regime, the TEE becomes
nonuniversal and decays. In crossover phases charge transfer correlations
between the opposite zigzag edges are present but fractional charges are not
well defined, or vice versa.

One can use a different but equivalent procedure to determine the phase
diagram. We verified that the same phase diagram can be obtained by analyzing
the presence of fractional charges and nonlocal correlations between the
opposite zigzag edges. At each grid point $(\Gamma,U,\delta N)$ in the
parameter space, we find the ground state and investigate whether the gap-edge
states display fractional charges and whether nonlocal correlations exist
between the opposite zigzag edges. By ultilizing this method, we have
successfully recovered the phase diagram shown in Fig. 3(d).

Figure 3: (a) Schematic to explain the geometric parameters of a ribbon. $W$
represents the ribbon width, which counts both $A$ and $B$ sites and is marked
by a solid line. $l_{\text{zig}},l_{\text{arm}}$, and $w$ are respectively,
zigzag length, arm length, and width of the Wilson loop. (b) Finite-size
scaling analysis of TEE in universal regime (blue circles, $(U,\Gamma,\delta
N)=(t,t,0)$) and strongly disordered regime (red triangles, $(U,\Gamma,\delta
N)=(t,15t,0)$). Ribbon sizes used in the analysis are
$(L,W,w,l_{\text{zig}},l_{\text{arm}})$:
$(200,172,12,171,68),\;(170,144,10,145,60),\;(130,112,7,111,46),\;(100,64,4,81,28),\;(80,68,4,69,24)$;
and $(50,32,2,41,12)$. At each data point, $N_{D}=50$ disorder realizations
were made. (c) TEE calculation for an undoped ZGNR $(\delta N=0)$: red squares
correspond to TEE with $U=0.5t$, black circles are with $U=t$, yellow hexagons
are with $U=1.5t$, and blue triangles are with $U=2t$. The values of TEE were
estimated by finite-size scaling. For each data point, the number of disorder
realizations is $N_{D}\sim 30$. The value of TEE is zero at $\Gamma/U=0$. (d)
Schematic phase diagram of finite-length ZGNRs obtained using the Hubbard
model. Regions with degraded color indicate crossover phases. (e), (f), (g)
Projections of phase diagrams onto 2D planes. Acronyms TO and COI (COII) stand
for topological order and crossover phase I (II), respectively. Here
$N_{\text{s}}=L\times W$ represents the total numbers of carbon sites of the
ribbon and $\delta N/N_{\text{s}}$ represents the doping concentration. The
degradation of color represents the variance of TEE (red means small
variance), and the white-color regions represent the trivial TEE ($\beta=0$).

## Topologically ordered phase

The universal region was investigated in Ref. [25], and therefore, we do not
describe this phase in detail here. However, we would like to mention some new
results. We elucidate the nature of correlations in topologically ordered
ZGNRs. A ZGNR is shown in Fig. 4(a). It consists of $8$ carbon lines labeled
$l=1,2,3,4,5,6,7,8$. In each pair of carbons lines $(1,8)$, $(2,7)$, $(3,6)$,
and $(4,5)$, an increase/decrease in the occupation number of one line is
correlated with a decrease/increase in that of the other line (see, for
example, lines $1$ and $8$ in Fig. 4(c)). It is not only the zigzag edges that
are correlated in this way, but also other carbon lines inside the ribbon that
are away from the edges. The corresponding site spins of the ribbon are shown
in Fig.4 (d). Mixed chiral gap states contribute to this effect (these gap-
edge states can decay slowly from the zigzag edges, unlike the fractional edge
states. A schematic picture of a mixed chiral state is shown in Fig. 4(b)).
Changes in the occupation numbers $\delta n_{i,\uparrow}$ and $\delta
n_{i,\downarrow}$ of an edge often coincide at nearly the same values as $k$,
which labels the site position along the ribbon direction. This effect can
lead to $n_{i,\uparrow}\approx n_{i,\downarrow}$ of the occupation numbers in
the presence of disorder, resulting in $s_{i}\approx 0$, i.e., the appearance
of spin-charge separation around a site on one of the edges [25].

Figure 4: (a) A ZGNR comprises lines of carbon atoms, labeled $1,2,3,\ldots$.
Red arrows connect two opposite lines that support charge transfers. (b)
Schematic of the probability densities of chiral and mixed chiral states. A
mixed state (black) comprises A-(red) and B-(blue) probability components
while a chiral state has only one component. The depicted mixed state in the
figure exhibits the possibility of overlapping A- and B-components; however,
it is also worth noting that they can be non-overlapping. (c) Change in the
spin-up occupation number $\delta
n_{i,\uparrow}=n_{i,\uparrow}-n^{0}_{i,\uparrow}$ of a ribbon size
$(L,W)=(120,8)$ with $(U,\Gamma,\delta N)=(2t,0.2t,0)$ (here
$n^{0}_{i,\uparrow}$ denotes the occupation numbers in the absence of
disorder). Occupation numbers are computed by single disorder realization
$N_{D}=1$. The $x$-axis label $k$ denotes sites along each carbon line
$l=1,2,\ldots,8$ (see Fig.2(a)). Blue arrows signify where charge transfer
takes place. Black dashed lines indicate where the occupation number changes.
Similar results were also found for spin-down occupation numbers. (d) Site
spins are displayed.

The following points should be also noted. The results in Fig.3(c) show that
the variance of $\beta$ decreases in the singular limit $\Gamma/U\rightarrow
0$. (Additional numerical results confirm this conclusion.) This result is
consistent with the previous finding that fractional charge of a midgap state
becomes accurate in the weak disorder regime and in the thermodynamic limit
(see Ref.[30]). In the opposite limit $\Gamma/U\gg 1$, the value of the TEE is
non-universal and decreases with increasing $U$ (see Fig. 3(c)). In addition,
the functional dependence of DOS on $E$ in the universal region is given by an
exponentially suppressed function, a linear function [25], or something in
between. The actual shape of the DOS is determined by the competition between
the strength of disorder and the on-site repulsion [28]; for example, the DOS
is linear for $(U,\Gamma,\delta N)=(2t,0.5t,0),$ but it is exponentially
suppressed in the weak disorder limit.

## Crossover phase I

We describe in detail the properties of undoped ZGNRs in the COI phase, where
$U\gtrsim t$ and $U\gtrsim\Gamma$ (the on-site repulsion $U$ is the dominant
energy in this phase). The TO phase gradually changes into the COI phase as
$U$ increases, as illustrated in Fig. 3(e). In this phase, $\beta$ is nonzero,
but its variance is significant, as shown in Fig 5(a). This phase has the
following properties: (i) The disorder-induced change in the edge occupation
numbers $\delta n_{i,\uparrow}=1/2$ for one type of spin $\sigma$ is entirely
transferred to the opposite edge, i.e., the zigzag edges are correlated in a
nontrivial manner (see Fig. 5(c) for $\Gamma/U=0.17$). However, the site
positions $k$ on the opposite zigzag edges, where changes in $\delta
n_{i,\uparrow}$ and $\delta n_{i,\downarrow}$ occur , do not coincide. (In
contrast, these positions are correlated at the nearly same values of $k$ in
the TO phase, as we mentioned before.) We believe that these edge transfer
correlations between zigzag edges change the ground state entanglement pattern
and yield a nonzero fluctuating TEE. The edge charge transfer correlations
become weaker when the disorder is stronger (see Fig. 5(d) for $\Gamma/t=4$),
leading to a smaller value of $\beta$ (see Fig. 3(c)). (ii) Although the
zigzag edge changes are fractional, $\delta n_{i,\uparrow}=1/2$ (Figs. 5(b)
and (c)), the A- and B-probability densities of the mixed chiral states
responsible for this feature overlap, see Fig. 4(c). Thus, fractional charges
are ill-defined. (iii) For spin-charge separation to be present, charge
transfers for both spins must occur at the same values $k$. These effects are
not observed in the COI phase. Note that, the condition
$S_{z}=\frac{1}{2}(n_{i,\uparrow}-n_{i,\downarrow})=0$ at site $i$ is not
sufficient for spin-charge separation. To fulfill the conditions, well-defined
fractional charges must exist.

Figure 5: (a) Distribution of $\beta$ values $(U,\Gamma,\delta N)=(3t,4.5t,0)$
with $N_{D}=134$. The dashed line indicates the mean value. (b)
$n_{i,\uparrow(\downarrow)}$ is calculated in a ribbon size $(L,W)=(120,8)$
with $(U,\Gamma,\delta N)=(3t,0.5t,0)$. The numbers $(1)$ and $(8)$ indicate
the upper and lower zigzag edges, as shown in Fig. 4(a). The $x$-axis label
$k$ denotes sites along each carbon line. Blue arrows signify where charge
transfer takes place. (In the lower graph, $n_{i,\uparrow(\downarrow)}$ are
displayed for the non-disorder case in the same condition as the upper graph.)
(c) Change in the spin-up (down) occupation number $\delta
n_{i,\uparrow(\downarrow)}=n_{i,\uparrow(\downarrow)}-n^{0}_{i,\uparrow(\downarrow)}$
of a ribbon size $(L,W)=(120,8)$ with $(U,\Gamma,\delta
N,N_{D})=(3t,0.5t,0,1)$. Edge site spins $S_{z}$ are also shown. (d) $\delta
n_{i,\uparrow}$ is displayed for a ribbon size $(L,W)=(120,8)$ with
$(U,\Gamma,\delta N,N_{D})=(3t,4t,0,1)$. In this case the correlations
disappear with increased disorder.

## Crossover phase II

For undoped ZGNRs, there is another CO phase for $\Gamma\gg U$ but
$U/t\lesssim 1$. We call this phase COII where the disorder strength $\Gamma$
is the dominant energy. As $\Gamma$ increases, the TO phase undergoes a
gradual transition into the COII phase, as demonstrated in Figure 3(e).
Concurrently, the gap is progressively filled with states, as depicted in the
upper graph of Figure 6(a). Similar to the COI phase (Fig. 5(a)) $\beta$ is
finite with a significant variance. But there are no charge-transfer
correlations between the zigzag edges (see the lower graph of Figure 6(a)).
However, some fractional charges may exist, see Fig. 6(b). This is consistent
with the following obtained results: (i) Some changes in the edge occupation
number are $\delta n_{i,\sigma}\approx\pm 1/2$. (ii) There are gap-edge states
with $q_{A}\approx 1/2$. (Here $q_{A}=\sum_{i\in A}|\psi_{i\sigma}(E)|^{2}$,
where $\psi_{i\sigma}(E)$ is the HF eigenstate with energy $E$, see Ref.[25].
The probability densities are summed over all sites of the A sublattice.)
However, the variance of $q_{A}$ in the energy interval $[E-\delta E,E+\delta
E]$ is large because $q_{A}$ varies substantially in this interval, as shown
in Fig. 6(a). (But the disorder averaged mean charge value of the states in
this energy interval is $e^{-}/2$.) Despite this, a fractional charge of a
state in the interval $[E-\delta E,E+\delta E]$ near $E=0$ does not overlap
significantly with the probability densities of other fractional and non-
fractional states in the same energy interval, provided that $\delta E$ is
small (for $U\sim t$ and $\Gamma\sim t$ this happens when $\delta E\sim
0.01t$). We believe that the interplay between localization and on-site
repulsion is responsible for this effect. However, since a gap is absent, the
fractional charges are less stable compared to the TO phase. We checked that
several states in the same value of $\delta E$ overlap in the absence of on-
site repulsion.

Thus far, we investigated the undoped case. Upon doping, the disorder-free
ZGNRs exhibit edge spin density waves instead of edge ferromagnetism of
undoped ribbons. If disorder is added to a doped ZGNR, the spin waves become
distorted [30]: There is a topological phase transition from modulated
ferromagnetic edges at zero doping to distorted spin-wave edges at finite
doping. Our results indicated that, when doping is substantial, this phase is
also a COII phase. The dependence of the mean value of $q_{A}$ of the states
in the mid-gap peak on the number of doped electrons is shown in Fig. 6(d)
(the DOS shows a sharp peak at the mid-gap energy, see Fig. 6(c)). At a low
doping concentration, the disorder-averaged value of $q_{A}$ is close to
$0.5$. The states in the mid-gap peak display well-defined fractional charges,
as we discussed below Fig.2 (note that the width of the midgap peak is $\delta
E\sim 0.005t$). As the doping concentration increases further, $q_{A}$
significantly deviates from 0.5, and simultaneously, the DOS mid-gap peak
starts to decrease [30]. These findings imply that even though fractional
charges can be found, their number decreases with an increase in doping. The
gradual change in $q_{A}$ as a function of $\delta N/N_{\text{s}}$indicates
that the transition from the phase of the distorted ferromagnetic edge to the
phase of the distorted edge spin-wave is not sharp. Figure 6(e) shows how
$\beta$ decreased with an increase in $\delta N/N_{\text{s}}$. For a large
$\delta N/N_{\text{s}}$, it is computationally demanding to calculate $\beta$
because the correlation length is expected to be longer (see Fig. 2(c)) in
comparison to that of undoped ZGNRs.

Figure 6: (a) $E$-$q_{A}$ diagram for the half-filling case $(U,\Gamma,\delta
N)=(t,5t,0)$. The geometry of the ribbon is $(L,W)=(100,8)$, and the number of
disorder realizations is $N_{D}=120$. Red (blue) dots represent states of spin
up (down). There are a lot of states that have $q_{A}\approx 1/2$ but
fractional charges barely exist. In lower graph, $\delta
n_{i,\uparrow(\downarrow)}$ is diplayed for a ribbon size $(L,W)=(120,8)$ with
$(U,\Gamma,\delta N)=(t,5t,0)$. (b) Plot of the probability density of several
mixed chiral gap-edge states are shown with $q_{A}\approx 1/2$ ($N_{D}=4$).
Green and blue represent densities on two different sublattices. All
parameters are similar to those of Fig. 6(a). Blue arrows point to a gap-edge
states that display fractional charges. The geometry of the ribbon is
$(L,W)=(100,8)$ (c) DOS of an undoped weakly disordered interacting ZGNRs is
shown as a dotted black line, which displays a soft gap [23] (left figure).
The mid-gap peak of the DOS in the low doping limit [29] is displayed as the
dotted red line. The mid-gap peak value of the DOS changes as a function of
doping concentration $\delta N/N_{\text{s}}$ (right figure, taken from Ref.
[30]) . Here, the size of the ribbon is $(L,W)=(101,8)$. The parameters are
$(\Gamma,U)=(0.01t,t)$. (d) Mean value of $q_{A}$ of the mid-gap states shown
in Fig. 6(c) is plotted as a function doping concentration. The number of
disorder realization is $N_{D}=200$. (e) Average beta is plotted for several
values of doping concentration. Here, we used $(U,\Gamma)=(t,t)$. The number
of disorder realization is $N_{D}=10$ for each point. Ribbon sizes used in the
analysis are $(L,W,w,l_{\text{zig}},l_{\text{arm}})$:
$(170,144,10,145,60),\;(130,112,7,111,46),\;(100,64,4,81,28),\;(80,68,4,69,24)$,
and $(50,32,2,41,12)$.

## Strongly disordered and strongly repulsive phases

We discuss the strongly disordered phase in region $\Gamma/U\gg 1$ (see Fig.
3(e)). The topological order is destroyed once the disorder strength reaches a
sufficiently large value (e.g., $\beta=0$ at $(U,\Gamma,\delta N)=(t,15t,0)$).
In this region the edge charge-transfer correlations and charge
fractionalization are not well-defined, which implies that TEE is zero. In
Fig. 7(a), site occupation numbers in weak ($\Gamma=0.03t$) and strong
($\Gamma=15t$) disorder regimes are shown side by side to highlight the
difference, where the ones in strong disorder regime highly fluctuate from
site to site. Edge magnetization is zero almost everywhere (see Fig. 7(b)).
The occupation numbers display sharp values of $n_{i,\sigma}=1$ at some sites
(see Fig. 7(b)) (they were also present in the DMRG calculations of disordered
ribbons, see Ref.[30]).

Figure 7: (a) Site spins and site occupation numbers for spin- up and -down
are shown for a ribbon size $(L,W)=(100,64)$ with $(U,\Gamma\delta
N)=(t,0.03t,0)$ and $(U,\Gamma\delta N)=(t,15t,0)$. Both cases are computed by
single disorder realization $N_{D}=1$. (b) Edge occupation numbers for spin-
up, spin-down, and edge site spins in strongly disordered regime:
$(U,\Gamma,\delta N)=(t,15t,0)$. Here, the geometry of the ribbon is
$(L,W)=(100,8)$ and these quantities are obtained by single disorder
realization $N_{D}=1$. (c) $E$ vs. $q_{A}$ graph with $(L,W)=(100,8)$,
$(U,\Gamma,\delta N)=(7t,0.5t,0)$ ($N_{D}=1$). Red (blue) dots are states of
spin up (down). Disorder-free DOS for $(U,\Gamma,\delta N)=(7t,0,0)$ shown in
the inset has a large gap, $\Delta\approx 6t$. Solid red and blue lines
represent $q_{A}$ values of spin-up and -down states for $\Gamma=0$. They are
zigzag gap-edge states localized either on the upper or lower zigzag edges.
(d) Probability density of several states near gap edges $\pm\Delta/2$
(horizontal thin black lines) for single disorder realization $N_{D}=1$ with
$q_{A}\approx 0.5$, $(U,\Gamma,\delta N)=(7t,0.5t,0)$. Here, the geometry of
the ribbon is $(L,W)=(100,8)$. (e) DOS with $(L,W)=(100,64)$ and
$(U,\Gamma,\delta N)=(4t,4t,0)$ ($N_{D}=100$ realizations were made). The
computed TEE value is $\beta=0$. (f) Zigzag edge occupation numbers
$n_{i,\sigma}$, $n_{i}=n_{i,\uparrow}+n_{i,\downarrow}$ and site spins $S_{z}$
are shown in the strong repulsive regime $(U,\Gamma,\delta N)=(5t,0.01t,0)$
with a ribbon size of $(L,W)=(120,8)$ and a disorder realization $N_{D}=1$.
(g) The same quantities are computed in $(U,\Gamma,\delta N)=(5t,0.1t,0)$ with
the same size of the ribbon with (f). At $\Gamma=0.1t$, there are more sites
that have $n_{i}\approx 2$ than those of $\Gamma=0.01t$.

Another phase, that is, the strongly repulsive phase ($U\gg t$) with no
fractional charges and zigzag edge correlations, is shown in Fig. 3(e). In
this case, $\beta\approx 0$. The $q_{A}$–$E$ diagram in Fig. 7(c) displays the
nonperturbative nature of disorder in this regime: the values $q_{A}$ are
scattered between $0$ and $1$ in the limit $\Gamma\rightarrow 0$, whereas they
are restricted to the four solid lines at $\Gamma=0$. Also, the $(E,q_{A})$
distribution indicates that in a strongly repulsive regime, even with the
presence of disorder, there is still a large energy gap. There are no states
with $q_{A}\approx 1/2$ near the mid-gap energy, as shown in Fig. 7(c). The A
and B components of the wave function of the states with $q_{A}\approx 1/2$
near the gap edges $\pm\Delta/2$ overlap (see Fig. 7(d)). For a stronger
disorder (larger values of $\Gamma$), the gap is filled with states such that
the DOS is finite at $E=0$, as shown in Fig. 7(e).

The main physics of this phase is illustrated by investigating the zigzag edge
structure: the occupation numbers are $n_{i,\sigma}=1$ or $0$ so charge
transfers are one ($\delta n_{i,\sigma}=\pm 1$) in the strongly repulsive
phase (the total site occupation number of each site is $n_{i}\approx 2$, $1$,
or $0$ despite a strong $U$). No transfer of fractional charges was observed
between zigzag edges. This is because mixed chiral gap edge states are not
present. Note that the edge magnetization displays sharp domain walls, as
indicated in Figs. 7(f)-(g).

## Summary and Discussion

We computed the phase diagram of zigzag graphene nanoribbons as a function of
the on-site repulsion $U$, doping $\delta N$, and disorder strength $\Gamma$.
We identified the universal, crossover, strongly disordered, and strongly
repulsive phases. Each phase of the phase diagram was defined by the TEE value
and its variance. We also investigated how the values of the TEE are related
to the following physical properties: the presence of charge fractionalization
and the edge charge transfer correlations between the opposite zigzag edges.
When both properties are present, in addition to correlations leading to spin-
charge separation, $\beta$ was universal. If only one of these properties is
present, $\beta$ was nonuniversal and its variance was significant. However,
when both were absent, $\beta$ was approximately zero. In addition, we found a
strongly repulsive phase with zero TEE in large on-site repulsion and weak
disorder limits. Its ground state contains abrupt kinks in zigzag edge
magnetizations without charge fractionalization, which is a consequence of the
singular perturbative nature of the disorder potential. There is another phase
with zero TEE, i.e., the strongly disordered phase in regime $\Gamma\gg U$. In
this phase, the edge site occupation numbers fluctuate highly from site to
site, and antiferromagnetic coupling between the two edges are nearly
destroyed. Each phase of the phase diagram has a different zigzag-edge
structure.

We also investigated the effect of the interplay between localization and on-
site repulsion on the charge quantization. In low-doped and/or weakly
disordered ZGNRs this interplay contributes to the spatial separation of
quasi-degenerate gap-edge states and protects the charge fractionalization
against quantum fluctuations. Even in the presence of moderately strong
disorder charge fractionalization is not completely destroyed.

We briefly discuss some experimental implications. It would be interesting to
observe the presence of nonlocal charge transfers between the zigzag edges of
the COI phase. This can be investigated by measuring correlations between the
edge site occupation numbers using a scanning tunneling microscope [42]. In
the COII phase, fractional $e^{-}/2$ edge charges are present; however,
unusual transport and magnetic susceptibility properties are not expected
because spin-charge separation is not present. (In contrast, the TO phase is
expected to display unusual transport and magnetic susceptibility because of
spin-charge separation [43, 25].) Similarly, scanning tunneling microscopy can
be used to verify the predicted edge occupation numbers of the strongly
repulsive and strongly disordered phases. In addition, investigation of
tunneling between zigzag edges, as in fractional quantum Hall bar systems [44,
45], may be fruitful.

## References

  * [1] Leinaas, J. M. & Myrheim, J.  On the theory of identical particles.  _Nuovo Cimento B_ 37, 1–23, DOI: 10.1007/BF02727953 (1977).
  * [2] Wilczek, F.  Quantum mechanics of fractional-spin particles.  _Phys. Rev. Lett._ 49, 957–959, DOI: 10.1103/PhysRevLett.49.957 (1982).
  * [3] Arovas, D., Schrieffer, J. R. & Wilczek, F.  Fractional statistics and the quantum hall effect.  _Phys. Rev. Lett._ 53, 722–723, DOI: 10.1103/PhysRevLett.53.722 (1984).
  * [4] Nakamura, J., Liang, S., Gardner, G. C. & Manfra, M. J.  Direct observation of anyonic braiding statistics.  _Nature Physics_ 16, 931–936, DOI: https://doi.org/10.1038/s41567-020-1019-1 (2020).
  * [5] Bartolomei, H. _et al._ Fractional statistics in anyon collisions.  _Science_ 368, 173–177, DOI: 10.1126/science.aaz5601 (2020).
  * [6] Girvin, S. M. & Yang, K.  _Modern condensed matter physics_ (Cambridge University Press, Cambridge, 2019).
  * [7] Pachos, J. K.  _Introduction to Topological Quantum Computation_ (Cambridge University Press, Cambridge, 2012).
  * [8] Wen, X.-G.  Colloquium: Zoo of quantum-topological phases of matter.  _Rev. Mod. Phys._ 89, 041004, DOI: 10.1103/RevModPhys.89.041004 (2017).
  * [9] Laughlin, R. B.  Anomalous quantum hall effect: An incompressible quantum fluid with fractionally charged excitations.  _Phys. Rev. Lett._ 50, 1395–1398, DOI: 10.1103/PhysRevLett.50.1395 (1983).
  * [10] S.-R. Eric Yang.  _Topologically Ordered Zigzag Nanoribbon_ (World Scientific, Singapore, 2023).
  * [11] Kitaev, A. & Preskill, J.  Topological entanglement entropy.  _Phys. Rev. Lett._ 96, 110404, DOI: 10.1103/PhysRevLett.96.110404 (2006).
  * [12] Levin, M. & Wen, X.-G.  Detecting topological order in a ground state wave function.  _Phys. Rev. Lett._ 96, 110405, DOI: 10.1103/PhysRevLett.96.110405 (2006).
  * [13] Li, H. & Haldane, F. D. M.  Entanglement spectrum as a generalization of entanglement entropy: Identification of topological order in non-abelian fractional quantum hall effect states.  _Phys. Rev. Lett._ 101, 010504, DOI: 10.1103/PhysRevLett.101.010504 (2008).
  * [14] Altshuler, B.  Inductory anderson localization.  In _Advanced Workshop on Anderson Localization, Nonlinearity and Turbulence: a Cross-Fertilization_ (International Centre for Theoretical Physics, 2010).
  * [15] Girvin, S. M.  The quantum hall effect: novel excitations and broken symmetries.  In _Aspects topologiques de la physique en basse dimension. Topological aspects of low dimensional systems_ , 53–175 (Springer, 1999).
  * [16] Fujita, M., Wakabayashi, K., Nakada, K. & Kusakabe, K.  Peculiar localized state at zigzag graphite edge.  _J. Phys. Soc. Jpn._ 65, 1920–1923, DOI: 10.1143/JPSJ.65.1920 (1996).
  * [17] Brey, L. & Fertig, H. A.  Electronic states of graphene nanoribbons studied with the dirac equation.  _Phys. Rev. B_ 73, 235411, DOI: 10.1103/PhysRevB.73.235411 (2006).
  * [18] Yang, L., Park, C.-H., Son, Y.-W., Cohen, M. L. & Louie, S. G.  Quasiparticle energies and band gaps in graphene nanoribbons.  _Phys. Rev. Lett._ 99, 186801, DOI: 10.1103/PhysRevLett.99.186801 (2007).
  * [19] Pisani, L., Chan, J. A., Montanari, B. & Harrison, N. M.  Electronic structure and magnetic properties of graphitic ribbons.  _Phys. Rev. B_ 75, 064418, DOI: 10.1103/PhysRevB.75.064418 (2007).
  * [20] Ruffieux, P. _et al._ On-surface synthesis of graphene nanoribbons with zigzag edge topology.  _Nature_ 531, 489–492, DOI: 10.1038/nature17151 (2016).
  * [21] Kolmer, M. _et al._ Rational synthesis of atomically precise graphene nanoribbons directly on metal oxide surfaces.  _Science_ 369, 571–575, DOI: 10.1126/science.abb8880 (2020).
  * [22] Brey, L., Seneor, P. & Tejeda, A. (eds.) _Graphene Nanoribbons_.  2053-2563 (IOP Publishing, 2019).
  * [23] Jeong, Y. H., S.-R. Eric Yang & Cha, M.-C.  Soliton fractional charge of disordered graphene nanoribbon.  _Journal of Physics: Condensed Matter_ 31, 265601, DOI: 10.1088/1361-648X/ab146b (2019).
  * [24] S.-R. Eric Yang.  Soliton fractional charges in graphene nanoribbon and polyacetylene: similarities and differences.  _Nanomaterials_ 9, 885, DOI: 10.3390/nano9060885 (2019).
  * [25] S.-R. Eric Yang, Cha, M.-C., Lee, H. J. & Kim, Y. H.  Topologically ordered zigzag nanoribbon: $e/2$ fractional edge charge, spin-charge separation, and ground-state degeneracy.  _Phys. Rev. Research_ 2, 033109, DOI: 10.1103/PhysRevResearch.2.033109 (2020).
  * [26] Dobrosavljevic, V., Trivedi, N. & Valles, J. M., Jr.  _Conductor-Insulator Quantum Phase Transitions_ (Oxford University Press, 2012).
  * [27] Belitz, D. & Kirkpatrick, T. R.  The anderson-mott transition.  _Rev. Mod. Phys._ 66, 261–380, DOI: 10.1103/RevModPhys.66.261 (1994).
  * [28] Byczuk, K., Hofstetter, W. & Vollhardt, D.  Competition between anderson localization and antiferromagnetism in correlated lattice fermion systems with disorder.  _Phys. Rev. Lett._ 102, 146403, DOI: 10.1103/PhysRevLett.102.146403 (2009).
  * [29] Kim, Y. H., Lee, H. J. & S.-R. Eric Yang.  Topological entanglement entropy of interacting disordered zigzag graphene ribbons.  _Phys. Rev. B_ 103, 115151, DOI: 10.1103/PhysRevB.103.115151 (2021).
  * [30] Kim, Y. H., Lee, H. J., Lee, H.-Y. & S.-R. Eric Yang.  New disordered anyon phase of doped graphene zigzag nanoribbon.  _Scientific Reports_ 12, 14551, DOI: 10.1038/s41598-022-18731-6 (2022).
  * [31] Éfros, A. L. & Shklovskii, B. I.  Coulomb gap and low temperature conductivity of disordered systems.  _Journal of Physics C: Solid State Physics_ 8, L49–L51, DOI: 10.1088/0022-3719/8/4/003 (1975).
  * [32] Lima, L. R. F., Pinheiro, F. A., Capaz, R. B., Lewenkopf, C. H. & Mucciolo, E. R.  Effects of disorder range and electronic energy on the perfect transmission in graphene nanoribbons.  _Phys. Rev. B_ 86, 205111, DOI: 10.1103/PhysRevB.86.205111 (2012).
  * [33] Canright, G. S., Girvin, S. M. & Brass, A.  Superconductive pairing of fermions and semions in two dimensions.  _Phys. Rev. Lett._ 63, 2295–2298, DOI: 10.1103/PhysRevLett.63.2295 (1989).
  * [34] Heeger, A. J., Kivelson, S., Schrieffer, J. R. & Su, W. P.  Solitons in conducting polymers.  _Rev. Mod. Phys._ 60, 781–850, DOI: 10.1103/RevModPhys.60.781 (1988).
  * [35] Zeng, B., Chen, X., Zhou, D.-L. & Wen, X.-G.  _Quantum Information Meets Quantum Matter_ (Springer New York, Springer, 2019).
  * [36] Stauber, T. _et al._ Interacting electrons in graphene: Fermi velocity renormalization and optical response.  _Phys. Rev. Lett._ 118, 266801, DOI: 10.1103/PhysRevLett.118.266801 (2017).
  * [37] Castro Neto, A. H., Guinea, F., Peres, N. M. R., Novoselov, K. S. & Geim, A. K.  The electronic properties of graphene.  _Rev. Mod. Phys._ 81, 109–162, DOI: 10.1103/RevModPhys.81.109 (2009).
  * [38] S.-R. Eric Yang, MacDonald, A. H. & Huckestein, B.  Interactions, localization, and the integer quantum hall effect.  _Phys. Rev. Lett._ 74, 3229–3232, DOI: 10.1103/PhysRevLett.74.3229 (1995).
  * [39] S.-R. Eric Yang & MacDonald, A. H.  Coulomb gaps in a strong magnetic field.  _Phys. Rev. Lett._ 70, 4110–4113, DOI: 10.1103/PhysRevLett.70.4110 (1993).
  * [40] Peschel, I.  Calculation of reduced density matrices from correlation functions.  _Journal of Physics A: Mathematical and General_ 36, L205, DOI: 10.1088/0305-4470/36/14/101 (2003).
  * [41] Jiang, H.-C., Wang, Z. & Balents, L.  Identifying topological order by entanglement entropy.  _Nature Phys_ 8, 902–905, DOI: 10.1038/nphys2465 (2012).
  * [42] Andrei, E. Y., Li, G. & Du, X.  Electronic properties of graphene: a perspective from scanning tunneling microscopy and magnetotransport.  _Reports on Progress in Physics_ 75, 056501, DOI: 10.1088/0034-4885/75/5/056501 (2012).
  * [43] Chung, T.-C., Moraes, F., Flood, J. D. & Heeger, A. J.  Solitons at high density in trans-${(\mathrm{CH})}_{x}$: Collective transport by mobile, spinless charged solitons.  _Phys. Rev. B_ 29, 2341–2343, DOI: 10.1103/PhysRevB.29.2341 (1984).
  * [44] Kang, W., Stormer, H., Pfeiffer, L., Baldwin, K. & West, K.  Tunnelling between the edges of two lateral quantum hall systems.  _Nature_ 403, 59–61, DOI: https://doi.org/10.1038/47436 (2000).
  * [45] Yang, I., Kang, W., Baldwin, K. W., Pfeiffer, L. N. & West, K. W.  Cascade of quantum phase transitions in tunnel-coupled edge states.  _Phys. Rev. Lett._ 92, 056802, DOI: 10.1103/PhysRevLett.92.056802 (2004).

## Acknowledgements

S.R.E.Y. was supported by the Basic Science Research Program of the National
Research Foundation of Korea (NRF), funded by the Ministry of Science and ICT
(MSIT) NRF-2021R1F1A1047759. Two other grants are also acknowledged: BK21 FOUR
(Fostering Outstanding Universities for Research) and KISTI Supercomputing
Center with supercomputing resources including technical support
KSC-2022-CRE-0345.

## Data availability

On reasonable request, the corresponding author will provide all relevant data
in this paper.

## Code availability

On reasonable request, the corresponding author will provide all numerical
codes in this paper.

## Author contributions

L.H.A., Y.H.K. and I.H.L. performed the HF calculations. S.R.E.Y. conceived
the project and supervised the study. All authors contributed to the writing
of the manuscript.

## Competing interests

The authors declare no competing interests.