# A Zero-Shot Adaptive Quadcopter Controller

Dingqi Zhang1, Antonio Loquercio2, Xiangyu Wu1, Ashish Kumar2, Jitendra
Malik2, Mark W. Mueller1 The authors are with the 1High Performance Robotics
Lab, Dept. of Mechanical Engineering, and the 2Dept. of Electrical Engineering
and Computer Science, University of California Berkeley, {dingqi, loquercio,
wuxiangyu, ashish_kumar, malik<EMAIL_ADDRESS>
###### Abstract

This paper proposes a universal adaptive controller for quadcopters, which can
be deployed zero-shot to quadcopters of very different mass, arm lengths and
motor constants, and also shows rapid adaptation to unknown disturbances
during runtime. The core algorithmic idea is to learn a single policy that can
adapt online _at test time_ not only to the disturbances applied to the drone,
but also to the robot dynamics and hardware in the same framework. We achieve
this by training a neural network to estimate a latent representation of the
robot and environment parameters, which is used to condition the behaviour of
the controller, also represented as a neural network. We train both networks
exclusively in simulation with the goal of flying the quadcopters to goal
positions and avoiding crashes to the ground. We directly deploy the same
controller trained in the simulation without any modifications on two
quadcopters with differences in mass, inertia, and maximum motor speed of up
to 4 times. In addition, we show rapid adaptation to sudden and large
disturbances (up to 35.7%) in the mass and inertia of the quadcopters. We
perform an extensive evaluation in both simulation and the physical world,
where we outperform a state-of-the-art learning-based adaptive controller and
a traditional PID controller specifically tuned to each platform individually.
Video results can be found at https://dz298.github.io/universal-drone-
controller/.

Figure 1: Demonstration of our adaptive controller on two quadrotors with
widely varying mass, arm length, and motor constants for the task of tracking
a straight and circular trajectory. In the middle of the trajectory, we add a
payload unknown to the drone of approximately 30% of the robot’s mass. Our
end-to-end controller is able to quickly estimate and react to the
disturbance. In both the demonstrations described above and all the simulation
and real-world experiments presented here, we use a single control policy
across different drones and tasks, which is deployed without any modifications
or tuning. Video results at https://dz298.github.io/universal-drone-
controller/

##  I Introduction

In this paper, we show the first universal controller for quadcopters.
Quadcopters are inherently unstable and require high-frequency control of up
to 500hz. Failure to adapt to perturbations within fractions of a second can
potentially lead to a crash. Consequently, a universal controller should not
only be able to handle diverse perturbations, it should do so very fast. In
this work, we learn a single controller which is capable of controlling a
variety of quadcopters with differences in mass, armlength, maximum motor
speed, etc. of up to 4 times within 0.8 second, without any modification or
fine-tuning for the specific quadcopters. In addition, our controller can
rapidly adapt to unknown disturbances in the mass and inertia changes of the
quadcopters. Our work builds on three main insights. The first is that the
dynamics parameters are observable from the history of states and actions.
This has long been known in the community: several works have shown how to
estimate parameters online with filters [1, 2] or with neural networks [3].
Yet, prior work overlooked the possibility of using these estimates to
condition the controller’s behaviour. The second insight is that, at test
time, we do not need an estimate of the parameters in some "ground truth"
sense. What matters is that the estimate leads to the "right" action, which
our end-to-end training procedure optimizes. This simplifies the estimation
problem and avoids identifiability issues, e.g., identical effects on
unobservable or unmatched uncertainties [4]. The third and final insight is
that a system can adapt to previously unseen disturbances as long as their
effect on the platform dynamics is in the convex hull of the training
disturbances (e.g., a motor losing efficiency has a similar effect to adding a
payload below the motor).

To operationalize this, we follow the approach initially proposed for legged
robots [5]. However, while this work performs online adaptation to terrains,
we use their approach to adapt to a diverse set of quadcopter bodies and
perturbations. Specifically, we estimate a latent representation of the
quadcopter’s body from a history of sensor observations and actions, which
conditions the behaviour of the controller. The space of quadcopter properties
in which we train is extremely diverse (Table I). This diversity enables
adaptation to sudden changes in environment conditions, e.g., a payload, for
which it was not explicitly trained. Our approach frees drone designers from
the estimation and tuning process required any time something changes, or the
risk that parameter changes unwittingly cause the control behaviour to
significantly change, potentially endangering the system. Furthermore, it
naturally lends itself for use in off-the-shelf autopilot systems, allowing
users who might otherwise lack the modelling abilities to control a custom
vehicle, simply by plugging in the autopilot and not requiring any parameter
tuning.

Figure 2: We show the training at the top and the deployment architecture of
our system. We train in two phases. In the first phase, we train a base policy
$\pi$ which takes the current state $x_{t}$, and the intrinsics vector $z_{t}$
which is a compressed version of the environment parameters $e_{t}$ generated
by the module $\mu$. Since we cannot deploy this policy in the real world
because we do not observe the environment parameters $e_{t}$, we learn an
adaptation module which takes the sensor history and action history, and
directly predicts the intrinsics vector $z_{t}$. This is done in phase two in
simulation using supervised learning. We can finally deploy the base policy
$\pi$ which takes as input the current state $x_{t}$ and the intrinsics vector
$\hat{z_{t}}$ predicted by the adaptation module $\phi$.

Our approach is related to existing methods in adaptive control. However, our
work shifts the meaning of adaptation to a different paradigm. While adaptive
control is generally concerned with estimating and counteracting disparities
between observations and a _reference_ model, our approach does not have this
notion. This key difference enables adaptation to a much wider range of
dynamics and disturbances. In addition, it waives the engineering and tuning
required by prior work when modifying the vehicle. Meta-learning is also
related to the context of our problem in providing fast online adaptation.
Although they have been demonstrated on real quadcopters with impressive
results on disturbance rejection [6, 7], they require real world learning
samples to adapt. When adapting to a much wider scale in our problem setting,
this approach could be potentially dangerous for a dynamic system like a
quadcopter, since failure to adapt quickly leads to catastrophic failure (a
crash). The closest related to our work are industrial systems like Beta-
Flight [8] and PX4 [9]. While they work well across a variety of vehicles,
they still require significant in-flight tuning for each robot type and are
generally tailored to human pilots.

##  II Related Work

The design of high-performance adaptive controllers for aerospace systems has
been a top priority for researchers and industry for more than 50 years. While
the overall goal remained the same over the decades, approaches greatly
evolved. In the following, we review both traditional and learning-based
solutions to adaptive control. While each method has its advantages and
limitations, they are mainly designed to handle model uncertainties and
disturbances. We are interested in exploring a much wider range of
adaptations, with substantial differences between platforms.

###  II-A Traditional Adaptive Control

One of the initial contributions in this space is the model reference adaptive
controller (MRAC), an extension of the well-known MIT-rule [10]. The empirical
success of this method sparked great interest in the aerospace community,
which led to the development of both practical tools and theoretical
foundations [11, 12]. From the many methods developed, one of the most popular
is $\mathcal{L}1$ adaptive control [13, 4]. The main reason behind its success
is the ability to provide rapid adaptation to model uncertainties and
disturbances with theoretical guarantees under (possibly restrictive)
assumptions. Its high-level working principle consists of estimating the
differences between the nominal (as predicted by the reference model) and
observed state transitions. Such differences are then compensated by
allocating a control authority proportional to the disturbance, effectively
driving the system to its reference behaviour. Applications of $\mathcal{L}1$
adaptive controllers span several types of aerial vehicles, from multi-rotors
to fixed wings [14, 15].

Mostly related to this work is the application of $\mathcal{L}1$ adaptive
control to quadcopters [16]. However, the performance of the classic
$\mathcal{L}1$ formulation degrades whenever the observed transitions differ
greatly from the (usually linear) reference model, which can happen due to
aerodynamic effects or large payloads. Therefore, recent work has combined
$\mathcal{L}1$ adaptive control with nonlinear online optimization [17, 18,
19]. While these methods achieved impressive results, they still require
explicit knowledge of (reference) system parameters, such as inertia, mass,
and motor characteristics. In addition, they generally require platform-
specific tuning to get the desired behaviour. Other approaches to adaptive
control on quadcopters include differential flatness [20] and nonlinear
dynamic inversion [21]. These methods have shown rapid adaptation to
aerodynamic effects and model uncertainties. However, they cannot cope with
large variations in the quadcopter’s dynamics, since the underlying
assumptions on locally linear disturbances are generally not fulfilled.

###  II-B Learning-Based Adaptive Control

Recent data-driven controllers have shown promising results for quadcopter
stabilization [22, 23], or waypoint tracking flight [24, 25]. As their classic
counterparts, learning-based controllers also allow for adaptation to
disturbances and model mismatches. One possibility to do so consists of
learning a model from the data and using the model to adapt the controller
[26, 27, 6, 28]. However, this has the limitation that the models are
difficult to carefully identify due to under-actuation of the platform and
sensing noise. This motivated model-free methods that, like ours, learn an
end-to-end adaptive policy [29]. Meta-learning has also been proposed to
augment the performance of the model-based controller for fast online
adaptation to wind [7] or suspended payloads [6]. These methods achieved
impressive results in disturbance rejection. However, they are still tailored
to a specific platform type, and lack an explicit mechanism for adaptation to
drastic changes in the drone’s model and actuation. Our work aims to fill this
gap and create a single control policy capable of flying vehicles with vastly
varying physical characteristics and under large disturbances.

##  III Learning a Universal Adaptive Drone Controller

We learn a universal controller to fly a quadcopter to a target position. This
controller takes a history of platform states and commands from a platform-
independent high level controller as inputs, and outputs individual motor
speed. Our resulting policy can control quadcopters with very different design
and hardware characteristics, and is robust to disturbances unseen at training
time, such as swing payloads and malfunctioned motors.

To achieve this, we use RMA [5]. However, we re-purpose the adaption module,
originally used to estimate parameters external to the robot, e.g. friction,
to estimate the robot’s intrinsic parameters (e.g. its mass, inertia, motor
constant, etc.) In the following, we recapitulate the method of RMA for
completeness. Our policy consists of a base policy $\pi$ and an adaptation
module $\phi$. The base policy $\pi$ takes the current state
$x_{t}\in\mathbb{R}^{23}$ and an intrinsics vector $z_{t}\in\mathbb{R}^{6}$ as
input and outputs the target motor speed $a_{t}\in\mathbb{R}^{4}$ for all
individual motors. The intrinsics vector $z_{t}$ allows the base policy to
adapt to variations in drone parameters, payloads, and disturbances such as
wind. Since we cannot directly measure $z_{t}$ in the real world, we instead
estimate it via the adaptation module $\phi$, which uses the discrepancy
between the commanded actions and the measured sensor readings to estimate it
online during deployment. More concretely,

$\displaystyle\hat{z_{t}}$
$\displaystyle=\phi\big{(}x_{t-k:t-1},a_{t-k:t-1}\big{)}$ (1) $\displaystyle
a_{t}$ $\displaystyle=\pi(x_{t},\hat{z_{t}})$ (2)

###  III-A Base Policy

We learn the base policy $\pi$ in simulation using model-free RL. The base
policy takes the current state $x_{t}$ and the intrinsics vector $z_{t}$,
which is a compressed version of the environment parameters $e_{t}$ containing
drone parameters, payload, etc. We use the network $\mu$ to compress $e_{t}$
to $z_{t}$. This gives us:

$\displaystyle z_{t}$ $\displaystyle=\mu(e_{t})$ (3) $\displaystyle a_{t}$
$\displaystyle=\pi(x_{t},z_{t})$ (4)

We implement $\mu$ and $\pi$ as MLPs and jointly train the base policy $\pi$
and the environmental factor encoder $\mu$ end to end using model-free
reinforcement learning. RL maximizes the following expected return of the
policy $\pi$:

$\displaystyle J(\pi)=\mathbb{E}_{\tau\sim
p(\tau|\pi)}\Bigg{[}\sum_{t=0}^{T-1}\gamma^{t}r_{t}\Bigg{]},$ (5)

where $\tau=\\{(x_{0},a_{0},r_{0}),(x_{1},a_{1},r_{1})...\\}$ is the
trajectory of the agent when executing the policy $\pi$, and $p(\tau|\pi)$
represents the likelihood of the trajectory under $\pi$.

#### RL Reward

The following reward function encourages the agent to hover at a goal position
and penalizes it for crash and oscillating motions. Let us denote the
acceleration in the $z$-axis i.e. the mass-normalized thrust as $\mathbf{c}$,
the angular velocity as $\bm{\omega}$, all in the quadcopter’s body frame. The
commanded mass-normalized thrust is defined as $\bm{c}_{des}$, commanded
angular velocity as $\bm{\omega}_{des}$, all given by the high-level
controller. We additionally define the motor speed as $\bm{m}$ and the
simulation step in the training as $\bm{\delta t}$. The reward at time $t$ is
defined as the sum of the following quantities:

  1. 1.

Angular Velocity Tracking : $-\|\bm{\omega}^{t}-\bm{\omega}_{des}^{t}\|$

  2. 2.

Z Acceleration Tracking: $-\|\mathbf{c}^{t}-\mathbf{c}_{des}^{t}\|$

  3. 3.

Output Command Smoothness: $-\|\bm{m}^{t}-\bm{m}^{t-1}\|$

  4. 4.

Survive: $\bm{\delta t}$

The scaling factor of each reward term is $0.01$, $0.02$, $0.0002$, $1$
respectively.

Parameters | Training Range | Testing Range
---|---|---
Mass (kg) | [0.142, 0.950] | [0.114, 1.140]
Arm length (m) | [0.046, 0.200] | [[0.037, 0.240]
Mass moment of inertia | [7.42e-5 , 5.60e-3] | [5.94e-5 , 6.72e-3]
around $x$, $y$ (kg$\cdot$m2)
Mass moment of inertia | [1.20e-4, 8.80e-3] | [9.60e-5, 1.06e-2]
around $z$ (kg$\cdot$m2)
Propeller constant $\kappa$ (m) | [0.0041, 0.0168] | [0.0033, 0.0201]
Payload (% of Mass) | [10, 50] | [5, 60]
Payload Location | [-10, 10] | [-10, 10]
from Center of Mass
(% of Arm length)
Motor Constant | [1.15e-7, 7.64e-6] | [9.16e-8, 9.17e-6]
Body drag coefficient | [0, 0.62] | [0, 0.74]
Max. motor speed (rad/s) | [707, 4895] | [566, 5874]
TABLE I:  Ranges of the drone and environmental parameters. Parameters without
units labelled are dimensionless quantities.

###  III-B Adaptation Module

During deployment, we do not have access to the vector $e_{t}$ and hence we
cannot compute the intrinsics vector $z$. Instead, we will use the sensor and
action history to directly estimate $z_{t}$ as proposed in [5]. We call this
module the adaptation module $\phi$, and we will train this in simulation
itself $a_{t}=\pi(x_{t},\hat{z_{t}})$. We can train this module in simulation
using supervised learning because we have access to both the ground truth
intrinsics $z_{t}$, and the sensor history and previous actions. We minimize
the mean squared error loss $\lVert z-\hat{z}\rVert^{2}$.

###  III-C Deployment

We directly deploy the base policy $\pi$ which uses the current state and the
intrinsics vector $\hat{z}$ predicted by the adaptation module $\phi$. We do
not calibrate or finetune our policy on any drones and use the same policy
without any modifications as the controller of the different drones under
different payload and conditions.

##  IV Experimental Setup

Simulation Environment. We use the Flightmare simulator [30] for training and
testing our control policies. We implement a custom high-level controller to
generates high-level commands at the level of body rates and collective
thrust. It is designed as a cascaded linear acceleration controller (with
desired acceleration mimicking a spring-mass-damper system with natural
frequency 2rad/s and damping ratio 0.7. The desired acceleration is then
converted to a desired total thrust and the desired thrust direction, and the
body rates are computed from this as proportional to the attitude error angle,
with a time constant of 0.2s. We define the task as hovering at a pre-defined
location. This high-level controller’s inputs are the platform’s state
(position, rotation, angular, and linear velocities) and the goal location.
The policy outputs individual motor speed commands, and we model the motors’
response using a first-order system. Each RL episode lasts for a maximum of 5s
of simulated time, with early termination if the quadcopter height drops below
2cm. The control frequency of the policy is $500$Hz, and the simulation step
is 2ms. We additionally implement an observation latency of 10ms.

Quadcopter and Environment Randomization. All our training and testing ranges
are listed in Table I. Of these, $e_{t}$ includes mass, arm length, propeller
torque to thrust ratio $\kappa$, motor constant, inertia (9dim), body drag
coefficients (3dim), maximum motor rotation and payload mass, which results in
an 18-dim vector. We randomize each of these parameters at the end of each
episode. In addition, at a randomly sampled time during each episode, the
parameters including mass, inertia, and the Center of Mass are also
randomized. The latter is used to mimic sudden variations in the quadcopter
parameters due to a sudden disturbance caused by a payload or wind.

Hardware Details. For all our real-world experiments we use two quadcopters,
which differ in mass by a factor of 4.5, and in size by a factor of 2.9. The
first one, which we name _largequad_ has a mass of 792g, a size of 16.6cm in
arm length, a thrust-to-weight ratio of 3.50, a diagonal inertia matrix of
[0.0047, 0.005, 0.0074]kg$\cdot$m2 (as expressed in the z-up body-fixed
frame), and a maximum motor speed of 943rad/s. The second one, _miniquad_ ,
has a mass of 177g, a size of 5.8cm in arm length, a thrust-to-weight ratio of
3.45, a diagonal inertia matrix of [92.7e-6, 92.7e-6, 158.57e-6]kg$\cdot$m2,
and a maximum motor speed of 3916rad/s. A motion capture system running at
$200$Hz provides estimates of the drone position and orientation, which is
followed by a Kalman filter to reduce noise and estimate linear velocity. An
onboard rate gyroscope measures the angular rotation of the robot, which is
low-pass filtered to reduce noise and remove outliers. The deployed policy
outputs motor speed commands, which are subsequently tracked by off-the-shelf
electronic speed controllers. We use as high-level a PD controller which takes
as input the goal position and outputs the mass normalized collective trust
and the body rates.

Network Architecture and Training Procedure. The base policy is a 3-layer
Multi-layer perceptron (MLP) with 256-dim hidden layers. This takes the drone
state and the vector of intrinsics as input to produce motor speeds. The
environment factor encoder is a 2-layer MLP with 128-dim hidden layers. The
policy and the value function share the same factor encoding layer. The
adaptation module projects the latest 400 state-action pairs into a 128-dim
representation. Then, a 3-layer 1-D CNN convolves the representation across
time. The flattened CNN output is linearly projected to estimate $z_{t}$. We
train the base policy and the environment encoder using PPO [31] for 100M
steps. We use the reward outlined in Section III. Policy training takes
approximately 2 hours on an ordinary desktop machine with 1 GPU. We finally
train the adaptation module with supervised learning by rolling out the
student policy. We train with the ADAM optimizer to minimize MSE loss. We run
the optimization process for 10M steps, training on data collected over the
last 1M steps. Training the adaptation module takes approximately 20 minutes.

##  V Results

Figure 3: Upper Left: Large and small quadcopters mounted with an inertia board. For the large quadcopter, we mount a wrench of 20.5cm and 140g. For the small quadcopter, a wrench of 14.5cm and 30g. Upper Right: Large and small quadcopters with a rigid payload. For the former, we add a load of 180g ($\approx 25\%$ of its weight). For the latter, we add a load of 50g, which corresponds to 35.7% of its weight. Lower Left: Large and small quadcopters with random pushes/pulls. Lower Right: Large and small quadcopters with a swing payload, which is the wrench in the inertia board test. Despite such large disturbances unknown to the control policy, our approach can always stabilize the quadcopter. Videos at https://dz298.github.io/universal-drone-controller/ | Vehicle | Method | Height | Ang. Vel. | Thrust | Success
---|---|---|---|---|---|---
| Err. (m) | Err. (rad/s) | Err. (m/s2) | Rate
Free Hover | small | PID | 0.06 | 0.51 | 0.57 | 100
LTF [32] | 0.09 | 0.78 | 0.78 | 80
Ours | 0.05 | 0.30 | 0.30 | 100
large | PID | 0.01 | 0.12 | 0.28 | 100
LTF [32] | 0.03 | 1.21 | 1.86 | 80
Ours | 0.03 | 0.19 | 0.36 | 100
Inertial Board | small | PID | - | - | - | 0
LTF [32] | - | - | - | 0
Ours | 0.08 | 1.14 | 1.09 | 100
large | PID | 0.31 | 0.37 | 1.46 | 100
LTF [32] | - | - | - | 0
Ours | 0.08 | 0.24 | 1.35 | 100
Payload | small | PID | 0.40 | 1.14 | 2.20 | 80
LTF [32] | 0.10 | 0.98 | 2.14 | 100
Ours | 0.05 | 0.88 | 0.51 | 100
large | PID | 0.19 | 1.14 | 1.21 | 100
LTF [32] | 0.06 | 1.01 | 4.2 | 100
| Ours | 0.04 | 0.34 | 1.19 | 100
TABLE II:  Real-World Results: We compare the performance of our controller to
two baselines: _LTF_ and a PID controller. The comparison is run on three
tasks for large and small quadcopters. Free Hover: hover at the goal position
without any disturbances. Inertia Board: hover at the goal position under an
unknown mass and inertia disturbance. Payload: hover at the goal position
under an unknown mass disturbance. Metrics are averaged over 5 experiments.

###  V-A Real World Deployment

We test our approach in the physical world and compare its performance to two
baselines: _LTF_ , which is a learning-based robust controller trained with an
error integration as one of inputs [32]; and a platform-specific PID
controller. The PID controller has access to the platform’s mass and inertia,
and it has been specifically tuned to the platform with in-flight tests. In
contrast, our approach has no knowledge whatsoever of the physical
characteristics of the system and requires no calibration or real-world fine
tuning. We compare the task of stabilization to a predefined set point without
any disturbances (free hover), or under an unknown mass and/or inertia
disturbance (Figure 3). We compare the three approaches under four metrics:
(i) the average height error to the goal point, and the average tracking error
of the (ii) angular velocity and (iii) mass normalized thrust of the high-
level controller’s commands, and the success rate. We define a failure if the
human operator had to intervene to avoid the quadcopter from crashing. The
results of these experiments are reported in table II. Our approach
significantly outperforms both baselines in all metrics under all
disturbances. In particular, our approach achieves a 100% success rate when
asymmetric disturbances are applied to the system, while the two baselines
both experience at least one total failure on either of the two platforms. The
PID baseline and our approach performs similarly in free hovering experiment,
with the PID controller slightly outperforming ours on _largequad_. The latter
difference in performance is justified, since the PID controller is
specifically tuned for each quadcopter, but ours does not have knowledge of
the system dynamics and hardware.

###  V-B Simulation Results

| Success | Height | Ang. Vel. | Thrust
---|---|---|---|---
| Rate | Err. (m) | Err. (rad/s) | Err. (m/s2)
Robust [33, 34] | 18% | 0.44 | 1.32 | 1.93
SysID [35] | 38% | 0.19 | 1.10 | 1.38
LTF [32] | 59% | 0.35 | 0.97 | 1.68
$\mathcal{L}1$ [17] | 59% | 0.17 | 0.92 | 1.60
Ours | 66% | 0.09 | 0.94 | 1.57
Expert (Phase I) | 69% | 0.09 | 0.91 | 1.29
TABLE III:  Simulation Testing Results: We compare our method with four baselines on the task of stabilization: _Robust_ , _SysID_ , _LTF_ , and $\mathcal{L}1$. We also list the results of our method in Phase I training, which has access to all ground-truth system parameters and can be regarded as the expert. The test ranges are defined in Table I. Metrics are averages over 100 experiments.  External Forces | Success | Height | Ang. Vel. | Thrust
---|---|---|---|---
Rate | Err. (m) | Err. (rad/s) | Err. (m/s2)
Robust [33, 34] | 5% | 0.39 | 2.05 | 1.51
SysID [35] | 2% | 0.22 | 1.38 | 1.23
LTF [32] | 32% | 0.23 | 1.53 | 0.99
$\mathcal{L}1$ [17] | 42% | 0.30 | 1.46 | 1.04
Ours | 49% | 0.09 | 1.40 | 0.85
Partially Failing Motors |  |  |
Robust [33, 34] | 1% | 0.53 | 1.94 | 1.59
SysID [35] | 14% | 0.33 | 1.25 | 1.23
LTF [32] | 21% | 0.38 | 1.57 | 1.35
$\mathcal{L}1$ [17] | 33% | 0.32 | 1.41 | 1.06
Ours | 38% | 0.26 | 1.36 | 1.01
TABLE IV:  Simulation Testing Results, Out-of-Distribution Disturbances: We
evaluate the performance of our method and all baselines on two types of
disturbances unseen at training time. External Forces: We apply a random force
of magnitude uniformly sampled between 0 and 50% of the weight and with
direction uniformly sampled on a cube. Partially Failing Motors: To simulate a
motor losing efficiency, we multiply the output of a randomly sampled motor’s
thrust force to a random number between 0 and 1. The duration of each
disturbance is random between the entire length of the episode (on and off
with 2% probability at every time stamp).

Finally, we compare our approach with a set of baselines in the simulation. We
select four baselines from prior work: _Robust_ , which consists of a policy
trained without access to environment factors or body parameters [34, 33];
_SysID_ , which directly predicts the ground-truth parameters $e_{t}$ [35]
instead of the low-dimensional intrinsics vector; _LTF_ , which essentially is
a robust policy with an error integration as additional inputs [32]; and
$\mathcal{L}1$, a model-based adaptive controller which estimates and
compensates the difference between the nominal and observed states to achieve
adaptive control [13, 4, 17]. We keep the same architecture and
hyperparameters as ours for all learning-based baselines. Similarly to real-
world experiments, we evaluate on the task of stabilization. The testing
ranges are listed in Table I. We rank the methods according to the success
rate, the average height error, and the tracking performance. At the beginning
of each experiment, the quadcopter is spawned with a random orientation, a
random position in x-y of [-1,1] and [0.5, 2.5] in z, and a random velocity on
each axis of [-1,1]. The experiment is considered successful if the end height
of the quadcopter is within 0.3m from the goal height. The results of these
experiments are reported in Table III. Given the very large amount of
quadcopter variations, the _Robust_ baseline trained without access to
environment parameters has the lowest success rate and largest tracking error.
This is because it is forced to learn a single conservative flying which can
fly all quadcopters under varying disturbances. Indeed, it either crashes to
the ground or flies outside the flying region. The baseline _LTF_ is _Robust_
with an additional observation of error integration. This additional input
helps it achieve higher success rate. However, similar to the _Robust_
baseline, _LTF_ fails in tracking the goal states. On the contrary, with
access to the environment parameters, the flight performance strongly
increases.

Figure 4: We analyze the change in behaviour of our policy as we incrementally
add total 290g payloads to our quadcopter. We plot all 6 components of the
intrinsics vector $\hat{z_{t}}$ predicted by the adaptation module. We see
that changes in intrinsics are strongly correlated with disturbances applied
to the quadcopter, indicating that the added payloads have been detected by
the adaptation module. When the payload is added, the quadcopter first sinks
and then recovers to the normal motion. The plotted components of the
intrinsics vector change in response to the disturbance, from which we know
the adaptation period takes around 2s.

Still, explicitly regressing the environment vector $e_{t}$, instead of its
low-dimensional representation $z_{t}$, scarifies the success rate since it is
trying to solve a harder estimation problem which is not necessary to solve
the task at hand. Since the tracking errors are computed only for successful
runs, the _SysID_ baseline achieves a slightly lower tracking error in one of
the metrics. $\mathcal{L}1$ is the baseline with the strongest performance in
our simulation experiments. However, it relies on the explicit knowledge of a
reference model which we choose as the median value of all parameters in the
testing range in Table I, while our method with other baselines has no prior
knowledge of the model. However, $\mathcal{L}1$ is not implemented in our
real-world experiments as a baseline because it is highly sensitive to latency
and requires significant engineering work. In contrast, our approach can
handle latency very well by simulating it during training. We also evaluate
the task of stabilization on held-out environmental disturbances, as
exemplified by random external forces and partial failure of a motor. The
results of these experiments are reported in Table IV. Our method outperforms
all baselines in both out-of-distribution disturbances cases. Compared to the
strongest baseline $\mathcal{L}1$, our method reduces the failure rate by up
to 1.14 times, the average height error by up to 67.7%, and the average
command tracking error by up to 22.4%.

###  V-C Adaptation Analysis

We analyze the estimated intrinsics vector $\hat{z_{t}}$ for adaptation on
incremental payloads. We incrementally add payloads in total of 290g to the
quadcopter in hovering. The quadcopter adapts successfully to payloads and
stabilizes itself at the target position. We plot all components of the
estimated intrinsics vector from the adaptation module during the experiment
in Figure 4. We find that whenever a payload is applied, there is a change in
all intrinsics components, indicating that the disturbance has been detected
by the adaptation module. Then the quadcopter recovers from the disturbance
and all components of the intrinsics vector converges to different values as
they are before the payload is applied. From the change of the plotted
components of intrinsics vector in response to the disturbance, we see that it
takes around 2s for our controller to detect the disturbance, estimate the
intrinsics vector, adapts from the disturbance, and finally stablizes. The
adaptation period of our controller is much faster than that of methods with
online estimation of ground truth model parameters such as [2], where it takes
10s to 15s in-flight to obtain an accurate model estimate.

##  VI Conclusion

In this work, we show how a single policy can control quadcopters of totally
different morphology, mass, and actuation. We successfully transfer a method
initially developed for legged locomotion in adapting terrains to quadcopters
in adapting a diverse set of quadcopter bodies and disturbances. Without any
additional tuning or modification, the single policy trained only in
simulation can be deployed zero-shot to real quadcopters of very different
design and hardware characteristics while showing rapid adaptation to unknown
disturbances at the same time. Adaptive control has achieved great success in
adaptation to model uncertainties and disturbances, but our work offers a
fresh view on the meaning of adaptation by completely bypassing the notion of
a _reference_ model that most adaptive control methods use. This feature
enables our method to adapt to a much wider range of robots and disturbances.

##  VII Acknowledgement

This work was supported by the DARPA Machine Common Sense program, Hong Kong
Center for Logistics Robotics, and the Graduate Division Block Grant of the
Dept. of Mechanical Engineering, UC Berkeley. The experimental testbed at the
HiPeRLab is the result of contributions of many people, a full list of which
can be found at hiperlab.berkeley.edu/members/. This work was also partially
supported by the ONR MURI award N00014-21-1-2801.

## References

  * [1] J. Svacha, J. Paulos, G. Loianno, and V. Kumar, “Imu-based inertia estimation for a quadrotor using newton-euler dynamics,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 3, pp. 3861–3867, 2020.
  * [2] V. Wüest, V. Kumar, and G. Loianno, “Online estimation of geometric and inertia parameters for multirotor aerial vehicles,” in _2019 International Conference on Robotics and Automation (ICRA)_. IEEE, 2019, pp. 1884–1890.
  * [3] M. Forgione and D. Piga, “Continuous-time system identification with neural networks: Model structures and fitting criteria,” _European Journal of Control_ , vol. 59, pp. 69–81, 2021.
  * [4] N. Hovakimyan and C. Cao, _L1 adaptive control theory: Guaranteed robustness with fast adaptation_. SIAM, 2010.
  * [5] A. Kumar, Z. Fu, D. Pathak, and J. Malik, “Rma: Rapid motor adaptation for legged robots,” _RSS: Robotics Science and Systems_ , 2021.
  * [6] S. Belkhale, R. Li, G. Kahn, R. McAllister, R. Calandra, and S. Levine, “Model-based meta-reinforcement learning for flight with suspended payloads,” _IEEE Robotics and Automation Letters_ , vol. 6, no. 2, pp. 1471–1478, 2021.
  * [7] M. O’Connell, G. Shi, X. Shi, K. Azizzadenesheli, A. Anandkumar, Y. Yue, and S.-J. Chung, “Neural-fly enables rapid learning for agile flight in strong winds,” _Science Robotics_ , vol. 7, no. 66, p. eabm6597, 2022.
  * [8] “Betaflight: Open source flight controller firmware,” https://betaflight.com/.
  * [9] “Open source autopilot for drones: Px4 autopilot,” https://px4.io/.
  * [10] I. M. Mareels, B. D. Anderson, R. R. Bitmead, M. Bodson, and S. S. Sastry, “Revisiting the mit rule for adaptive control,” _IFAC Proceedings Volumes_ , vol. 20, no. 2, pp. 161–166, 1987.
  * [11] K. J. Åström and B. Wittenmark, _Adaptive control_. Courier Corporation, 2013.
  * [12] E. Lavretsky and K. A. Wise, “Robust adaptive control,” in _Robust and adaptive control_. Springer, 2013, pp. 317–353.
  * [13] C. Cao and N. Hovakimyan, “Design and analysis of a novel l1 adaptive control architecture with guaranteed transient performance,” _IEEE Transactions on Automatic Control_ , vol. 53, no. 2, pp. 586–591, 2008.
  * [14] S. Mallikarjunan, B. Nesbitt, E. Kharisov, E. Xargay, N. Hovakimyan, and C. Cao, “L1 adaptive controller for attitude control of multirotors,” in _AIAA guidance, navigation, and control conference_ , 2012, p. 4831.
  * [15] I. Gregory, C. Cao, E. Xargay, N. Hovakimyan, and X. Zou, “L1 adaptive control design for nasa airstar flight test vehicle,” in _AIAA guidance, navigation, and control conference_ , 2009, p. 5738.
  * [16] M. Schreier, “Modeling and adaptive control of a quadrotor,” in _2012 IEEE international conference on mechatronics and automation_. IEEE, 2012, pp. 383–390.
  * [17] D. Hanover, P. Foehn, S. Sun, E. Kaufmann, and D. Scaramuzza, “Performance, precision, and payloads: Adaptive nonlinear mpc for quadrotors,” _IEEE Robotics and Automation Letters_ , vol. 7, no. 2, pp. 690–697, 2021.
  * [18] J. Pravitra, K. A. Ackerman, C. Cao, N. Hovakimyan, and E. A. Theodorou, “L1-adaptive mppi architecture for robust and agile control of multirotors,” in _2020 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2020, pp. 7661–7666.
  * [19] K. Pereida and A. P. Schoellig, “Adaptive model predictive control for high-accuracy trajectory tracking in changing conditions,” in _2018 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2018, pp. 7831–7837.
  * [20] M. Faessler, A. Franchi, and D. Scaramuzza, “Differential flatness of quadrotor dynamics subject to rotor drag for accurate tracking of high-speed trajectories,” _IEEE Robotics and Automation Letters_ , vol. 3, no. 2, pp. 620–626, 2017.
  * [21] E. J. Smeur, Q. Chu, and G. C. De Croon, “Adaptive incremental nonlinear dynamic inversion for attitude control of micro air vehicles,” _Journal of Guidance, Control, and Dynamics_ , vol. 39, no. 3, pp. 450–461, 2016.
  * [22] J. Hwangbo, I. Sa, R. Siegwart, and M. Hutter, “Control of a quadrotor with reinforcement learning,” _IEEE Robotics and Automation Letters_ , vol. 2, no. 4, pp. 2096–2103, 2017.
  * [23] W. Koch, R. Mancuso, R. West, and A. Bestavros, “Reinforcement learning for uav attitude control,” _ACM Transactions on Cyber-Physical Systems_ , vol. 3, no. 2, pp. 1–21, 2019.
  * [24] Y. Song, M. Steinweg, E. Kaufmann, and D. Scaramuzza, “Autonomous drone racing with deep reinforcement learning,” in _2021 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2021, pp. 1205–1212.
  * [25] E. Kaufmann, L. Bauersfeld, and D. Scaramuzza, “A benchmark comparison of learned control policies for agile quadrotor flight,” _arXiv preprint arXiv:2202.10796_ , 2022.
  * [26] N. O. Lambert, D. S. Drew, J. Yaconelli, S. Levine, R. Calandra, and K. S. Pister, “Low-level control of a quadrotor with deep model-based reinforcement learning,” _IEEE Robotics and Automation Letters_ , vol. 4, no. 4, pp. 4224–4230, 2019.
  * [27] G. Shi, X. Shi, M. O’Connell, R. Yu, K. Azizzadenesheli, A. Anandkumar, Y. Yue, and S.-J. Chung, “Neural lander: Stable drone landing control using learned dynamics,” in _2019 International Conference on Robotics and Automation (ICRA)_. IEEE, 2019, pp. 9784–9790.
  * [28] G. Torrente, E. Kaufmann, P. Föhn, and D. Scaramuzza, “Data-driven mpc for quadrotors,” _IEEE Robotics and Automation Letters_ , vol. 6, no. 2, pp. 3769–3776, 2021.
  * [29] C.-H. Pi, W.-Y. Ye, and S. Cheng, “Robust quadrotor control through reinforcement learning with disturbance compensation,” _Applied Sciences_ , vol. 11, no. 7, p. 3257, 2021.
  * [30] Y. Song, S. Naji, E. Kaufmann, A. Loquercio, and D. Scaramuzza, “Flightmare: A flexible quadrotor simulator,” in _Proceedings of the 2020 Conference on Robot Learning_ , 2021, pp. 1147–1157.
  * [31] J. Schulman, F. Wolski, P. Dhariwal, A. Radford, and O. Klimov, “Proximal policy optimization algorithms,” _arXiv preprint arXiv:1707.06347_ , 2017\.
  * [32] J. Xu, T. Du, M. Foshey, B. Li, B. Zhu, A. Schulz, and W. Matusik, “Learning to fly: computational controller design for hybrid uavs with reinforcement learning,” _ACM Transactions on Graphics (TOG)_ , vol. 38, no. 4, pp. 1–12, 2019.
  * [33] X. B. Peng, M. Andrychowicz, W. Zaremba, and P. Abbeel, “Sim-to-real transfer of robotic control with dynamics randomization,” in _2018 IEEE international conference on robotics and automation (ICRA)_. IEEE, 2018, pp. 3803–3810.
  * [34] J. Tobin, R. Fong, A. Ray, J. Schneider, W. Zaremba, and P. Abbeel, “Domain randomization for transferring deep neural networks from simulation to the real world,” in _2017 IEEE/RSJ international conference on intelligent robots and systems (IROS)_. IEEE, 2017, pp. 23–30.
  * [35] W. Yu, J. Tan, C. K. Liu, and G. Turk, “Preparing for the unknown: Learning a universal policy with online system identification,” in _Robotics: Science and Systems XIII, Massachusetts Institute of Technology, Cambridge, Massachusetts, USA, July 12-16, 2017_ , N. M. Amato, S. S. Srinivasa, N. Ayanian, and S. Kuindersma, Eds., 2017. [Online]. Available: http://www.roboticsproceedings.org/rss13/p48.html


# Discovery of the first anti-glitch event in the rotation-powered pulsar PSR
B0540-69

Youli Tuo Institut für Astronomie und Astrophysik, Kepler Center for Astro and
Particle Physics, Eberhard Karls Universität Tübingen, Sand 1, 72076 Tübingen,
Germany  Muhammed Miraç Serim Institut für Astronomie und Astrophysik, Kepler
Center for Astro and Particle Physics, Eberhard Karls Universität Tübingen,
Sand 1, 72076 Tübingen, Germany  Marco Antonelli CNRS/in2p3, Laboratoire de
Physique Corpusculaire, 14050 Caen, France  Lorenzo Ducci Institut für
Astronomie und Astrophysik, Kepler Center for Astro and Particle Physics,
Eberhard Karls Universität Tübingen, Sand 1, 72076 Tübingen, Germany  ISDC
Data Center for Astrophysics, Université de Genève, 16 chemin d’Écogia, 1290
Versoix, Switzerland  Armin Vahdat Institut für Astronomie und Astrophysik,
Kepler Center for Astro and Particle Physics, Eberhard Karls Universität
Tübingen, Sand 1, 72076 Tübingen, Germany  Mingyu Ge Key Laboratory of
Particle Astrophysics, Institute of High Energy Physics, Chinese Academy of
Sciences, 19B Yuquan Road, Beijing 100049, China  University of Chinese
Academy of Sciences, Chinese Academy of Sciences, Beijing 100049, People’s
Republic of China  Andrea Santangelo Institut für Astronomie und Astrophysik,
Kepler Center for Astro and Particle Physics, Eberhard Karls Universität
Tübingen, Sand 1, 72076 Tübingen, Germany  Fei Xie Guangxi Key Laboratory for
Relativistic Astrophysics, School of Physical Science and Technology, Guangxi
University, Nanning 530004, China

###### Abstract

Using data from the Neutron star Interior Composition ExploreR (NICER)
observatory, we identified a permanent spin frequency decrease of
$\Delta\nu=-(1.04\pm 0.07)\times 10^{-7}\,\mathrm{Hz}$ around MJD 60132 in the
rotation-powered pulsar PSR B0540-69, which exhibits a periodic signal at a
frequency of $\nu\sim 19.6\,\mathrm{Hz}$. This points to an anti-glitch event,
a sudden decrease of the pulsar’s rotational frequency without any major
alteration in the pulse profile or any significant increase of the pulsed
flux. Additionally, no burst activity was observed in association with the
anti-glitch. To date, observations of the few known anti-glitches have been
made in magnetars or accreting pulsars. This is the first anti-glitch detected
in a rotation-powered pulsar. Given its radiatively quiet nature, this anti-
glitch is possibly of internal origin. Therefore, we tentatively frame this
event within a proposed mechanism for anti-glitches where the partial
‘evaporation’ of the superfluid component leads to an increase of the normal
component’s moment of inertia and a decrease of the superfluid one.

Neutron stars(1108) – Pulsars(1306) – Rotation powered pulsars(1408) – Compact
objects(288)

††software:  astropy (Robitaille et al., 2013), TEMPO2 (Hobbs et al., 2006),
Stingray (Huppenkothen et al., 2019), TAT-pulsar (Tuo et al., 2022).

##  1 Introduction

Isolated pulsars are neutron stars (NSs) characterized by stable rotational
dynamics, with periods that span from milliseconds to several tens of seconds
(Manchester et al., 2016; Enoto et al., 2019). They exhibit long-term spin-
down due to angular momentum loss via electromagnetic radiation, together with
possible contributions from relativistic particle outflows which form the
pulsar wind nebula, and the emission of gravitational waves (Hobbs et al.,
2010; Aasi et al., 2015). Rotation-powered pulsars (RPPs) are a subset of
isolated pulsars characterized by emissions primarily driven by magnetic
braking. In contrast, another subset of isolated pulsars, magnetars, power
their emission largely from the decay of their strong magnetic fields (e.g.,
Kaspi & Beloborodov, 2017; Mereghetti et al., 2015). However, the distinction
between RPPs and magnetars is blurred because some peculiar high B-field RPPs
– objects with intermediate properties between the two classes – are known (Ng
& Kaspi, 2011; Harding, 2013; Borghese & Esposito, 2023).

The spin-down of isolated neutron stars is not completely regular and
predictable, and two kinds of timing instabilities are often observed:
glitches and timing noise (D’Alessandro, 1996; Çerri-Serim et al., 2019).
Glitches, commonly observed in RPPs, are sudden increases in the rotation
frequency of a pulsar often followed by an exponential relaxation (see Haskell
& Melatos, 2015; Antonelli et al., 2022; Antonopoulou et al., 2022; Zhou et
al., 2022, for recent reviews on the glitch phenomenon). Conversely, there are
few observed events where the overall glitch contribution to the pulsar’s
rotation frequency is negative, referred to as ‘anti-glitch’ or ‘spin-down
glitch’.

So far, most observed glitches are spin-ups that are instantaneous to the
timing data’s accuracy, interrupting the otherwise regular spin-down of a
canonical pulsar or a magnetar (e.g., Espinoza et al., 2011; Yu et al., 2013;
Serim et al., 2017; Basu et al., 2022). On the other hand, anti-glitch events
have been observed only in magnetars or in accreting pulsars. Archibald et al.
(2013) discovered a sudden decrease in the rotational frequency in magnetar 1E
2259+586. Moreover, in binary pulsars where accretion from the companion star
primarily governs the long-term spin evolution, anti-glitch events have been
observed in NGC 300 ULX-1 (Ray et al., 2019). The same internal mechanism
invoked for glitches observed in spinning down pulsars – the angular momentum
transfer from the internal superfluid component to the normal one (Anderson &
Itoh, 1975) – can result in an ‘anti-glitch’ in accreting neutron stars that
are spinning up (Ducci et al., 2015; Howitt & Melatos, 2022; Antonelli et al.,
2022). However, even after the seminal work of Anderson & Itoh (1975), it is
not obvious which internal mechanism could lead to an anti-glitch in spinning
down RPPs. Apart from this theoretical difficulty, the absence of radiatively
quiet anti-glitches in RPPs sparks a debate about whether RPPs can abruptly
spin down without a change in their external braking torque. This work
presents the first anti-glitch in an RPP, the pulsar PSR B0540-69.

The young rotation-powered pulsar PSR B0540-69 (also known as PSR J0540-6919)
was discovered in the early 1980s by the Einstein X-ray Observatory (Seward et
al., 1984). It is located in the Large Magellanic Cloud, surrounded by a
pulsar wind nebula, and has a characteristic age of $\sim$1700 years (Petre et
al., 2007). With a spin period of $\sim$50, PSR B0540-69 shows typical
characteristics of young pulsars, particularly in its timing behaviour (it is
sometimes referred to as one of the ‘Crab’s twins’). In 2011, PSR B0540-69
experienced a sudden change in spin-down rate (Marshall et al., 2015),
followed by the enhancement of the flux of pulsar wind nebula by 32% (Ge et
al., 2019). The braking index evolution of PSR B0540-69 has been extensively
studied (Marshall et al., 2016; Wang et al., 2020), and several spin-up
glitches were discovered in Rossi X-ray Timing Explorer data (Ferdman et al.,
2015).

Our paper is structured as follows. Section 2 details the X-ray observations
and data reduction. Section 3.1 presents the timing analysis of PSR B0540-69,
including the glitch model solution. The analysis of flux variations is
presented in Section 3.2. Finally, in Section 4, we discuss and tentatively
interpret our findings in the light of the anti-glitch models of Yim et al.
(2024), Kantor et al. (2016) and Kantor & Gusakov (2014).

##  2 Observation and data reduction

PSR B0540-69 was observed using NICER. NICER’s X-ray Timing Instrument (XTI)
is equipped with 56 X-ray ‘concentrator’ optics (XRC) paired with silicon
drift detectors (SDD), each accompanied by a Focal Plane Module to collect
X-ray photons (Gendreau et al., 2016). NICER/XTI offers excellent timing
stability, maintaining a root-mean-square (RMS) of 100 nanoseconds relative to
Universal Time (Enoto et al., 2021). This precision is particularly
advantageous for pulsar timing monitoring and other timing analyses.

Our study encompasses the entire publicly available NICER/XTI dataset for PSR
B0540-69, spanning from March 14, 2018, to December 25, 2023 comprising a
total of 145 observations. We perform the standard Level 2 data generation
processes by using the NICER pipeline tool nicerl2 implemented in HEADAS
(v6.31.1). Two instrumental parameters are applied to exclude the electron
precipitation-type flare: an overshoot rate of ¡20 counts/s, and a Cutoff
Rigidity (COR_SAX) of ¿1.5 GeV. The barycentric correction for each photon is
applied using the barycor tool implemented in HEADAS. The solar ephemeris
JPL.DE430 is utilized. The precise celestial coordinates of PSR B0540-69, with
a right ascension of $05^{\mathrm{h}}40^{\mathrm{m}}10.84^{\mathrm{s}}$ and a
declination of $-69^{\circ}19^{\prime}54.2^{\prime\prime}$ according to SIMBAD
Astronomical Database 111https://simbad.u-strasbg.fr/simbad/ is employed
(Brown et al., 2021).

Following the recommended GTI filtering criteria, the good time interval (GTI)
was empty for some observations near the glitch epoch (ObsID: 6020010132,
6020010133, 6020010134, and 6020010139). Therefore, we used the Level 1 data
for their timing analysis. This approach ensured maximal data utilization, as
data outside the valid GTIs, despite their substantial background and particle
events, do not significantly impact the timing behavior and pulse profile
shape due to their random distribution (see Section 3.1 for details). Our
analysis includes photons in the energy range of 0.5–12 keV for the timing and
pulsed flux measurements.

##  3 Results

###  3.1 Timing

We performed the timing analysis for the whole public NICER archive. Here, we
only present data around the anti-glitch epoch, from MJD 59901 to MJD 60303.
The time of arrival (ToA) for each NICER observation is obtained via a phase-
coherent analysis. As usual, the ToAs and their uncertainties are calculated
by cross-correlating each pulse profile with a long-term cumulative profile
template (Taylor, 1992; Huppenkothen et al., 2019). Before the glitch epoch,
the long-term evolution of the ToAs in the phase domain can be approximated by
a truncated Taylor series,

$\Phi(t)=\Phi_{0}+\nu(t-t_{0})+\dfrac{\dot{\nu}}{2}(t-t_{0})^{2}+\dfrac{\ddot{\nu}}{6}(t-t_{0})^{3},$
(1)

where $\nu$, $\dot{\nu}$, and $\ddot{\nu}$ are the spin frequency, frequency
derivative, and second derivative of frequency at a reference time $t_{0}$
before the glitch, respectively. The timing model before the glitch is
sufficient to describe the spin-down trend of the pulsar, with an RMS residual
of 737.259. The values for each parameter are listed in Table 1.

As shown in the top panel of Figure 1, at around MJD 60132 the ToAs exhibited
a sudden deviation to the timing model, which indicates a glitch event. To fit
the ToAs, the assumed model for the post-glitch phase residues is (cf.
equation A7)

$\begin{split}&\Delta\Phi(t)=\Delta\nu\,\delta
t+\frac{\Delta\dot{\nu}}{2}\delta t^{2}+\Delta\nu_{d}\,\tau\left[1-e^{-\delta
t/\tau}\right]\\\ &\delta t=t-t_{g}>0\,,\end{split}$ (2)

where $t_{g}$ is the time when the glitch/anti-glitch occurs, $\Delta\nu$ and
$\Delta\dot{\nu}$ are the permanent changes in frequency and frequency
derivative. The above expression is valid for $\delta t\geq 0$, while for
$\delta t<0$ we assume that the spin-down is well described by the phase model
in (1). An exponential post-glitch recovery term is often required in glitch
modelling, where $\tau$ is the recovery time scale when a transient frequency
increment $\Delta\nu_{d}$ decays exponentially to zero (see, e.g., Wong et
al., 2001; Wang et al., 2012; Yu et al., 2013).

We use Markov Chain Monte Carlo (MCMC) techniques, as implemented in the emcee
package (Foreman-Mackey et al., 2013), to fit ToAs and determine the optimal
parameters for the timing model in equation (2). Details on the priors are
given in Appendix A.3. We choose a Gaussian likelihood – see, e.g., equation
(15) of Montoli et al. (2020) for the explicit likelihood’s formula – and
employ 16 walkers for the MCMC. We find that 320,000 steps for each walker are
enough for the statistical convergence of our Bayesian fit of (1) and (2). The
posterior probability distributions and the marginal distributions for each
parameter are shown in Figure 2.

Unfortunately, our Bayesian analysis indicates that the exponential terms in
equation (2) cannot be reliably constrained. This ambiguity arises from the
circumstance that the decay timescale is extremely long, leading to a
reduction to a linear term that contributes to the permanent increment in
frequency, or a case of rapid recovery over a short timescale, where
$\tau=1.06_{-0.79}^{+1.22}\,\mathrm{day}$). The F-test indicates that the
inclusion of the frequency derivative increment is favoured over a model that
incorporates only a single term for frequency increment. This conclusion is
supported by an F-statistic value of 17.6146 and a probability of $9.4\times
10^{-5}$.

The inferred parameters of (2), along with the timing ephemeris for the spin-
down model in (1), are reported in Table 1. The uncertainty associated with
each parameter is represented by the 68% credible interval derived from the
posterior distribution. The RMS of timing residual after the overall fit is
$\SI{829.8}{\micro\second}$. The negative value of
$\Delta\nu=-1.04_{-0.07}^{+0.07}\times 10^{-7}\,\mathrm{Hz}$ indicates the
anti-glitch nature of this event, see the discussion in Section 4.1.

Parameters | Values
---|---
R.A. (J2000) | $05^{\mathrm{h}}40^{\mathrm{m}}10.84^{\mathrm{s}}$
Decl. (J2000) | $-69^{\circ}19^{\prime}54.2^{\prime\prime}$
$\nu$ (Hz) | 19.636085141(2)
$\dot{\nu}$ ($\times 10^{-10}\,\mathrm{Hz}\cdot\mathrm{s}^{-1}$) | $-2.521868(3)$
$\ddot{\nu}$ ($\times 10^{-21}\,\mathrm{Hz}\cdot\mathrm{s}^{-2}$) | $4.6(1)$
Epoch (MJD) | 60041.21699
Valid Range (MJD) | 59901–60318
Ephemeris | JPL-DE430
$\Delta\nu$ ($\times 10^{-7}\,\mathrm{Hz}$) | $-1.042_{-0.074}^{+0.076}$
$\Delta\dot{\nu}$ ($\times 10^{-15}\,\mathrm{Hz}\cdot\mathrm{s}^{-1}$) | $-7.4_{-6.1}^{+6.2}$
$t_{g}$ (MJD) | $60132.158_{-4.633}^{+5.224}$
$|\Delta\nu/\nu|$ ($\times 10^{-9}$) | $5.306_{-0.037}^{+0.038}$
RMS residual () | 829.8
Table 1:  The inferred spin and glitch parameters. Errors refer to the 68%
confidence interval of the posterior distributions.  Figure 1: Timing results
and the pulsed flux evolution of the radiatively quiet anti-glitch of PSR
B0540-69. The inferred epoch $t_{g}=\mathrm{MJD}\,60132$ is marked with a blue
dashed line, and the shaded band around it represents its uncertainty. (a) The
phase evolution, based on the pre-glitch polynomial ephemeris. The pulsar’s
phase exhibits a linear deviation from the pre-glitch solution. (b) The
residuals are fitted with the glitch model, which incorporates the terms for
permanent frequency and frequency derivative increments. (c) The pulsed flux
variation around the glitch based on NICER observations in the 0.5–12 keV
range. Here, black circles represent individual observations, while red
squares denote averages compiled over time.  Figure 2:  Corner plot for the
posterior probability distributions of the glitch model parameters. In the
two-dimensional projections, the contours for each pair of parameters denote
the 1, 2, and 3-$\sigma$ credible intervals. The marginal distributions for
each parameter are shown along the diagonal: the median values of each
parameter are shown, alongside the 68% credible intervals. The labels
$\delta\nu$ and $\delta\dot{\nu}$ indicate that we shift the posterior
distributions for $\nu$ and $\dot{\nu}$ by their median posterior values,
which are exactly the ones reported in Table 1 for $\nu$ and $\dot{\nu}$.
Figure 3:  Cumulative pulse profiles before and after the glitch. (a) The pre-
glitch pulse profile is represented by the blue line, while the red line
indicates the post-glitch pulse profile. The grey area highlights the off-
pulse range chosen as the background for normalization, within which each
pulse profile has been background subtracted and normalized by dividing by the
peak counts. (b) The residuals between the pre- and post-glitch pulse
profiles. The error bars show 1 $\sigma$ uncertainties.

###  3.2 Flux

Looking for burst events or flux variations associated with the anti-glitch is
an important clue for determining its, possibly internal, nature. To this end,
we have examined the pulsed flux, burst activities, and variations in the
pulse profile properties of PSR B0540-69.

Our analysis of pulsed flux relies on the NICER observations in the energy
band 0.5-12 keV. The background, identified as the un-pulsed phase, is
selected for each pulse profile and is shown as the gray area in Figure 3. The
pulsed fluxes are normalized by the total number of Focal Plane Modules (FPMs)
active during each NICER observation. These normalized results are displayed
in the bottom panel of Figure 1. The red squares in the bottom panel represent
the averaged pulsed flux over time. We find no significant variation in the
pulsed flux observed during the epoch of the anti-glitch event. On the other
hand, we report the cumulative pulse profiles observed before and after the
glitch in Figure 3. Again, we find no significant changes in the pulse
profiles before and after the anti-glitch.

We also employ the burst searching techniques described in Cai et al. (2021)
for NICER light curves across several time scales, 1, 0.1, 1, and 10. The
energy band focused on is 2–8 keV to limit the effects of possible particle
background in the high energy band and the visible-light effect in the low
energy band. This method is used to scrutinize NICER observations to identify
X-ray bursts associated with the glitch event. We employ a time window to
examine each time bin of the light curve, aiming to identify outliers
exceeding background fluctuations with $3\,\sigma$ significance. The
background fluctuations are based on Poisson statistics, determined by
averaging the flux level in the intervals before and after the intervals of
interest. There is no significant burst event discovered in the filtered NICER
light curve. We also search the gamma-ray burst catalogue in the General
Coordinates Network Circulars222https://gcn.nasa.gov/circulars. No gamma-ray
bursts were associated with PSR B0540-69 during the period of glitch activity.
The findings suggest that the anti-glitch we discovered in PSR B0540-69
exhibits a radiatively quiet nature.

To summarize, we classify this event as an asymptotic anti-glitch,
characterized by a negative post-glitch permanent frequency increment (see
Figure 4). The glitch is radiatively quiet, exhibiting no significant
variations in the shape of the pulse profile or pulsed flux, and no burst
activities were observed in association with the glitch.

##  4 discussion

Glitches typically exhibit fractional amplitudes of $\Delta\nu/\nu\gtrsim
10^{-9}$, coinciding with a decrease in $\dot{\nu}$. On the other hand, RPPs
occasionally undergo ‘micro-glitches’ with smaller amplitudes
($\Delta\nu/\nu\lesssim 10^{-10}$, Chukwude & Urama 2010; Zhou et al. 2022),
manifesting as either positive or negative shifts in $\dot{\nu}$, albeit often
evading detection (Vahdat et al., 2024). This could contribute to the observed
rarity of anti-glitches in RPP.

Seven magnetars have been found exhibiting an anti-glitch or a spin-up glitch
followed by a rapid spin-down resulting in a net spin-down effect: 1E 2259+586
(Younes et al., 2020; Archibald et al., 2013; Içdem et al., 2012), 4U 0142+61
(Archibald et al., 2017; Gavriil et al., 2011), 1E 1841-045 (Dib & Kaspi,
2014), SGR J1935+2154 (Ge et al., 2024; Younes et al., 2023), PSR J1846-0258
(Livingstone et al., 2010; Şaşmaz Muş et al., 2014), SGR J1900+14 (Woods et
al., 1999), and PSR J1119-6127 (Dai et al., 2018; Archibald et al., 2017). Due
to the sparse observations around the glitch epoch, it is crucial to determine
whether the glitch observed in PSR B0540-69 represents a true anti-glitch, or
if the initial spin-up process was missed, leading the pulsar to over-recover
to a net spin-down state, as in Figure 4. Thus, we extend our analysis further
by applying a limiting case of the phenomenological model used by Montoli et
al. (2020) to our source, aiming to constrain the timing solution.

###  4.1 Glitch or anti-glitch?

To assess whether the glitch of PSR B0540-69 is a regular spin-up glitch or an
anti-glitch, we adopt a phenomenological model for the residuals based on a
reduction of the one used for Vela’s 2016 glitch (Montoli et al., 2020), see
Appendix A. For $\delta t>0$ the phase residual model can be written as (cf.
equation A4)

$\displaystyle\Delta\Phi(t)=\Delta\nu_{\infty}\,\delta
t+\frac{\Delta\nu_{0}-\Delta\nu_{\infty}}{\lambda}\left(1-e^{-\delta
t\,\lambda}\right)$ (3)

where $\Delta\nu_{\infty}$ is the asymptotic frequency offset and
$\Delta\nu_{0}$ is the immediate post-glitch frequency jump. The values of
$\Delta\nu_{\infty}$ and $\Delta\nu_{0}$ define different cases. For example,
for spin-up glitches we can have that:

  1. 1.

$0<\Delta\nu_{0}<\Delta\nu_{\infty}$: A regular spin-up glitch with a fast
initial jump and a following delayed spin-up, similar to what observed in the
Crab (Shaw et al., 2018; Ge et al., 2020).

  2. 2.

$0<\Delta\nu_{\infty}<\Delta\nu_{0}$: An overshooting glitch (Ashton et al.,
2019; Pizzochero et al., 2020). If the fast initial transient is not resolved,
this resembles a regular spin-up glitch with a certain amount of recovery.

The anti-glitch case is analogous, just the sign of the amplitudes is
reversed, see Figure 4. Therefore, we now explore the parameter space for
$\Delta\nu_{\infty}$ and $\Delta\nu_{0}$ in the phenomenological model,
allowing them to be negative, to aptly describe the anti-glitch behavior. To
this end, a similar MCMC approach is adopted, with the priors for each
parameter detailed in Appendix A.3. We find that $\Delta\nu_{0}$ – the post-
glitch frequency jump extrapolated at $t=t_{g}$ – is broadly consistent with
zero, $\Delta\nu_{0}=4.2_{-92.6}^{+115.4}\times 10^{-9}\,\mathrm{Hz}$, and it
is accompanied by a recovery timescale of $5.7_{-1.9}^{+5.0}\,\mathrm{days}$.
Regarding the asymptotic frequency offset, the analysis points at an anti-
glitch with $\Delta\nu_{\infty}=-1.08_{-0.09}^{+0.07}\times
10^{-7}\,\mathrm{Hz}$. Given the large uncertainty on $\Delta\nu_{0}$, the
Bayesian analysis suggests that this event is an _asymptotic_ anti-glitch,
which could be preceded by a quickly recovered spin-up glitch or even by an
undershoot, as shown in Figure 4. Unfortunately, continuous observation around
the glitch epoch $t_{g}$ was unavailable, so we can not distinguish between
these possibilities. In any case, it is worth noting that the recovery
timescale for this event is shorter than that of the majority of other
glitches observed in RPPs, at least those characterised by a single
exponential recovery timescale (Yu et al., 2013). There are only a few cases
in the Crab Pulsar where the recovery timescales are as short as a few days
(Wong et al., 2001), and the post-glitch frequency does not fully recover to
the pre-glitch solution. However, in our case, the small fractional change of
$\Delta\dot{\nu}/\dot{\nu}\sim 2.7\times 10^{-5}$ is not sufficient to resolve
an over-recovery scenario for the spin-up glitch in such a short recovery
timescale. In conclusion, the likelihood that the glitch event observed in PSR
B0540-69 represents a spin-up glitch appears to be marginal.

###  4.2 Anti-glitch model

Our timing analysis of PSR B0540-69 reveals a radiatively quiet anti-glitch in
an RPP. To our knowledge, this is the first study to report such an
occurrence.

There are several proposed models for the anti-glitches in magnetars, but are
not suitable for RPPs. For example, García & Ranea-Sandoval (2015) and
Mastrano et al. (2015) interpreted anti-glitches as a result of the decay of
its internal toroidal magnetic field component, which transforms a stable
prolate configuration into an unstable one. Unlike magnetars, RPPs do not
possess sufficiently strong magnetic fields to form a prolate configuration;
instead, the oblate configuration of RPPs (Baym & Pines, 1971) tends towards a
more spherical configuration. Thus, when a glitch occurs in RPPs, the moment
of inertia decreases, leading to a spin-up glitch (Ruderman, 1969; Baym et
al., 1969).

Tong (2014) posited that the particle wind contributes to the observed flux
enhancement in 1E 2259+586 and exerts a net spin-down torque on the NS.
Alternatively, Huang & Geng (2014) speculate that anti-glitches may result
from collisions of a pulsar with a small body, where gravitational potential
energy is released either in a short hard X-ray burst or in a soft X-ray
afterglow. However, our results on PSR B0540-69 do not indicate any
significant radiative change that can be associated with the anti-glitch
event.

Recently, Yim et al. (2024) established a toy model based on ejecta from the
magnetosphere that can lead to glitch or anti-glitch events. According to this
model, some mass $M_{0}$ can be expelled from the magnetic pole but may get
trapped as it travels outside through the magnetosphere. According to their
analysis (refer to Figures 3 and 4 therein), the ejecta results in an anti-
glitch depending on two conditions:

the magnetic inclination angle from which ejecta are expelled must be large,
and the ejecta must be expelled to a significant distance, at least $10R_{0}$,
where $R_{0}$ is the NS radius. For an NS like PSR B0540-69, the co-rotation
radius $R_{\mathrm{co}}$ is

$\frac{R_{\mathrm{co}}}{R_{0}}\approx
170\left[\frac{M}{1.4\,M_{\odot}}\right]^{\frac{1}{3}}\\!\left[\frac{\nu}{1\,\mathrm{Hz}}\right]^{-\frac{2}{3}}\\!\left[\frac{R_{0}}{10\,\mathrm{km}}\right]^{-1}\\!\\!\\!\approx
23,$ (4)

which is at the edge of the anti-glitch condition. This elucidates both the
possibility and the rarity of anti-glitch events, indicating that, for such an
event to occur, the ejecta must be expelled to a distance nearing the co-
rotation radius. However, this toy model posits that the motion of the trapped
ejecta induces free precession, leading to modulations in the pulse profile,
polarization, and timing. Since we have not observed variations in the X-ray
band, additional observational evidence is required to establish this scenario
as a dominant mechanism underlying the anti-glitch phenomenon.

Regarding the possible internal origin, Kantor & Gusakov (2014) developed an
anti-glitch model that extends the standard scenario of pulsar glitches.
Glitches in RPPs are interpreted as a momentum transfer between the superfluid
and the normal components in the inner crust or outer core, depending on where
the quantized vortices can pin (Antonelli et al., 2022). Kantor & Gusakov
(2014) extended this framework by considering that the superfluid fraction
also depends on the superfluid current (it increases with decreasing velocity
between the normal and superfluid components). This results in a mass
redistribution between the superfluid and normal part of an NS when the
velocity lag between the two components changes after the quantized vortices
unpin. Following the notation of Kantor & Gusakov (2014), the total angular
momentum conservation reads

$I_{\mathrm{c0}}\Omega_{\mathrm{c0}}+I_{\mathrm{s0}}\Omega_{\mathrm{s0}}=I_{\mathrm{c1}}\Omega_{c1}+I_{\mathrm{s1}}\Omega_{\mathrm{s1}}\,,$
(5)

where $I$ denotes the moment of inertia, $\Omega$ represents the angular
velocity, and the subscripts ‘c’ and ‘s’ refer to the normal (i.e., everything
that is non-superfluid) and superfluid components, respectively, with ‘0’
indicating pre-glitch and ‘1’ post-glitch parameters. Taking into account the
changes in the moments of inertia, $I_{\mathrm{c}}^{\prime}$ and
$I_{\mathrm{s}}^{\prime}$ due to mass transfusion between the components
($I_{\mathrm{c}}^{\prime}+I_{\mathrm{s}}^{\prime}=0$), the change in angular
velocity of the normal component after the glitch is

$\delta\Omega_{c}=-\frac{I_{s0}-I^{\prime}_{c}\Delta\Omega_{0}}{I_{c0}+I^{\prime}_{c}\Delta\Omega_{0}}\delta\Omega_{s}\,,$
(6)

where $\Delta\Omega_{0}=\Omega_{\mathrm{s0}}-\Omega_{\mathrm{c0}}>0$ is the
angular velocity lag before the unpinning and
$I_{s,c}^{\prime}=dI_{s,c}/d\Delta\Omega$ is the rate of change of the partial
moment of inertia with the lag $\Delta\Omega=\Omega_{s}-\Omega_{c}$. This
gives an anti-glitch ($\delta\Omega_{c}<0$) when
$|I^{\prime}_{s}|\Delta\Omega_{0}>I_{s0}$, a condition that can be met
depending on two physical parameters: the stellar temperature and the
rotational lag $\Delta\Omega$. The higher the internal temperature ($T\gtrsim
10^{7}\,\mathrm{K}$) and the lag before vortex unpinning ($\Delta\Omega\gtrsim
1\,$rad/s), the more likely is to meet the anti-glitch condition (see Figures
2 and 3 in Kantor & Gusakov 2014). PSR B0540-69 is arguably one of the
youngest pulsars and is a possible candidate to satisfy the anti-glitch
criterion. First, the internal temperature in young pulsars like PSR B0540-69
is expected to be relatively high, $T\gtrsim 10^{8}$K. Moreover, the increase
in the spin-down rate in 2011 (Marshall et al., 2015; Ge et al., 2019) might
have boosted the lag $\Delta\Omega$ in PSR B0540-69. It is possible to provide
an upper limit to the increase of $\Delta\Omega$ since the last reported
glitch, around MJD 52927 (Ferdman et al., 2015). The long-term spin-down rate
changed from $-1.87\times 10^{-10}\,\mathrm{Hz/s}$ to $-2.52\times
10^{-10}\,\mathrm{Hz/s}$ after a spin-down rate change at MJD 55866 (Ge et
al., 2019). Therefore, the maximum variation of the lag achieved during this
time interval is about $\delta\Omega\approx 0.88\,$rad/s, assuming no
undetected glitches have occurred and the superfluid’s angular velocity was
conserved between the last reported glitch and our anti-glitch event (i.e.,
vortices were always pinned).

This value represents an upper limit to how much $\Delta\Omega$ could have
changed from MJD 52927 to $t_{g}\approx$ MJD 60134. Its order of magnitude is
consistent with the typical values of $\Delta\Omega_{0}\sim 2$ rad/s needed to
trigger an anti-glitch for a star whose internal temperature is $T\sim 10^{8}$
K. This seems to indicate that triggering an anti-glitch is indeed possible,
considering that it is unlikely that the previous glitch at MJD 52927 would
have emptied the entire superfluid’s momentum reservoir.

This anti-glitch model is likely better suited for magnetars, as it requires a
strong toroidal magnetic field to pin vortices in the outer core (Alpar,
2017): vortex unpinning in the outer core would result in an anti-glitch
because the critical temperature for superfluid neutrons is significantly
lower there than in the inner crust, allowing for a more efficient matter
transfusion from the two components.

Thus, even though PSR B0540-69 is not classified as a magnetar, its internal
toroidal component would have to be much stronger than its poloidal field, but
this is not necessarily a problem. The Bayesian analysis of Vela’s 2016 glitch
also provides a first observational clue for pinning in the outer core due to
a strong magnetic field (Montoli et al., 2020), a scenario often invoked in
pulsar glitch models (e.g., Gügercinoğlu & Alpar, 2014; Sourie & Chamel,
2020). Furthermore, magnetohydrodynamic simulations reveal that it is possible
to develop stable internal toroidal magnetic fields that are relatively high
with respect to the poloidal component (Ciolfi & Rezzolla, 2013; Ciolfi, 2014;
Lander, 2014). Therefore, given the current observational results, we cannot
rule out the anti-glitch scenario of (Kantor & Gusakov, 2014) for PSR
B0540-69.

Finally, we note that Kantor et al. (2016) also proposed a scenario for the
evolution of rapidly rotating NSs in low-mass X-ray binaries. In such objects,
a resonant interaction of normal r-modes with superfluid inertial modes
produces new stable regions in the $T$-$\nu$ plane (again, $T$ is the internal
temperature). The spin evolution of these hot and rapidly rotating NSs moves
around the boundary between stable and unstable regions, periodically crossing
it. The periodic excitation of r-modes implies the possibility of anti-
glitches in these objects, whose frequency jumps unfold over a time scale that
varies, depending on the system, from hours to months. Although this model can
be applied to spinning down pulsars, PSR B0540-69 is located in a different
region of the $T$-$\nu$ plane, where r-modes are expected to be stable.
However, future observations of anti-glitches in PSR B0540-69 may reveal the
extent of their periodicity. In fact, while for glitches (or anti-glitches) in
the standard vortex scenario, it is natural to expect no periodicity (vortex
unpinning is a multi-threshold phenomenon, see the discussion in Antonelli et
al. 2022), the r-mode scenario relies on a single and global instability
threshold, resulting in a periodic behaviour.

In conclusion, our Bayesian analysis points to an asymptotic anti-glitch event
in an RPP, marking the first occurrence of such a phenomenon observed beyond
magnetars and accreting pulsars. The radiatively quiet nature associated with
this glitch event suggests that internal processes within the neutron star are
responsible for this anti-glitch. To date, the possible physical mechanisms
that can lead to anti-glitches are based on magnetar properties, as they were
the only isolated neutron stars that exhibited anti-glitches. We speculate
that the anti-glitch of PSR B0540-69 can be explained through the transfer of
mass from the superfluid to the normal component, as suggested by (Kantor &
Gusakov, 2014) under the condition that the internal toroidal field of the
pulsar is significantly higher than its dipole fields, a possibility that may
also be met in the Vela pulsar (Gügercinoğlu & Alpar, 2014; Montoli et al.,
2020). Further observations on timing and multi-wavelength monitoring,
together with a more systematic analysis of proposed models for anti-glitches,
are needed to understand these phenomena in RPPs.

The authors would like to express our sincere gratitude to the referee for
their insightful comments. This work is supported by the National Key R&D
Program of China (2021YFA0718500) from the Minister of Science and Technology
of China (MOST). The authors thank the support from the National Natural
Science Foundation of China (Grants 12373051, 12333007) and the International
Partnership Program of the Chinese Academy of Sciences (Grant No.
113111KYSB20190020). This work is supported by National Key R&D Program of
China (2023YFE0117200) and National Natural Science Foundation of China (Grant
12373041). We acknowledge the use of data from the _Neutron star Interior
Composition Explorer_ (NICER) mission, a project led by the NASA/Goddard Space
Flight Center. This research also uses data obtained via the High Energy
Astrophysics Science Archive Research Center Online Service, provided by the
NASA/Goddard Space Flight Center. The authors acknowledge NASA’s Astrophysics
Data System (ADS) Bibliographic Services and the arXiv repository. We are
grateful for the insightful discussions with Dr. Abdujappar Rusul, Dr. Ce Cai,
and Prof. Lin Lin.

##  Appendix A Phenomenological model for the glitch/anti-glitch residuals

For completeness, we discuss the parameter space of the phase residue model in
equation (3), and show that it can be recovered as a limit of the one used for
the fit of Vela’s 2016 glitch (Ashton et al., 2019; Montoli et al., 2020).

###  A.1 Five-parameter model for overshooting glitches and anti-glitches

Assume that $\nu_{pre}(t)$ and $\nu_{post}(t)$ describe the pulsar evolution
before and after the glitch/anti-glitch epoch, that is conveniently split at
$t_{g}=0$. The rotation frequency $\nu(t)$ and its residue $\Delta\nu(t)$ for
any time $t$ are

$\nu(t)=\theta(t)\,\nu_{post}(t)+\theta(-t)\,\nu_{pre}(t)=\nu_{pre}(t)+\Delta\nu(t)\qquad\qquad\Delta\nu(t)=\theta(t)\left(\nu_{post}(t)-\nu_{pre}(t)\right)\,,$
(A1)

where $\theta(t)$ is the Heaviside step function. Now, we have to choose a
suitable phenomenological form for $\Delta\nu(t)$. The one used for the
Bayesian fit of Vela’s 2016 glitch is (cf. equation 4 of Montoli et al., 2020)

$\begin{split}&\Delta\nu(t)\,=\,\theta(t)\,\Delta\nu_{\infty}\left[1-\omega\,e^{-t\lambda_{+}}-(1-\omega)\,e^{-t\lambda_{-}}\right]\\\
&\Delta\Phi(t)=\theta(t)\int_{0}^{t}\Delta\nu(t^{\prime})dt^{\prime}\,=\theta(t)\,\Delta\nu_{\infty}\left[t+\dfrac{\omega}{\lambda_{+}}\left(e^{-t\lambda_{+}}-1\right)+\dfrac{1-\omega}{\lambda_{-}}\left(e^{-t\lambda_{-}}-1\right)\right]\end{split}$
(A2)

where $\omega$ is a dimensionless parameter that sets the overall behaviour of
the glitch (see Figure 4), $0<\lambda_{-}<\lambda_{+}$ are the inverse of two
recovery timescales (the short one given by $\tau_{s}=1/\lambda_{+}$ and a
longer one, $\tau_{l}=1/\lambda_{-}$) and $\Delta\nu_{\infty}$ is the
asymptotic glitch/anti-glitch amplitude. We have ‘asymptotic spin-up’ events
for $\Delta\nu_{\infty}>0$, and ‘asymptotic spin-down’ events for
$\Delta\nu_{\infty}<0$. The model in equation (A2) allows fitting situations
where the normal component overshoots333 Given a glitch triggered at $t=0$,
assume that there a is a time $t_{\rm max}>0$ such that
$\Delta\dot{\nu}(t_{\rm max})=0$. A spin-up glitch ‘overshoots’ if
$\Delta{\nu}(t_{\rm max})$ is a maximum, $\Delta\ddot{\nu}(t_{\rm max})$¡0.
Conversely, an anti-glitch ‘undershoots’ if $\Delta{\nu}(t_{\rm max})$ is a
minimum, $\Delta\ddot{\nu}(t_{\rm max})$¡0. , which is realised when (cf.
equation 5 of Montoli et al., 2020)

$\frac{t_{\rm
max}}{\tau_{l}-\tau_{s}}\,=\,\log\left(\frac{\tau_{l}}{\tau_{s}}\right)+\log\left(\frac{\omega}{\omega-1}\right)\,>0\,.$
(A3)

The time $t_{\rm max}$ is a real positive number only when $\omega<0$ or
$\omega>1$. All the possible qualitative behaviours of (A2) are sketched in
Figure 4:

  1. 1.

$\omega>1$: Overshooting glitch ($\Delta\nu_{\infty}>0$) or undershooting
anti-glitch ($\Delta\nu_{\infty}<0$).

  2. 2.

$0<\omega<1$: Two-timescale spin-up/spin-down, depending on the sign of
$\Delta\nu_{\infty}$. The case $\Delta\nu_{\infty}>0$ resembles the ‘delayed
spin-up’ behaviour observed during the spin-up phase of Crab’s 2017 glitch
(Shaw et al., 2018).

  3. 3.

$\omega<0$: In this case, the glitch/anti-glitch distinction depends on
whether one looks at the short-timescale or the long-timescale behaviour.
Since this is the most exotic possibility, in Figure 4 we provide two examples
for this case (3a, $-1<\omega<0$ and 3b, $\omega<-1$), that are not
qualitatively different.

When (A2) is used in a fit procedure, we have to consider that the
glitch/anti-glitch epoch $t_{g}$ is an extra parameter that should be inferred
together with the other parameters: $\omega\in\mathbb{R}$,
$t_{g}\in\mathbb{R}$, $\Delta\nu_{\infty}\in\mathbb{R}$,
$0<\lambda_{-}<\lambda_{+}$ (or, equivalently, $0<\tau_{s}<\tau_{l}$), for a
total of 5 real parameters to be extracted from the timing data.

Figure 4:  Sketch of the possible realizations of model (A2) (left panel,
$\tau_{s}=\tau_{l}/7$) and model (A4) (right panel, $\tau_{s}=0$). The blue
curves represent events that are ‘asymptotic spin-up’ events
($\Delta\nu_{\infty}>0$), and the red ones are ‘asymptotic spin-down’ events
($\Delta\nu_{\infty}<0$). The solid curves refer to case 1: $\omega=2.5$
(overshooting glitch/undershooting anti-glitch). The dash-dotted curves refer
to case 2: $\omega=0.5$ (two-timescale spin-up/spin-down). Case 3 is
represented by the dotted (3a, $\omega=-0.5$) and dashed (3b, $\omega=-2.5$)
curves. The left panel contains the same curves in the limit
$\tau_{s}=1/\lambda_{+}=0$.

###  A.2 Unresolved short timescale: four-parameter model

Assume that the timescale $\tau_{s}$ can not be resolved with the available
timing data. The limit $\tau_{s}\ll\tau_{l}$ gives:

$\Delta\nu(t)\,=\,\theta(t)\,\Delta\nu_{\infty}\left[1-(1-\omega)\,e^{-t\lambda_{-}}\right]\qquad\quad\Delta\Phi(t)=\theta(t)\Delta\nu_{\infty}\left[t+\dfrac{1-\omega}{\lambda_{-}}\left(e^{-t\lambda_{-}}-1\right)\right]$
(A4)

Hence, we have an instantaneous jump of amplitude
$\Delta\nu_{0}=\omega\Delta\nu_{\infty}$ at $t=t_{g}$, that then relaxes to
$\Delta\nu_{\infty}$ for $t\gg\tau_{l}$. The model (A4) is formally equivalent
to the one of Baym et al. (1969), see Antonelli et al. (2022) for a revised
discussion of the original model. This can be seen by introducing the ‘healing
parameter’ $Q$,

$\Delta\nu(t)\,=\,\theta(t)\,\Delta\nu_{0}\left[1-Q\left(1-e^{-t\lambda_{-}}\right)\right]\qquad\quad
Q=(\Delta\nu_{0}-\Delta\nu_{\infty})/\Delta\nu_{0}=(\omega-1)/\omega\,.$ (A5)

Typically, $0<Q<1$ is observed (Crawford & Demiański, 2003): $Q=0$ is a step-
like glitch with no detectable relaxation, and $Q=1$ is a glitch with complete
relaxation. Since the condition $0<Q<1$ is equivalent to $\omega>1$, the
simplest model to fit glitches formally corresponds to the $\tau_{s}=0$ limit
of an agnostic model for overshooting glitches. Moreover, as long as (A5) is
regarded as an agnostic model, we do not necessarily have to stick to the
parameter space of typical spin-up glitches, and we can assume that
$Q\in\mathbb{R}$. Despite its simplicity, the model (A4) is already quite rich
as it can represent all the possibilities sketched in Figure 4. Hence, it
provides a reasonably simple prescription to fit situations where one is
unsure about the event’s glitch/anti-glitch nature, as in our case.

###  A.3 Agnostic model for the Bayesian fit

The Bayesian fit is performed by considering all the timing data available
before and after the (unknown) glitch epoch $t_{g}$. For the frequency
residues model, we adopt the form (A4) but we also extend it to account for
the possibility of a permanent change of the frequency derivative:

$\Delta\nu(t)=\theta(\delta t)\left[\Delta\nu+\Delta\dot{\nu}\,\delta
t+\Delta\nu_{d}\,e^{-\delta t/\tau}\right]\,,$ (A6)

cf. with equation (A5) by using $\Delta\nu=\Delta\nu_{\infty}$ (the
‘permanent’ frequency jump is just the ‘asymptotic’ frequency jump),
$\Delta\nu_{d}=\Delta\nu_{0}-\Delta\nu_{\infty}$,
$Q=\Delta\nu_{d}/\Delta\nu_{0}$, $\tau=1/\lambda_{-}$. Therefore, the full
model that we have to fit is $\nu(t)=\nu_{pre}(t)+\Delta\nu(t)$. Since we have
to fit the ToA, in practice we fit the corresponding phase model, namely
$\Phi(t)=\Phi_{pre}(t)+\Delta\Phi(t)$, where $\Phi_{pre}(t)$ is exactly given
in (1) and

$\Delta\Phi(t)=\theta(\delta t)\left[\Delta\nu\,\delta
t+\Delta\nu_{d}\,\tau\left(1-e^{-\delta
t/\tau}\right)+\frac{\Delta\dot{\nu}}{2}\delta t^{2}\right]\,,$ (A7)

for a total of 8 real parameters (i.e., $\nu$, $\dot{\nu}$, $\ddot{\nu}$,
$t_{g}$, $\Delta\nu$, $\Delta\dot{\nu}$, $\Delta\nu_{d}$, $\tau$) to be
extracted from the timing data. For these parameters, we take flat prior
distributions within reasonably large intervals:

$\displaystyle\nu\in[19.63,19.64]\,\text{Hz}\qquad\dot{\nu}\in[-2.6,-2.5]\times
10^{-10}\,\text{Hz/s}\qquad\ddot{\nu}\in[4.3\times 10^{-24},4.3\times
10^{-19}]\,\text{Hz/s}^{2}\qquad$ $\displaystyle\Delta{\nu}\in[-1,1]\times
10^{-4}\,\text{Hz}\qquad\Delta{\dot{\nu}}\in[-1,1]\times
10^{-12}\,\text{Hz/s}^{2}\qquad\Delta{\nu}_{d}\in[-1,1]\times
10^{-4}\,\text{Hz}\qquad$ (A8) $\displaystyle\tau\in[0,200]\,\text{day}\qquad
t_{g}\in[60120,60150]\,\text{MJD}.$

## References

  * Aasi et al. (2015) Aasi, J., Abbott, B., Abbott, R., et al. 2015, Physical Review D, 91, 022004
  * Alpar (2017) Alpar, M. A. 2017, Journal of Astrophysics and Astronomy, 38, 44, doi: 10.1007/s12036-017-9473-6
  * Anderson & Itoh (1975) Anderson, P., & Itoh, N. 1975, Nature, 256, 25
  * Antonelli et al. (2022) Antonelli, M., Montoli, A., & Pizzochero, P. M. 2022, Insights into the physics of neutron star interiors from pulsar glitches (World Scientific), 219–281, doi: 10.1142/9789811220944_0007
  * Antonopoulou et al. (2022) Antonopoulou, D., Haskell, B., & Espinoza, C. M. 2022, Reports on Progress in Physics
  * Archibald et al. (2017) Archibald, R., Kaspi, V., Scholz, P., et al. 2017, The Astrophysical Journal, 834, 163
  * Archibald et al. (2013) Archibald, R. F., Kaspi, V., Ng, C.-Y., et al. 2013, Nature, 497, 591
  * Ashton et al. (2019) Ashton, G., Lasky, P. D., Graber, V., & Palfreyman, J. 2019, Nature Astronomy, 3, 1143
  * Basu et al. (2022) Basu, A., Shaw, B., Antonopoulou, D., et al. 2022, Monthly Notices of the Royal Astronomical Society, 510, 4049
  * Baym et al. (1969) Baym, G., Pethick, C., Pines, D., & Ruderman, M. 1969, Nature, 224, 872, doi: 10.1038/224872a0
  * Baym & Pines (1971) Baym, G., & Pines, D. 1971, Annals of Physics, 66, 816
  * Borghese & Esposito (2023) Borghese, A., & Esposito, P. 2023, in Handbook of X-ray and Gamma-ray Astrophysics (eds. C. Bambi, 146, doi: 10.1007/978-981-16-4544-0_102-1
  * Brown et al. (2021) Brown, A. G., Vallenari, A., Prusti, T., et al. 2021, Astronomy & Astrophysics, 649, A1
  * Cai et al. (2021) Cai, C., Xiong, S., Li, C., et al. 2021, Monthly Notices of the Royal Astronomical Society, 508, 3910
  * Çerri-Serim et al. (2019) Çerri-Serim, D., Serim, M. M., Şahiner, Ş., Inam, S. C., & Baykal, A. 2019, MNRAS, 485, 2, doi: 10.1093/mnras/sty3213
  * Chukwude & Urama (2010) Chukwude, A. E., & Urama, J. O. 2010, Monthly Notices of the Royal Astronomical Society, 406, 1907
  * Ciolfi (2014) Ciolfi, R. 2014, Astronomische Nachrichten, 335, 624, doi: 10.1002/asna.201412083
  * Ciolfi & Rezzolla (2013) Ciolfi, R., & Rezzolla, L. 2013, MNRAS, 435, L43, doi: 10.1093/mnrasl/slt092
  * Crawford & Demiański (2003) Crawford, F., & Demiański, M. 2003, The Astrophysical Journal, 595, 1052
  * Dai et al. (2018) Dai, S., Johnston, S., Weltevrede, P., et al. 2018, Monthly Notices of the Royal Astronomical Society, 480, 3584
  * D’Alessandro (1996) D’Alessandro, F. 1996, Ap&SS, 246, 73, doi: 10.1007/BF00637401
  * Dib & Kaspi (2014) Dib, R., & Kaspi, V. M. 2014, The Astrophysical Journal, 784, 37
  * Ducci et al. (2015) Ducci, L., Pizzochero, P., Doroshenko, V., et al. 2015, Astronomy & Astrophysics, 578, A52
  * Enoto et al. (2019) Enoto, T., Kisaka, S., & Shibata, S. 2019, Reports on Progress in Physics, 82, 106901
  * Enoto et al. (2021) Enoto, T., Terasawa, T., Kisaka, S., et al. 2021, Science, 372, 187
  * Espinoza et al. (2011) Espinoza, C. M., Lyne, A. G., Stappers, B. W., & Kramer, M. 2011, Monthly Notices of the Royal Astronomical Society, 414, 1679
  * Ferdman et al. (2015) Ferdman, R. D., Archibald, R. F., & Kaspi, V. M. 2015, The Astrophysical Journal, 812, 95
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, Publications of the Astronomical Society of the Pacific, 125, 306
  * García & Ranea-Sandoval (2015) García, F., & Ranea-Sandoval, I. F. 2015, Monthly Notices of the Royal Astronomical Society: Letters, 449, L73
  * Gavriil et al. (2011) Gavriil, F. P., Dib, R., & Kaspi, V. M. 2011, The Astrophysical Journal, 736, 138
  * Ge et al. (2019) Ge, M., Lu, F., Yan, L., et al. 2019, Nature Astronomy, 3, 1122
  * Ge et al. (2020) Ge, M., Zhang, S., Lu, F., et al. 2020, The Astrophysical Journal, 896, 55
  * Ge et al. (2024) Ge, M.-Y., Yang, Y.-P., Lu, F.-J., et al. 2024, Research in Astronomy and Astrophysics, 24, 015016
  * Gendreau et al. (2016) Gendreau, K. C., Arzoumanian, Z., Adkins, P. W., et al. 2016, in Space telescopes and instrumentation 2016: Ultraviolet to gamma ray, Vol. 9905, SPIE, 420–435
  * Gügercinoğlu & Alpar (2014) Gügercinoğlu, E., & Alpar, M. A. 2014, ApJ, 788, L11, doi: 10.1088/2041-8205/788/1/L11
  * Harding (2013) Harding, A. K. 2013, Frontiers of Physics, 8, 679, doi: 10.1007/s11467-013-0285-0
  * Haskell & Melatos (2015) Haskell, B., & Melatos, A. 2015, International Journal of Modern Physics D, 24, 1530008
  * Hobbs et al. (2006) Hobbs, G., Edwards, R., & Manchester, R. 2006, Monthly Notices of the Royal Astronomical Society, 369, 655
  * Hobbs et al. (2010) Hobbs, G., Archibald, A., Arzoumanian, Z., et al. 2010, Classical and Quantum Gravity, 27, 084013
  * Howitt & Melatos (2022) Howitt, G., & Melatos, A. 2022, MNRAS, 514, 863, doi: 10.1093/mnras/stac1358
  * Huang & Geng (2014) Huang, Y., & Geng, J. 2014, The Astrophysical Journal Letters, 782, L20
  * Huppenkothen et al. (2019) Huppenkothen, D., Bachetti, M., Stevens, A. L., et al. 2019, The Astrophysical Journal, 881, 39
  * Içdem et al. (2012) Içdem, B., Baykal, A., & İnam, S. Ç. 2012, Monthly Notices of the Royal Astronomical Society, 419, 3109
  * Kantor & Gusakov (2014) Kantor, E., & Gusakov, M. 2014, The Astrophysical Journal Letters, 797, L4
  * Kantor et al. (2016) Kantor, E. M., Gusakov, M. E., & Chugunov, A. I. 2016, MNRAS, 455, 739, doi: 10.1093/mnras/stv2352
  * Kaspi & Beloborodov (2017) Kaspi, V. M., & Beloborodov, A. M. 2017, Annual Review of Astronomy and Astrophysics, 55, 261
  * Lander (2014) Lander, S. K. 2014, Monthly Notices of the Royal Astronomical Society, 437, 424
  * Livingstone et al. (2010) Livingstone, M. A., Kaspi, V. M., & Gavriil, F. P. 2010, The Astrophysical Journal, 710, 1710
  * Manchester et al. (2016) Manchester, R., Hobbs, G., Teoh, A., & Hobbs, M. 2016, VizieR Online Data Catalog, B
  * Marshall et al. (2015) Marshall, F., Guillemot, L., Harding, A., Martin, P., & Smith, D. 2015, The Astrophysical journal letters, 807, L27
  * Marshall et al. (2016) Marshall, F. E., Guillemot, L., Harding, A., Martin, P., & Smith, D. A. 2016, The Astrophysical journal letters, 827, L39
  * Mastrano et al. (2015) Mastrano, A., Suvorov, A. G., & Melatos, A. 2015, Monthly Notices of the Royal Astronomical Society, 453, 522
  * Mereghetti et al. (2015) Mereghetti, S., Pons, J. A., & Melatos, A. 2015, Space Sci. Rev., 191, 315, doi: 10.1007/s11214-015-0146-y
  * Montoli et al. (2020) Montoli, A., Antonelli, M., Magistrelli, F., & Pizzochero, P. 2020, Astronomy & Astrophysics, 642, A223
  * Ng & Kaspi (2011) Ng, C. Y., & Kaspi, V. M. 2011, in American Institute of Physics Conference Series, Vol. 1379, AstroPhysics of Neutron Stars 2010: A Conference in Honor of M. Ali Alpar, ed. E. Göğüş, T. Belloni, & Ü. Ertan (AIP), 60–69, doi: 10.1063/1.3629486
  * Petre et al. (2007) Petre, R., Hwang, U., Holt, S., Safi-Harb, S., & Williams, R. 2007, The Astrophysical Journal, 662, 988
  * Pizzochero et al. (2020) Pizzochero, P. M., Montoli, A., & Antonelli, M. 2020, A&A, 636, A101, doi: 10.1051/0004-6361/201937019
  * Ray et al. (2019) Ray, P. S., Guillot, S., Ho, W. C., et al. 2019, The Astrophysical Journal, 879, 130
  * Robitaille et al. (2013) Robitaille, T. P., Tollerud, E. J., Greenfield, P., et al. 2013, Astronomy & Astrophysics, 558, A33
  * Ruderman (1969) Ruderman, M. 1969, Nature, 223, 597, doi: 10.1038/223597b0
  * Şaşmaz Muş et al. (2014) Şaşmaz Muş, S., Aydın, B., & Göğüş, E. 2014, Monthly Notices of the Royal Astronomical Society, 440, 2916
  * Serim et al. (2017) Serim, M. M., Şahiner, Ş., Çerri-Serim, D., Inam, S. C., & Baykal, A. 2017, MNRAS, 471, 4982, doi: 10.1093/mnras/stx1771
  * Seward et al. (1984) Seward, F., Harnden Jr, F., & Helfand, D. J. 1984, Astrophysical Journal, Part 2-Letters to the Editor (ISSN 0004-637X), vol. 287, Dec. 1, 1984, p. L19-L22. Research supported by the Alfred P. Sloan Foundation and the Danish Space Board., 287, L19
  * Shaw et al. (2018) Shaw, B., Lyne, A. G., Stappers, B. W., et al. 2018, MNRAS, 478, 3832, doi: 10.1093/mnras/sty1294
  * Sourie & Chamel (2020) Sourie, A., & Chamel, N. 2020, MNRAS, 493, L98, doi: 10.1093/mnrasl/slaa015
  * Taylor (1992) Taylor, J. H. 1992, Philosophical Transactions of the Royal Society of London. Series A: Physical and Engineering Sciences, 341, 117
  * Tong (2014) Tong, H. 2014, The Astrophysical Journal, 784, 86
  * Tuo et al. (2022) Tuo, Y., Li, X., Ge, M., et al. 2022, The Astrophysical Journal Supplement Series, 259, 14
  * Vahdat et al. (2024) Vahdat, A., Posselt, B., Pavlov, G. G., et al. 2024, ApJ, 963, 138, doi: 10.3847/1538-4357/ad20ca
  * Wang et al. (2012) Wang, J., Wang, N., Tong, H., & Yuan, J. 2012, Ap&SS, 340, 307, doi: 10.1007/s10509-012-1058-x
  * Wang et al. (2020) Wang, L., Ge, M., Wang, J., et al. 2020, Monthly Notices of the Royal Astronomical Society, 494, 1865
  * Wong et al. (2001) Wong, T., Backer, D., & Lyne, A. 2001, The Astrophysical Journal, 548, 447
  * Woods et al. (1999) Woods, P. M., Kouveliotou, C., Van Paradijs, J., et al. 1999, The Astrophysical Journal, 524, L55
  * Yim et al. (2024) Yim, G., Gao, Y., Kang, Y., Shao, L., & Xu, R. 2024, Monthly Notices of the Royal Astronomical Society, 527, 2379
  * Younes et al. (2020) Younes, G., Ray, P. S., Baring, M. G., et al. 2020, The Astrophysical Journal Letters, 896, L42
  * Younes et al. (2023) Younes, G., Baring, M., Harding, A., et al. 2023, Nature Astronomy, 7, 339
  * Yu et al. (2013) Yu, M., Manchester, R., Hobbs, G., et al. 2013, Monthly Notices of the Royal Astronomical Society, 429, 688
  * Zhou et al. (2022) Zhou, S., Gügercinoğlu, E., Yuan, J., Ge, M., & Yu, C. 2022, Universe, 8, 641