††thanks: These two authors contributed equally††thanks: These two authors
contributed equally

# Time-dependent inertia of self-propelled particles: the Langevin rocket

Alexander R. Sprenger  Institut für Theoretische Physik II: Weiche Materie,
Heinrich-Heine-Universität Düsseldorf, D-40225 Düsseldorf, Germany  Soudeh
Jahanshahi  Institut für Theoretische Physik II: Weiche Materie, Heinrich-
Heine-Universität Düsseldorf, D-40225 Düsseldorf, Germany  Alexei V. Ivlev
Max-Planck-Institut für Extraterrestrische Physik, 85748 Garching, Germany
Hartmut Löwen  Institut für Theoretische Physik II: Weiche Materie, Heinrich-
Heine-Universität Düsseldorf, D-40225 Düsseldorf, Germany

###### Abstract

Many self-propelled objects are large enough to exhibit inertial effects but
still suffer from environmental fluctuations. The corresponding basic
equations of motion are governed by active Langevin dynamics which involve
inertia, friction and stochastic noise for both the translational and
orientational degrees of freedom coupled via the self-propulsion along the
particle orientation. In this paper, we generalize the active Langevin model
to time-dependent parameters and explicitly discuss the effect of time-
dependent inertia for achiral and chiral particles. Realizations of this
situation are manifold ranging from minirockets which are self-propelled by
burning their own mass, dust particles in plasma which lose mass by
evaporating material to walkers with expiring activity. Here we present
analytical solutions for several dynamical correlation functions such as the
mean-square displacement and the orientational and velocity autocorrelation
functions. If the parameters exhibit a slow power-law in time, we obtain
anomalous superdiffusion with a non-trivial dynamical exponent. Finally we
constitute the “Langevin rocket” model by including orientational fluctuations
in the traditional Tsiolkovsky rocket equation. We calculate the mean reach of
the Langevin rocket and discuss different mass ejection strategies to maximize
it. Our results can be tested in experiments on macroscopic robotic or living
particles or in self-propelled mesoscopic objects moving in media of low
viscosity such as complex plasma.

##  I Introduction

The non-equilibrium dynamics of active Brownian particles –also referred to as
microswimmers– are typically described in the overdamped limit, where inertial
effects are sufficiently small relative to viscous ones Howse  _et al._
(2007); Romanczuk  _et al._ (2012); Elgeti  _et al._ (2015); Bechinger  _et
al._ (2016). This is an excellent approximation for micron-sized self-
propelled particles swimming in a viscous Newtonian liquid such as water
Purcell (1977) at low Reynolds number. The standard model of a single active
Brownian particle Howse  _et al._ (2007); ten Hagen  _et al._ (2011);
Bechinger  _et al._ (2016) involves a translational and an orientational
degree of freedom and includes Stokesian friction and fluctuations. These
degrees of freedom are coupled via the self-propulsion along the particle
orientation which is modeled in a simple averaged way by an internal velocity,
sometimes referred to as the particle activity.

However, inertial effects become relevant for larger particle sizes or the
motion in gaseous media of lower viscosity. Though highly relevant for
swimming and flying organisms as well as for autonomous machines (e.g., flying
insect-drones, marine robots) Klotsa (2019) mesoscale active matter at
intermediate Reynolds number has been much less studied. Aiming at a simple
description of a single particle first, one basic model is that of active
Langevin motion Scholz  _et al._ (2018a); Callegari and Volpe (2019); Um  _et
al._ (2019); Löwen (2020): it generalizes the common overdamped model of
active Brownian motion Howse  _et al._ (2007); ten Hagen  _et al._ (2011);
Bechinger  _et al._ (2016) towards underdamped dynamics by including finite
particle mass and moment of inertia in the equations of motion Enculescu and
Stark (2011); Ghosh  _et al._ (2015); Joyeux and Bertin (2016); Manacorda and
Puglisi (2017); Takatori and Brady (2017); Mokhtari  _et al._ (2017); Das  _et
al._ (2019); Sandoval (2020); Gutierrez-Martinez and Sandoval (2020). The
inertial self-propelled particles may therefore be called “microflyers”
(rather than “microswimmers”), sometimes they are also termed “runners”,
“walkers” or “hoppers” Ebeling  _et al._ (1999); Weber  _et al._ (2013);
Parisi  _et al._ (2018). Examples for inanimate inertial self-propelled
particles modeled by active Langevin dynamics are manifold, they include a
complex plasma consisting of mesoscopic dust particles in a weakly ionized gas
Morfill and Ivlev (2009); Sütterlin  _et al._ (2009); Couëdel  _et al._
(2010); Chaudhuri  _et al._ (2011); Ivlev  _et al._ (2015); Nosenko  _et al._
(2020), vibration-driven granular particles Narayan  _et al._ (2007); Kudrolli
_et al._ (2008); Deseigne  _et al._ (2010); Giomi  _et al._ (2013); Weber  _et
al._ (2013); Klotsa  _et al._ (2015); Patterson  _et al._ (2017); Junot  _et
al._ (2017); Ramaswamy (2017); Deblais  _et al._ (2018); Dauchot and Démery
(2019), autorotating seeds and fruits Rabault  _et al._ (2019); Fauli  _et
al._ (2019), camphor surfers Leoni  _et al._ (2020), hexbug crawlersLeoni  _et
al._ (2020), trapped aerosols Di Leonardo  _et al._ (2007) and mini-robots
Rubenstein  _et al._ (2014); Fujiwara  _et al._ (2014); Tolba  _et al._
(2015); Zhakypov  _et al._ (2019); Yang  _et al._ (2019). Moreover, there are
numerous examples of animals moving at intermediate Reynolds number such as
swimming organisms (nematodes, brine shrimps, whirligig beetles etc) Klotsa
(2019); Turner  and flying insects and birds Toner and Tu (1995, 1998);
Chiappini (2008); Bartussek and Lehmann (2016); Mukundarajan  _et al._ (2016);
Bartussek and Lehmann (2018); Attanasi  _et al._ (2014).

In this paper, we extend the active Langevin model to time-dependent
parameters such as time-dependent inertia, self-propulsion and friction. This
is a situation frequently encountered in nature and realizable in laboratory
experiments on artificial self-propelled objects. Let us mention some
examples: scallops move their shells and accelerate by jet propulsion.
Therefore, they become smaller in the course of the motion such that their
moment of inertia and their friction coefficients become time-dependent.
Moving animals typically have a finite energy reservoir Ebeling  _et al._
(1999) implying that their self-propulsion velocity is getting slower as a
function of time. The maneuverability of animal motion is provided by changes
in the body shape Vandenberg and Rayner (1995); Qiao  _et al._ (2014), which
implies a change in the moment of inertia at fixed total mass. Likewise, in
the inanimate world, minirockets, which are propelled by ejecting mass, are
getting lighter as a function of time Krasheninnikov  _et al._ ; Nosenko  _et
al._ (2010). Similarly, inflated toy balloons Müller and Strehlow (1999);
Mangan and Destrade (2015) are self-propelled by jet propulsion and strongly
subject to random fluctuations in their orientation, their body size shrinks
as a function of time and so does the mass, the moment of inertia, the
friction coefficient and the self-propulsion speed. Granulate hoppers equipped
with an internal vibration motor (”hexbugs”) Dauchot and Démery (2019), will
consume energy such that the self-propulsion speed will slowly expire and fade
away as a function of time. Robots that pick up or release objects possess a
mass variation Eke and Mao (2002), too, and a time-dependent mass can bring
about time-dependent friction coefficients Panchal and Jangid (2008). Last but
not least, any prescribed time-dependence can be programmed artificially at
wish for man-made robots, artificial walkers and microswimmer: the self-
propulsion speed can be made time-dependent by exposing particles to external
optical fields Lozano  _et al._ (2019), the noise strength can be steered by
external fields Fernandez-Rodriguez  _et al._ (2020); Sprenger  _et al._
(2020), the damping by the solvent viscosity Embs  _et al._ (2000); Stickler
_et al._ (2018) or both by the external vibration amplitude and frequency
Couder and Fort (2006); Klotsa  _et al._ (2015); Valani  _et al._ (2019);
Scholz  _et al._ (2018a).

For the active Langevin model with prescribed time-dependent parameters, we
present here analytical solutions for several dynamical correlation functions
such as the orientational and velocity autocorrelation function, the mean
displacements, the mean-square displacement and the delay function. Our
results are as follows: First, we constitute a model which we refer to as the
Langevin rocket. In doing so, we combine orientational fluctuations and mass
loss described by the traditional Tsiolkovsky rocket equation Tsiolkovsky
(1903). We calculate the mean reach of the Langevin rocket and discuss
different mass ejection strategies to maximize it. For increasing rotational
noise, the optimal strategy to achieve a maximal reach changes discontinuously
from a complete mass ejection extended over a long time to a instantaneous
ejection of a mass fraction. Second, we compare different set-ups of time-
dependent inertia such as directed and isotropic mass ejection and isotropic
shape changes with constant mass. Last, we study the case of slow
(”adiabatic”) variation of system parameters. In particular, for a change in
the system parameters described by a power law in time we predict a
superdiffusive anomalous diffusion involving a mean-square displacement
$\propto t^{\alpha}$ which scales as a power law in time $t$ with a non-
trivial exponent $\alpha$ Metzler  _et al._ (1999); Golestanian (2009);
Hoefling and Franosch (2013); Bleibel  _et al._ (2014); Morin  _et al._
(2017); Taloni  _et al._ (2017); Charalambous  _et al._ (2019); Chaki and
Chakrabarti (2019); Bodrova  _et al._ (2016). In particular, we discuss chiral
particles exposed to a torque which exhibit circling motion. This generalizes
earlier work for overdamped systems van Teeffelen and Löwen (2008); Babel  _et
al._ (2014); Olsen (2020). Our predictions can be tested in various
experimental set-ups ranging from macroscopic vibrated granular matter, robots
or living systems to self-propelled micron-particles that are flying in a
gaseous medium or in a plasma.

The paper is organized as follows: in Sec. II, we introduce the theoretical
model for active Langevin motion describing an inertial particle. In Section
III we recapitulate the case of time-independent self-propulsion, inertia,
damping and fluctuations found earlier Ghosh  _et al._ (2015); Scholz  _et
al._ (2018a) but include here also new results such as an analytical
expression for the time-resolved mean trajectory and mean-square-displacement.
In Sec. IV, we demonstrate how time-dependent parameters change the dynamics
of the system: in particular, we introduce the Langevin rocket model and study
slow temporal variations. Finally, we conclude in Sec. V.

##  II Basic model and different set-ups

Figure 1:  Self-propelled inertial particle with center position
$\bm{\mathrm{R}}(t)$ at time $t$ moving with its center in the two-dimensional
$xy$-plane. The particle position is indicated as $\bm{\mathrm{R}}(t)$ (black
arrow). Moreover the particle possesses an orientational degree of freedom
which is characterized by a unit vector
$\bm{\hat{\mathbf{n}}}=(\cos\phi,\sin\phi)$ with $\phi$ denoting the angle
relative to the $x$-axis. The particle self-propels along its orientation with
the velocity $v_{0}\bm{\hat{\mathbf{n}}}$ (red arrow). It may also experience
a torque $M$ along the $z$-axis leading to rotational motion as indicated by
the blue arrow. The translational motion is further influenced by a
translational friction $\xi$ and the noise strength $D$ (as indicated by the
light red horizontal double arrow) while the rotational motion is influenced
by a rotational friction $\xi_{r}$ and the orientational noise strength
$D_{r}$ (as indicated by the light blue curved double arrow).  Figure 2:
Schematic illustration of the different special set-ups for an active inertial
particle. The particle is shown as a dark-gray sphere and its inertia is
characterized by the particle mass $m$ and its moment of inertia $J$. (a)
Time-independent inertia with constant $m$ and $J$ as a reference situation
(gray background). (b) Directed mass ejection: Per unit time, the mass
$\dot{m}(t)$ is ejected centrally with a velocity $-u\bm{\hat{\mathbf{n}}}(t)$
along the particle orientation $\bm{\hat{\mathbf{n}}}(t)$ which leads to a
change $-u\bm{\hat{\mathbf{n}}}\dot{m}(t)$ in the translational momentum of
the particle and a time-dependent particle mass $m(t)$ but a constant moment
of inertia $J$ (red background). (c) Isotropic mass evaporation. Here the
translational and the angular momentum of the particle are both conserved but
the particle mass $m(t)$ and moment of inertia $J(t)$ are time-dependent
(green background). (d) Isotropic change in the particle shape. Here again the
linear and the angular momentum of the particle are both conserved, the
particle mass $m$ is constant but the moment of inertia $J(t)$ is time-
dependent (blue background).

In this section, we define the basic model of underdamped Langevin motion for
a self-propelled particle with time-dependent inertia. We consider a self-
propelled inertial particle with a center-of-mass coordinate
$\bm{\mathrm{R}}(t)$ at time $t$ moving with its center in the two-dimensional
$xy$-plane, see Figure 1 for a sketch. The particle is polar such that it
possesses an orientational degree of freedom characterized by a unit vector
$\bm{\hat{\mathbf{n}}}(t)=(\cos\phi(t),\sin\phi(t))$ where $\phi(t)$ is the
angle relative to the $x$-axis. The particle self-propels along its
orientation with the self-propulsion velocity $v_{0}\bm{\hat{\mathbf{n}}}$,
also indicated in Figure 1. It may also additionally be exposed to an external
or internal torque $M$ along the $z$-axis leading to an angular velocity as
shown by the blue arrow in Figure 1. As the particle has inertia in both
translation and rotation, its configuration fully specified by its center-of-
mass coordinate $\bm{\mathrm{R}}(t)$, its center-of-mass velocity
$\dot{\bm{\mathrm{R}}}(t)=d\bm{\mathrm{R}}(t)/dt$, its orientational angle
$\phi(t)$ and its angular velocity $\dot{\phi}(t)$.

While previous work Scholz  _et al._ (2018a); Löwen (2020); Sandoval (2020)
has considered constant particle mass and moment of inertia, here we
generalize the model towards time-dependent parameters with a particular focus
on a time-dependent particle mass $m(t)$ and a time-dependent moment of
inertia $J(t)$, which we define with respect to the center-of-mass to describe
the rotation around the $z$-axis. It turns out that the corresponding
equations of motion need to be discussed with care as they depend on the
physical origin of the change in inertia. In order to do this systematically
step by step, we first consider four different set-ups which are outlined in
Figure 2 and which are actually realizable in nature. We then give the most
general model equation which accommodates all these set-ups as special cases.

###  II.1 Time-independent inertia

First of all, as a reference, the special case of time-independent inertia is
considered. This set-up is sketched in Figure 2(a) (gray background). The
particle has a constant mass $m$ and a constant moment of inertia $J$. In this
case, the Langevin equation of motion reads as

$\displaystyle
m\,\ddot{\bm{\mathrm{R}}}(t)=\xi\,v_{0}\,\bm{\hat{\mathbf{n}}}(t)-\xi\,\dot{\bm{\mathrm{R}}}(t)+\xi\sqrt{2D}\,\bm{\mathrm{f}}_{\rm
st}(t),$ (1) $\displaystyle
J\,\ddot{\phi}(t)=M-\xi_{r}\,\dot{\phi}(t)+\xi_{r}\sqrt{2D_{r}}\,\tau_{\rm
st}(t).$ (2)

As far as the translational dynamics is concerned, there is a frictional
damping force $-\xi\,\dot{\bm{\mathrm{R}}}(t)$ and a self-propelling effective
force along the particle orientation $\xi\,v_{0}\,\bm{\hat{\mathbf{n}}}(t)$
which gives rise to the particle self-propulsion velocity $v_{0}$ ten Hagen
_et al._ (2015). The latter does not stem from mass ejection but is of another
origin, such as diffusiophoresis or photophoresis. This self-propulison force
couples the orientational and translational degrees of freedom. Furthermore
there is a stochastic force (“noise”) $\xi\sqrt{2D}\,\bm{\mathrm{f}}_{\rm
st}(t)$ where the effective translational diffusion coefficient $D$ quantifies
the noise strength. We describe the stochastic term $\bm{\mathrm{f}}_{\rm
st}(t)$ as zero-mean Gaussian white noise with unit variance

$\overline{\bm{\mathrm{f}}_{\rm st}(t)\otimes\bm{\mathrm{f}}_{\rm
st}(t^{\prime})}=\delta(t-t^{\prime})\bm{\mathbb{I}},$ (3)

where $\,\overline{\dotsb}\,$ indicates a noise average and $\bm{\mathbb{I}}$
is the unit matrix. Likewise, the rotational dynamics in Eq. (2) involves a
frictional torque $-\xi_{r}\,\dot{\phi}$ and an imposed torque $M$ plus the
stochastic torque $\xi_{r}\sqrt{2D_{r}}\,\tau_{\rm st}(t)$ where the effective
rotational diffusion coefficient $D_{r}$ now quantifies the rotational noise
strength and the Gaussian noise $\tau_{\rm st}(t)$ has zero-mean and unit
variance

$\overline{\tau_{\rm st}(t)\tau_{\rm st}(t^{\prime})}=\delta(t-t^{\prime}).$
(4)

One of the best experimental realization of active Langevin motion (see Eqs.
(1) and (2)) can be found in self-propelled granular particles. These
particles are capable of transferring the energy of a vibrating surface or an
internal motor to translational or rotational motion. Asymmetry in the
particle design causes them to jump forward or to rotate when lifted from the
ground. From a recent experiment on these active granular particles Scholz
_et al._ (2018a), we list exemplary orders of magnitude for our model
parameters $m=1\leavevmode\nobreak\ \text{g}$, $J=10\leavevmode\nobreak\
\text{g}\,\text{mm}^{2}$, $\xi=10\leavevmode\nobreak\ \text{g}/\text{s}$,
$\xi_{r}=100\leavevmode\nobreak\ \text{g}\leavevmode\nobreak\
\text{mm}^{2}/\text{s}$, $D=100\leavevmode\nobreak\ \text{mm}^{2}/\text{s}$,
$D_{r}=1/s$, $v_{0}=10-100\leavevmode\nobreak\ \text{mm}/\text{s}$ and
$M=10^{-7}\leavevmode\nobreak\ \text{N}\,\text{m}$. We shall revisit this
standard situation again in Section III. In the absence of any inertial
effects, i.e. when $m=J=0$, the equations of motion are overdamped and lead to
the standard picture of active Brownian motion Howse  _et al._ (2007); ten
Hagen  _et al._ (2011); Bechinger  _et al._ (2016).

###  II.2 Directed mass ejection

A rocket is self-propelled by directed mass ejection, so it establishes a
fundamental set-up of time-dependent inertia. In the typical geometry assumed
here and shown in Figure 2(b) (red background) the direction of the mass
ejection is centrally outwards opposite to the particle orientation. For
simplicity, the mass ejection occurs with a constant velocity $u$ relative to
the moving rocket ($u>0$) and the outlet coincides with the center of mass as
indicated by a wedge in Fig. 2(b). The general case where the ejection occurs
not from the center but from a point distant to the center leads to additional
terms which complicate the analysis and is left for future studies.

We assume, however, here for more generality that the rocket also has an
internal motor, which leads to an additional self-propulsion of velocity
$v_{0}$. In typical descriptions of macroscopic rockets translational and
rotational fluctuations are ignored. While this is a reasonable assumption for
macroscopic rockets, it breaks down for mini-rockets. The characteristic
equations of motion for a self-propelled particle with directed mass ejection
are

$\displaystyle\frac{d}{dt}\Big{(}m(t)\,\dot{\bm{\mathrm{R}}}(t)\Big{)}=$
$\displaystyle\xi\,v_{0}\,\bm{\hat{\mathbf{n}}}(t)-\xi\,\dot{\bm{\mathrm{R}}}(t)+\xi\sqrt{2D}\,\bm{\mathrm{f}}_{\rm
st}(t)$
$\displaystyle-\dot{m}(t)(u\,\bm{\hat{\mathbf{n}}}(t)-\dot{\bm{\mathrm{R}}}(t)),$
(5)

and the orientational equation of motion is given by (2).

In discussing the basic physics of Eq. (5), we use Newton’s second postulate
which states that the total change in translational momentum is the total
force which is in this case the sum of friction, translational stochastic and
self-propulsion forces. But even in the force-free case the ejected mass
carries away the momentum
$\dot{m}(t)(u\,\bm{\hat{\mathbf{n}}}(t)-\dot{\bm{\mathrm{R}}}(t))$ per unit
time which needs to be included in the balance of (5) with a minus sign due to
the conservation of total momentum, see also Sommerfeld (1952); Thomson
(1966); Plastino and Muzzio (1992). This constitutes in fact the thrust force
which accelerates the rocket. It is important to note here that the special
case of the traditional Tsiolkovsky rocket equation is obtained as a special
limit of no fluctuations, no frictions, no additional self-propulsion, no
external torque and a vanishing initial angular velocity, i.e., for
$D=D_{r}=\xi=\xi_{r}=v_{0}=M=\dot{\phi}(t=0)=0$ Tsiolkovsky (1903).

Since we assume that the outlet/tank of the particle coincides with the
center-of-mass, the moment of inertia is not affected by the mass ejection and
remains constant. Hence the orientational motion is identical to the case of
time-independent inertia. Clearly, via the mass ejection, the two equations
(5) and (2) are coupled.

Realizations of the rocket-like self-propelled objects can in principle be
found for self-propelled Janus particles in a complex plasma which are laser-
heated such that they evaporate mass in a certain direction Krasheninnikov
_et al._ ; Nosenko  _et al._ (2010, 2020) or even for inflated toy balloons
Müller and Strehlow (1999); Mangan and Destrade (2015) or active granular
particles equipped with compressed air tanks. For the latter, we would expect
the particle loss mass at a rate of approximately
$\dot{m}=-1\leavevmode\nobreak\ \text{g}/\text{s}$ by exhausting air at a
velocity $u=100\leavevmode\nobreak\ \text{mm}/\text{s}$. Initial mass and
moment of inertia are $m_{0}=10\leavevmode\nobreak\ \text{g}$ and
$J=100\leavevmode\nobreak\ \text{g}\,\text{mm}^{2}$. The remaining parameters
are of the order of $\xi=10\leavevmode\nobreak\ \text{g}/\text{s}$,
$\xi_{r}=100\leavevmode\nobreak\ \text{g}\,\text{mm}^{2}/\text{s}$,
$D=100\leavevmode\nobreak\ \text{mm}^{2}/\text{s}$, $D_{r}=1/s$,
$v_{0}=10-100\leavevmode\nobreak\ \text{mm}/\text{s}$ and
$M=10^{-7}\leavevmode\nobreak\ \text{N}\,\text{m}$. We finally remark that
there is some overdamped counterpart of rocket-like motion in the osmoto-
phoresis of semipermeable vesicles Anderson (1983) where the ejection of
molecules out of vesicle body leads to self-propulsion driven by the osmotic
pressure difference Nardi  _et al._ (1999) and for Janus-particles and
nanorockets driven by reactive momentum-transfer Eloul  _et al._ (2020); Li
_et al._ (2016).

###  II.3 Isotropic mass evaporation

A different situation occurs if the mass ejection is not directed but
isotropic as sketched in Figure 2(c) (green background). Imagine a particle
coated with an isotropic layer that evaporates likewise in all directions, as
realizable in dusty plasmas Krasheninnikov  _et al._ ; Nosenko  _et al._
(2010, 2020). In this case, the ejected mass does only carry away the
translational momentum given by $-\dot{m}(t)\,\dot{\bm{\mathrm{R}}}(t)$ such
that the translational equations of motion for this case coincide with Eq. (5)
for $u=0$. However, the mass ejection is radial only in the body frame but for
a rotating particle the angular momentum ${\dot{J}}(t)\,{\dot{\phi}}(t)$ is
taken away in the laboratory frame even in the absence of any torque.
Therefore the orientational equation of motion reads as (2) with $J$ replaced
by $J(t)$, as follows

$J(t)\,\ddot{\phi}(t)=M-\xi_{r}\,\dot{\phi}(t)+\xi_{r}\sqrt{2D_{r}}\,\tau_{\rm
st}(t).$ (6)

This set-up could be realized in experiments by placing a leaking water tank
or evaporating material on an active granular particle. The order of magnitude
of the parameters might be $m_{0}=10\leavevmode\nobreak\ \text{g}$,
$\dot{m}=-1\leavevmode\nobreak\ \text{g}/\text{s}$,
$J_{0}=100\leavevmode\nobreak\ \text{g}\,\text{mm}^{2}$,
$\xi=10\leavevmode\nobreak\ \text{g}/\text{s}$,
$\xi_{r}=100\leavevmode\nobreak\ \text{g}\,\text{mm}^{2}/\text{s}$,
$D=100\leavevmode\nobreak\ \text{mm}^{2}/\text{s}$, $D_{r}=1/s$,
$v_{0}=10-100\leavevmode\nobreak\ \text{mm}/\text{s}$ and
$M=10^{-7}\leavevmode\nobreak\ \text{N}\,\text{m}$. Finally we remark that the
inverse situation of mass adsorption can be treated in a similar way with a
positive sign of $\dot{m}(t)$.

###  II.4 Isotropic shape change

The pirouette of figure skating is an example of a fourth situation where the
total mass $m$ of the body is time-independent but the moment of inertia does
change due to a shape change of the body. In this special case, sketched in
Figure 2(d) (blue background), the shape change does not carry away angular
momentum but the total angular momentum is conserved. Consequently while the
translational equation of motion is identical with Eq. (1), the orientational
equation of motion is given by

$\frac{d}{dt}\Big{(}J(t)\,\dot{\phi}(t)\Big{)}=M-\xi_{r}\,\dot{\phi}(t)+\xi_{r}\sqrt{2D_{r}}\,\tau_{\rm
st}(t).$ (7)

Lastly this model could describe an active granular particle with a stretched
elastic material attached to it. In that way, the initial moment of inertia
could be increased by an order of magnitude $J_{0}=100\leavevmode\nobreak\
\text{g}\,\text{mm}^{2}$, relaxing over a few seconds to its equilibrium shape
with $\dot{J}=-1\leavevmode\nobreak\ \text{g}\,\text{mm}^{2}/\text{s}$. The
order of magnitude of the other parameters might be $m=1\leavevmode\nobreak\
\text{g}$, $\xi=10\leavevmode\nobreak\ \text{g}/\text{s}$,
$\xi_{r}=100\leavevmode\nobreak\ \text{g}\,\text{mm}^{2}/\text{s}$,
$D=100\leavevmode\nobreak\ \text{mm}^{2}/\text{s}$, $D_{r}=1/s$,
$v_{0}=10-100\leavevmode\nobreak\ \text{mm}/\text{s}$ and
$M=10^{-7}\leavevmode\nobreak\ \text{N}\,\text{m}$.

###  II.5 General model

The lesson to be learned from the previous examples is that the equations of
motion depend on the imposed set-up of mass change. In order to proceed in a
general way, we now present a general framework of equations of motion which
accommodates all previous special cases. To define this model as general as
possible, we also assume an effective time-dependent self-propulsion speed
$v_{0}(t)$, a time-dependent internal torque $M(t)$, a time-dependent
translational $\xi(t)$ and rotational friction coefficient $\xi_{r}(t)$, as
well as, time-dependent translational $D(t)$ and rotational diffusion
coefficient $D_{r}(t)$ and a time-dependent mass ejection velocity $u(t)$.

We now consider the following general Langevin equations governing the
translational and the rotational motion for a self-propelled particle

$\displaystyle\frac{d}{dt}\Big{(}m(t)\,\dot{\bm{\mathrm{R}}}(t)\Big{)}=$
$\displaystyle\xi(t)\big{(}v_{0}(t)\,\bm{\hat{\mathbf{n}}}(t)-\dot{\bm{\mathrm{R}}}(t)+\sqrt{2D(t)}\,\bm{\mathrm{f}}_{\rm
st}(t)\big{)}$
$\displaystyle-\dot{m}(t)\big{(}u(t)\,\bm{\hat{\mathbf{n}}}(t)-\dot{\bm{\mathrm{R}}}(t)\big{)},$
(8) $\displaystyle\frac{d}{dt}\Big{(}J(t)\,\dot{\phi}(t)\Big{)}=$
$\displaystyle
M(t)-\xi_{r}(t)\,\dot{\phi}(t)+\xi_{r}(t)\sqrt{2D_{r}(t)}\,\tau_{\rm st}(t)$
$\displaystyle+\nu\,\dot{J}(t)\,\dot{\phi}(t).$ (9)

Clearly, all situations discussed so far and shown in Figure 2 are obtained
from these equations as special cases: of course, 2(a) is the special limit
where the parameters $m$, $J$, $\xi$, $\xi_{r}$ $D$, $D_{r}$ $v_{0}$ and $M$
are constant. The rocket set-up in Figure 2(b) coincides with the general
equations (8) and (9) when the parameters $J$, $\xi$, $\xi_{r}$ $D$, $D_{r}$
$v_{0}$, $M$, the relative velocity $u$ are constant. The isotropic mass
evaporation (Figure 2(c)) is contained in (8) and (9) when the parameters
$\xi$, $\xi_{r}$ $D$, $D_{r}$, $v_{0}$, $M$ are constant, the relative
velocity vanishes $u=0$ and $\nu=1$. Finally the equations for an isotropic
shape change (Figure 2(d)) follow when in (8) and (9) the parameters $m$,
$\xi$, $\xi_{r}$ $D$, $D_{r}$, $v_{0}$, $M$ are constant and $\nu=0$.

At this stage we remark that more realistic situations can also be
accommodated into the general equations (8) and (9). These include, for
example, a rocket where the outlet of the mass ejection does not coincide with
the center-of-mass or where the ejection direction is not parallel to the
particle orientation Rankin and Rosenhead (1949).

From a mathematical point of view the equations of motion (8) and (9) are
stochastic differential equations with Gaussian noise. The rotational equation
(9) is linear so that the distribution of the angle and angular velocity is
Gaussian for any time. We give the corresponding general solutions of (8) and
(9) in the Appendix A.

##  III  Time-independent inertia

We now turn to the special case of time-independent parameters defined by
equations (1) and (2). These equations of motion were studied before in
references Scholz  _et al._ (2018a); Löwen (2020); Sandoval (2020). Here we
summarize essential known results but also provide new analytical results for
the full time resolved mean displacement, velocity correlation function and
mean-square displacement. In doing so, we first consider the noise-free case
and then include effects of noise.

###  III.1 Results for vanishing noise

For given initial orientations $\phi_{0}=\phi(0)$ and angular velocities
$\dot{\phi}_{0}=\dot{\phi}(0)$ at time $t=0$, the deterministic solution of
the general orientational equation of motion (2) in the absence of noise is

$\phi(t)=\phi_{0}+\omega
t+\frac{\dot{\phi}_{0}-\omega}{\gamma_{r}}\left(1-e^{-\gamma_{r}t}\right).$
(10)

with the spinning frequency $\omega=M/\xi_{r}$ and rotational damping rate
$\gamma_{r}=\xi_{r}/J$. Plugging this solution into the noise-free
translational equation (1), we obtain for given initial positions
$\bm{\mathrm{R}}_{0}=\bm{\mathrm{R}}(0)$ and velocities
$\dot{\bm{\mathrm{R}}}_{0}=\dot{\bm{\mathrm{R}}}(0)$ at time $t=0$ the
particle velocity

$\displaystyle\dot{\bm{\mathrm{R}}}(t)=$
$\displaystyle\dot{\bm{\mathrm{R}}}_{0}e^{-\gamma
t}+v_{0}\bm{\mathrm{\hat{P}}}\Big{[}\tilde{\gamma}\,\big{(}i\theta\big{)}^{\tilde{\gamma}+i\tilde{\omega}}e^{i(\phi_{0}+\theta)}$
(11)
$\displaystyle\times\Gamma\big{(}-(\tilde{\gamma}+i\tilde{\omega}),i\theta
e^{-\gamma_{r}t},i\theta\big{)}\Big{]}e^{-\gamma t},$

where we introduced the translational damping rate $\gamma=\xi/m$ and the
notations $\tilde{\gamma}=\gamma/\gamma_{r}$,
$\tilde{\omega}=\omega/\gamma_{r}$,
$\theta=(\dot{\phi}_{0}-\omega)/\gamma_{r}$. Moreover, $\Gamma(s,x_{1},x_{2})$
denotes the generalized Gamma function Paris (2010),

$\Gamma(s,x_{1},x_{2})=\int_{x_{1}}^{x_{2}}d{t}\,t^{s-1}e^{-t},$ (12)

and the operator $\bm{\mathrm{\hat{P}}}$ formally transforms a complex number
$z$ into its two-dimensional vector $(\operatorname{Re}z,\operatorname{Im}z)$
in the complex plane. This results in the particle trajectory

$\displaystyle\bm{\mathrm{R}}(t)=$
$\displaystyle\bm{\mathrm{R}}_{0}+\frac{\dot{\bm{\mathrm{R}}}_{0}}{\gamma}\big{(}1-e^{-\gamma
t}\big{)}$ (13)
$\displaystyle+\frac{v_{0}}{\gamma_{r}}\bm{\mathrm{\hat{P}}}\Big{[}\big{(}i\theta\big{)}^{i\tilde{\omega}}e^{i(\phi_{0}+\theta)}\Gamma\big{(}-i\tilde{\omega},i\theta
e^{-\gamma_{r}t},i\theta\big{)}\Big{]}$
$\displaystyle-\frac{v_{0}}{\gamma_{r}}\bm{\mathrm{\hat{P}}}\Big{[}$
$\displaystyle\big{(}i\theta\big{)}^{\tilde{\gamma}+i\tilde{\omega}}e^{i(\phi_{0}+\theta)}\Gamma\big{(}-(\tilde{\gamma}+i\tilde{\omega}),i\theta
e^{-\gamma_{r}t},i\theta\big{)}\Big{]}e^{-\gamma t},$

In the limit of long times, the angular velocity reaches the spinning
frequency, $\lim_{t\to\infty}{\dot{\phi}}(t)=\omega$ so that the particle is
rotating with this frequency around a circle of radius

$r=\frac{v_{0}}{\omega}\sqrt{\frac{\gamma^{2}}{\gamma^{2}+\omega^{2}}},$ (14)

centered at the position

$\bm{\mathrm{R}}_{c}=\bm{\mathrm{R}}_{0}+\frac{\dot{\bm{\mathrm{R}}}_{0}}{\gamma}+\frac{v_{0}}{\gamma_{r}}\bm{\mathrm{\hat{P}}}\Big{[}\big{(}i\theta\big{)}^{i\tilde{\omega}}e^{i(\phi_{0}+\theta)}\Gamma\big{(}-i\tilde{\omega},0,i\theta\big{)}\Big{]}.$
(15)

Clearly, the spinning frequency $\omega$ does not depend on any inertia.
However, the circle radius $r$ depends on the mass $m$ via the translational
damping rate $\gamma$ due to the centrifugal force but is independent on the
moment of inertia $J$. The center of the circle depends on
$\bm{\mathrm{R}}_{0}$, $\dot{\bm{\mathrm{R}}}_{0}$, $\phi_{0}$,
$\dot{\phi}_{0}$, demonstrating that for vanishing noise even the long-time
limit may depend on the initial conditions. Finally, in the overdamped limit
of vanishing inertia, the results reduce to that of Brownian circle swimmers
van Teeffelen and Löwen (2008); Kurzthaler and Franosch (2017).

###  III.2 Effect of Brownian noise

Subjected to Brownian noise, the particle will relax to a steady state after a
long time forgetting about its initial conditions $\bm{\mathrm{R}}_{0}$,
$\dot{\bm{\mathrm{R}}}_{0}$, $\phi_{0}$, $\dot{\phi}_{0}$ at time $t=0$. The
static and dynamical correlation in the steady state can be calculated as a
time average over a very long time window which we shall denote with brackets
$\langle\dots\rangle$. In the sequel we shall consider several of such
dynamical correlations. In the steady state, one can also calculate
conditional averages. For example one can built dynamical averages in the
steady state after a lag time under the condition that the particle’s position
and orientation is prescribed at an initial time. We shall compile analytical
results for the different correlation functions first and show also examples
for numerical evaluations of the analytical formula.

####  III.2.1 Velocity correlation function

First we introduce the translational velocity correlation function Hansen and
McDonald (1990),

$Z(t)=\langle\dot{\bm{\mathrm{R}}}(t)\cdot\dot{\bm{\mathrm{R}}}(0)\rangle,$
(16)

where $t$ now denotes a lag time and $\dot{\bm{\mathrm{R}}}(0)$ is taken from
the velocity distribution in the steady state. We remark that the latter was
computed recently for small inertia Herrera and Sandoval (2021) and for the
formally equivalent model of an overdamped particle in a harmonic potential
Malakar  _et al._ (2020). The velocity distribution is non-Gaussian (i.e. non
Maxwellian) and its second moment
$Z(0)=\langle\dot{\bm{\mathrm{R}}}(0)\cdot\dot{\bm{\mathrm{R}}}(0)\rangle$
which is proportional to the the mean kinetic energy is known analytically
Scholz  _et al._ (2018a) as

$Z(0)=2D\gamma+v_{0}^{2}\operatorname{Re}\Big{[}\tilde{\gamma}e^{\tilde{D}_{r}}\tilde{D}_{r}^{-\Omega_{+}}\,\Gamma(\Omega_{+},0,\tilde{D}_{r})\Big{]},$
(17)

where we introduced $\tilde{D}_{r}=D_{r}/\gamma_{r}$ and
$\Omega_{\pm}=\big{(}D_{r}\pm(\gamma+i\omega)\big{)}/\gamma_{r}$. For an
active inertial particle considered here, we have obtained the analytical
result

$Z(t)=2D\gamma e^{-\gamma
t}+\frac{v_{0}}{2}\big{(}\langle\dot{\bm{\mathrm{R}}}(t)\cdot\bm{\hat{\mathbf{n}}}(0)\rangle+\langle\dot{\bm{\mathrm{R}}}(0)\cdot\bm{\hat{\mathbf{n}}}(t)\rangle\big{)}$
(18)

with

$\displaystyle\langle\dot{\bm{\mathrm{R}}}(t)\cdot\bm{\hat{\mathbf{n}}}(0)\rangle=$
$\displaystyle
v_{0}\operatorname{Re}\Big{[}\tilde{\gamma}e^{\tilde{D}_{r}}\Big{(}\tilde{D}_{r}^{-\Omega_{-}}\Gamma(\Omega_{-},\tilde{D}_{r}e^{-\gamma_{r}t},\tilde{D}_{r})$
$\displaystyle\qquad\,\,\,\,+\tilde{D}_{r}^{-\Omega_{+}}\Gamma(\Omega_{+},0,\tilde{D}_{r})\Big{)}e^{-\gamma
t}\Big{]},$ (19)

and

$\langle\dot{\bm{\mathrm{R}}}(0)\cdot\bm{\hat{\mathbf{n}}}(t)\rangle=v_{0}\operatorname{Re}\Big{[}\tilde{\gamma}e^{\tilde{D}_{r}}\tilde{D}_{r}^{-\Omega_{+}}\Gamma(\Omega_{+},0,\tilde{D}_{r}e^{-\gamma_{r}t})e^{\gamma
t}\Big{]}.$ (20)

which implies that the long-time behavior of $Z(t)$ is exponential in time.

####  III.2.2 Orientation correlation function

Similarly the dynamical orientational correlation function
$C(t)=\langle\bm{\hat{\mathbf{n}}}(t)\cdot\bm{\hat{\mathbf{n}}}(0)\rangle$ in
the steady state can be expressed analytically as a double exponential as

$\displaystyle C(t)=\cos(\omega
t)\,e^{-D_{r}(t-\gamma_{r}^{-1}(1-e^{-\gamma_{r}t}))},$ (21)

which was found previously in another context by Ghosh and coworkers Ghosh
_et al._ (2015) for $\omega=0$ and for general $\omega$ in Ref. Scholz  _et
al._ (2018a). Again it decays exponentially in time for long times. A
characteristic orientational persistence time $\tau_{p}$ can be determined as

$\tau_{p}=\int_{0}^{\infty}\\!C(t)dt=\frac{1}{D_{r}}\operatorname{Re}\Big{[}\tilde{D}_{r}e^{\tilde{D}_{r}}\tilde{D}_{r}^{-\Omega}\,\Gamma(\Omega,0,\tilde{D}_{r})\Big{]}.$
(22)

with $\Omega=(D_{r}-i\omega)/\gamma_{r}$. For vanishing inertia we recover the
known result of the persistence time $\tau_{p}=D_{r}/(D_{r}^{2}+\omega^{2})$
van Teeffelen and Löwen (2008); Kurzthaler and Franosch (2017) which
simplifies further to the standard result $\tau_{p}=1/D_{r}$ for linear
swimmer Howse  _et al._ (2007).

####  III.2.3  Mean displacement

Figure 3:  Mean displacement $\langle\Delta\bm{\mathrm{R}}(t)\rangle$ in the
$xy$-plane for a chiral particle with initial orientation along the $x$-axis
for different moment of inertia, $J=0.1\leavevmode\nobreak\ \xi_{r}/D_{r}$
(orange), $J=1\leavevmode\nobreak\ \xi_{r}/D_{r}$ (red) and
$J=10\leavevmode\nobreak\ \xi_{r}/D_{r}$ (purple). Lengths are given in units
of $l_{p}=v_{0}/D_{r}$. The parameters are $\omega=4\leavevmode\nobreak\
D_{r}$, $m=0.1\leavevmode\nobreak\ \xi/D_{r}$. The starting point at $t=0$ is
denoted by a black dot. The spira mirabilis of the overdamped limit is plotted
on the left (black) for comparison.

Next, we address the mean displacement
$\langle\Delta\bm{\mathrm{R}}(t)\rangle=\langle\bm{\mathrm{R}}(t)-\bm{\mathrm{R}}_{0}\rangle$
of the particle in the steady state as a function of time $t$. The average is
now taken in the steady state but under the condition that for the initial
time $t=0$, the position $\bm{\mathrm{R}}(0)=\bm{\mathrm{R}}_{0}$ and the
orientation $\hat{\bm{\mathrm{n}}}(0)$ (embodied in $\phi(0)=\phi_{0}$) are
prescribed. Since the particle velocities and the orientations are correlated
in the steady state, the average over the translational velocity
$\langle\dot{\bm{\mathrm{R}}}(0)\rangle$ is not vanishing due to the
prescribed orientation $\hat{\bm{\mathrm{n}}}(0)$. This average is given by

$\langle\dot{\bm{\mathrm{R}}}(0)\rangle=v_{0}\bm{\mathrm{\hat{P}}}\Big{[}\tilde{\gamma}e^{\tilde{D}_{r}}\tilde{D}_{r}^{-\Omega_{+}}\Gamma(\Omega_{+},0,\tilde{D}_{r})e^{i\phi_{0}}\Big{]}.$
(23)

We obtain for the mean displacement

$\displaystyle\langle\Delta\bm{\mathrm{R}}(t)\rangle=$
$\displaystyle\frac{\langle\dot{\bm{\mathrm{R}}}(0)\rangle}{\gamma}\Big{(}1-e^{-\gamma
t}\Big{)}$ (24)
$\displaystyle+\frac{v_{0}}{D_{r}}\bm{\mathrm{\hat{P}}}\Big{[}\tilde{D}_{r}e^{\tilde{D}_{r}}\Big{(}\tilde{D}_{r}^{-\Omega}\Gamma(\Omega,\tilde{D}_{r}e^{-\gamma_{r}t},\tilde{D}_{r})$
$\displaystyle+\tilde{D}_{r}^{-\Omega_{-}}\Gamma(\Omega_{-},\tilde{D}_{r}e^{-\gamma_{r}t},\tilde{D}_{r})e^{-\gamma
t}\Big{)}e^{i\phi_{0}}\Big{]}.$

For short times $t$, the particle moves proceeds on average ballistically
(i.e. linearly in time) with

$\langle\Delta\bm{\mathrm{R}}(t)\rangle=\langle\dot{\bm{\mathrm{R}}}(0)\rangle
t+\mathcal{O}(t^{2}).$ (25)

Then the rotational noise decorrelates the current orientation from the
initial orientation and the mean displacement saturates to a finite
persistence length
$\bm{\mathrm{L}}_{p}=\lim_{t\to\infty}\langle\Delta\bm{\mathrm{R}}(t)\rangle$
given as

$\bm{\mathrm{L}}_{p}=\frac{\langle\dot{\bm{\mathrm{R}}}(0)\rangle}{\gamma}+\frac{v_{0}}{D_{r}}\bm{\mathrm{\hat{P}}}\Big{[}\tilde{D}_{r}e^{\tilde{D}_{r}}\tilde{D}_{r}^{-\Omega}\,\Gamma(\Omega,0,\tilde{D}_{r})e^{i\phi_{0}}\Big{]}.$
(26)

In case of a vanishing spinning frequency ($\omega=0$), the persistence length
simplifies to
$\bm{\mathrm{L}}_{p}=\langle\dot{\bm{\mathrm{R}}}(0)\rangle/\gamma+v_{0}\tau_{p}\bm{\hat{\mathbf{n}}}_{0}$
with $\tau_{p}$ given by (22). In the overdamped limit we obtain the standard
results of the persistence length for linear microswimmers
$\bm{\mathrm{L}}_{p}=v_{0}\bm{\hat{\mathbf{n}}}_{0}/D_{r}$ ($\omega=0$) Howse
_et al._ (2007). Moreover, for overdamped circle swimmer, the full time-
resolved mean displacement given by (24) simplifies to a spira mirabilis van
Teeffelen and Löwen (2008); Kümmel  _et al._ (2013). The presence of inertia
will distort the ideal spira mirabilis give rise to a more complex mean
trajectory. This is shown in Figure 3 where three shapes of the mean
trajectory for increasing moment of inertia $J$ are compared to the overdamped
case. Increasing $J$ reduces effectively the role of fluctuations such that
there are more turns until the particle reaches half of the distance to its
final fixpoint. Even though $\bm{\hat{\mathbf{n}}}(0)$ is oriented towards the
positive $x$-axis in all cases, the inertial mean trajectory first
”oversteers” the initial orientation due to the velocity average, an effect
which we shall elaborate and quantify further in subchapter 5.

####  III.2.4 Mean-square displacement

The full time-resolved mean-square-displacement (MSD) can be calculated as

$\langle\Delta\bm{\mathrm{R}}^{2}(t)\rangle=4D_{L}t+\frac{2}{\gamma^{2}}\left(Z(t)-Z(0)\right)+2\frac{v_{0}^{2}}{\gamma_{r}^{2}}F(t),$
(27)

with the long-time diffusion coefficient

$D_{L}=D+\frac{v_{0}^{2}}{2D_{r}}\operatorname{Re}\Big{[}\tilde{D}_{r}e^{\tilde{D}_{r}}\tilde{D}_{r}^{-\Omega}\,\Gamma(\Omega,0,\tilde{D}_{r})\Big{]},$
(28)

and the function

$\displaystyle F(t)=\operatorname{Re}\Bigg{[}$
$\displaystyle\frac{e^{\tilde{D}_{r}}}{\Omega^{2}}\Bigg{(}{}_{2}F_{2}{\left[\genfrac{.}{.}{0.0pt}{}{\Omega{\mathchar
44\relax}\mskip 8.0mu\Omega}{\Omega+1{\mathchar 44\relax}\mskip
8.0mu\Omega+1};-\tilde{D}_{r}\right]}$ (29)
$\displaystyle-{}_{2}F_{2}{\left[\genfrac{.}{.}{0.0pt}{}{\Omega{\mathchar
44\relax}\mskip 8.0mu\Omega}{\Omega+1{\mathchar 44\relax}\mskip
8.0mu\Omega+1};-\tilde{D}_{r}e^{-\gamma_{r}t}\right]}e^{-\gamma_{r}\Omega
t}\Bigg{)}\Bigg{]},$

where ${}_{p}F_{q}$ represents the generalized hypergeometric function Askey
and Daalhuis (2010).

Figure 4:  Mean-square-displacement as a function of time on a double-
logarithmic plot. (a) For an achiral particle with fixed mass $m$ and varied
moment of inertia $J$. (b) For a chiral particle with fixed mass $m$ and
varied moment of inertia $J$. The fixed parameters are
$m=0.1\leavevmode\nobreak\ \xi/D_{r}$, $D=0$ and the moment of inertia is
$J=0.1\leavevmode\nobreak\ \xi_{r}/D_{r}$ (orange), $J=1\leavevmode\nobreak\
\xi_{r}/D_{r}$ (red) or $J=10\leavevmode\nobreak\ \xi_{r}/D_{r}$ (purple). (c)
For an achiral particle with varied $m$ and fixed $J$. (d) For a chiral
particle with varied $m$ and fixed $J$. Here, the fixed parameters are
$J=0.1\leavevmode\nobreak\ \xi_{r}/D_{r}$, $D=0$ and the mass is
$m=0.1\leavevmode\nobreak\ \xi/D_{r}$ (orange), $m=1\leavevmode\nobreak\
\xi/D_{r}$ (red) or $m=10\leavevmode\nobreak\ \xi/D_{r}$ (purple).

Figures 4(a)-4(d) compare the temporal behavior of the mean-square
displacement of an achiral particle to that of a chiral particle for different
mases and moments of inertia $J$. All curves exhibit the characteristic
crossover from a short-time ballistic behavior

$\langle\Delta\bm{\mathrm{R}}^{2}(t)\rangle=Z(0)\,t^{2}+\mathcal{O}(t^{3}),$
(30)

to the long-time diffusive behavior governed by

$\langle\Delta\bm{\mathrm{R}}^{2}(t)\rangle\sim 4D_{L}t.$ (31)

In the limit of small $J$, the short-time ballistic dynamics is

$\lim_{J\to
0}\langle\Delta\bm{\mathrm{R}}^{2}(t)\rangle=\Big{(}2D\gamma+v_{0}^{2}\,\frac{\gamma(\gamma+D_{r})}{(\gamma+D_{r})^{2}+\omega^{2}}\Big{)}t^{2}+\mathcal{O}(t^{3}),$
(32)

while for large $J$, we have

$\lim_{J\to\infty}\langle\Delta\bm{\mathrm{R}}^{2}(t)\rangle=\Big{(}2D\gamma+v_{0}^{2}\,\frac{\gamma^{2}}{\gamma^{2}+\omega^{2}}\Big{)}t^{2}+\mathcal{O}(t^{3}).$
(33)

In general, the long-time diffusion coefficient $D_{L}$ (see Eq. (28)) can be
represented as

$D_{L}=D+\frac{v_{0}^{2}}{2}\tau_{\text{p}},$ (34)

where the first term in Eq. (22) captures the diffusive behavior of a passive
particle and the second is consistent with the standard picture of a typical
jump length of $v_{0}\tau_{p}$ and a typical jump time of $\tau_{p}$, similar
to the overdamped expression of microswimmers when $\omega=0$ Howse  _et al._
(2007). It was emphasized in Ref. Scholz  _et al._ (2018a) that $D_{L}$
depends on the moment of inertia $J$ but not explicitly on the mass $m$.

In case of small moments of inertia, the long-time diffusion coefficient of
the circle flyer asymptotically goes to Scholz  _et al._ (2018a)

$D_{L}=D+\frac{v_{0}^{2}}{2}\frac{D_{r}}{D_{r}^{2}+\omega^{2}}\left(1+\frac{D_{r}}{\xi_{r}}J\right)+\mathcal{O}(J^{2}),$
(35)

which grows dominantly proportional to moment of inertia. The asymptotic
behavior of the long-time diffusion coefficient for large moments of inertia
is Scholz  _et al._ (2018a),

$D_{L}\sim\begin{cases}v_{0}^{2}\sqrt{\frac{\pi}{8D_{r}\xi_{r}}}\sqrt{J},&(\omega=0)\\\
D,&(\omega\neq 0)\end{cases}$ (36)

As the moment of inertia grows for $\omega\not=0$, the activity-induced part
of the long-time diffusion coefficient goes asymptotically to zero (see Eq.
(36)) since diffusion is hampered by systematic circling motion, i.e., the
particle is get trapped in a circular path due to its huge moment of inertia.

Figures 4(a) and 4(c) show data for an achiral swimmer with different moments
of inertia $J$ and different masses. The short-time ballistic prefactor is
pretty independent of $J$ but decreases with increasing $m$. The latter trend
follows from the fact that for large $m$ the particle cannot accelerate
towards its self-propulsion velocity $v_{0}$. Conversely, the long-time
diffusivity is also increasing with $J$ according to (35) and (36) as the
persistence in orientation increases with $J$ but it is independent of $m$.
For a chiral particle, shown in Figures 4(b) and 4(d), the MSD exhibits
wiggles due to the circling.

An immediate consequence of (35) and (36) is that the long-time diffusivity
behaves non-monotonic in $J$. Explicit data are presented in Figure 5 which
illustrates the non-monotonic dependence of $D_{L}$ on the moment of inertia
$J$ for different spinning frequencies $\omega$ for the special case $D=0$.
There is an intermediate maximum in $D_{L}$ which is indicated in Figure 5 by
a red point. This peak could be exploited for an optimal exploration of an
unknown territory by adapting the moment of inertia accordingly. The
associated optimal moment of inertia is plotted as a function of the spinning
frequency $\omega$ in the inset of Figure 5.

Figure 5:  Long-time diffusion coefficient $D_{L}$ as a function of the moment
of inertia $J$ for different circling frequencies $\omega=10D_{r}$,
$\omega=1D_{r}$, $\omega=0.1D_{r}$ and $\omega=0$. The translational diffusion
coefficient was set to zero, $D=0$. In the inset, the global maximum point
$J_{\text{max}}$ for a given circling frequency $\omega$ is depicted. The
corresponding maximal value $D_{L}(J_{\text{max}})$ is shown as a red dot in
the main figure.

####  III.2.5  Delay function

Opposed to the overdamped case, the velocity of an inertial particle does not
coincide with its self-propulsion direction and in Ref. Scholz  _et al._
(2018a) a dynamical correlation function, referred to as delay function
$d(t)$, was introduced to quantify the delay between the velocity and
orientation dynamics

$d(t)=\langle\dot{\bm{\mathrm{R}}}(t)\cdot\bm{\hat{\mathbf{n}}}(0)\rangle-\langle\dot{\bm{\mathrm{R}}}(0)\cdot\bm{\hat{\mathbf{n}}}(t)\rangle.$
(37)

The ”mixed” difference ensures that this function is trivially zero in the
overdamped limit but when non-zero its sign contains valuable information
about the delay process between $\bm{\hat{\mathbf{n}}}(t)$ and
$\dot{\bm{\mathrm{R}}}(t)$. If for example $d(t)$ is positive, this means that
- on average - first the particle orientation changes and then the velocity
will follow that change after a time $t$. A positive $d(t)$ is the standard
behavior exploited by the oversteering of racing cars which is also expected
for achiral particles. The full analytical result for $d(t)$ directly follows
from (19) and (20) and was given in Ref. Scholz  _et al._ (2018a). Most
notably, for an achiral particle, $d(t)$ has a positive peak after a typical
delay time, while for a chiral particle, $d(t)$ oscillates due to the
systematic change in orientation. The latter oscillation was recently observed
in macroscopic whirligig beetles swimming at the water surface Turner .

Here we also provide analytical limits of small and large moments of inertia
$J$. For small $J$ we get

$d(t)=2v_{0}A(t)\left(1+\frac{D_{r}}{\xi_{r}}J\right)+\mathcal{O}(J^{2}),$
(38)

with

$\displaystyle A(t)=$
$\displaystyle\frac{\gamma\,D_{r}\,(\gamma^{2}-D_{r}^{2}-\omega^{2})(\cos(\omega
t)\,e^{-D_{r}t}-e^{-\gamma
t})}{\big{(}(\gamma+D_{r})^{2}+\omega^{2}\big{)}\big{(}(\gamma-
D_{r})^{2}+\omega^{2}\big{)}}$
$\displaystyle+\frac{\gamma\,\omega\,(\gamma^{2}+D_{r}^{2}+\omega^{2})\sin(\omega
t)\,e^{-D_{r}t}}{\big{(}(\gamma+D_{r})^{2}+\omega^{2}\big{)}\big{(}(\gamma-
D_{r})^{2}+\omega^{2}\big{)}}.$ (39)

Since $A(t)$ is positive for small times $t$, the delay effect is enhanced for
increasing $J$. Moreover the oscillatory behavior of chiral particle can be
seen here directly. In the opposite limit of large moment of inertia, the
inertial delay approaches

$\lim_{J\to\infty}d(t)=2v_{0}\frac{\gamma\,\omega}{\gamma^{2}+\omega^{2}}\sin(\omega
t),$ (40)

independent of the rotational diffusion coefficient $D_{r}$ which documents
again the oscillatory behavior for chiral particles.

##  IV  Time-dependent inertia

Here we study the effect of time-dependent inertia on the Langevin motion of
underdamped particle. We first introduce a reduced Langevin rocket model in
which the mass of the particle gets burned to accelerate the particle giving
rise to a time-dependent mass and propulsion speed. Then we compare the four
different set-ups introduced in chapter II. Last, we consider the limiting
case of slowly varying parameters with respect to time.

###  IV.1  Langevin rocket

We define the “Langevin rocket” model by including orientational fluctuations
in the traditional Tsiolkovsky rocket equation Tsiolkovsky (1903). The effect
of noise on rocket motion has been considered previously, see e.g. Srivastava
_et al._ (2012), but a simple basis reference model for that is missing. We
therefore simplify the general Eqs. (5) and (2) for directed mass ejection and
assume vanishing moment of inertia, torque and translational diffusion ($J=0$,
$M=0$ and $D=0$) The Langevin rocket dynamics for a prescribed $m(t)$ is then
given by

$\displaystyle
m(t)\,\ddot{\bm{\mathrm{R}}}(t)+\xi\,\dot{\bm{\mathrm{R}}}(t)=-u\,\dot{m}(t)\,\bm{\mathrm{n}}(t),$
(41) $\displaystyle\dot{\phi}(t)=\sqrt{2D_{r}}\,\tau_{\rm st}(t).$ (42)

This set of equations approaches the ideal Tsiolkovsky rocket equation,
$m(t)\,\ddot{\bm{\mathrm{R}}}(t)=-u\,\dot{m}(t)\,\bm{\mathrm{n}}_{0}$, in the
limit of vanishing damping $(\xi=0)$ and noise ($D_{r}=0$) Tsiolkovsky (1903).

For the sake of simplicity, we assume that the rocket is ejecting mass at a
constant rate $(m_{\infty}-m_{0})/\Delta t$ where $m_{0}$ denotes the initial
mass, $m_{\infty}$ the final rest mass of the rocket and $\Delta t$ is the
total burn time. The ejection process happens in the window $0<t<\Delta t$
such that the time dependent particle mass is

$m(t)=m_{0}+(m_{\infty}-m_{0})\frac{\min(t,\Delta t)}{\Delta t}.$ (43)

In the following, we discuss the average reach of the rocket (i.e. its mean
displacement) as a function time. In particular, we investigate the final
reach for long times as a function of the burn time $\Delta t$ and the
propellant mass fraction

$\zeta=\frac{m_{0}-m_{\infty}}{m_{0}}.$ (44)

####  IV.1.1 Results for vanishing noise

In the absence of rotational noise, the displacement of the rocket for a
vanishing initial velocity at $t=0$ is

$\displaystyle\Delta\bm{\mathrm{R}}(t)=u\frac{\min(t,\Delta
t)}{S_{1}+1}\bm{\hat{\mathbf{n}}}_{0}-\frac{u}{\gamma_{0}}\frac{m(t)}{m_{0}}\frac{1-\big{(}\textstyle\frac{m(t)}{m_{0}}\big{)}^{S_{1}}}{S_{1}+1}\bm{\hat{\mathbf{n}}}_{0}$
(45)
$\displaystyle+\frac{u}{\gamma_{0}}\frac{m(t)}{m_{0}}\frac{1-\big{(}\textstyle\frac{m(t)}{m_{0}}\big{)}^{S_{1}}}{S_{1}}\Bigg{(}1-e^{-\gamma_{\infty}(\max(t,\Delta
t)-\Delta t)}\Bigg{)}\bm{\hat{\mathbf{n}}}_{0},$

with the initial damping rate $\gamma_{0}=\xi/m_{0}$, the final damping rate
$\gamma_{\infty}=\xi/m_{\infty}$ and the reduced burn time
$S_{1}=\gamma_{0}\Delta t/\zeta$.

For short times, the rocket exhibits an acceleration by ejecting mass such
that the displacement scales with $t^{2}$,

$\Delta\bm{\mathrm{R}}(t)=\frac{u\zeta}{2\Delta
t}\bm{\hat{\mathbf{n}}}_{0}\,t^{2}+\mathcal{O}\left(t^{3}\right).$ (46)

After the burn time $\Delta t$, the rocket reaches its maximal velocity which
is subsequently exponentially damped with the final damping rate
$\gamma_{\infty}$ until the rocket comes to a standstill. The total long-time
displacement $\Delta
R_{\infty}=\lim_{t\to\infty}\lvert\Delta\bm{\mathrm{R}}(t)\rvert$ is given by

$\Delta R_{\infty}=\frac{u\Delta
t}{S_{1}+1}+\frac{u}{\gamma_{0}}\frac{(1-\zeta)\big{(}1-(1-\zeta)^{S_{1}}\big{)}}{S_{1}(S_{1}+1)}.$
(47)

In Figure 6, we show the long-time displacement $\Delta R_{\infty}$ as a
function of the propellant mass fraction $\zeta$ for different burn times
$\Delta t$. For long burn times $\Delta t\gg 1/\gamma_{0}$, the ultimate
displacement increases linearly with the propellant mass fraction $\Delta
R_{\infty}\sim u\,\zeta/\gamma_{0}$. The rocket reaches the longest distance
when $\zeta_{\text{max}}\sim 1$, a situation which can be called complete
extended mass ejection. Interestingly, however, there is a qualitative
different behavior for burn times which are comparable or smaller than the
characteristic damping time $1/\gamma_{0}$ where the displacements behave non-
monotonic in the mass fraction $\zeta$. This can be intuitively understood as
follows: for small mass fractions, more ejection means more propulsion and
acceleration such that $\Delta R_{\infty}$ increases with $\zeta$. Conversely
for $\zeta$ close to 1, the rocket becomes very light after the burn time and
therefore very quickly stops within an extremely short damping time
$1/\gamma_{\infty}$ which reduces its reach relative to a situation of smaller
$\zeta$. Consequently the optimal value $\zeta_{\text{max}}$ for the mass
ratio for which the reach is maximal is smaller than $1$. These corresponding
optimal mass ratios are marked by red points in Figure 6 and plotted as a
function of the reduced burn time in the inset. For decreasing burn times the
optimal mass ratio $\zeta_{\text{max}}$ exhibits a bifurcation-like behavior
from complete mass ejection to a finite fraction with the special limit of
$\zeta_{\text{max}}\sim 1-e^{-1}\approx 0.63$ as $\Delta t$ approaches zero.

Figure 6:  Maximal reach $\Delta R_{\infty}$ as a function of the propellant
mass fraction $\zeta$ for several burn times $\Delta t=10^{-2}/\gamma_{0}$,
$\Delta t=1/\gamma_{0}$ and $\Delta t=10^{2}/\gamma_{0}$. The inset depicts
the optimal propellant mass fraction $\zeta_{\text{max}}$ for a given burn
time $\Delta t$. The corresponding maximal value
$R_{\infty}(\zeta_{\text{max}})$ is shown as a red dot in the main figure.

The special limit of $\Delta t\ll 1/\gamma_{0}$ deserves some more attention.
In this case of fractional instantaneous mass ejection, the particle ejects
only a fraction of its propellant to gain momentum very quickly. But it keeps
a rest mass in order to still proceed during the subsequent damping time. In
this limit we obtain

$\Delta\bm{\mathrm{R}}(t)=-\frac{u}{\gamma_{0}}(1-\zeta)\ln\big{(}1-\zeta\big{)}\left(1-e^{-\gamma_{\infty}t}\right)\bm{\hat{\mathbf{n}}}_{0},$
(48)

which scales for $\Delta t\ll t\ll 1/\gamma_{0}$ linearly in time

$\Delta\bm{\mathrm{R}}(t)=-u\,\ln\big{(}1-\zeta\big{)}\,\bm{\hat{\mathbf{n}}}_{0}\,t+\mathcal{O}\left(t^{2}\right),$
(49)

For long times, $t\gg 1/\gamma_{0}$, we obtain

$\Delta R_{\infty}=-\frac{u}{\gamma_{0}}(1-\zeta)\ln(1-\zeta).$ (50)

We finally remark that one can consider a full optimization problem with
respect to both burn time $\Delta t$ and the mass fraction $\zeta$ by posing
the question: what is the maximal reach of the rocket if the burn time $\Delta
t$ and the mass fraction $\zeta$ can be varied freely and independently? The
answer in the fluctuation-free case is simple: the best strategy is to burn
all mass $\zeta_{\text{max}}\to 1$ and do this over a very long time $\Delta
t\to\infty$. Then one achieves the maximum

$\max(\Delta R_{\infty})=\frac{u}{\gamma_{0}},$ (51)

in the upper right corner of Figure 6. In other terms, the strategy of
complete extended mass ejection always outperforms that of an fractional
instantaneous mass ejection. This simple answer will change if orientational
noise is included, a case which we shall address next.

####  IV.1.2 Noise-averaged mean reach and noise-induced transition between
two mass ejection strategies

Figure 7:  (a) The optimal mean reach $\max\big{(}\overline{\Delta
R_{\infty}}\big{)}$ maximized with respect to the propellant mass fraction
$\zeta$ and the burn time $\Delta t$ as a function of the rotational noise
strength $D_{r}/\gamma_{0}$. (b) Optimal burn time $\Delta t_{\text{max}}$ as
a function of the rotational noise strength $D_{r}/\gamma_{0}$. (c) Optimal
propellant mass fraction $\zeta_{\text{max}}$ as a function of the rotational
noise strength $D_{r}/\gamma_{0}$. The transition from the complete extended
mass ejection strategy to that of the fractional instantaneous mass ejection
is marked by vertical black lines at $D_{r,\,\text{crit}}\approx
0.72\,\gamma_{0}$ in all three figures.

In case of finite rotational noise $(D_{r}>0$), we obtain for the noise-
averaged displacement of the Langevin rocket the analytical result

$\displaystyle\overline{\Delta\bm{\mathrm{R}}(t)}=\frac{u}{D_{r}}\frac{1}{S_{1}+1}\left(1-e^{-D_{r}\min(t,\Delta
t)}\right)\bm{\hat{\mathbf{n}}}_{0}$ (52)
$\displaystyle+\frac{u}{D_{r}}\operatorname{Re}\bigg{[}e^{S_{2}}\big{(}-S_{2}\big{)}^{S_{1}+1}\Gamma\Big{(}-S_{1},-S_{2}\big{(}\textstyle\frac{m(t)}{m_{0}}\big{)},-S_{2}\Big{)}\bigg{]}\bm{\hat{\mathbf{n}}}_{0}$
$\displaystyle\times\bigg{(}\frac{\big{(}\textstyle\frac{m(t)}{m_{0}}\big{)}^{S_{1}+1}}{S_{1}+1}-\frac{\big{(}\textstyle\frac{m(t)}{m_{0}}\big{)}^{S_{1}+1}}{S_{1}}\left(1-e^{-\gamma_{\infty}(\max(t,\Delta
t)-\Delta t)}\right)\bigg{)},$

with $S_{2}=D_{r}\Delta t/\zeta$ proportional to the rotational noise.
Orientational fluctuations do not contribute to the short-time behavior as
witnessed by the fact that in this limit the mean displacement coincides with
the noise-free acceleration behavior of Eq. (46). For long times, on the other
hand, the mean reach of the Langevin rocket is

$\displaystyle\overline{\Delta
R_{\infty}}=\frac{u}{D_{r}}\frac{1}{S_{1}+1}\left(1-e^{-D_{r}\Delta
t}\right)+\frac{u}{D_{r}}\frac{(1-\zeta)^{S_{1}+1}}{S_{1}(S_{1}+1)}$
$\displaystyle\times\operatorname{Re}\Big{[}e^{S_{2}}\big{(}-S_{2}\big{)}^{S_{1}+1}\Gamma\big{(}-S_{1},-S_{2}(1-\zeta),-S_{2}\big{)}\Big{]}.$
(53)

Returning to the previous optimization problem, we now maximize the mean reach
as a function of burn time $\Delta t$ and mass fraction $\zeta$ for fixed
prescribed noise strength $D_{r}/\gamma_{0}$. In Figure 7(a), the resulting
maximal reach $\max(\overline{\Delta R_{\infty}})$ is shown for varied noise
strength $D_{r}/\gamma_{0}$ in units of its universal noise-free limit
$u/\gamma_{0}$ of complete extended mass ejection. The associated optimal burn
time $\Delta t_{\text{max}}$ and optimal mass fraction $\zeta_{\text{max}}$
are also presented (see Figure 7(b) and 7(c)). If rotational noise is
increased, the complete extended mass ejection is still the best strategy, but
it is optimal to burn the full mass over a finite burn time. This strategy
defies best the ultimate orientational decorrelation which reduces the mean
reach. In the opposite limit of very large orientational noise, the best
strategy is to get momentum quickly by ejecting apart of the mass and use it
to proceed further within the characteristic damping time. If one would eject
the mass completely, the system is overdamped after the burn time and will
stop immediately lacking the additional benefit of the inertia. Hence the
fractional instantaneous mass ejection is the optimal strategy. Interestingly,
there is a sharp noise-induced discontinuous transition between the two
strategies for an intermediate finite value

$\displaystyle D_{r,\text{crit}}\approx 0.72\,\gamma_{0},$ (54)

of the orientational noise. The latter is signaled by a sharp jump in the
optimal burn time from $1.39\,\gamma_{0}$ to $0$ (see Figure 7(b)). The
optimal propellant mass fraction jumps from 1 to the universal value of
$1-e^{-1}\approx 0.63$ (see Figure 7(c)) and can thus be viewed as the ”order
parameter” of the transition.

###  IV.2  Comparison between the different set-ups

We now compare the different set-ups for time-dependent inertia as discussed
in Section II in more detail (see again Figure 2). In the case of directed
mass ejection or isotropic mass evaporation (Figure 2(b) and 2(c)), we assume
a mass loss exponentially in time $t$ as

$m(t)=m_{\infty}+(m_{0}-m_{\infty})e^{-\gamma_{m}t},$ (55)

where $m_{0}$ is the initial mass, $m_{\infty}$ the rest mass which remains
after the fuel is burned and $\gamma_{m}$ the mass decay rate. As outlined in
Appendix B, an exponential mass loss occurs in particular for a rocket which
ejects gas molecules at constant speed from a tank under isothermal and
isochoric conditions. In this case the exponential mass reduction follows from
the reduction of the gas density in the tank. Accordingly we also assume an
exponential decrease in the moment of inertia

$J(t)=J_{\infty}+(J_{0}-J_{\infty})e^{-\gamma_{J}t},$ (56)

where $J_{0}$ is the initial and $J_{\infty}$ the final moment of inertia and
$\gamma_{J}$ the decay rate of the moment of inertia. For the isotropic shape
change (Figure 2(d)) the mass is assumed to be constant and only an
exponential loss in the moment of inertia is prescribed.

| time-independent inertia | directed mass ejection | isotropic mass ejection | isotropic shape change
---|---|---|---|---
$m(t)$ | $m_{0}$ | $m_{\infty}+(m_{0}-m_{\infty})e^{-\gamma_{m}t}$ | $m_{\infty}+(m_{0}-m_{\infty})e^{-\gamma_{m}t}$ | $m_{0}$
$\gamma_{m}$ | - | 0.1 $D_{r}$ | 0.1 $D_{r}$ | -
$m_{\infty}/m_{0}$ | - | 0.1 | 0.1 | -
$u$ | - | 1 $v_{0}$ | - | -
$J(t)$ | $J_{0}$ | $J_{0}$ | $J_{\infty}+(J_{0}-J_{\infty})e^{-\gamma_{J}t}$ | $J_{\infty}+(J_{0}-J_{\infty})e^{-\gamma_{J}t}$
$\gamma_{J}$ | - | - | 0.1 $D_{r}$ | 0.1 $D_{r}$
$J_{\infty}/J_{0}$ | - | - | 0.1 | 0.1
Table 1:  Simulation parameter for the different set-ups.

The protocol is as follows. At time $t=0$, we start from a steady state
achieved for constant parameters and then initiate the mass loss and moment of
inertia change (or in general arbitrary time-dependencies). For the different
dynamical correlation functions we correlate the system configuration after a
time $t$ with the steady state condition at time $t=0$ (over which we perform
the average). For the different dynamical correlation functions we correlate
the steady state condition at time $t=0$ over which we perform the average
with the system configuration after a time $t$. Under these conditions, we
obtain general analytical results for arbitrary time-dependencies. Since the
system is relaxing or ”aging”, the two-point correlation functions now depend
explicitly on two times $t_{1}$, $t_{2}$, not just on the time difference as
in the steady state.

For $t_{1}<t_{2}$, the orientational correlation function
$C(t_{1},t_{2})=\langle\bm{\hat{\mathbf{n}}}(t_{1})\cdot\bm{\hat{\mathbf{n}}}(t_{2})\rangle$
is given by

$C(t_{1},t_{2})=\cos\big{(}\mu(t_{1},t_{2})\big{)}\,e^{-\textstyle\frac{1}{2}\sigma(t_{1},t_{2})},$
(57)

with the mean angle difference

$\mu(t_{1},t_{2})=\int_{t_{1}}^{t_{2}}d{t^{\prime\prime}}\int_{-\infty}^{t^{\prime\prime}}d{t^{\prime}}\,\frac{\xi_{r}(t^{\prime})}{J(t^{\prime})}\,\omega(t^{\prime})e^{-\Gamma_{r}(t^{\prime},t^{\prime\prime})},$
(58)

the corresponding variance

$\displaystyle\sigma(t_{1},t_{2})=4\int_{t_{1}}^{t_{2}}d{t^{\prime\prime\prime}}\int_{t_{1}}^{t^{\prime\prime\prime}}d{t^{\prime\prime}}$
(59)
$\displaystyle\quad\times\bigg{(}\int_{-\infty}^{t^{\prime\prime}}d{t^{\prime}}\bigg{(}\frac{\xi_{r}(t^{\prime})}{J(t^{\prime})}\bigg{)}^{2}D_{r}(t^{\prime})e^{-2\Gamma_{r}(t^{\prime},t^{\prime\prime})}\bigg{)}e^{-\Gamma_{r}(t^{\prime\prime},t^{\prime\prime\prime})}.$

and the rotational damping function

$\Gamma_{r}(t_{1},t_{2})=\int_{t_{1}}^{t_{2}}dt^{\prime}\,\frac{\xi_{r}(t^{\prime})}{J(t^{\prime})}+(1-\nu)\ln\left(\frac{J(t_{2})}{J(t_{1})}\right).$
(60)

Here, $\nu=0$ in the case of isotropic shape change and $\nu=1$ in the case of
isotropic mass evaporation.

Similarly, the velocity correlation function
$Z(t_{1},t_{2})=\langle\dot{\bm{\mathrm{R}}}(t_{1})\cdot\dot{\bm{\mathrm{R}}}(t_{2})\rangle$
for $t_{1}<t_{2}$ is

$\displaystyle
Z(t_{1},t_{2})=4\int_{-\infty}^{t_{1}}d{t^{\prime}}\bigg{(}\frac{\xi(t^{\prime})}{m(t^{\prime})}\bigg{)}^{2}D(t^{\prime})e^{-\Gamma(t^{\prime},t_{1})}e^{-\Gamma(t^{\prime},t_{2})}$
(61)
$\displaystyle+\int_{-\infty}^{t_{1}}d{t^{\prime}}\int_{-\infty}^{t_{2}}d{t^{\prime\prime}}a(t^{\prime})a(t^{\prime\prime})\langle\bm{\hat{\mathbf{n}}}(t^{\prime})\cdot\bm{\hat{\mathbf{n}}}(t^{\prime\prime})\rangle
e^{-\Gamma(t^{\prime},t_{1})}e^{-\Gamma(t^{\prime\prime},t_{2})},$

with the acceleration

$a(t)=\frac{\xi(t)}{m(t)}v_{0}(t)-\frac{\dot{m}(t)}{m(t)}u(t),$ (62)

and the translational damping function

$\Gamma(t_{1},t_{2})=\int_{t_{1}}^{t_{2}}dt^{\prime}\,\frac{\xi(t^{\prime})}{m(t^{\prime})}.$
(63)

For the delay function
$d(t_{1},t_{2})=\langle\dot{\bm{\mathrm{R}}}(t_{2})\cdot\bm{\hat{\mathbf{n}}}(t_{1})\rangle-\langle\dot{\bm{\mathrm{R}}}(t_{1})\cdot\bm{\hat{\mathbf{n}}}(t_{2})\rangle$,
we obtain

$\displaystyle d(t_{1},t_{2})=$
$\displaystyle\int_{-\infty}^{t_{2}}d{t^{\prime}}a(t^{\prime})\,\langle\bm{\hat{\mathbf{n}}}(t^{\prime})\cdot\bm{\hat{\mathbf{n}}}(t_{1})\rangle\,e^{-\Gamma(t^{\prime},t_{2})}$
(64)
$\displaystyle-\int_{-\infty}^{t_{1}}d{t^{\prime}}a(t^{\prime})\,\langle\bm{\hat{\mathbf{n}}}(t^{\prime})\cdot\bm{\hat{\mathbf{n}}}(t_{2})\rangle\,e^{-\Gamma(t^{\prime},t_{1})}.$

The general expression for the mean displacement
$\langle\Delta\bm{\mathrm{R}}(t_{1},t_{2})\rangle=\langle\bm{\mathrm{R}}(t_{2})-\bm{\mathrm{R}}(t_{1})\rangle$
is

$\langle\Delta\bm{\mathrm{R}}(t_{1},t_{2})\rangle=\int_{t_{1}}^{t_{2}}d{t^{\prime}}\int_{-\infty}^{t^{\prime}}d{t^{\prime\prime}}a(t^{\prime\prime})\,\langle\bm{\hat{\mathbf{n}}}(t^{\prime\prime})|\bm{\hat{\mathbf{n}}}(t_{1})\rangle
e^{-\Gamma(t^{\prime\prime},t^{\prime})},$ (65)

where the conditional average

$\displaystyle\langle\bm{\hat{\mathbf{n}}}(t_{2})|\bm{\hat{\mathbf{n}}}(t_{1})\rangle$
(66)
$\displaystyle=\begin{cases}\bm{\hat{\mathbf{P}}}\Big{[}e^{-\textstyle\frac{1}{2}\sigma(t_{2},t_{1})+i(\phi_{1}+\mu(t_{2},t_{1}))}\Big{]},&\text{for}\,t_{2}>t_{1},\\\
\bm{\hat{\mathbf{P}}}\Big{[}e^{-\textstyle\frac{1}{2}\sigma(t_{1},t_{2})+i(\phi_{1}+\mu(t_{1},t_{2}))}\Big{]},&\text{for}\,t_{2}<t_{1},\end{cases}$

denotes the mean orientation under the condition that the particle has the
angle $\phi(t_{1})=\phi_{1}$ at time $t_{1}$.

Last, the mean-square displacement
$\langle\Delta\bm{\mathrm{R}}^{2}(t_{1},t_{2})\rangle=\langle(\bm{\mathrm{R}}(t_{2})-\bm{\mathrm{R}}(t_{1}))^{2}\rangle$
is

$\langle\Delta\bm{\mathrm{R}}^{2}(t_{1},t_{2})\rangle=\int_{t_{1}}^{t_{2}}d{t^{\prime}}\int_{t_{1}}^{t_{2}}d{t^{\prime\prime}}Z(t^{\prime},t^{\prime\prime}).$
(67)

For time-independent parameters, we recover the results discussed in section
III. In particular, we have $C(t_{1},t_{2})=C(\lvert t_{1}-t_{2}\rvert)$ (see
Eq. (21)), $Z(t_{1},t_{2})=Z(\lvert t_{1}-t_{2}\rvert)$ (see Eq. (18)),
$d(t_{1},t_{2})=d(\lvert t_{1}-t_{2}\rvert)$ (see Eq. (37)),
$\langle\Delta\bm{\mathrm{R}}(t_{1},t_{2})\rangle=\langle\Delta\bm{\mathrm{R}}(\lvert
t_{1}-t_{2}\rvert)\rangle$ (see Eq. (24)) and
$\langle\Delta\bm{\mathrm{R}}^{2}(t_{1},t_{2})\rangle=\langle\Delta\bm{\mathrm{R}}^{2}(\lvert
t_{1}-t_{2}\rvert)\rangle$ (see Eq. (27)).

Figure 8:  Comparison of the different special set-ups for an achiral active
particle ($\omega=0$) with time-dependent inertia: (a) orientation
autocorrelation function $C(0,t)$, (b) velocity autocorrelation function
$Z(0,t)$, (c) delay function $d(0,t)$, (d) mean displacement along the initial
orientation
$\langle\Delta\bm{\mathrm{R}}(0,t)\rangle\cdot\bm{\hat{\mathbf{n}}}_{0}$, (e)
mean-square displacement $\langle\Delta\bm{\mathrm{R}}^{2}(0,t)\rangle$ and
(f) the corresponding scaling exponent $\alpha(0,t)$ for time-independent
inertia (dashed black), directed mass ejection (red), isotropic mass
evaporation (green) and an isotropic change in the particle shape (blue).
Velocities are given in units of $v_{0}$, times in $1/D_{r}$ and lengths in
$l_{p}=v_{0}/D_{r}$. The time-dependencies of the mass $m(t)$ and the moment
of inertia $J(t)$ for the different set-ups are summarized in Table 1. The
remaining parameters are $D=0$, $\gamma_{0}=\xi/m_{0}=0.1D_{r}$ and
$\gamma_{r,0}=\xi_{r}/J_{0}=0.1D_{r}$.  Figure 9:  Same as in Figure 8 for a
chiral particle with a spinning frequency of $\omega=0.1D_{r}$.

Numerical data for the special case of an exponential mass loss (see Eq. (55))
and/or an exponential decay of the moment of inertia (see Eq. (56)) as
summarized in Table 1 are presented in Figures 8 and 9. Figure 8 is for an
achiral particle and Figure 9 for a chiral particle. The case of time-
independent inertia (with the parameters at time $t=0$) is shown as a
reference, too. Equations (8) and (9) were discretized to perform Brownian
dynamics simulations. For these simulations, we chose the time step $\Delta
t=10^{-2}/D_{r}$ and we performed $10^{6}$ realizations to calculate the
respective ensemble averages.

We first discuss the case of an achiral particle. For isotropic shape change,
the orientational correlation function $C(0,t)$ decorrelates faster (see
Figure 8(a)), since the rotational noise is amplified during the decay of the
moment of inertia. The velocity autocorrelation $Z(0,t)$ as well as the delay
function $d(0,t)$ decorrelate faster if the particle actually looses mass (see
Figure 8(b) and 8(c)). For the particle with directed mass ejection, we see an
increase in the velocity autocorrelation for short times and a more pronounced
peak in the delay function due to the additional acceleration which enhances
the particle velocity. The mean displacement along the initial displacement
$\langle\Delta\bm{\mathrm{R}}(0,t)\rangle\cdot\bm{\hat{\mathbf{n}}}_{0}$ is
displayed in Figure 8(d). Although the particle with directed mass ejection is
the fastest for short times, it gets overtaken for long times by the particle
with time-independent inertia. Last we discuss the mean-square displacement
$\langle\Delta\bm{\mathrm{R}}^{2}(0,t)\rangle$. Beside additional acceleration
for the particle with directed mass ejection for short times, the long time
diffusivity is identical to the case of time-independent inertia. In contrast,
the cases of isotropic mass evaporation and isotropic shape end up with a
decreased long-time diffusion coefficient (see Figure 8(e)) due to a smaller
persistence. The differences between the set-ups become clearer by considering
the logarithmic derivative of the mean-square displacement

$\alpha(t_{1},t_{2})=\frac{d\ln\big{(}\langle\Delta\bm{\mathrm{R}}^{2}(t_{1},t_{2})\rangle\big{)}}{d\ln(t_{2})},$
(68)

If the mean-square displacement follows a power law
$\langle\Delta\bm{\mathrm{R}}^{2}(t_{1},t_{2})\rangle\sim(t_{2}-t_{1})^{\alpha}$,
$\alpha(t_{1}-t_{2})$ is equal to the power-law exponent $\alpha$. This
scaling exponent is shown in Figure 8(f). All set-ups start in a ballistic
regime ($\alpha=2$) for short times and end up in a diffusive regime
($\alpha=1$) for long times. Again for the particle with directed mass
ejection we observe faster motion for short times indicated by a super-
ballistic scaling $\alpha>2$ due to the acceleration. For times greater than
the inverse decay rate of the moment of inertia $1/\gamma_{J}$, the particle
with isotropic mass evaporation and an isotropic shape change behave sub-
diffusively with $\alpha<1$ since their effective diffusivity decreases.

Now we turn to the case of a chiral particle. First of all, even for constant
parameters, the presence of the torque $M$ yields systematic oscillations in
the orientation and velocity autocorrelations, and also in the delay function
(see Figure 9(a)-9(c)). Indeed such oscillations in the delay function have
been found recently in data for whirligig beetles Turner . Turning to the
time-dependent cases, similar to the pirouette of figure skating, the particle
with an isotropic shape contraction is spinning with a higher frequency during
the decay of the moment of inertia. This is visible in the orientational and
velocity autocorrelation functions and the delay function (see Figure
9(a)-9(c)). Also, when the particle loses mass, the oscillation becomes more
pronounced since the particle can adapt more easily to orientation changes. In
contrast to the achiral case, the long-time behavior of the mean-square
displacement increases for the time-dependent set-ups when the moment of
inertia $J(t)$ decreases (see Figure 9(e)) in line with the trend discussed
previously in Figure 5. This is marked by a peak in the scaling exponent for
times larger than $1/\gamma_{J}$ (see Figure 9(f)).

###  IV.3 Adiabatic approximation for slow variations

When the parameters (such as mass $m(t)$, moment of inertia $J(t)$, friction
coefficients $\xi(t)$ and $\xi_{r}(t)$, noise strengths $D(t)$ and $D_{r}(t)$,
self-propulsion velocity $v_{0}(t)$) change very slowly in time, i.e. much
slower than any other time-scale inherent in the model, the system can be
analyzed using the adiabatic approximation. In other words, one can take the
expressions for the dynamical correlation function with constant parameters
(as discussed in section III) and insert into these expressions the slowly
varying time-dependent parameters. This approximation becomes exact if the two
time scales (largest system time scale and fastest time scale governing the
change of all parameters) are separated completely.

Let us elaborate on the adiabatic approximation for the mean-square
displacement (MSD) by considering an achiral active particle. Corresponding
analytical expressions for the MSD in the two limits of small and high moments
of inertia $J$ are given by (35) and (36) respectively. Using the long-time
limit (31) we obtain within the adiabatic approximation for large $J$

$\langle\Delta\bm{\mathrm{R}}^{2}(t)\rangle\sim
4\Bigg{(}D(t)+\frac{v_{0}^{2}(t)}{4}\sqrt{\frac{2\pi
J(t)}{D_{r}(t)\,\xi_{r}(t)}}\Bigg{)}t,$ (69)

when the moment of inertia $J$ becomes sufficiently large and

$\langle\Delta\bm{\mathrm{R}}^{2}(t)\rangle\sim
4\Bigg{(}D(t)+\frac{v_{0}^{2}(t)}{2D_{r}(t)}+\frac{v_{0}^{2}(t)J(t)}{2\xi_{r}(t)}\Bigg{)}t,$
(70)

in the case of a small moment of inertia $J$. Let us now assume a slow power
law in time for the moment of inertia, the self-propulsion, the rotational
friction, and the diffusion coefficients

$v_{0}(t)\sim t^{\beta},\quad J(t)\sim t^{\delta},\quad\xi_{r}(t)\sim
t^{\varepsilon},\quad D_{r}(t)\sim t^{\eta},$ (71)

with prescribed dynamical exponents $\beta$, $\delta$, $\epsilon$ and $\eta$.
Plugging this into the expressions (69) and (70) we obtain a power law for the
long-time MSD of the active particle

$\langle\Delta\bm{\mathrm{R}}^{2}(t)\rangle\sim t^{\alpha},$ (72)

with

$\alpha=\max(1,1+2\beta-\frac{1}{2}\left(\varepsilon-\delta+\eta\right)),$
(73)

for large $J$ and

$\alpha=\max(1,1+2\beta-\min(\varepsilon-\delta,\eta)),$ (74)

for small $J$. If $\alpha>1$, adiabatic term is dominating overwhelming the
standard diffusion such that the particle exhibits anomalous super-diffusion.
If $\alpha=1$, the full MSD is dominated by the translational diffusion. We
finally remark that simpler scaling laws were obtained earlier in the
overdamped limit Babel  _et al._ (2014).

##  V  Conclusions

To conclude, we have investigated the dynamics of an inertia-dominated
Brownian particle, referred to as active Langevin dynamics. Dynamical
correlations within a simple model were calculated for a single ”microflyer”
which is simultaneously subjected to self-propulsion, inertia, damping and
fluctuations and analytical results known for the overdamped limit of
microswimmers were generalized to the inertial situation. In particular, we
considered the case of time-dependent inertia. Furthermore we identified a
basic Langevin model for a rocket-like particle self-propelled by the ejection
of mass for which we calculated its mean reach and found a noise-induced
discontinuous transition in the optimal propulsion strategy for reaching the
furthest distance. The case of chiral particles referred to as circle-flyers
was included. One characteristic dynamical correlation absent in the
overdamped case concerns the inertial delay between the orientation variations
and the subsequent changes in the velocity direction. For achiral particles
with vanishing spinning frequency, the inertial delay decays to zero after a
profound peak at a typical delay time. Conversely, for chiral particles, the
inertial delay correlation may oscillate between positive and negative values.
Finally we have also addressed the limiting ”adiabatic” case of very slow
inertia variation and have highlighted that a microflyer can undergo anomalous
diffusion if the parameters are varying as a power law in time.

Future work should generalize the present model to external potentials like
optical fields, disorder and confinement Bianchi  _et al._ (2016); Caprini
_et al._ (2019); Dauchot and Démery (2019); Militaru  _et al._ (2020); Breoni
_et al._ (2020) and to motion in noninertial rotating frames Löwen (2019);
Zheng and Löwen (2020). Furthermore, anisotropic particles which show out-of-
plane orientations and positions relevant for active complex plasmas Landau
and Lifshitz (1976) should be considered in the future. In this case, the
equations of motion are getting more complex involving friction and inertia
tensors significantly more complicated than in the overdamped limit Wittkowski
and Löwen (2012); Kraft  _et al._ (2013). Next, the ”rocket-like” particles
studied here should be realized in experiments, the most promising way seems
to be dust particles in the plasma with evaporating mass. Moreover it would be
interesting to study collective effects of inertia-dominated active particles
such as motility-induced-phase-separation Suma  _et al._ (2014); Scholz  _et
al._ (2018b); Petrelli  _et al._ (2018); Mayya  _et al._ (2019); Mandal  _et
al._ (2019); Caprini and Marconi (2020) or pattern formation in general Arold
and Schmiedeberg (2020). Finally it would be interesting to generalize the
more coarse-grained Ornstein-Uhlenbeck model for inertial active particles
Caprini and Marini Bettolo Marconi (2021); Nguyen  _et al._ to the situation
of time-dependent parameters.

##  VI Acknowledgments

We thank Christian Scholz, Ian Williams and Anton Ldov for helpful
discussions. This work is supported by the German Research Foundation through
grants LO 418/23-1 and IV 20/3-1.

##  Appendix A General solution

For an analytical solution of the equations of motion we first consider the
rotational part (see Eq. (9)) For $\phi_{0}=\phi(t=0)$ and
$\dot{\phi}_{0}=\dot{\phi}(t=0)$, the solution of Eq. (9) is

$\displaystyle\dot{\phi}(t)$
$\displaystyle=\dot{\phi}_{0}\,e^{-\Gamma_{r}(0,t)}$ (75)
$\displaystyle+\int_{0}^{t}d{t^{\prime}}\frac{\xi_{r}(t^{\prime})}{J(t^{\prime})}\,\omega(t^{\prime})\,e^{-\Gamma_{r}(t^{\prime},t)}$
$\displaystyle+\int_{0}^{t}d{t^{\prime}}\frac{\xi_{r}(t^{\prime})}{J(t^{\prime})}\sqrt{2D_{r}(t^{\prime})}\tau_{\rm
st}(t^{\prime})\,e^{-\Gamma_{r}(t^{\prime},t)},$

and thus

$\displaystyle\phi(t)$
$\displaystyle=\phi_{0}+\int_{0}^{t}d{t^{\prime}}\dot{\phi}_{0}\,e^{-\Gamma_{r}(0,t^{\prime})}$
(76)
$\displaystyle+\int_{0}^{t}d{t^{\prime}}\int_{0}^{t^{\prime}}d{t^{\prime\prime}}\frac{\xi_{r}(t^{\prime\prime})}{J(t^{\prime\prime})}\,\omega(t^{\prime\prime})e^{-\Gamma_{r}(t^{\prime\prime},t^{\prime})}$
$\displaystyle+\int_{0}^{t}d{t^{\prime}}\int_{0}^{t^{\prime}}d{t^{\prime\prime}}\frac{\xi_{r}(t^{\prime\prime})}{J(t^{\prime\prime})}\sqrt{2D_{r}(t^{\prime\prime})}\tau_{\rm
st}(t^{\prime\prime})e^{-\Gamma_{r}(t^{\prime\prime},t^{\prime})},$

where

$\Gamma_{r}(t_{1},t_{2})=\int_{t_{1}}^{t_{2}}dt^{\prime}\,\frac{\xi_{r}(t^{\prime})}{J(t^{\prime})}+(1-\nu)\ln\left(\frac{J(t_{2})}{J(t_{1})}\right).$
(77)

The translational equation of motion yields for the particle velocity

$\displaystyle\dot{\bm{\mathrm{R}}}(t)$
$\displaystyle=\dot{\bm{\mathrm{R}}}_{0}\,e^{-\Gamma(0,t)}$ (78)
$\displaystyle+\int_{0}^{t}d{t^{\prime}}a(t^{\prime})\,\bm{\hat{\mathbf{n}}}(t^{\prime})\,e^{-\Gamma(t^{\prime},t)}$
$\displaystyle+\int_{0}^{t}d{t^{\prime}}\frac{\xi(t^{\prime})}{m(t^{\prime})}\sqrt{2D(t^{\prime})}\bm{\mathrm{f}}_{\rm
st}(t^{\prime})\,e^{-\Gamma(t^{\prime},t)}.$

Hence, the center-of-mass position is calculated as

$\displaystyle\bm{\mathrm{R}}(t)$
$\displaystyle=\bm{\mathrm{R}}_{0}+\int_{0}^{t}d{t^{\prime}}\dot{\bm{\mathrm{R}}}_{0}\,e^{-\Gamma(0,t^{\prime})}$
(79)
$\displaystyle+\int_{0}^{t}d{t^{\prime}}\int_{0}^{t^{\prime}}d{t^{\prime\prime}}a(t^{\prime\prime})\,\bm{\hat{\mathbf{n}}}(t^{\prime\prime})e^{-\Gamma(t^{\prime\prime},t^{\prime})}$
$\displaystyle+\int_{0}^{t}d{t^{\prime}}\int_{0}^{t^{\prime}}d{t^{\prime\prime}}\frac{\xi(t^{\prime\prime})}{m(t^{\prime\prime})}\sqrt{2D(t^{\prime\prime})}\bm{\mathrm{f}}_{\rm
st}(t^{\prime\prime})e^{-\Gamma(t^{\prime\prime},t^{\prime})},$

where

$a(t)=\frac{\xi(t)}{m(t)}v_{0}(t)-\frac{\dot{m}(t)}{m(t)}u(t).$ (80)

and

$\Gamma(t_{1},t_{2})=\int_{t_{1}}^{t_{2}}dt^{\prime}\,\frac{\xi(t^{\prime})}{m(t^{\prime})}.$
(81)

Here $\bm{\mathrm{R}}_{0}$ and $\dot{\bm{\mathrm{R}}}_{0}$ are the initial
position and velocity at time $t=0$.

##  Appendix B Exponential mass loss

In an isothermal environment of temperature $T$, the mass loss through a small
leak of cross section $S$ in the rocket tank of volume $V$ in quasi-
equilibrium is governed by

$\dot{m}_{\text{fuel}}(t)=-\frac{1}{6}\frac{S}{V}\sqrt{\frac{3k_{B}T}{m_{\text{mol}}}}\,m_{\text{fuel}}(t)=-\gamma_{m}\,m_{\text{fuel}}(t),$
(82)

where $m_{\text{mol}}$ is the mass of the ejected molecules and $k_{B}$ is the
Boltzmann constant. Equation (82) implies an exponential decay of the rocket
fuel, i.e., $m_{\text{fuel}}(t)=m_{\text{fuel}}(0)\,e^{-\gamma_{m}t}$ with the
mass decay rate $\gamma_{m}$ and thus motivates Eq. (55).

## References

  * Howse  _et al._ (2007) J. R. Howse, R. A. L. Jones, A. J. Ryan, T. Gough, R. Vafabakhsh, and R. Golestanian, Phys. Rev. Lett. 99, 048102 (2007).
  * Romanczuk  _et al._ (2012) P. Romanczuk, M. Bär, W. Ebeling, B. Lindner, and L. Schimansky-Geier, Eur. Phys. J. Spec. Top. 202, 1 (2012).
  * Elgeti  _et al._ (2015) J. Elgeti, R. G. Winkler, and G. Gompper, Rep. Prog. Phys. 78, 056601 (2015).
  * Bechinger  _et al._ (2016) C. Bechinger, R. Di Leonardo, H. Löwen, C. Reichhardt, G. Volpe, and G. Volpe, Reviews of Modern Physics 88, 045006 (2016).
  * Purcell (1977) E. M. Purcell, Am. J. Phys. 45, 3 (1977).
  * ten Hagen  _et al._ (2011) B. ten Hagen, S. van Teeffelen, and H. Löwen, J. Phys. Condens. Matter 23, 194119 (2011).
  * Klotsa (2019) D. Klotsa, Soft Matter 15, 8946 (2019).
  * Scholz  _et al._ (2018a) C. Scholz, S. Jahanshahi, A. Ldov, and H. Löwen, Nat. Commun. 9, 5156 (2018a).
  * Callegari and Volpe (2019) A. Callegari and G. Volpe,  _Flowing Matter_ (Springer International Publishing, 2019) pp. 211–238.
  * Um  _et al._ (2019) J. Um, T. Song, and J.-H. Jeon, Front. Phys. 7, 143 (2019).
  * Löwen (2020) H. Löwen, The Journal of Chemical Physics 152, 040901 (2020).
  * Enculescu and Stark (2011) M. Enculescu and H. Stark, Phys. Rev. Lett. 107, 058301 (2011).
  * Ghosh  _et al._ (2015) P. K. Ghosh, Y. Li, G. Marchegiani, and F. Marchesoni, J. Chem. Phys. 143, 211101 (2015).
  * Joyeux and Bertin (2016) M. Joyeux and E. Bertin, Phys. Rev. E 93, 032605 (2016).
  * Manacorda and Puglisi (2017) A. Manacorda and A. Puglisi, Phys. Rev. Lett. 119, 208003 (2017).
  * Takatori and Brady (2017) S. C. Takatori and J. F. Brady, Phys. Rev. Fluids 2, 094305 (2017).
  * Mokhtari  _et al._ (2017) Z. Mokhtari, T. Aspelmeier, and A. Zippelius, EPL 120, 14001 (2017).
  * Das  _et al._ (2019) S. Das, G. Gompper, and R. G. Winkler, Sci. Rep. 9, 6608 (2019).
  * Sandoval (2020) M. Sandoval, Phys. Rev. E 101, 012606 (2020).
  * Gutierrez-Martinez and Sandoval (2020) L. L. Gutierrez-Martinez and M. Sandoval, The Journal of Chemical Physics 153, 044906 (2020).
  * Ebeling  _et al._ (1999) W. Ebeling, F. Schweitzer, and B. Tilch, Bio Systems 49 1, 17 (1999).
  * Weber  _et al._ (2013) C. A. Weber, T. Hanke, J. Deseigne, S. Léonard, O. Dauchot, E. Frey, and H. Chaté, Phys. Rev. Lett. 110, 208001 (2013).
  * Parisi  _et al._ (2018) D. R. Parisi, R. Cruz Hidalgo, and I. Zuriguel, Sci. Rep. 8, 9133 (2018).
  * Morfill and Ivlev (2009) G. E. Morfill and A. V. Ivlev, Rev. Mod. Phys. 81, 1353 (2009).
  * Sütterlin  _et al._ (2009) K. R. Sütterlin, A. Wysocki, A. V. Ivlev, C. Räth, H. M. Thomas, M. Rubin-Zuzic, W. J. Goedheer, V. E. Fortov, A. M. Lipaev, V. I. Molotkov, O. F. Petrov, G. E. Morfill, and H. Löwen, Phys. Rev. Lett. 102, 085003 (2009).
  * Couëdel  _et al._ (2010) L. Couëdel, V. Nosenko, A. V. Ivlev, S. K. Zhdanov, H. M. Thomas, and G. E. Morfill, Phys. Rev. Lett. 104, 195001 (2010).
  * Chaudhuri  _et al._ (2011) M. Chaudhuri, A. V. Ivlev, S. A. Khrapak, H. M. Thomas, and G. E. Morfill, Soft Matter 7, 1287 (2011).
  * Ivlev  _et al._ (2015) A. V. Ivlev, J. Bartnick, M. Heinen, C.-R. Du, V. Nosenko, and H. Löwen, Phys. Rev. X 5, 011035 (2015).
  * Nosenko  _et al._ (2020) V. Nosenko, F. Luoni, A. Kaouk, M. Rubin-Zuzic, and H. Thomas, Phys. Rev. Research 2, 033226 (2020).
  * Narayan  _et al._ (2007) V. Narayan, S. Ramaswamy, and N. Menon, Science 317, 105 (2007).
  * Kudrolli  _et al._ (2008) A. Kudrolli, G. Lumay, D. Volfson, and L. S. Tsimring, Phys. Rev. Lett. 100, 058001 (2008).
  * Deseigne  _et al._ (2010) J. Deseigne, O. Dauchot, and H. Chaté, Phys. Rev. Lett. 105, 098001 (2010).
  * Giomi  _et al._ (2013) L. Giomi, N. Hawley-Weld, and L. Mahadevan, Proc. Royal Soc. A 469, 20120637 (2013).
  * Klotsa  _et al._ (2015) D. Klotsa, K. A. Baldwin, R. J. A. Hill, R. M. Bowley, and M. R. Swift, Phys. Rev. Lett. 115, 248102 (2015).
  * Patterson  _et al._ (2017) G. A. Patterson, P. I. Fierens, F. Sangiuliano Jimka, P. G. König, A. Garcimartín, I. Zuriguel, L. A. Pugnaloni, and D. R. Parisi, Phys. Rev. Lett. 119, 248301 (2017).
  * Junot  _et al._ (2017) G. Junot, G. Briand, R. Ledesma-Alonso, and O. Dauchot, Phys. Rev. Lett. 119, 028002 (2017).
  * Ramaswamy (2017) S. Ramaswamy, J. Stat. Mech.: Theory Exp. 2017, 054002 (2017).
  * Deblais  _et al._ (2018) A. Deblais, T. Barois, T. Guerin, P. H. Delville, R. Vaudaine, J. S. Lintuvuori, J. F. Boudet, J. C. Baret, and H. Kellay, Phys. Rev. Lett. 120, 188002 (2018).
  * Dauchot and Démery (2019) O. Dauchot and V. Démery, Phys. Rev. Lett. 122, 068002 (2019).
  * Rabault  _et al._ (2019) J. Rabault, R. A. Fauli, and A. Carlson, Phys. Rev. Lett. 122, 024501 (2019).
  * Fauli  _et al._ (2019) R. A. Fauli, J. Rabault, and A. Carlson, Phys. Rev. E 100, 013108 (2019).
  * Leoni  _et al._ (2020) M. Leoni, M. Paoluzzi, S. Eldeen, A. Estrada, L. Nguyen, M. Alexandrescu, K. Sherb, and W. W. Ahmed, “Surfing and crawling macroscopic active particles under hard confinement – inertial dynamics,” (2020), arXiv:2005.10768.
  * Di Leonardo  _et al._ (2007) R. Di Leonardo, G. Ruocco, J. Leach, M. J. Padgett, A. J. Wright, J. M. Girkin, D. R. Burnham, and D. McGloin, Phys. Rev. Lett. 99, 010601 (2007).
  * Rubenstein  _et al._ (2014) M. Rubenstein, A. Cornejo, and R. Nagpal, Science 345, 795 (2014).
  * Fujiwara  _et al._ (2014) R. Fujiwara, T. Kano, and A. Ishiguro, Advanced Robotics 28, 639 (2014).
  * Tolba  _et al._ (2015) S. Tolba, R. Ammar, and S. Rajasekaran, IEEE Symposium on Computers and Communication (ISCC) , 1007 (2015).
  * Zhakypov  _et al._ (2019) Z. Zhakypov, K. Mori, K. Hosoda, and J. Paik, Nature 571, 381 (2019).
  * Yang  _et al._ (2019) X. Yang, C. Ren, K. Cheng, and H. P. Zhang, “Robust boundary flow in chiral active fluid,” (2019), arXiv:1910.14341.
  * (49) M. Turner, Private discussion.
  * Toner and Tu (1995) J. Toner and Y. Tu, Phys. Rev. Lett. 75, 4326 (1995).
  * Toner and Tu (1998) J. Toner and Y. Tu, Phys. Rev. E 58, 4828 (1998).
  * Chiappini (2008) E. Chiappini,  _Encyclopedia of Entomology_ (Springer Netherlands, 2008) pp. 152–154.
  * Bartussek and Lehmann (2016) J. Bartussek and F. O. Lehmann, Royal Soc. Open Sci. 3, 150562 (2016).
  * Mukundarajan  _et al._ (2016) H. Mukundarajan, T. C. Bardon, D. H. Kim, and M. Prakash, J. Exp. Biol. 219, 752 (2016).
  * Bartussek and Lehmann (2018) J. Bartussek and F. O. Lehmann, J. Royal Soc. Interface 15, 20180408 (2018).
  * Attanasi  _et al._ (2014) A. Attanasi, A. Cavagna, L. Del Castello, I. Giardina, T. S. Grigera, A. Jelic, S. Melillo, L. Parisi, O. Pohl, E. Shen, and M. Viale, Nat. Phys. 10, 692 (2014).
  * Vandenberg and Rayner (1995) C. Vandenberg and J. Rayner, J. Exp. Biol. 198, 1655 (1995).
  * Qiao  _et al._ (2014) M. Qiao, B. Brown, and D. L. Jindrich, J. Exp. Biol. 217, 432 (2014).
  * (59) S. I. Krasheninnikov, A. Y. Pigarov, R. D. Smirnov, and T. K. Soboleva, Contributions to Plasma Physics 50, 410.
  * Nosenko  _et al._ (2010) V. Nosenko, A. V. Ivlev, and G. E. Morfill, Physics of Plasmas 17, 123705 (2010).
  * Müller and Strehlow (1999) I. Müller and P. Strehlow, Physikalische Blätter 55, 37 (1999).
  * Mangan and Destrade (2015) R. Mangan and M. Destrade, International Journal of Non-Linear Mechanics 68, 52 (2015).
  * Eke and Mao (2002) F. O. Eke and T. C. Mao, IJMEE 30, 123 (2002).
  * Panchal and Jangid (2008) V. Panchal and R. Jangid, Nucl. Eng. Des. 238, 1304 (2008).
  * Lozano  _et al._ (2019) C. Lozano, B. Liebchen, B. ten Hagen, C. Bechinger, and H. Löwen, Soft Matter 15, 5185 (2019).
  * Fernandez-Rodriguez  _et al._ (2020) M. A. Fernandez-Rodriguez, F. Grillo, L. Alvarez, M. Rathlef, I. Buttinoni, G. Volpe, and L. Isa, Nature Communications 11, 4223 (2020).
  * Sprenger  _et al._ (2020) A. R. Sprenger, M. A. Fernandez-Rodriguez, L. Alvarez, L. Isa, R. Wittkowski, and H. Löwen, Langmuir 36, 7066 (2020).
  * Embs  _et al._ (2000) J. Embs, H. W. Müller, C. Wagner, K. Knorr, and M. Lücke, Phys. Rev. E 61, R2196 (2000).
  * Stickler  _et al._ (2018) B. A. Stickler, B. Schrinski, and K. Hornberger, Phys. Rev. Lett. 121, 040401 (2018).
  * Couder and Fort (2006) Y. Couder and E. Fort, Phys. Rev. Lett. 97, 154101 (2006).
  * Valani  _et al._ (2019) R. N. Valani, A. C. Slim, and T. Simula, Phys. Rev. Lett. 123, 024503 (2019).
  * Tsiolkovsky (1903) K. E. Tsiolkovsky,  _The Exploration of Cosmic Space by Means of Reaction Devices_ , Vol. 5 (The Science Review, 1903).
  * Metzler  _et al._ (1999) R. Metzler, E. Barkai, and J. Klafter, Phys. Rev. Lett. 82, 3563 (1999).
  * Golestanian (2009) R. Golestanian, Phys. Rev. Lett. 102, 188305 (2009).
  * Hoefling and Franosch (2013) F. Hoefling and T. Franosch, Rep. Prog. Phys. 76, 046602 (2013).
  * Bleibel  _et al._ (2014) J. Bleibel, A. Dominguez, F. Gunther, J. Harting, and M. Oettel, Soft Matter 10, 2945 (2014).
  * Morin  _et al._ (2017) A. Morin, D. L. Cardozo, V. Chikkadi, and D. Bartolo, Phys. Rev. E 96, 042611 (2017).
  * Taloni  _et al._ (2017) A. Taloni, O. Flomenbom, R. Castaneda-Priego, and F. Marchesoni, Soft Matter 13, 1096 (2017).
  * Charalambous  _et al._ (2019) C. Charalambous, M. Á. García-March, G. Muñoz-Gil, P. R. Grzybowski, and M. Lewenstein, “Control of anomalous diffusion of a Bose polaron,” (2019), arXiv:1910.01571.
  * Chaki and Chakrabarti (2019) S. Chaki and R. Chakrabarti, J. Chem. Phys. 150, 094902 (2019).
  * Bodrova  _et al._ (2016) A. S. Bodrova, A. V. Chechkin, A. G. Cherstvy, H. Safdari, I. M. Sokolov, and R. Metzler, “Underdamped scaled brownian motion: (non-)existence of the overdamped limit in anomalous diffusion,”  (2016).
  * van Teeffelen and Löwen (2008) S. van Teeffelen and H. Löwen, Phys. Rev. E 78, 020101 (2008).
  * Babel  _et al._ (2014) S. Babel, B. ten Hagen, and H. Löwen, J. Stat. Mech.: Theory Exp. 2014, P02011 (2014).
  * Olsen (2020) K. S. Olsen, “Perturbative effective diffusivity of microswimmers in the presence of oscillating torques,” (2020), arXiv:2012.05415.
  * ten Hagen  _et al._ (2015) B. ten Hagen, R. Wittkowski, D. Takagi, F. Kümmel, C. Bechinger, and H. Löwen, J. Phys. Condens. Matter 27, 194110 (2015).
  * Sommerfeld (1952) A. Sommerfeld,  _Mechanics_ , Lectures on Theoretical Physics, Vol. Vol. 1 (New York, 1952) p. p. 28.
  * Thomson (1966) W. T. Thomson, AIAA Journal 4, 766 (1966).
  * Plastino and Muzzio (1992) A. R. Plastino and J. C. Muzzio, Celestial Mechanics and Dynamical Astronomy 53, 227 (1992).
  * Anderson (1983) J. L. Anderson, Phys. Fluids 26, 2871 (1983).
  * Nardi  _et al._ (1999) J. Nardi, R. Bruinsma, and E. Sackmann, Phys. Rev. Lett. 82, 5168 (1999).
  * Eloul  _et al._ (2020) S. Eloul, W. C. K. Poon, O. Farago, and D. Frenkel, Phys. Rev. Lett. 124, 188001 (2020).
  * Li  _et al._ (2016) J. Li, I. Rozen, and J. Wang, ACS Nano 10, 5619 (2016).
  * Rankin and Rosenhead (1949) R. A. Rankin and L. Rosenhead, Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences 241, 457 (1949).
  * Paris (2010) R. B. Paris, “Incomplete gamma function,” in  _NIST Handbook of Mathematical Functions_ (Cambridge University Press, 2010).
  * Kurzthaler and Franosch (2017) C. Kurzthaler and T. Franosch, Soft Matter 13, 6396 (2017).
  * Hansen and McDonald (1990) J.-P. Hansen and I. R. McDonald,  _Theory of Simple Liquids_ (Elsevier, 1990).
  * Herrera and Sandoval (2021) P. Herrera and M. Sandoval, Phys. Rev. E 103, 012601 (2021).
  * Malakar  _et al._ (2020) K. Malakar, A. Das, A. Kundu, K. V. Kumar, and A. Dhar, Phys. Rev. E 101, 022610 (2020).
  * Kümmel  _et al._ (2013) F. Kümmel, B. ten Hagen, R. Wittkowski, I. Buttinoni, R. Eichhorn, G. Volpe, H. Löwen, and C. Bechinger, Physical Review Letters 110, 198302 (2013).
  * Askey and Daalhuis (2010) R. A. Askey and A. B. O. Daalhuis, “Generalized hypergeometric function,” in  _NIST Handbook of Mathematical Functions_ (Cambridge University Press, 2010).
  * Srivastava  _et al._ (2012) N. Srivastava, P. T. Tkacik, and R. G. Keanini, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 468, 3965 (2012).
  * Bianchi  _et al._ (2016) S. Bianchi, R. Pruner, G. Vizsnyiczai, C. Maggi, and R. di Leonardo, Scientific Reports 6, 27681 (2016).
  * Caprini  _et al._ (2019) L. Caprini, U. Marini Bettolo Marconi, A. Puglisi, and A. Vulpiani, Journal of Chemical Physics 150, 024902 (2019).
  * Militaru  _et al._ (2020) A. Militaru, M. Innerbichler, M. Frimmer, F. Tebbenjohanns, L. Novotny, and C. Dellago, “Escape dynamics of active particles in multistable potentials,” (2020), arXiv:2012.04478.
  * Breoni  _et al._ (2020) D. Breoni, M. Schmiedeberg, and H. Löwen, Phys. Rev. E 102, 062604 (2020).
  * Löwen (2019) H. Löwen, Phys. Rev. E 99, 062608 (2019).
  * Zheng and Löwen (2020) Y. Zheng and H. Löwen, Phys. Rev. Research 2, 023079 (2020).
  * Landau and Lifshitz (1976) L. D. Landau and E. M. Lifshitz,  _Mechanics_ , Course of Theoretical Physics, Vol. Vol. 1 (Pergamon Press, 1976).
  * Wittkowski and Löwen (2012) R. Wittkowski and H. Löwen, Phys. Rev. E 85, 021406 (2012).
  * Kraft  _et al._ (2013) D. J. Kraft, R. Wittkowski, B. ten Hagen, K. V. Edmond, D. J. Pine, and H. Löwen, Phys. Rev. E 88, 050301 (2013).
  * Suma  _et al._ (2014) A. Suma, G. Gonnella, D. Marenduzzo, and E. Orlandini, EPL (Europhysics Letters) 108, 56004 (2014).
  * Scholz  _et al._ (2018b) C. Scholz, M. Engel, and T. Pöschel, Nat. Commun. 9, 931 (2018b).
  * Petrelli  _et al._ (2018) I. Petrelli, P. Digregorio, L. F. Cugliandolo, G. Gonnella, and A. Suma, The European Physical Journal E 41, 128 (2018).
  * Mayya  _et al._ (2019) S. Mayya, G. Notomista, D. Shell, S. Hutchinson, and M. Egerstedt, “Surfing and crawling macroscopic active particles under hard confinement – inertial dynamics,” (2019), arXiv:1902.10662.
  * Mandal  _et al._ (2019) S. Mandal, B. Liebchen, and H. Löwen, Phys. Rev. Lett. 123, 228001 (2019).
  * Caprini and Marconi (2020) L. Caprini and U. M. B. Marconi, “Spatial velocity correlations in inertial systems of active brownian particles,” (2020), arXiv:2012.11916.
  * Arold and Schmiedeberg (2020) D. Arold and M. Schmiedeberg, Journal of Physics: Condensed Matter 32, 315403 (2020).
  * Caprini and Marini Bettolo Marconi (2021) L. Caprini and U. Marini Bettolo Marconi, The Journal of Chemical Physics 154, 024902 (2021).
  * (119) P. Nguyen, R. Wittmann, and H. Löwen, “Active Ornstein-Uhlenbeck model for self-propelled particles with inertia,” In preparation.