# Harmonically trapped inertial run-and-tumble particle in one dimension

Debraj Dutta  S. N. Bose National Centre for Basic Sciences, Kolkata 700106,
India  Anupam Kundu  International Centre for Theoretical Sciences, Bengaluru
560089, India  Sanjib Sabhapandit  Raman Research Institute, Bengaluru 560080,
India  Urna Basu  S. N. Bose National Centre for Basic Sciences, Kolkata
700106, India

###### Abstract

We study the nonequilibrium stationary state of a one-dimensional inertial
run-and-tumble particle (IRTP) trapped in a harmonic potential. We find that
the presence of inertia leads to two distinct dynamical scenarios, namely,
overdamped and underdamped, characterized by the relative strength of the
viscous and the trap time-scales. We also find that inertial nature of the
active dynamics leads to the particle being confined in specific regions of
the phase plane in the overdamped and underdamped cases, which we compute
analytically. Moreover, the interplay of the inertial and active time-scales
gives rise to several sub-regimes, which are characterized by very different
behaviour of position and velocity fluctuations of the IRTP. In particular, in
the underdamped regime, both the position and velocity undergoes transitions
from a novel multi-peaked structure in the strongly active limit to a single
peaked Gaussian-like distribution in the passive limit. On the other hand, in
the overdamped scenario, the position distribution shows a transition from a
U-shape to a dome-shape, as activity is decreased. Interestingly, the velocity
distribution in the overdamped scenario shows two transitions—from a single-
peaked shape with an algebraic divergence at the origin in the strongly active
regime to a double peaked one in the moderately active regime to a dome-shaped
one in the passive regime.

##  I Introduction

Active particles perform self-propelled motion by consuming energy at an
individual level and thereby breaking detailed balance. Active motion is
observed in nature over multiple length scales ranging from bacterial
motility, cell migration, and insect movements, to bird flocks and fish
schools at macroscopic scales Cavagna and Giardina (2014); Feinerman  _et al._
(2018); Pavlov and Kasumyan (2000); Mukundarajan  _et al._ (2016). Recent
years have also seen a huge advancement in the synthesis of artificial active
agents such as nanobots, Janus colloids to granular particles and self-
propelled robots Jiang and Granick (2012); Bechinger  _et al._ (2016); Walther
and Muller (2013); Nelson  _et al._ (2010). Most commonly, the motions of
individual active agents are described theoretically using overdamped Langevin
equations, such as Active Brownian process (ABP) Basu  _et al._ (2018); Howse
_et al._ (2007), Run-and-Tumble process (RTP) Malakar  _et al._ (2018);
Tailleur and Cates (2008) and Active Ornstein-Uhlenbeck process (AOUP) Bonilla
(2019); Martin  _et al._ (2021). In such a description, the inertial effects
are ignored assuming that the velocity relaxes over a time-scale much smaller
than the observation time-scale.

The typical relaxation time of the velocity of an object of linear size $\ell$
and density $\rho$ scales as $\rho\ell^{2}$. For example, for a micron-sized
particle, the typical relaxation time of the velocity is around $100$ ns,
while the typical temporal resolution of experiments is usually above
microseconds, making the overdamped Langevin description fairly accurate. On
the other hand, for a macro-sized particle, say, with $\ell\sim 1$ mm, the
relaxation time is about $0.1$ s, necessitating the inclusion of inertia in
the description of its motion. Therefore, to describe the motion of macro-
sized active particles, one needs to consider the underdamped Langevin
equations.

Recent years have seen an increased effort to explore the effect of inertia on
active particle motion, both at individual and collective levels. It was found
that the presence of inertia drastically changes the collective properties of
active matter. For example, it changes the nature of motility induced phase
transition Mandal  _et al._ (2019); Su  _et al._ (2021), hinders
crystallization De Karmakar and Ganesh (2020); Liao  _et al._ (2021), promotes
hexatic ordering in homogeneous phases Negro  _et al._ (2022) and, in general,
reduces spatial velocity correlations Marconi  _et al._ (2021); Caprini  _et
al._ (2021); Caprini and Marconi (2021). It has also been shown that the
presence of inertia makes an active reservoir behave more like an equilibrium
one Khali  _et al._ (2024).

The inclusion of inertia also leads to remarkable effects on the statistical
properties of individual active particles. It has been shown that adding
inertia to ABPs impacts their behaviour both in free space as well as in the
presence of external potential, manifest in the change of kinetic temperature
and swim pressure Patel and Chaudhuri (2023, ); Sprenger  _et al._ (2023).
Moreover, presence of rotational inertia also affects the long-time
diffusivity for an ABP Scholz  _et al._ (2018); Breoni  _et al._ (2020); Lisin
_et al._ (2022). The effect of inertia on the large deviation function of a
harmonically trapped RTP in the strongly passive regime has recently been
studied Smith and Farago (2022). The dynamics of inertial AOUP has also been
investigated recently both for free particles and harmonically trapped
particles Caprini and Marini Bettolo Marconi (2021), as well as in the
presence of underdamped self-propulsion force Sprenger  _et al._ (2023). The
behaviour of inertial ABP and AOUP have also been studied in the presence of a
magnetic field Muhsin and Sahoo (2022); Obreque  _et al._ , coulomb friction
Antonov  _et al._ and in viscoelastic media Muhsin and Sahoo (2023); Adersh
_et al._ . Since AOUP is a linear Gaussian process, the corresponding position
and velocity distributions are Gaussian, completely characterized by the
correlation matrix. For non-Gaussian processes like ABP and RTP, a
comprehensive understanding of the effect of inertia on the full position and
velocity distributions is still lacking.

In this paper, we address this question in the context of an inertial run-and-
tumble particle in the presence of a harmonic potential. Run-and-tumble
particle, which is one of the simplest and most-studied model of active
motion, describes the overdamped motion of a particle along an internal
direction which itself changes stochastically Malakar  _et al._ (2018);
Tailleur and Cates (2008). Inertia induces a longer delay for the particles to
change their speed, increasing the persistence of trajectories and particle
mobility. In the presence of a trap, the particle reaches a nonequilibrium
stationary state. We find that the interplay between inertia and self-
propulsion leads to non-trivial features in the stationary position and
velocity distributions, that we characterize analytically.

The paper is organized as follows. In the next section we introduce the the
inertial Run-and-Tumble particle (IRTP) model and provide a summary of our
main results. We then investigate the behaviour of the position and velocity
moments in Sec. III to get an idea about the typical position and velocity
fluctuations. In Sec. IV and Sec. V, we provide a detailed discussion of the
stationary position and velocity distributions in the context of different
regimes and sub-regimes arising from the interplay of different time-scales.
Finally, we conclude in Sec. VI with some open questions.

Figure 1: Phase space trajectory of a damped harmonic oscillator for the
overdamped (a) and underdamped (b) scenarios, starting from
$(x_{0},v_{0})=(1,0)$ [indicated by dark green squares]. The insets show the
corresponding time-evolution of $x(t)$ and $v(t)$ separately. In both the
cases, the particle eventually reaches $(0,0)$ [indicated by red discs]. Here
we have used $\mu=1.0,m=1.0,a_{0}=1.0$ with $\gamma=3.0$ for the overdamped
case and $\gamma=0.8$ for the underdamped case.

##  II Model and summary of the results

We consider the one-dimensional motion of an inertial run-and-tumble particle
(IRTP) in a harmonic potential $U(x)=\mu x^{2}/2$. The position $x(t)$ and the
velocity $v(t)$ of the particle evolve according to the Langevin equations,

$\displaystyle\dot{x}(t)$ $\displaystyle=v(t),$ (1) $\displaystyle
m\dot{v}(t)$ $\displaystyle=-\mu x(t)-\gamma v(t)+a_{0}\sigma(t),$ (2)

where $\gamma$ denotes the dissipation coefficient. The activity of the
particle originates from the self-propulsion force $a_{0}\sigma(t)$, where
$a_{0}$ is a constant and $\sigma(t)$ is a dichotomous noise, alternating
between $\pm 1$ with a constant rate $\tau_{a}^{-1}$. Note that, between two
such successive ‘tumbling’ events, i.e., when $\sigma$ remains constant, the
particle undergoes a deterministic damped motion in a harmonic trap centered
at $a_{0}\sigma/\mu$. The tumbling results in a switch in the position of the
trap between $\pm a_{0}/\mu$.

In general, there can also be an additional Gaussian white noise with a
strength $\sqrt{2\gamma k_{B}T}$ in Eq. (2), arising from thermal
fluctuations, which can be neglected in most physical situations. Moreover,
mathematically, the effect of the white noise can be incorporated
straightforwardly as a convolution with the $T=0$ solution. Hence, it suffices
to consider the $T=0$ case. We are interested in the position and velocity
fluctuations of this IRTP. In particular, we investigate how the presence of
the finite mass $m$ affects the marginal position and velocity distributions.

Figure 2: Typical phase trajectory for the IRTP for the overdamped (a) and
underdamped (b) scenarios. The orange and green solid lines correspond to
trajectories with $\sigma=\pm 1$ respectively and the dotted vertical lines
indicate the alternating positions $\pm a_{0}/\mu$ of the trap minima.
Additionally, the dashed vertical line marks the boundary of the support of
the joint distribution $P_{\text{st}}(x,v)$ in the phase plane. The
corresponding time evolution of $x(t)$ and $v(t)$ are also shown separately
for the overdamped (c) and underdamped (d) cases. Here we have used
$\tau_{a}=1,\mu=1,m=1$ with $\gamma=1$ for the overdamped case and
$\gamma=2.5$ for the underdamped case.

In the absence of the active noise, i.e., for $a_{0}=0$, equations Eqs.
(1)-(2) describe the deterministic motion of a particle in a harmonic trap
centered at the origin, where the particle eventually comes to rest at the
minimum of the potential. In the $\mu\to 0$ limit, the viscous time-scale
$\tau_{1}=m/\gamma$ controls the relaxation of the velocity. On the other
hand, when $m\to 0$, the time-scale $\tau_{2}=\gamma/\mu$ controls the
relaxation of the particle position in the harmonic trap. For non-zero $m$ and
$\mu$, the presence of the two time-scales $\tau_{1}$ and $\tau_{2}$ give rise
to two scenarios, namely,

  * •

Underdamped case: $\gamma^{2}<4m\mu$, i.e., $\tau_{2}<4\tau_{1}$,

  * •

Overdamped case: $\gamma^{2}\geq 4m\mu$, i.e., $\tau_{2}\geq 4\tau_{1}$,

where the particle exhibits qualitatively different kinds of motion [see Fig.
1]. In the presence of the stochastic active force, the particle shows a
rather interesting random motion [see Fig. 2 for typical trajectories]. In
fact, in this case, an additional time-scale $\tau_{a}$ appears, which
characterize the activity of the system. The relative strength of $\tau_{a}$
with respect to $\tau_{1}$ and $\tau_{2}$ and $\tau_{a}$ gives rise to a range
of sub-regimes for the IRTP dynamics, in both the overdamped and underdamped
scenarios. In this paper we characterize the position and velocity
fluctuations of the IRTP in all these regimes.

Figure 3: Plot for the joint stationary distribution $P_{\text{st}}(x,v)$ in
the strongly overdamped ($\tau_{1}\ll\tau_{2}$) [(a), (b), (c)] and
underdamped ($\tau_{2}\ll\tau_{1}$) [(d), (e)] scenarios, obtained using
numerical simulations. In both cases, the distributions show a shape
transition with changing activity (from left to right). The overdamped case is
characterized by three qualitatively distinct shapes of the joint distribution
$P_{\text{st}}(x,v)$ in three different activity regimes. In contrast, the
underdamped case is characterized by two activity regimes, resulting in two
qualitatively different behaviours of $P_{\text{st}}(x,v)$. The analytically
obtained bound for the joint distribution $P_{\text{st}}(x,v)$ in the (a)
overdamped [see Eqs. (55)-(56)] and (d) underdamped [see Eqs. (44)-(43)] cases
are also drawn on the $x-v$ plane with black-dashed line. Here, we have used
parameter values $m=0.1,\mu=1.0,a_{0}=1.0,\gamma=10.0$ and
$\tau_{a}=100.0,1.0,0.001$ (left to right) for the overdamped case and
$m=1.0,\mu=1.0,\gamma=1.0,a_{0}=1.0$ and $\tau_{a}=2.5,0.1$ (left to right)
for the underdamped case.

Note that, in the viscous limit $\tau_{1}\to 0$, Eq. (2) reduces to the well-
studied ordinary run-and-tumble particle (RTP) Malakar  _et al._ (2018) with
self-propulsion speed $\bar{v}_{0}$ in a harmonic trap of strength
$\bar{\mu}$, where,

$\displaystyle\bar{v}_{0}=a_{0}/\gamma,\leavevmode\nobreak\
\leavevmode\nobreak\ \text{and}\quad\bar{\mu}=\mu/\gamma.$ (3)

In this case, the stationary state position distribution, bounded in the
regime $-\bar{v}_{0}/\bar{\mu}\leq x\leq\bar{v}_{0}/\bar{\mu}$, undergoes a
shape transition at a critical activity $\tau_{a}^{\star}=\gamma/\mu$ Dhar
_et al._ (2019). The presence of inertia is expected to affect both position
and velocity fluctuations, which we investigate here.

Strongly Overdamped ($\tau_{1}\ll\tau_{2}$):
• Position bound $|x|\leq X_{\text{OD}}$
[see Eq. (12)] • Velocity bound $|v|\leq V_{\text{OD}}$
[see Eq. (13) ] |  Strongly active regime ($\tau_{2}\ll\tau_{a}$) [see Fig. 3(a)]
$\circ$ Position distribution: U-shaped with divergences at $x=\pm X_{\text{OD}}$. $\circ$ Velocity distribution: Single algebraic divergence at $v=0$. |  Moderately active regime ($\tau_{1}\ll\tau_{a}\ll\tau_{2}$) [see Fig. 3(b)]
$\circ$ Position distribution: Dome shaped with peak at $x=0$. $\circ$ Velocity distribution: Double peaked with minimum at $v=0$. |  Strongly passive regime ($\tau_{a}\ll\tau_{1}$) [see Fig. 3(c)]
$\circ$ Position distribution: Gaussian around $x=0$ $\circ$ Velocity
distribution: Gaussian around $v=0$
---|---|---|---
Strongly Underdamped ($\tau_{2}\ll\tau_{1}$):
• Position bound $|x|\leq X_{\text{UD}}$
[see Eq. (14)] • Velocity bound $|v|\leq V_{\text{UD}}$
[see Eq. (15)] |  Strongly active regime ($\tau_{1}\ll\tau_{a}$) [see Fig. 3(d)]
$\circ$ Position distribution: Multi-peaked with preliminary peaks at $x=\pm a_{0}/\mu$. $\circ$ Velocity distribution: Multi-peaked with preliminary peak at $v=0$. |  Strongly passive regime ($\tau_{a}\ll\tau_{1}$) [see Fig. 3(e)]
$\circ$ Position distribution: Gaussian around $x=0$ $\circ$ Velocity
distribution: Gaussian around $v=0$
Table 1: Summary of the different regimes emerging from the interplay of the
inertial time-scale $\tau_{1}=m/\gamma$, the viscous time-scale
$\tau_{2}=\gamma/\mu$ and the activity time-scale $\tau_{a}$.

We first discuss the overdamped case ($4\tau_{1}\leq\tau_{2}$). The linear
nature of the Langevin equations (2) allows us to solve it exactly, albeit
formally [see Appendix A for the details]. In the overdamped scenario, the
formal solution reads,

$\displaystyle x(t)=$
$\displaystyle\frac{2\bar{v}_{0}}{\lambda}\int_{0}^{t}ds\leavevmode\nobreak\
\sigma(t-s)\exp\Big{[}-\frac{s}{2\tau_{1}}\Big{]}\sinh\frac{\lambda
s}{2\tau_{1}},$ (4) $\displaystyle v(t)=$
$\displaystyle\frac{\bar{v}_{0}}{\tau_{1}}\int_{0}^{t}ds\,\sigma(t-s)\exp\Big{[}-\frac{s}{2\tau_{1}}\Big{]}\Big{\\{}\cosh\frac{\lambda
s}{2\tau_{1}}$ (5)
$\displaystyle\qquad\qquad\qquad\qquad-\frac{1}{\lambda}\sinh\frac{\lambda
s}{2\tau_{1}}\Big{\\}},$ (6)

where we have defined,

$\displaystyle\lambda=\sqrt{1-\frac{4\tau_{1}}{\tau_{2}}}\quad\text{with}\quad\tau_{1}=\frac{m}{\gamma},\leavevmode\nobreak\
\tau_{2}=\frac{\gamma}{\mu}.$ (7)

Here we assume, without any loss of generality, that at $t=0$ the particle is
at rest at the origin, i.e., $x(0)=0,v(0)=0$. Evidently, in this overdamped
regime, the combination of the two time-scales $\tau_{1}$ and $\tau_{2}$ gives
rise to the emergent time-scales $\tau_{1}/(1\pm\lambda)$. However, for the
strongly overdamped scenario $\tau_{1}\ll\tau_{2}$, these time-scales reduce
to the intrinsic ones $\tau_{1}$ and $\tau_{2}$, since in this limit,
$1+\lambda\simeq 2$ and $1-\lambda\simeq 2\tau_{1}/\tau_{2}$. The interplay of
the three time-scales, $\tau_{1},\tau_{2}$ and ${\tau_{a}}$ leads to the
emergence of three regimes of activity, namely, (i) the strongly active regime
($\tau_{1},\tau_{2}\ll\tau_{a}$), (ii) the moderately active
regime,($\tau_{1}\ll\tau_{a}\ll\tau_{2}$) and (iii) the passive regime
$\tau_{a}\ll\tau_{1},\tau_{2}$ [see Table 1]. We analyse the behaviour of the
position and velocity distributions in these regimes in Sec V.

Next, we discuss the underdamped oscillatory scenario ($\tau_{2}<4\tau_{1}$),
where the solution of Eq. (2) is given by,

$\displaystyle x(t)=$
$\displaystyle\frac{2\bar{v}_{0}}{\Lambda}\int_{0}^{t}ds\leavevmode\nobreak\
\sigma(t-s)\exp\Big{[}-\frac{s}{2\tau_{1}}\Big{]}\sin\frac{\Lambda
s}{2\tau_{1}},$ (8) $\displaystyle v(t)=$
$\displaystyle\frac{\bar{v}_{0}}{\tau_{1}}\int_{0}^{t}ds\,\sigma(t-s)\exp\Big{[}-\frac{s}{2\tau_{1}}\Big{]}\Big{\\{}\cos\frac{\Lambda
s}{2\tau_{1}}$ (9)
$\displaystyle\qquad\qquad\qquad\qquad-\frac{1}{\Lambda}\sin\frac{\Lambda
s}{2\tau_{1}}\Big{\\}},$ (10)

with,

$\displaystyle\Lambda=\sqrt{\frac{4\tau_{1}}{\tau_{2}}-1}.$ (11)

Here, once again, we have assumed $x(0)=0$ and $v(0)=0$.

Clearly, from Eqs. (8)-(10), the time-scale $\tau_{2}$, appearing through
$\Lambda$, only governs the oscillation period, and the relaxation is solely
determined by the time-scale $\tau_{1}$. Therefore, the interplay of this
time-scale $\tau_{1}$ and the active time-scale $\tau_{a}$ gives rise to only
two regimes $\tau_{1}\ll\tau_{a}$ (strongly active) and $\tau_{a}\ll\tau_{1}$
(strongly passive) [see Table 1].

The bounded nature of the active noise $\sigma(t)$ ensures that the stationary
distribution $P_{\text{st}}(x,v)$ is supported in a finite region of the $x-v$
phase plane. However, the supports are different in the overdamped and
underdamped cases. Figure 2(a) shows a typical trajectory of the IRTP in the
overdamped scenario which indicates that the particle position, in this case,
is always confined between the two possible trap centres $\pm a_{0}/\mu$. In
fact, we compute the boundary of the supporting region in the phase space [see
Fig. 14 and Eqs. (55)-(56)]. In particular, we show that the marginal position
and velocity distributions are supported in the regions,

$\displaystyle|x|\leq X_{\text{OD}}$ $\displaystyle\equiv\frac{a_{0}}{\mu},$
(12) $\displaystyle|v|\leq V_{\text{OD}}$
$\displaystyle\equiv\frac{2a_{0}}{\sqrt{m\mu}}\exp\Big{[}-\frac{1}{\lambda}\text{sech}^{-1}\frac{2\sqrt{m\mu}}{\gamma}\Big{]},$
(13)

respectively, where $\lambda$ is defined in Eq. (7). On the other hand, in the
underdamped case, the inertia allows the particle to move past the trap
centres, thereby extending the accessible region in the phase space [see Fig.
2(b)]. We show that, in this case, the bound depends on the mass of the
particle [see Fig. 6 and Eqs. (44)-(43)]. In particular, we find that the
marginal distributions are supported in the region,

$\displaystyle|x|\leq X_{\text{UD}}$
$\displaystyle\equiv\frac{a_{0}}{\mu}\coth\frac{\pi}{2\Lambda},$ (14)
$\displaystyle|v|\leq V_{\text{UD}}$
$\displaystyle\equiv\frac{a_{0}}{\sqrt{m\mu}}\Big{[}1+\coth\frac{\pi}{2\Lambda}\Big{]}\exp\bigg{[}-\frac{1}{\Lambda}\tan^{-1}\Lambda\bigg{]},$
(15)

where $\Lambda$ is defined in Eq. (11). Figure 2 illustrates typical
trajectories of the particle on the $x-v$ plane as well as time evolution of
$x(t)$ and $v(t)$, in both the scenarios.

The stationary distribution $P_{\text{st}}(x,v)$ in the underdamped and
overdamped scenarios, obtained from numerical simulations, are shown in Fig.
3. In both the cases, the distributions undergo intricate shape transitions,
from a multi-peaked structure in the strongly active regime to a single-peaked
one in the passive regime. However, the multi-peaked structures appearing in
the strongly active regime are very different in the overdamped and the
underdamped scenarios. To analytically characterize these shape transitions,
we focus on the marginal position and velocity distributions separately.

In the strongly active regime, the distributions can be obtained by
considering the deterministic motion of the particles between successive
tumbling events. Nevertheless, the difference in the nature of the
deterministic trajectories in the overdamped and underdamped cases leads to
drastically different shapes of the distributions. In the former case, the
peaks in the marginal position distribution always appear at the edges of the
support $\pm X_{\text{OD}}$ whereas in the latter case, multiple peaks appear
inside the support. On the other hand, the marginal distribution of velocity
in the overdamped case has a single peak at $v=0$ in contrast to a multiple
peaked structure in the underdamped case. These peaks smoothen out as the
activity is decreased and eventually, in the passive limit, both the position
and velocity distributions converge to the Boltzmann measure
$P_{\text{eq}}(x,v)\propto e^{-\beta_{\text{eff}}(U(x)+mv^{2}/2)}$ with
$\beta_{\text{eff}}=2\gamma/(\tau_{a}a_{0}^{2})$.

For the overdamped case, the transition from the strongly active to the
passive regime goes through an intermediate, moderately active, regime as
shown in Fig. 3 (a)-(c) [also see Table 1]. In particular, we compute the
position and velocity distributions analytically in the strongly overdamped
limit ($\gamma^{2}\gg m\mu$). The position distribution undergoes a transition
from a U-shape in the strongly active regime to a dome-shape in the moderately
active regime at a critical activity, $\tau_{a}^{\star}(m)\equiv
1/\kappa=\tau_{a}^{\star}[1-2\tau_{1}/\tau_{2}+O(\tau_{1}^{2})]$, where
$\kappa$ is defined in Eq. (65). Eventually, in the strongly passive regime,
the position distribution converges to the Gaussian. On the other hand, we
find the velocity distribution in the active regime as,

$\displaystyle P_{\text{st}}(v)=\frac{2^{\alpha}}{\nu
B(\alpha,\alpha)}\Big{(}\frac{|v|}{\nu}\Big{)}^{\alpha-1}\Big{(}1-\frac{|v|}{2\nu}\Big{)}^{\alpha},\text{with}\leavevmode\nobreak\
\leavevmode\nobreak\ \alpha=\frac{1}{\tau_{a}\kappa},$ (16)

where $B(u,u)$ denotes the beta function (DLMF, , (5.12.1)) and $\nu$ is
defined in Eq. (65). Clearly, the velocity distribution also undergoes a shape
transition at $\alpha=1$, i.e., $\tau_{a}^{\star}(m)=1/\kappa$, from a
distribution exhibiting algebraic divergence at $v=0$ with exponent
$1-1/(\tau_{a}\kappa)$, in the strongly active regime
$\tau_{a}^{\star}(m)\ll\tau_{a}$, to a bimodal distribution with a minimum at
$v=0$ in the moderately active regime. Interestingly, the velocity
distribution undergoes an additional shape transition to a dome-shaped
distribution, as the activity is decreased further, eventually becoming
Gaussian in the strongly passive regime [see Fig. 19]. Table 1 briefly
summarises the results in the different regimes and subregimes.

In the following sections, we analyze the position and velocity distributions
of the IRTP in the overdamped and underdamped regimes separately. We start
with computing the moments and correlations of position and velocity in the
next section.

##  III Moments

To get an idea about the nature of the position and velocity fluctuations, we
first investigate the behaviour of the moments of $x$ and $v$ in the
stationary state. Let $P_{\sigma}(x,v,t)$ denote the joint probability density
for the particle to have position $x$ and velocity $v$, along with self-
propulsion force direction $\sigma$, at time $t$. It is straightforward to see
that the set of Fokker-Planck equations describing the evolution of
$P_{\sigma}(x,v,t)$ are given by,

$\displaystyle\frac{\partial P_{\sigma}}{\partial t}$ $\displaystyle=$
$\displaystyle-\frac{\partial}{\partial
x}(vP_{\sigma})+\frac{1}{m}\frac{\partial}{\partial v}\big{[}(\gamma v+\mu
x-a_{0}\sigma)P_{\sigma}\big{]}$ (17)
$\displaystyle-(P_{\sigma}-P_{-\sigma})/\tau_{a}.$

It is difficult to solve this Fokker-Planck equation directly, even in the
stationary state. However, as we show below, the stationary correlations of
the form $\langle x^{k}v^{n}\rangle$, for any arbitrary values of $k$ and $n$,
can be computed recursively.

To this end, we define the quantities
$P(x,v,t)=\sum_{\sigma}P_{\sigma}(x,v,t)$ and $Q(x,v,t)=\sum_{\sigma}\sigma
P_{\sigma}(x,v,t)$. Recasting the Fokker Planck equations (17) in terms of $P$
and $Q$, and taking the stationary limit, $\partial_{t}P=0$ and
$\partial_{t}Q=0$ we get,

$\displaystyle\frac{a_{0}}{m}\frac{\partial Q}{\partial v}$
$\displaystyle=-v\frac{\partial P}{\partial
x}+\frac{\gamma}{m}\frac{\partial}{\partial
v}(vP)+\frac{\mu}{m}x\frac{\partial P}{\partial v},$ (18)
$\displaystyle\frac{a_{0}}{m}\frac{\partial P}{\partial v}$
$\displaystyle=-v\frac{\partial Q}{\partial
x}+\frac{\gamma}{m}\frac{\partial}{\partial
v}(vQ)+\frac{\mu}{m}x\frac{\partial Q}{\partial v}-\frac{2}{\tau_{a}}Q.$ (19)

We further define the correlation functions,

$\displaystyle M(k,n)$ $\displaystyle=$
$\displaystyle\intop_{-\infty}^{\infty}dx\intop_{-\infty}^{\infty}dv\,x^{k}v^{n}P(x,v),$
(20) $\displaystyle L(k,n)$ $\displaystyle=$
$\displaystyle\intop_{-\infty}^{\infty}dx\intop_{-\infty}^{\infty}dv\,x^{k}v^{n}Q(x,v),$
(21)

where $k$ and $n$ are non-negative integers. Clearly,
$M(k,n)\equiv\left\langle x^{k}v^{n}\right\rangle$ refer to the correlations
we are interested in. Multiplying both sides of Eqs. (18) and (19) with
$x^{k}v^{n}$, and integrating over $x$ and $v$, we get a set of coupled
recursion relations,

$\displaystyle\mu nM(k+1,n-1)=mkM(k-1,n+1)-\gamma nM(k,n)$ (22)
$\displaystyle+a_{0}nL(k,n-1),\leavevmode\nobreak\ \leavevmode\nobreak\ $ (23)
$\displaystyle\mu nL(k+1,n-1)=mkL(k-1,n+1)-\gamma nL(k,n)$ (24)
$\displaystyle+a_{0}nM(k,n-1)-2mL(k,n)/\tau_{a}.\leavevmode\nobreak\
\leavevmode\nobreak\ $ (25)

Note that, to get the above equations we have used the fact that both $P(x,v)$
and $Q(x,v)$ vanish at $x\to\pm\infty$ and $v\to\pm\infty$.

By definition, $M(k,n)=0$ for negative values of $k<0$ and $n<0$ and
normalization of $P(x,v)$ gives the boundary condition $M(0,0)=1$. We further
note that the Eqs. (18)-(19) are invariant under the transformation
$x\rightarrow-x,\,v\rightarrow-v,\,\sigma\rightarrow-\sigma$, which ensures,

$\displaystyle M(k,n)=0,$ $\displaystyle\qquad\forall(k+n)\in\text{odd},$ (26)
$\displaystyle L(k,n)=0,$ $\displaystyle\qquad\forall(k+n)\in\text{even}.$
(27)

Substituting $n=0$ in Eqs. (23)-(25), we obtain two more constraints on
$M(k,n)$ and $L(k,n)$,

$M(k,1)=0,\quad
L(k,0)=\frac{k\tau_{a}}{2}L(k-1,1)\quad\text{for}\leavevmode\nobreak\ k\geq
0.$ (28)

The relations Eqs. (23)-(28), along with the normalization condition,
$M(0,0)=1$, allow the computations of the correlations $M(k,n)$ and $L(k,n)$
recursively, as discussed below.

Figure 4: Schematic diagram illustrating the recursive connections between the
correlation function : $M(k,n)$ (red disks) and $L(k,n)$ (green squares) drawn
on the $k-n$ lattice.

Figure 4 shows the $k-n$ lattice where the $d$-th diagonal, given by $k+n=d$,
contains $L(k,n)$ or $M(k,n)$ for odd and even $d$, respectively. We start
with the first diagonal $k+n=1$ which contains the entries $L(1,0)$ and
$L(0,1)$. These quantities satisfy

$\displaystyle\mu L(1,0)=-\gamma L(0,1)+a_{0}M(0,0)-2mL(0,1)/\tau_{a},$ (29)

which is obtained by setting $k=1,n=0$ in Eq. (25), along with the boundary
condition $M(0,0)=1$. Moreover, we have, from Eq. (28)

$\displaystyle\mu L(1,0)=\frac{\tau_{a}}{2}L(0,1).$ (30)

The above equations can be solved immediately to obtain,

$\displaystyle L(0,1)$ $\displaystyle=$
$\displaystyle\frac{2a_{0}\tau_{a}}{(\mu\tau_{a}^{2}+2\gamma\tau_{a}+4m)},$
(31) $\displaystyle L(1,0)$ $\displaystyle=$
$\displaystyle\frac{a_{0}\tau_{a}^{2}}{(\mu\tau_{a}^{2}+2\gamma\tau_{a}+4m)}.$
(32)

Next, we move to the second diagonal $k+n=2$ which contains the entries
$M(2,0),M(1,1)$ and $M(0,2)$. Setting $k=1,n=1$ and $k=0,n=2$ in Eq. (23) we
find the equations satisfied by these quantities, which are given by,

$\displaystyle\mu M(2,0)$ $\displaystyle=mM(0,2)-\gamma M(1,1)+a_{0}L(1,0),$
(33) $\displaystyle\mu M(1,1)$ $\displaystyle=-\gamma M(0,2)+a_{0}L(0,1).$
(34)

Additionally, we have $M(1,1)=0$, from Eq. (28). These three equations can be
solved to obtain the correlations appearing on the second diagonal,

$\displaystyle\langle x^{2}\rangle$ $\displaystyle\equiv$ $\displaystyle
M(2,0)=\frac{a_{0}^{2}\tau_{a}(\gamma\tau_{a}+2m)}{\gamma\mu(\mu\tau_{a}^{2}+2\gamma\tau_{a}+4m)},$
(35) $\displaystyle\langle v^{2}\rangle$ $\displaystyle\equiv$ $\displaystyle
M(0,2)=\frac{2a_{0}^{2}\tau_{a}}{\gamma(\mu\tau_{a}^{2}+2\gamma\tau_{a}+4m)},$
(36)

where we have used Eq. (32) obtained in the previous step. One can proceed in
a recursive manner to compute the higher order correlations. In general, the
$d$-th diagonal contains $d+1$ entries, which satisfy $d$ equations, obtained
from Eq. (23)-(25), along with an additional constraint from Eq. (28). These
$(d+1)$ coupled equations can be solved to obtain the correlations $M(k,d-k)$
or $L(k,d-k)$.

It is worthwhile to look at the stationary variance of the position and
velocity, given by Eqs. (36). Interestingly, the position variance $\langle
x^{2}\rangle$ shows a non-monotonic behaviour as the activity $\tau_{a}$
increases in the regime $\tau_{2}<\tau_{1}$. As the activity $\tau_{a}$ is
increased, the variance first increases, attaining a maximum value at
$\tau_{a}=2m/(\sqrt{m\mu}-\gamma)$. The variance then decreases with
$\tau_{a}$, eventually becoming $a_{0}^{2}/\mu^{2}$ asymptotically in the
limit $\tau_{a}\to\infty$. This is in contrast with the ordinary RTP, i.e.,
the $m=0$ case, where the variance of the position increases monotonically
with $\tau_{a}$. On the other hand, the velocity variance $\langle
v^{2}\rangle$ always shows a non-monotonic behaviour with $\tau_{a}$ with a
maximum at $\tau_{a}=\sqrt{4m/(\mu)}$.

It is also useful to look at the time evolution of the position and velocity
correlations. To this end, we compute the full time-dependent correlations
$\langle x^{2}(t)\rangle$, $\langle v^{2}(t)\rangle$ and $\langle
x(t)v(t)\rangle$ from the formal solution of the Langevin equations (78)-(80).
The details of this calculation is provided in Appendix B, here we only
mention the notable features. At short-times, i.e., for
$t\ll\min(\tau_{1},\tau_{2},\tau_{a})$, we have,

$\displaystyle\langle x^{2}(t)\rangle$
$\displaystyle=\big{(}\frac{a_{0}}{2m}\big{)}^{2}t^{4}+O(t^{5}),$ (37)
$\displaystyle\langle x(t)v(t)\rangle$
$\displaystyle=\frac{1}{2}\big{(}\frac{a_{0}}{m}\big{)}^{2}t^{3}+O(t^{4}),$
(38) $\displaystyle\langle v^{2}(t)\rangle$
$\displaystyle=\big{(}\frac{a_{0}}{m}\big{)}^{2}t^{2}+O(t^{3}).$ (39)

On the other hand, at intermediate times, the correlations show oscillatory
behaviour in the underdamped case, which is absent in the overdamped case.
Figure 5 shows the time evolution of the second order correlations, $\langle
x^{2}(t)\rangle$, $\langle x(t)v(t)\rangle$ and $\langle v^{2}(t)\rangle$
obtained from numerical simulations, using Euler discretization, along with
the corresponding analytical expressions Eqs. (83)-(91) for the underdamped
and the overdamped cases.

Figure 5: Time evolution of the moments in the (a) underdamped and (b)
overdamped scenarios. The symbols indicate the data obtained from numerical
simulations while the solid lines indicate the analytical predictions [see
Eqs. (83)-(91)]. The insets illustrate the anomalous growth at short-times
[see Eqs. (37)-(39)] using the same data plotted in log-scales. Here we have
used: $\tau_{a}=1.0,\mu=1.0,m=1.0,\gamma=0.5$ and $a_{0}=1.0$ for the
underdamped case, and $\tau_{a}=1,\mu=1,m=1.0,\gamma=5.0$ and $a_{0}=1.0$ for
the overdamped case.

##  IV Underdamped Case

In this section, we study the marginal position and velocity distributions for
the underdamped scenario $\gamma^{2}<4m\mu$, i.e. when $4\tau_{1}>\tau_{2}$.
In this case, between successive tumbling events, where $\sigma$ remains
constant, the particle undergoes a deterministic underdamped motion in a
harmonic trap centred at $a_{0}\sigma/\mu$ [see Eqs. (8)-(10)]. Each tumbling
event results in a sudden change of the position of the trap minimum, i.e.,
$a_{0}\sigma/\mu\to-a_{0}\sigma/\mu$, with the resulting motion now directed
towards the new minimum [see Fig. 2]. It is evident from Eqs. (8)-(10), this
underdamped run-and-tumble motion is characterized by the inertial time-scale
$\tau_{1}=m/\gamma$ associated with the relaxation of the particle in a
harmonic trap, and the active time-scale $\tau_{a}$ associated with the
switching of the trap position. The interplay between these two time-scales
gives rise to certain novel features for the joint stationary distribution
$P_{\text{st}}(x,v)$ [see Fig. 3 (d),(e)]. To understand this unique behaviour
in the underdamped regime, we discuss the marginal position and velocity
distributions separately.

Figure 6: (a) Underdamped scenario: Typical displacement maximizing trajectory
of the particle shows that the particle settles into a limit-cycle (black
dashed line) described by the trajectories Eqs. (43)-(44). Here we have used
$\mu=1.0,a_{0}=1.0,\gamma=0.5$ and $m=1.5$. (b) Trajectories [see Eqs.
(43)-(44)] bounding the joint stationary distribution $P_{\text{s}t}(x,v)$ for
different values of mass $m$. Here we have used $\mu=1.0,a_{0}=1.0$ and
$\gamma=0.5$.

We start by computing the finite supports of the distributions which arise due
to the bounded nature of the active force $\sigma(t)$. To compute these bounds
we consider the trajectory which maximizes the displacement of the particle.
The displacement of an underdamped particle from its initial position $x_{0}$
is maximum when the particle reaches the first turning point, i.e. where its
velocity changes sign for the first time. Thus, for an IRTP, the displacement
will be maximized if the tumbling events $\sigma\to-\sigma$ coincide with the
velocity $v$ changing its sign. Hence, the bound on the position fluctuation
can be obtained by considering such a trajectory [see Fig. 6(a)].

Let $x_{\ell}$ denote the particle position at the $\ell$-th turning point of
this trajectory. By definition, the $\ell$-th turning point coincides with the
$\ell$-th tumbling event, $\sigma_{\ell}\to-\sigma_{\ell-1}$. Considering the
harmonic nature of the motion between successive tumbling events, it is
straightforward to show that $x_{\ell}$ satisfies the recursion relation [see
Appendix C for the details],

$\displaystyle x_{\ell+1}=cx_{\ell}+\frac{a_{0}}{\mu}(1-c)\sigma_{\ell},$ (40)

where $c=-e^{-\frac{\pi}{\Lambda}}$ and $\Lambda=\sqrt{4m\mu/\gamma^{2}-1}$
[see Eq. (11)]. Without any loss of generality, we assume that the particle
starts at rest from position $|x_{0}|<a_{0}/\mu$ with propulsion direction
$\sigma_{0}$. Then, the solution of this recursion relation is given by,

$\displaystyle
x_{\ell}=c^{\ell}x_{0}+\frac{a_{0}\sigma_{0}}{\mu}\Big{[}\frac{1-c}{1+c}\Big{]}\big{[}(-1)^{\ell-1}+c^{\ell}\big{]},$
(41)

where we have used the fact that $\sigma_{\ell}=(-1)^{\ell}\sigma_{0}$. In
that case, it is straightforward to show that $|x_{\ell}|-|x_{\ell-1}|>0$ for
all $\ell>0$, i.e., the distance of the particle from the origin increases at
successive turning points. Since $c<1$, $x_{\ell}$ approaches a finite value
in the limit $\ell\to\infty$. Taking this limit, $\lim_{l\to\infty}x_{\ell}$,
gives the maximum distance the particle can reach, which by definition is the
bound for the position distribution quoted in Eq. (14). Note that, this bound
on the position distribution of an IRTP has also been obtained earlier by
considering the large deviation function Smith and Farago (2022).

Figure 7: Underdamped scenario: Plot of position distribution
$P_{\text{st}}(x)$ (a) and velocity distribution $P_{\text{st}}(v)$, obtained
from numerical simulations, for different values of activity $\tau_{a}$. Here
we have used $\gamma=1.0,\mu=1.0$ and $m=2.0$.

The bound on the velocity fluctuations can also be computed from this
trajectory, following a similar argument. We obtain this bound explicitly
[quoted in Eq. (15)] by computing the velocity of the particle where the
acceleration changes sign in the trajectory [see Appendix C for the details].
Note that, in the strongly underdamped limit $m\rightarrow\infty$, the
position bound $X_{\text{UD}}\sim\sqrt{m}$, i.e., position distribution
becomes unbounded, whereas the bound for the velocity distribution remains
finite, with $V_{\text{UD}}\to 4a_{0}/\gamma\pi$.

In fact, it is easy to see that the particle motion is bounded by the
trajectories which connect the points $(\pm X_{\text{UD}},0)$. This allows us
to compute the bound of the joint distribution by computing these trajectories
parametrically. Using Eqs. (78)-(80) we obtain these trajectories explicitly
by setting $\sigma=\pm 1$, $x_{0}=\mp X_{\text{UD}}$ and $v_{0}=0$ as the
initial conditions,

$\displaystyle X_{\text{UD}}(t)$
$\displaystyle=\pm\frac{a_{0}}{\mu}\mp\frac{a_{0}}{\mu}e^{-\frac{t}{2\tau_{1}}}(1+\coth\frac{\pi}{2\Lambda})$
(42) $\displaystyle\qquad\qquad\times\big{[}\frac{1}{\Lambda}\sin\frac{\Lambda
t}{2\tau_{1}}+\cos\frac{\Lambda t}{2\tau_{1}}\big{]},$ (43) $\displaystyle
V_{\text{UD}}(t)$
$\displaystyle=\frac{2\bar{v}_{0}}{\Lambda}e^{-\frac{t}{2\tau_{1}}}\sin\frac{\Lambda
t}{2\tau_{1}}(1+\coth\frac{\pi}{2\Lambda}),$ (44)

where $0\leq t\leq 2\tau_{1}\pi/\Lambda$. Figure 6 shows typical trajectories
describing the boundary of the joint stationary distribution in the
underdamped scenario for different values of particle mass. In the following,
we separately discuss the behaviour of the position and the velocity
distributions.

###  IV.1 Position Distribution

Figure 8: Underdamped strongly active regime: Typical trajectory of the IRTP
in the phase space. The inset shows the turning points in the trajectory
marked using symbols—disks for the position turning points; and squares for
the velocity turning points. Here we have used
$\gamma=1.0,\mu=1.0,a_{0}=1.0,m=4.0$ and $\tau_{a}=100$.

In this section, we investigate the behaviour of the marginal position
distribution of the IRTP in the underdamped case. The stationary position
distribution $P_{\text{st}}(x)$ undergoes a shape transition with change in
the activity $\tau_{a}$. For $\tau_{1}\ll\tau_{a}$, i.e., in the strongly
active regime, the position distribution shows a novel shape — it has two
strong peaks at $x=\pm a_{0}/\mu$, accompanied by several smaller peaks and
shoulder-like structures near the tails. These peaks smooth out as $\tau_{a}$
is decreased, eventually giving rise to a dome-shaped distribution in the
passive regime $\tau_{a}\ll\tau_{1}$. This shape transition is illustrated in
Fig. 7 (a) which shows $P_{\text{st}}(x)$, measured from numerical
simulations, for different values of activity $\tau_{a}$.

To understand the multi-peaked shape of the position distribution in the
strongly active regime, we again consider the trajectory of the IRTP. Let us
recall that, for a constant $\sigma$, i.e., between two successive tumbling
events, the IRTP undergoes a deterministic underdamped motion in a harmonic
trap centred at $a_{0}\sigma/\mu$. For $\tau_{1}\ll\tau_{a}$, the typical time
between two successive tumbling events is much larger compared to the time
required by the particle to reach the trap minimum. Thus, in the strongly
active regime, the trajectory of the IRTP consists of alternate segments where
an underdamped harmonic oscillator starting at rest from $x=a_{0}\sigma/\mu$
relaxes to $x=-a_{0}\sigma/\mu$, and stays there until the next tumbling event
[see Fig. 8 for a typical trajectory]. Moreover, before reaching the trap
minimum, the particle undergoes an oscillatory motion around it, with its
velocity changing sign at the turning points. Correspondingly, the particle
predominantly resides near these turning points, giving rise to the multi-
peaked structure of the distribution.

Figure 9: Underdamped strongly active regime: Plot of the velocity
distribution $P_{\text{st}}(v)$ obtained from numerical simulations of IRTP
(symbols) and the approximate deterministic dynamics (dashed black lines), for
different values of mass $m$. The inset compares the shoulder-like structures
observed in IRTP (symbols) with the same obtained by incorporating a uniformly
distributed waiting time distribution in the approximate deterministic
dynamics (solid lines). Here we have used parameter values
$\mu=1.0,\gamma=1.0,a_{0}=1.0$ and $\tau_{a}=100$.

To quantitatively estimate this position distribution, we assume that, in the
strongly active regime, the tumbling dynamics can be approximated by a
deterministic one, with flips taking place at a fixed interval $\tau_{a}$. We
measure the position distribution of this approximate deterministic dynamics
using numerical simulations, which is shown in Fig. 9 along with the same
obtained from the actual IRTP dynamics. It is apparent that the two show an
excellent agreement—the approximate dynamics correctly reproduces the multi-
peaked shape of the position distribution (black dashed lines) except for the
shoulder-like structures at the tails.

This approximate dynamics allows us to locate the position of the peaks,
which, as mentioned already, are nothing but the turning points of the
deterministic trajectory of the particle travelling between $\pm a_{0}/\mu$.
These can be obtained systematically by setting $v=0$ in the equation of the
trajectory starting with $x_{0}=\pm a_{0}/\mu,v_{0}=0$ with $\sigma=\mp 1$.
The details of this computation is provided in Appendix C, here we only quote
the main result—in the strongly active regime, the position distribution is
peaked around the points,

$\displaystyle
x^{*}_{n}=\pm\frac{a_{0}}{\mu}\Big{[}1-2(-1)^{n}e^{-\frac{n\pi}{\Lambda}}\Big{]},\leavevmode\nobreak\
\leavevmode\nobreak\ \text{with}\leavevmode\nobreak\ \leavevmode\nobreak\
n=1,2,3\dots$ (45)

where $\Lambda$ is defined in Eq. (11). Clearly, in the $n\to\infty$ limit,
$x^{*}$ converges to the primary peak positions $\pm a_{0}/\mu$.

Note that, this approximate dynamics ignores the effect of possible tumbling
events occurring while the particle is traversing between $\pm a_{0}/\mu$. The
effects of such events can be incorporated into the approximate dynamics by
allowing the waiting times to be drawn from a uniform distribution in
$[0,2\tau_{a}]$ which has the same mean $\tau_{a}$ as the original exponential
distribution in the IRTP dynamics. The inset in Fig. 9 shows that this
inclusion of minimal stochasticity reproduces the shoulder-like structures
qualitatively.

Figure 10: Underdamped strongly passive regime: Plot showing the stationary
position distribution $P_{\text{st}}(x)$ of the IRTP for different values of
activity $\tau_{a}$. The symbols are for simulation results and the solid
lines stand for the corrected Gaussian given by Eq. (48). Here we have used
$\mu=1.0,\gamma=3.0$ and $m=10$.

As $\tau_{a}$ is decreased (below $\tau_{1}$), the position distribution
eventually becomes dome-shaped with a single peak at the origin [see Fig. 7].
In the strongly passive limit, i.e., $\tau_{a}\ll\tau_{1}$, the dichotomous
noise $\sigma(t)$ behaves like a Gaussian white noise with strength
$a_{0}^{2}/(2\gamma)$ [see Santra  _et al._ (2021)] Correspondingly, Eq. (8)
implies that the typical position distribution also becomes a Gaussian with
width $a_{0}\sqrt{\tau_{a}/2\mu\gamma}$. For any finite $\tau_{a}$,
corrections to the Gaussian appear at the tails. In the following, we show
that the corrections to the Gaussian can be systematically obtained from the
higher moments of the position.

Figure 11: Velocity distribution in the underdamped case in the strongly
active regime: Plot of $P_{\text{st}}(v)$ obtained from numerical simulations
of IRTP (symbols) and the approximate deterministic dynamics (dashed black
lines), for different values of mass $m$. The inset compares the shoulder-like
structures observed in IRTP (symbols) with the same obtained by incorporating
a uniformly distributed waiting time distribution in the approximate
deterministic dynamics (solid lines). Here we have taken
$\mu=1.0,\gamma=1.0,a_{0}=1.0$ and $\tau_{a}=100$.

To this end, we recall that the cumulant generating function (CGF) of the
position distribution, which is symmetric around $x=0$, is given by,

$\displaystyle Q(q)=-\sum_{k=1}^{\infty}\frac{q^{2k}}{(2k)!}\langle
x^{2k}\rangle_{c},$ (46)

where $\langle x^{2k}\rangle_{c}$ denote the corresponding even cumulants. The
CGF, in turn, is related to the moment-generating function,
$\tilde{P}(q)=\langle e^{iqx}\rangle$ by,

$\displaystyle Q(q)=\log\tilde{P}(q).$ (47)

If all the stationary cumulants of the distribution are known, then the exact
stationary distribution can, in principle, be computed using Eqs. (46)-(47),
and taking the inverse Fourier transform of $\tilde{P}(q)$. Clearly, if one
considers only $k=1$ term in Eq. (46), up to the variance $\langle
x^{2}\rangle$, this leads to a Gaussian distribution. Corrections to this
Gaussian distribution can be systematically computed by truncating Eq. (46) by
keeping upto $k=N$ terms.

To this end, we compute the first $N$ non-zero cumulants of the position
distribution from the corresponding moments, which, in turn, are computed
using the recursive procedure illustrated in Sec. III. Using this truncated
CGF, we then calculate $\tilde{P}(q)$ from Eq. (47). Numerically inverting the
Fourier transform gives us the position distribution.

$\displaystyle
P_{\text{st}}(x)=\intop_{-\infty}^{\infty}\frac{dq}{2\pi}\leavevmode\nobreak\
\tilde{P}(q)\,\cos{qx}.$ (48)

Figure 10 shows the distribution obtained keeping $N=5$ terms, which shows an
excellent agreement with the data obtained from numerical simulations. Note
that, the $P_{\text{st}}(x)$ shows significant deviation from the Gaussian as
the activity $\tau_{a}$ is increased.

It should be noted that, this method works well only in the passive regime,
where the contributions from the higher order cumulants decrease
progressively. In the active regime, where the distribution has a multiple
peaks away from the origin, significant contributions arise from the higher
order cumulants, and this method is not expected to work.

###  IV.2 Velocity Distribution

Figure 12: Underdamped strongly passive regime: Plot showing the stationary
velocity distribution $P_{\text{st}}(v)$ of the IRTP for different values of
activity $\tau_{a}$. The symbols are for simulation results and the solid
lines stand for the corrected Gaussian given by Eq. (53). Here we have used
$\mu=1.0,\gamma=3.0$ and $m=10$.

In this section, we investigate the nature of the velocity fluctuations in the
underdamped scenario. The stationary velocity distribution also shows a shape
transition similar to the position distribution— for $\tau_{1}\ll\tau_{a}$,
i.e., in the strongly active regime, the distribution has a primary peak at
$v=0$, with numerous smaller peaks, and shoulder-like structures near the
tail. Decreasing activity $\tau_{a}$, is marked by the smoothing of the peaks,
till a dome-shaped distribution emerges in the passive regime
$\tau_{a}\ll\tau_{1}$. This shape transition is illustrated in Figure 7(b)
which shows $P_{\text{st}}(v)$, measured from numerical simulations, for
different values of activity $\tau_{a}$.

The novel multi-peaked structure of the velocity distribution in the strongly
active regime can also be understood using the trajectory-based picture,
similar to the position distribution. As already discussed, in this limiting
scenario, the typical time between two successive tumbling events is much
larger than the time required to travel between the points $(\pm
a_{0}/\mu,0)$. From Fig. 8, it is clear that the velocity also undergoes
oscillations before the particle eventually comes to rest at the trap minimum
$a_{0}\sigma/\mu$. Intuitively, one expects the velocity distribution to peak
around the points where the acceleration is minimum, i.e., where
$\dot{v}(t)=0$.

The deterministic dynamics approximation discussed in the context of position
distribution can also be used to estimate the velocity distribution
quantitatively. We measure the velocity distribution by numerically simulating
this approximate dynamics, which is shown in Fig. 11, along with the same
obtained from the actual IRTP dynamics. Clearly, the deterministic dynamics
reproduces the multi-peaked shape of the velocity distribution fairly
accurately (black dashed lines), except for the shoulder-like structures at
the tail.

As before, we use this deterministic dynamics to locate the peaks in the
velocity distribution. These can be obtained by setting $\dot{v}(t)=0$ in the
trajectory starting with $x_{0}=\pm a_{0}/\mu,v_{0}=0$ and for $\sigma=\mp 1$
[see appendix C], and are given by,

$\displaystyle v^{\star}_{n}$
$\displaystyle=\pm\frac{2a_{0}}{\sqrt{m\mu}}\exp\Big{[}-\frac{1}{\Lambda}\big{(}n\pi+\tan^{-1}\Lambda\big{)}\Big{]}$
(49) $\displaystyle\hskip 85.35826pt\leavevmode\nobreak\ \leavevmode\nobreak\
\text{with}\leavevmode\nobreak\ \leavevmode\nobreak\ n=0,1,2\dots$ (50)

Moreover, the shoulder-like structures near the tails in $P_{\text{st}}(v)$
can also be qualitatively reproduced by allowing tumbling events to occur
while the particle is traversing, by drawing the waiting times from the
uniform distribution in $[0,2\tau_{a}]$, similar to the position distribution.
This is illustrated in the inset of Fig. 11.

Figure 13: Overdamped scenario: Plot of the position distribution
$P_{\text{st}}(x)$ (a) and the velocity distribution $P_{\text{st}}(v)$,
obtained from numerical simulations, for different values of activity
$\tau_{a}$. Here we have used $\mu=0.5,m=0.5,\gamma=5$ and $a_{0}=1.0$.

We now turn to the strongly passive regime $\tau_{a}\ll\tau_{1}$. In this
regime, the typical velocity distribution approaches a Gaussian. Similar to
the position distribution, the corrections to the Gaussian, for any non-zero
$\tau_{a}$, can be systematically obtained by considering the contributions
from higher order cumulants of the velocity, in the CGF,

$\displaystyle Q_{v}(q)=-\sum_{k=1}^{\infty}\frac{q^{2k}}{(2k)!}\langle
v^{2k}\rangle_{c}.$ (51)

Consequently, the velocity distribution can be computed by taking an inverse
Fourier transform of the velocity moment-generating function, which is related
to the CGF by,

$\displaystyle Q_{v}(q)=\log\tilde{P}_{v}(q).$ (52)

To this end, we construct a truncated CGF by computing the first $N$ non-zero
cumulants of velocity from the corresponding moments. These moments, in turn,
are obtained using the recursive procedure illustrated in Sec. III.
Numerically inverting the Fourier transform gives us the velocity
distribution,

$\displaystyle
P_{\text{st}}(v)=\intop_{-\infty}^{\infty}\frac{dq}{2\pi}\tilde{P}_{v}(q)\cos
qv.$ (53)

Figure 12 compares the velocity distribution obtained using numerical
simulations, with the analytical prediction (53) with $N=5$, which shows an
excellent agreement.

##  V Overdamped Case

The overdamped scenario emerges for $\gamma^{2}\geq 4m\mu$, i.e., when
$4\tau_{1}\leq\tau_{2}$. In the fully overdamped scenario, which arises in the
$\tau_{1}\to 0$ limit, the IRTP dynamics reduces to that of an ordinary RTP
with self-propulsion speed $\bar{v}_{0}=a_{0}/\gamma$ in a trap of strength
$\bar{\mu}=\mu/\gamma$. The stationary position distribution, in this limit,
is bounded in the interval $[-a_{0}/\mu,a_{0}/\mu]$, and is given by Dhar  _et
al._ (2019),

$\displaystyle
P_{\text{st}}(x)=\frac{2^{1-2\gamma/(\tau_{a}\mu)}}{B(\frac{\gamma}{\tau_{a}\mu},\frac{\gamma}{\tau_{a}\mu})}\frac{\mu}{a_{0}}\Big{[}1-\Big{(}\frac{\mu
x}{a_{0}}\Big{)}^{2}\Big{]}^{\frac{\gamma}{\tau_{a}\mu}-1},$ (54)

where $B(u,u)$ is the beta function. This distribution undergoes a shape
transition at the critical activity $\tau_{a}^{\star}=\tau_{2}=\gamma/\mu$. In
the active regime $\tau_{2}<\tau_{a}$, the distribution has a U-shape with
algebraic divergences occurring at the two ends $x=\pm a_{0}/\mu$. On the
other hand, in the passive regime $\tau_{a}<\tau_{2}$, a dome-shaped
distribution emerges with $P_{\text{st}}(x)$ vanishing at the two ends.

Figure 14: Overdamped scenario: Typical phase trajectories bounding the joint
stationary distribution $P_{\text{s}t}(x,v)$ for different values of mass $m$
using analytic predictions from Eqs. (55)-(56). Here we have used
$\mu=1.0,\gamma=5.0$ and $a_{0}=1.0$.

In this section, we investigate how the presence of a finite mass $m$ affects
the position and velocity distributions in the overdamped regime of an IRTP.
We start by investigating how the finite support of the stationary
distribution changes in the overdamped scenario. A typical trajectory of the
overdamped IRTP is shown in Fig. 2. Clearly, in this case, the maximum
displacement of the particle corresponds to the trajectory with no tumbling
events which eventually takes the particle to the trap minimum
$a_{0}\sigma/\mu$, without allowing it to cross it. This implies that, in the
stationary state, the position of the particle is bounded in the region
$-a_{0}/\mu<x<a_{0}/\mu$. A more rigorous proof of this statement is provided
in Appendix C. This bound can also be obtained taking the limit
$m\to\gamma^{2}/(4\mu)$, in Eqs. (14), which gives Eq. (12). It is noteworthy
that the position bound in this case does not depend on the mass of the IRTP.

Figure 15: Overdamped scenario: Plot of the position distribution
$P_{\text{st}}(x)$ obtained from numerical simulations of the IRTP (symbols)
along with the approximate distribution given by Eq. (67) (solid line) for
different values of activity $\tau_{a}$. Here we have taken
$m=0.2,\mu=1.0,a_{0}=1.0$ and $\gamma=10$.

The bound in the phase-space in the overdamped case can be obtained following
an argument similar to the one used in the underdamped case. The particle
motion must be bounded by the trajectories connecting the phase points $(\pm
X_{\text{OD}},0)$ and the velocity bound is given by the maximum velocity
attained by the particle in these trajectories. By setting $\sigma=\pm
1,x_{0}=\mp X_{\text{OD}}$ and $v_{0}=0$ in Eq. (80), the bound on the
velocity fluctuations is explicitly obtained [see appendix C for details]
which is quoted in Eq. (13). Interestingly, unlike the position, the velocity
bound explicitly depends on the mass of the particle. The boundary for the
joint distribution on the $x-v$ plane is given by setting $\sigma=\pm
1,x_{0}=\mp X_{OD}$ and $v_{0}=0$ simultaneously in Eqs. (78)-(80),

$\displaystyle X_{\text{OD}}(t)$
$\displaystyle=\pm\frac{\bar{v}_{0}}{\tau_{2}}\Bigg{[}1-2e^{-{t}/{(2\tau_{1})}}\left(\cosh\frac{\lambda
t}{2\tau_{1}}+\frac{1}{\lambda}\sinh\frac{\lambda
t}{2\tau_{1}}\right)\Bigg{]},$ (55) $\displaystyle V_{\text{OD}}(t)$
$\displaystyle=\pm\frac{4\bar{v}_{0}}{\lambda}e^{-{t}/(2\tau_{1})}\sinh\frac{\lambda
t}{2\tau_{1}},$ (56)

where $0\leq t<\infty.$ Fig. 14 shows typical trajectories describing the
boundary of the joint stationary distribution in the overdamped scenario for
different values of particle mass.

We now turn to the marginal distributions in the overdamped case. Figure 13(a)
shows $P_{\text{st}}(x)$ measured from numerical simulations, for different
values of activity $\tau_{a}$. Clearly, the position distribution of the IRTP
shows qualitatively similar behaviour to the ordinary RTP—it undergoes a
transition from a U-shape to a dome shape as the activity is decreased.
However, the presence of a finite mass, i.e., non-zero $\tau_{1}$ changes the
distribution quantitatively, including the critical value of activity, which
now depends on the mass of the particle. Figure 13(b) illustrates the velocity
distribution $P_{\text{st}}(v)$ of the IRTP measured from numerical
simulations for different activity $\tau_{a}$. As evident, velocity
distribution, shows a very different behaviour than in the underdamped case
[see Fig. 7(b)]. In particular, it undergoes two distinct shape transitions
with decreasing activity, first from a unimodal distribution peaked at $v=0$
to a bimodal distribution with minimum at $v=0$, and then from the bimodal
distribution to a Gaussian-like distribution. To further understand these
behaviours in the overdamped regime, we consider the case of the strongly
overdamped scenario $\tau_{1}\ll\tau_{2}$.

From Eqs. (4)-(6), it is apparent that the particle motion in the overdamped
scenario is characterized by two inertial time-scales $\tau_{1}/(1+\lambda)$
and $\tau_{1}/(1-\lambda)$, where $\lambda$ is defined in Eq. (7), along with
the active time-scale $\tau_{a}$. In the strongly overdamped regime, i.e., for
$\gamma^{2}\gg 4m\mu$ the inertial time-scales reduce to $\tau_{1}=m/\gamma$
and $\tau_{2}=\gamma/\mu$ which are well separated. Consequently, three
regimes of activity emerge, namely, strongly active
$\tau_{a}\gg\tau_{2}\gg\tau_{1}$, moderately active regime
$\tau_{2}\gg\tau_{a}\gg\tau_{1}$ and the passive regime
$\tau_{2}\gg\tau_{1}\gg\tau_{a}$. In the following, we separately explore the
behaviour of the position and velocity distributions in these regimes.

Figure 16: Overdamped strongly passive regime: Plot showing the stationary
position distribution $P_{\text{st}}(x)$ for different values of activity
$\tau_{a}$. The symbols are for simulation results and the solid lines stand
for the corrected Gaussian given by Eq. (48). Here we have used $m=0.5,\mu=1$
and $\gamma=5$.

###  V.1 Position Distribution

Figure 17: Plot showing the dependence of the shape transition point
$\tau_{a}^{\star}$ of the stationary position distribution on the viscous
time-scale $\tau_{1}$ for different values of trap strength $\mu$ and for
$a_{0}=1.0,\gamma=10$.

In this section, we analytically characterize the position distribution in the
strongly overdamped regime $\tau_{1}\ll\tau_{2}$. We start from Eq. (4), the
formal solution for $x(t)$ in the overdamped case. It is straightforward to
see that Eq. (4) can be recast as,

$\displaystyle x(t)=x_{1}(t)-x_{2}(t),$ (57)

where, $\\{x_{1},x_{2}\\}$ are two ordinary RTPs, propelled by the same active
noise $\sigma(t)$ with strength $\bar{v}_{0}/\lambda$, but moving in two
harmonic traps of strengths $\mu_{1}=(1-\lambda)/2\tau_{1}$ and
$\mu_{2}=(1+\lambda)/2\tau_{1}$, respectively. Thus, $x_{1}(t)$ and $x_{2}(t)$
evolve according to the Langevin equations,

$\displaystyle\dot{x}_{1}(t)$
$\displaystyle=-\mu_{1}x_{1}(t)+\bar{a}_{0}\sigma(t),$ (58)
$\displaystyle\dot{x}_{2}(t)$
$\displaystyle=-\mu_{2}x_{2}(t)+\bar{a}_{0}\sigma(t).$ (59)

It is convenient to introduce the variable $z(t)\equiv x_{1}(t)+x_{2}(t)$, and
write the Langevin equations for $x(t)$ and $z(t)$ using Eq. (59),

$\displaystyle\dot{x}(t)=$ $\displaystyle-k_{+}x(t)+k_{-}z(t),$ (60)
$\displaystyle\dot{z}(t)=$ $\displaystyle-
k_{+}z(t)+k_{-}x(t)+2\bar{a}_{0}\sigma(t),$ (61)

where,

$\displaystyle k_{\pm}=\frac{1}{2}(\mu_{2}\pm\mu_{1}).$ (62)

Using the above equations, we compute the distribution of $x$ in the strongly
overdamped limit ($\tau_{2}\ll\tau_{1}$) perturbatively, to the leading order
in $\tau_{1}$. The details of this computation is provided in the Appendix D,
here we sketch the main steps. From Eq. (60), we can express $x$ in terms of
$z$ and $\dot{z}$,

$\displaystyle
x(t)=z(t)-2\Big{(}\frac{z(t)}{\tau_{2}}+\dot{z}\Big{)}\tau_{1}+O(\tau_{1}^{2}).$
(63)

Moreover, we find that, to this order, $z$ satisfies a Langevin equation,

$\displaystyle\dot{z}(t)\approx-\kappa z(t)+\nu\sigma(t),$ (64)

where,

$\displaystyle\kappa\approx\frac{1}{\tau_{2}}\Big{(}1+\frac{2\tau_{1}}{\tau_{2}}\Big{)},\leavevmode\nobreak\
\leavevmode\nobreak\ \text{and}\leavevmode\nobreak\ \leavevmode\nobreak\
\nu=\bar{v}_{0}\tau_{2}\kappa.$ (65)

Clearly, Eq. (64) describes the motion of an ordinary RTP in a harmonic trap
of strength $\kappa$ and self-propulsion speed $\nu$. The corresponding joint
probability distribution $P_{\sigma}(z)$ is known exactly Dhar  _et al._
(2019) which can be used to obtain the distribution of $x$ [see Appendix D for
the details],

$\displaystyle P_{\text{st}}(x)$
$\displaystyle\approx\frac{(\kappa/(2\nu))^{2\alpha}}{B(\alpha,\alpha)}\Bigg{[}\frac{\big{(}X_{-}-x\big{)}^{\alpha-1}}{\big{(}X_{+}+x\big{)}^{-\alpha}}\Theta\Big{(}X_{+}+x\Big{)}\Theta\Big{(}X_{-}-x\Big{)}$
(66)
$\displaystyle+\frac{\big{(}X_{-}+x\big{)}^{\alpha-1}}{\big{(}X_{+}-x\big{)}^{-\alpha}}\Theta\Big{(}X_{+}-x\Big{)}\Theta\Big{(}X_{-}+x\Big{)}\Bigg{]},$
(67)

where $\alpha=1/(\tau_{a}\kappa)$ and
$b=2\bar{v}_{0}\tau_{1}(1-2\tau_{1}/\tau_{2})$ and $\Theta(u)$ indicates the
Heaviside Theta function, and we have defined,

$\displaystyle
X_{\pm}=X_{\text{OD}}\left[1\pm\frac{2\tau_{1}}{\tau_{2}}\right].$ (68)

Figure 18: Overdamped active regime: Plot of the velocity distribution
$P_{\text{st}}(v)$ obtained from numerical simulation of the IRTP (symbols)
along with the prediction obtained from the effective model [see Eq. (16)]
(solid line) for different values of activity $\tau_{a}$. Here we have taken
$m=0.5,\mu=1.0,a_{0}=1.0$ and $\gamma=5$.

The above equation implies that, to this order, the position distribution is
bounded in the region $|x|<X_{+}$. However, as we have shown before, the exact
bound for the position distribution is given by
$X_{\text{OD}}=(X_{+}+X_{-})/2$. Hence, the approximation Eqs. (63)-(64) leads
to $O(\tau_{1})$ deviation from the exact bound of the position distribution.
Nevertheless, we expect Eq. (67) to accurately describe the position
distribution of the IRTP in this regime for $x$ away from $X_{\text{OD}}$. In
fact, Eq. (67) obviously predicts a shape transition for the distribution at
the critical activity,

$\displaystyle\tau_{a}^{*}=\kappa^{-1}=\tau_{2}\Big{[}1-2\frac{\tau_{1}}{\tau_{2}}+O(\tau_{1}^{2})\Big{]}.$
(69)

Hence, in the active regime $\tau_{a}>\tau_{a}^{*}$, the distribution diverges
at $x=\pm X_{-}$, while in the passive regime $\tau_{a}<\tau_{a}^{*}$ the
distribution vanishes at $x=\pm X_{-}$. Clearly, the presence of inertia
decreases the value of the critical activity.

The approximate distribution given by Eq. (67) is compared with the data from
numerical simulations in Fig. 15, for a small value of $\tau_{1}$ which shows
an excellent agreement. The position distribution indeed shows the shape
transition predicted by Eq. (67), from a U-shaped curve for
$\tau_{a}>\tau_{a}^{*}$ to a inverted U-shaped curve for
$\tau_{a}<\tau_{a}^{*}$. We also numerically estimate $\tau_{a}^{*}$ for
different values of $m$, i.e., $\tau_{1}$ and compare it with the prediction
Eq. (69) [see Fig. 17].

###  V.2 Velocity Distribution

In this section, we analytically characterize the position distribution in the
strongly overdamped regime $\tau_{1}\ll\tau_{2}$. It turns out that, unlike
the position distribution, the velocity distribution shows qualitatively
different behaviour in the strongly active, moderately active and passive
regimes.

Figure 19: Overdamped moderately active and passive regime: Plot showing
second velocity transition in the velocity distributions $P_{\text{st}}(v)$ by
increasing activity $\tau_{a}$ beyond $\tau_{2}$. Here we have taken
$m=0.5,\mu=0.5,a_{0}=1.0$ and $\gamma=5$.

To understand the behaviour of the velocity distribution in the strongly
active and the moderately active regimes, we use the same approximation
introduced in Sec. V.1. Taking a time derivative of (63), we get,

$\displaystyle v=-\kappa z+\bar{v}_{0}\sigma+O(\tau_{1}).$ (70)

where, $\kappa$ and $\nu$ are defined in Eq. (65). Keeping the leading order
terms, i.e., ignoring the $O(1)$ terms, we can compute the stationary
distribution of the velocity [see Appendix D for details], quoted in Eq. (16).

Clearly, to this order, the velocity distribution $P_{\text{st}}(v)$ is
bounded in the region $|v|<2\bar{v}_{0}$. Moreover, it predicts a shape
transition in the velocity distribution at the critical activity
$\tau_{a}^{\star}=\kappa^{-1}$ defined in Eq. (69). In the active regime
$\tau_{a}>\tau_{a}^{\star}$, the distribution diverges at $v=0$ algebraically
with exponent $1-1/(\tau_{a}\kappa)$, while in the passive regime
$\tau_{a}<\tau_{a}^{\star}$ the distribution becomes bimodal with a minimum at
$v=0$. Figure 18 compares the approximate velocity distribution given by Eq.
(16) with the data obtained from numerical simulations, which shows an
excellent agreement between the two for a small value of $\tau_{1}$, i.e., in
the strongly overdamped regime.

Figure 20: Overdamped strongly passive regime: Plot showing the velocity
distribution $P_{\text{st}}(v)$ for different values of activity $\tau_{a}$.
The symbols are for simulation results and the solid lines stand for the
corrected Gaussian given by Eq. (53). Here we have used $m=0.5,\mu=1.0$ and
$\gamma=5$.

As the activity is decreased further, the velocity distribution undergoes a
second shape transition from the bimodal distribution to a dome-shaped
distribution as shown in Fig. 19. However, since the approximation Eq. (70)
ignores $O(1)$ terms, it works only in the strongly and moderately active
regimes. We cannot use this approximation to obtain the velocity distribution
in the passive regime, i.e., $\tau_{a}\ll\tau_{1}$. The behaviour of the
velocity distribution in the passive regime can be obtained using the same
argument as before—for $\tau_{a}\ll\tau_{1}$, the Langevin equation (2)
reduces to that of a diffusive particle with an effective temperature
$\tau_{a}a_{0}^{2}/(2\gamma)$. Consequently, the typical velocity distribution
approaches a Gaussian. As already discussed in Sec. IV.2 for the underdamped
scenario, the corrections to the Gaussian can be systematically computed by
considering contributions from higher-order cumulants of the velocity. Figure
20 compares the velocity distribution obtained using numerical simulations,
with the analytical prediction (53) with $N=5$, which shows an excellent
agreement.

##  VI Conclusion

In this paper, we have studied the behaviour of an inertial run-and-tumble
particle of mass $m$ in the presence of a harmonic trap of strength $\mu$. The
active nature of the dynamics, characterized by the tumbling rate
$\tau_{a}^{-1}$ gives rise to a nonequilibrium stationary state which is
strikingly different than the equilibrium state of a passive inertial particle
in a harmonic trap. We show that, the interplay of the three inherent time-
scales present in the system, namely, the viscous time-scale
$\tau_{1}=m/\gamma$, trap time-scale $\tau_{2}=\gamma/\mu$ and the activity
time-scale $\tau_{a}$, gives rise to a wide range of novel features. In
particular, we find that, two qualitatively different scenarios emerge
depending on whether $\tau_{2}>4\tau_{1}$ or not, i.e., whether the
corresponding passive dynamics is overdamped or underdamped. We characterize
the behaviour of the IRTP by investigating the position and velocity
distributions in these two scenarios.

The marginal distributions in the underdamped scenario have certain unique
features which are very distinct from that of ordinary RTP. Firstly, the
distributions have a mass-dependent finite support which we compute
explicitly. Secondly, both the position and the velocity distributions undergo
shape transition from a multi-peaked structure in the strongly active regime
($\tau_{1}\ll\tau_{a}$) to a single-peaked one in the strongly passive regime
($\tau_{a}\ll\tau_{1}$). We characterize the multi-peaked distribution
analytically using a simplified approximate dynamics. On the other hand, the
marginal distributions in the passive regime result from fluctuations that are
typically Gaussian.

In the overdamped case, the presence of finite inertia also leads to various
interesting features. The position distribution, in this case, is
qualitatively similar to that of ordinary RTP, which has a mass-independent
bound and undergoes a transition from a U-shaped distribution to a dome-shaped
one as the activity is decreased. However, now the critical activity depends
on the mass of the particle. The velocity distribution, which has a mass-
dependent bound also undergoes a shape transition—from a distribution
exhibiting algebraic divergence at $v=0$ in the strongly active regime to a
bimodal distribution with minimum at $v=0$ in the moderately active regime.
Interestingly, the effect of inertia is manifest by the presence of a second
shape transition while moving from the moderately active regime to the
strongly passive regime where the distribution eventually converges to a
Gaussian. We analytically characterize these distributions in the strongly
overdamped limit $\tau_{1}\ll\tau_{2}$.

Here we have investigated the marginal distributions of the IRTP in the
presence of well-separated time scales. However, it remains an interesting
open question to find the distributions when the time-scales are comparable.
It would also be interesting to investigate the behaviour of an IRTP in higher
dimensions, which is particularly relevant experimentally. Moreover, it would
be intriguing to study the motion of an IRTP in the presence of a magnetic
field to investigate questions like how the non-Gaussian nature of the active
noise changes the behaviour compared to inertial AOUP Muhsin and Sahoo (2022)
and passive Brownian particles Jiménez-Aquino and Romero-Bastida (2006).
Finally, IRTP provides a simple set-up to investigate how inertia affects the
collective behaviour of active particles.

##  Appendix A Formal Solution of the Langevin Equation

In this Appendix, we provide the formal solution of the Langevin Eq. (2). To
solve these equations, which is a set of first-order linear coupled
inhomogeneous equations, we first recast it in the form of a matrix equation,

$\displaystyle\frac{dX}{dt}=MX+F,$ (71)

where,

$\displaystyle X=\left[{x\atop v}\right],\leavevmode\nobreak\
M=\left[{0\atop-\frac{\mu}{m}}\leavevmode\nobreak\
{1\atop-\frac{\gamma}{m}}\right],\leavevmode\nobreak\
F(t)=\left[{0\atop\frac{a_{0}}{m}\sigma(t)}\right].$ (72)

The formal solution to Eq. (71) is given by,

$\displaystyle X(t)=\Phi(t)X_{0}+\Phi(t)\int_{0}^{t}ds\leavevmode\nobreak\
\Phi^{-1}(s)F(s),$ (73)

where, $\Phi(t)=e^{Mt}$ and $X_{0}=[x_{0}\leavevmode\nobreak\
v_{0}]^{\text{T}}$ denotes the initial condition on position and velocity. It
is straightforward to compute $\Phi(t)$ explicitly, and it is given by,

$\displaystyle\Phi(t)=e^{-t/(2\tau_{1})}\begin{pmatrix}\displaystyle{\cosh\frac{\omega
t}{2\tau_{1}}\atop+\frac{1}{\omega}\sinh\frac{\omega
t}{2\tau_{1}}}&\frac{2\tau_{1}}{\omega}\sinh\frac{\omega
t}{2\tau_{1}}\\\\[15.00002pt] -\frac{2}{\omega\tau_{2}}\sinh\frac{\omega
t}{2\tau_{1}}&\displaystyle{\cosh\frac{\omega
t}{2\tau_{1}}\atop-\frac{1}{\omega}\sinh\frac{\omega
t}{2\tau_{1}}}\end{pmatrix},$ (74)

where, $\omega=\sqrt{1-4\tau_{1}/\tau_{2}}$. Moreover,

$\displaystyle\Phi^{-1}(s)F(s)=\frac{a_{0}\sigma(s)}{m}e^{s/(2\tau_{1})}\left[\begin{array}[]{c}-\frac{2\tau_{1}}{\omega}\sinh\frac{\omega
s}{2\tau_{1}}\\\ \cosh\frac{\omega
s}{2\tau_{1}}+\frac{1}{\omega}\sinh\frac{\omega
s}{2\tau_{1}}\end{array}\right].$ (77)

Substituting the above matrix in Eq. (73), we get the general solution,

$\displaystyle x(t)$
$\displaystyle=e^{-\frac{t}{2\tau_{1}}}\left[x_{0}(\cosh\frac{\omega
t}{2\tau_{1}}+\frac{1}{\omega}\sinh\frac{\omega
t}{2\tau_{1}})+\frac{2\tau_{1}v_{0}}{\omega}\sinh\frac{\omega
t}{2\tau_{1}}\right]+\frac{2\bar{v}_{0}}{\omega}\int_{0}^{t}ds\leavevmode\nobreak\
\sigma(s)e^{-\frac{1}{2\tau_{1}}(t-s)}\sinh\left[\frac{\omega}{2\tau_{1}}(t-s)\right],$
(78) $\displaystyle v(t)$
$\displaystyle=e^{-\frac{t}{2\tau_{1}}}\left[-\frac{2x_{0}}{\tau_{2}\omega}\sinh\frac{\omega
t}{2\tau_{1}}+v_{0}(\cosh\frac{\omega
t}{2\tau_{1}}-\frac{1}{\omega}\sinh\frac{\omega t}{2\tau_{1}})\right]$ (79)
$\displaystyle\hskip
120.0pt+\frac{\bar{v}_{0}}{\tau_{1}}\int_{0}^{t}ds\leavevmode\nobreak\
\sigma(s)e^{-\frac{1}{2\tau_{1}}(t-s)}\Bigg{[}\cosh\Big{[}\frac{\omega}{2\tau_{1}}(t-s)\Big{]}-\frac{1}{\omega}\sinh\Big{[}\frac{\omega}{2\tau_{1}}(t-s)\Big{]}\Bigg{]},$
(80)

where $\bar{v}_{0}=a_{0}/\gamma$ is defined in (3). For the overdamped case,
$4\tau_{1}\leq\tau_{2}$, and using initial conditions, $x=0$ and $v=0$ at
$t=0$, we recover Eqs. (4)-(6) by identifying $\omega=\lambda$. On the other
hand, for the underdamped case, $4\tau_{1}>\tau_{2}$, we substitute
$\omega=i\Lambda$ in Eqs. (78)-(80), which gives Eqs. (8)-(10) using the same
set initial conditions, $x=0$ and $v=0$ at $t=0$.

##  Appendix B Computation of moments from Langevin equation

In this section, we explicitly compute the exact time-dependent second-order
moments and correlations of position and velocity, using the Langevin Eq. (2).
To this end, we require the autocorrelation function for the dichotomous noise
$\sigma(t)$ Hänggi and Jung (1994), which reads,

$\displaystyle\langle\sigma(t_{1})\sigma(t_{2})\rangle=e^{-2|t_{1}-t_{2}|/\tau_{a}}$
(81)

Using the initial condition $x_{0}=v_{0}=0$ for simplicity, the correlations
can then be explicitly computed using Eqs. (78) and 80). This yields, for the
overdamped case,

$\displaystyle\langle x^{2}(t)\rangle$
$\displaystyle=\frac{\bar{v}_{0}^{2}\tau_{2}}{\lambda^{2}\Big{[}\big{(}\frac{4\tau_{1}}{\tau_{a}^{2}}+\frac{1}{\tau_{2}}\big{)}^{2}-\frac{4}{\tau_{a}^{2}}\Big{]}}\Bigg{[}e^{-t/\tau_{1}}\Big{[}\frac{2}{\tau_{a}}\big{(}\frac{2\tau_{1}}{\tau_{a}}+1\big{)}+\frac{1}{\tau_{2}}\Big{]}\Big{[}\frac{2\tau_{1}}{\tau_{2}}\big{(}\frac{4\tau_{1}}{\tau_{a}}-1\big{)}+\big{(}1-\frac{2\tau_{1}}{\tau_{a}}-\frac{2\tau_{1}}{\tau_{2}}\big{)}\cosh\frac{\lambda
t}{\tau_{1}}-\big{(}\frac{2\tau_{1}}{\tau_{a}}-1\big{)}$ (82)
$\displaystyle\lambda\sinh\frac{\lambda
t}{\tau_{1}}\Big{]}+\lambda^{2}\Big{(}\frac{2\tau_{1}}{\tau_{a}}+1\Big{)}\Big{(}\frac{4\tau_{1}}{\tau_{a}^{2}}-\frac{2}{\tau_{a}}+\frac{1}{\tau_{2}}\Big{)}-\frac{2\lambda}{\tau_{2}}e^{-t(4\tau_{1}/\tau_{a}+1)/(2\tau_{1})}\Big{[}\lambda\cosh\frac{\lambda
t}{2\tau_{1}}+\big{(}\frac{4\tau_{1}}{\tau_{a}}+1\big{)}\sinh\frac{\lambda
t}{2\tau_{1}}\Big{]}\Bigg{]},$ (83) $\displaystyle\langle v^{2}(t)\rangle$
$\displaystyle=\frac{8\bar{v}_{0}^{2}\tau_{1}/\tau_{2}}{\lambda^{2}\Big{[}\lambda^{2}-\big{(}\frac{4\tau_{1}}{\tau_{a}}+1\big{)}^{2}\Big{]}\Big{[}\lambda^{2}-\big{(}\frac{4\tau_{1}}{\tau_{a}}-1\big{)}^{2}\Big{]}}\Bigg{[}16\lambda\frac{\tau_{1}}{\tau_{a}}e^{-t(4\tau_{1}/\tau_{a}+1)/(2\tau_{1})}\Big{[}\frac{\tau_{2}}{\tau_{a}}\lambda\cosh\frac{\lambda
t}{2\tau_{1}}-\big{(}1+\frac{\tau_{2}}{\tau_{a}}\big{)}\sinh\frac{\lambda
t}{2\tau_{1}}\Big{]}$ (84)
$\displaystyle-\frac{\tau_{2}}{\tau_{a}}\lambda^{2}\Big{[}\lambda^{2}-\big{(}\frac{4\tau_{1}}{\tau_{2}}-1\big{)}^{2}\Big{]}-e^{-t/\tau_{1}}\Big{[}\lambda^{2}-\big{(}1+\frac{4\tau_{1}}{\tau_{a}}\big{)}^{2}\Big{]}\Big{[}\frac{4\tau_{1}}{\tau_{a}}-1+\big{(}1-\frac{\tau_{2}}{\tau_{a}}\big{)}\cosh\frac{\lambda
t}{\tau_{1}}+\frac{\tau_{2}}{\tau_{a}}\lambda\sinh\frac{\lambda
t}{\tau_{1}}\Big{]}\Bigg{]},$ (85) $\displaystyle\langle x(t)v(t)\rangle$
$\displaystyle=\frac{8\bar{v}_{0}^{2}\tau_{1}\sinh\frac{\lambda
t}{2\tau_{1}}}{\lambda^{2}\big{[}\lambda^{2}-(\frac{4\tau_{1}}{\tau_{a}}-1)^{2}\big{]}}\Bigg{[}e^{-t/\tau_{1}}\Big{[}\lambda\cosh\frac{\lambda
t}{2\tau_{1}}-\big{(}\frac{4\tau_{1}}{\tau_{a}}-1\big{)}\sinh\frac{\lambda
t}{2\tau_{1}}\Big{]}-\lambda
e^{-t(4\tau_{1}/\tau_{a}+1)/(2\tau_{1})}\Bigg{]},$ (86)

where $\lambda$ is defined in Eq. (7). Similarly, in the underdamped case, the
correlations take the following form,

$\displaystyle\langle x^{2}(t)\rangle$
$\displaystyle=\frac{\bar{v}_{0}^{2}\tau_{2}}{\Lambda^{2}\Big{[}\big{(}\frac{4\tau_{1}}{\tau_{a}^{2}}+\frac{1}{\tau_{2}}\big{)}^{2}-\frac{4}{\tau_{a}^{2}}\Big{]}}\Bigg{[}-e^{-t/\tau_{1}}\Big{[}\frac{2}{\tau_{a}}\big{(}\frac{2\tau_{1}}{\tau_{a}}+1\big{)}+\frac{1}{\tau_{2}}\Big{]}\Big{[}\frac{2\tau_{1}}{\tau_{2}}\big{(}\frac{4\tau_{1}}{\tau_{a}}-1\big{)}+\big{(}1-\frac{2\tau_{1}}{\tau_{a}}-\frac{2\tau_{1}}{\tau_{2}}\big{)}\cos\frac{\Lambda
t}{\tau_{1}}+\big{(}\frac{2\tau_{1}}{\tau_{a}}-1\big{)}$ (87)
$\displaystyle\Lambda\sin\frac{\Lambda
t}{\tau_{1}}\Big{]}+\Lambda^{2}\Big{(}\frac{2\tau_{1}}{\tau_{a}}+1\Big{)}\Big{(}\frac{4\tau_{1}}{\tau_{a}^{2}}-\frac{2}{\tau_{a}}+\frac{1}{\tau_{2}}\Big{)}-\frac{2\Lambda}{\tau_{2}}e^{-t(4\tau_{1}/\tau_{a}+1)/(2\tau_{1})}\Big{[}\Lambda\cos\frac{\Lambda
t}{2\tau_{1}}+\big{(}\frac{4\tau_{1}}{\tau_{a}}+1\big{)}\sin\frac{\Lambda
t}{2\tau_{1}}\Big{]}\Bigg{]},$ (88) $\displaystyle\langle v^{2}(t)\rangle$
$\displaystyle=\frac{8\bar{v}_{0}^{2}\tau_{1}/\tau_{2}}{\Lambda^{2}\Big{[}\Lambda^{2}+\big{(}\frac{4\tau_{1}}{\tau_{a}}+1\big{)}^{2}\Big{]}\Big{[}\Lambda^{2}+\big{(}\frac{4\tau_{1}}{\tau_{a}}-1\big{)}^{2}\Big{]}}\Bigg{[}16\Lambda\frac{\tau_{1}}{\tau_{a}}e^{-t(4\tau_{1}/\tau_{a}+1)/(2\tau_{1})}\Big{[}\frac{\tau_{2}}{\tau_{a}}\Lambda\cos\frac{\Lambda
t}{2\tau_{1}}-\big{(}1+\frac{\tau_{2}}{\tau_{a}}\big{)}\sin\frac{\Lambda
t}{2\tau_{1}}\Big{]}$ (89)
$\displaystyle+\frac{\tau_{2}}{\tau_{a}}\Lambda^{2}\Big{[}\Lambda^{2}+\big{(}\frac{4\tau_{1}}{\tau_{2}}-1\big{)}^{2}\Big{]}-e^{-t/\tau_{1}}\Big{[}\Lambda^{2}+\big{(}1+\frac{4\tau_{1}}{\tau_{a}}\big{)}^{2}\Big{]}\Big{[}\frac{4\tau_{1}}{\tau_{a}}-1+\big{(}1-\frac{\tau_{2}}{\tau_{a}}\big{)}\cos\frac{\Lambda
t}{\tau_{1}}-\frac{\tau_{2}}{\tau_{a}}\Lambda\sin\frac{\Lambda
t}{\tau_{1}}\Big{]}\Bigg{]},$ (90) $\displaystyle\langle x(t)v(t)\rangle$
$\displaystyle=\frac{8\bar{v}_{0}^{2}\tau_{1}\sin\frac{\Lambda
t}{2\tau_{1}}}{\Lambda^{2}\big{[}\Lambda^{2}+(\frac{4\tau_{1}}{\tau_{a}}-1)^{2}\big{]}}\Bigg{[}e^{-t/\tau_{1}}\Big{[}\big{(}\frac{4\tau_{1}}{\tau_{a}}-1\big{)}\sin\frac{\Lambda
t}{2\tau_{1}}-\Lambda\cos\frac{\Lambda t}{2\tau_{1}}\Big{]}+\Lambda
e^{-t(4\tau_{1}/\tau_{a}+1)/(2\tau_{1})}\Bigg{]}.$ (91)

where $\Lambda$ is defined in Eq. (11).

##  Appendix C Calcuation of stationary state distribution bound

In this section, we explicitly compute the bounds of the stationary state
distribution. As discussed in Sec. IV, in the underdamped case, this bound is
obtained from the trajectory which maximizes the displacement of the particle
— when tumbling events $\sigma\to-\sigma$ coincide with the velocity $v$
changing its sign. To construct such a trajectory, we first compute the
position of a particle at the turning points of its underdamped motion. To
this end, we evaluate the integrals in Eqs. (78) and (80) with the particle
starting at rest with some initial position $x_{0}$ and for a fixed propulsion
direction $\sigma$,

$\displaystyle x(t)$
$\displaystyle=\frac{a_{0}}{\mu}\sigma+(x_{0}-\frac{a_{0}}{\mu}\sigma)e^{-t/(2\tau_{1})}\Big{[}\cos\frac{\Lambda
t}{2\tau_{1}}+\frac{1}{\Lambda}\sin\frac{\Lambda t}{2\tau_{1}}\Big{]},$ (92)
$\displaystyle v(t)$
$\displaystyle=\frac{2e^{-t/(2\tau_{1})}}{\Lambda\tau_{2}}\big{(}\frac{a_{0}}{\mu}\sigma-
x_{0}\big{)}\sin\frac{\Lambda t}{2\tau_{1}},$ (93)

where $\Lambda$ is defined in Eq. (11). Then the time instants of the turning
points are computed by solving for the roots of the equation $v(t)=0$,

Figure 21: Underdamped scenario: Typical displacement maximizing trajectory of
IRTP with initial position ($x_{0}=0$) and velocity ($v_{0}=1$) chosen outside
the supported region in the $x-v$ plane. The particle eventually settles into
a limit cycle (black dashed line) described by the trajectories Eq. (44)-(43),
which gives the bound of joint stationary distribution $P_{\text{st}}(x,v)$.
Here we have used $\gamma=0.5,\mu=1.0,a_{0}=1.0$ and $m=1.5$.

$\displaystyle t_{n}=\frac{2\pi n\tau_{1}}{\Lambda},\quad n=1,2,3,\ldots$ (94)

Consequently, the position of the $n$-th turning point is given by
substituting $t_{n}$ in Eq. (92),

$\displaystyle
x_{n}=\frac{a_{0}}{\mu}\sigma+(-1)^{n}e^{-n\pi/\Lambda}(x_{0}-\frac{a_{0}}{\mu}\sigma).$
(95)

By definition, the tumbling event occurs at the first turning point $x_{1}$.
Consequently, we have the underdamped particle now starting at rest with
initial position $x_{1}$ and with self-propulsion direction $-\sigma$. The
next tumbling event again occurs at the first turning point of the next
trajectory, and so on — iteratively generating the full trajectory.

Let $x_{\ell}$ denote the particle position at the $\ell$-th tumbling event
$\sigma_{\ell}=-\sigma_{\ell-1}$ in the trajectory. Considering the iterative
nature of the trajectory, we obtain a recursive relation for $x_{\ell}$, by
substituting $x_{\ell-1}$ as the initial position and $\sigma_{\ell-1}$ as the
propulsion direction in Eq. (95),

$\displaystyle x_{\ell}=cx_{\ell-1}+\frac{a_{0}}{\mu}(1-c)\sigma_{\ell-1},$
(96)

where $c=-e^{-\pi/\Lambda}$ and $\Lambda$ is defined in Eq. (11). Depending on
the sign of $|x_{\ell}|-|x_{\ell-1}|$, the particle distance from the origin
can either increase or decrease at successive tumbling events. The two
scenarios for $|x_{0}|<a_{0}/\mu$ and $|x_{0}|>a_{0}/\mu$ are illustrated in
Figs. 6 and 21 respectively. Consequently, the bound can be computed from the
asymptotic value of $x_{\ell}$ in the limit $\ell\to\infty$. Assuming that the
particle starts at rest with some initial position $x_{0}$, and propulsion
direction $\sigma_{0}$, the position of the $\ell$-th turning point can be
computed iteratively, using Eq. (96), as quoted in Eq. (41). Then, the maximum
distance that the particle can reach is given by the limiting value
$\lim_{\ell\to\infty}x_{\ell}$, which we identify as the bound of the
stationary position distribution quoted in (14),

Alternately, the bound can also be obtained using the condition,
$|x_{l}|=|x_{l-1}|$, which gives the limit cycle describing the boundary of
the joint distribution $P_{\text{st}}(x,v)$. Then by setting
$x_{\ell}=-x_{\ell-1}$ in Eq. (96), we obtain the bound for the position
distribution explicitly.

Similarly, the bound for the velocity distribution is given by the maximum
attained velocity in this limit cycle, i.e. when the particle acceleration
changes sign for the first time. To locate this velocity turning point, we
first compute the time instants of these turning points by solving for the
roots of the equation $dv/dt=0$, using Eq. (93), which gives,

$\displaystyle t_{n}=\frac{2\tau_{1}}{\Lambda}(\tan^{-1}\Lambda+n\pi),\quad
n=0,1,2,\ldots.$ (97)

Then, the $n$-th turning point in the velocity is computed by evaluating the
velocity $v(t_{n})$ at time $t_{n}$

$\displaystyle
v_{n}=\frac{(-1)^{n}}{\sqrt{m/\mu}}\Big{(}\frac{a_{0}}{\mu}\sigma-
x_{0}\Big{)}e^{-\frac{1}{\Lambda}\big{[}n\pi+\tan^{-1}\Lambda\big{]}}$ (98)

It is then straightforward to calculate the velocity bound from the first
turning point $(n=0)$ for the particle starting from an initial position
$x_{0}=-\sigma X_{\text{UD}}$, which gives us the expression quoted in Eq.
(15).

Lastly, we note that these turning points also allow us to compute the
position of the peaks in the position and velocity distributions. By
substituting $x_{0}=\pm a_{0}/\mu$ in Eqs. (95) and (95), we arrive at the
expressions quoted in Eqs. (45) and (50) respectively.

We now turn to the overdamped scenario. As already discussed in Sec. V, in the
overdamped case, the displacement of the particle is maximum when no tumbling
event takes place. This is easily seen by examining the integrand in the
formal solution for $x(t)$, quoted in Eq. (78). Apart from $\sigma(s)$, the
integrand is a purely positive quantity, i.e., in a given time interval the
integral evaluates to the maximum value when no tumbling event takes place.
Performing the integral in Eqs. (78) and (80) with particle starting at with
position $x_{0}$ and with a fixed propulsion direction $\sigma$, we get for
the overdamped case,

$\displaystyle x(t)$
$\displaystyle=\frac{a_{0}}{\mu}\sigma+(x_{0}-\frac{a_{0}}{\mu}\sigma)e^{-t/(2\tau_{1})}\Big{[}\cosh\frac{\lambda
t}{2\tau_{1}}+\frac{1}{\lambda}\sin\frac{\lambda t}{2\tau_{1}}\Big{]}.$ (99)
$\displaystyle v(t)$
$\displaystyle=\frac{2e^{-t/(2\tau_{1})}}{\lambda\tau_{2}}\big{(}\frac{a_{0}}{\mu}\sigma-
x_{0}\big{)}\sinh\frac{\lambda t}{2\tau_{1}},$ (100)

where $\lambda$ is defined Eq. (7). Taking the limit, $t\to\infty$ in Eq.
(99), gives the maximum displacement of the particle for a given propulsion
direction, which we identify as the bound of the position distribution,
$X_{\text{OD}}$, in the overdamped case, quoted in Eq. (12). Consequently, the
joint distribution must be bounded by the trajectories connecting the point
$(\pm X_{\text{OD}},0)$, allowing us to compute the bound on the velocity
distribution.

In particular, the velocity extremum between the points $(\pm
X_{\text{OD}},0)$ marks the velocity bound in this case. We compute this bound
explicitly by first solving for the equation $dv/dt=0$, using Eq. (100), which
gives,

$\displaystyle
t=\frac{2\tau_{1}}{\lambda}\cosh^{-1}\frac{\gamma}{2\sqrt{m\mu}},$ (101)

and by computing the velocity at this time instant with initial position
$x_{0}=-\sigma X_{\text{OD}}$ in Eq. (100), which gives the expression quoted
in Eq. (13).

##  Appendix D Effective models in the strongly overdamped limit

In this Appendix, we explicitly compute the marginal position distribution in
the strongly overdamped limit. We start from Eq. (61), the Langevin equations
for the variables $x(t)=x_{1}(t)-x_{2}(t)$ and $z(t)=x_{1}(t)+x_{2}(t)$,

$\displaystyle\dot{x}(t)=$ $\displaystyle-k_{+}x(t)+k_{-}z(t),$ (102)
$\displaystyle\dot{z}(t)=$ $\displaystyle-
k_{+}z(t)+k_{-}x(t)+2\bar{a}_{0}\sigma(t).$ (103)

where, $k_{\pm}$ are defined in Eq. (62). In the strongly overdamped limit,
$\tau_{1}\ll\tau_{2}$, we treat $\tau_{1}$ as the small parameter and expand
$x(t)$ in orders of $\tau_{1}$ using Eq. (102), to get,

$\displaystyle
x=z-2\Big{(}\dot{x}+\frac{z}{\tau_{2}}\Big{)}\tau_{1}-2\Big{(}\frac{\tau_{1}}{\tau_{2}}\Big{)}^{2}\dot{z}+O(\tau_{1}^{3}).$
(104)

By taking time derivative of the above Eq. (104), we get,

$\displaystyle\dot{x}=\dot{z}-2\Big{(}\ddot{x}+\frac{z}{\tau_{2}}\Big{)}\tau_{1}-2\Big{(}\frac{\tau_{1}}{\tau_{2}}\Big{)}^{2}\ddot{z}+O(\tau_{1}^{3}).$
(105)

Substituting Eq. (105) in Eq. (104) and keeping up to the leading order term
in $\tau_{1}$, we get,

$\displaystyle
x\approx\Big{(}1-2\frac{\tau_{1}}{\tau_{2}}\Big{)}z-2\tau_{1}\dot{z}.$ (106)

Substituting $x(t)$ Eq. (106) in (103) we get an effective equation for
$z(t)$, quoted in Eq. (64). Clearly, $z(t)$ is a fully overdamped RTP moving
in a harmonic trap of stiffness $\kappa$ with propulsion speed $\nu$. The
joint distribution $P_{\sigma}(z)$ is known explicitly Dhar  _et al._ (2019),
and is given by,

$\displaystyle
P_{\sigma}(z)=\frac{4^{-\alpha}\kappa/\nu}{B(\alpha,\alpha)}\frac{[1-(\frac{\kappa
z}{\nu})^{2}]^{\alpha}}{\Big{(}1-\sigma\frac{\kappa z}{\nu}\Big{)}},$ (107)

where $\alpha=1/(\tau_{a}\kappa)$, and is bounded in the region
$-\nu/\kappa<z<\nu/\kappa$. The distribution of $x$ can be obtained using
(106), (64) and (107). Substituting Eq. (64) in Eq. (106) and keeping upto the
leading order term, we get,

$\displaystyle x(t)=az(t)-b\sigma(t),$ (108)

with,

$\displaystyle
a=\frac{\lambda\tau_{2}}{\tau_{2}-2\tau_{1}}\leavevmode\nobreak\
\leavevmode\nobreak\ \text{and}\leavevmode\nobreak\ \leavevmode\nobreak\
b=2\bar{v}_{0}\tau_{1}a.$ (109)

Thus the stationary distribution of $x$ is given by,

$\displaystyle
P_{\text{st}}(x)=\sum_{\sigma}\intop_{-\nu/\kappa}^{\nu/\kappa}dz\leavevmode\nobreak\
\delta(x-az+b\sigma)P_{\sigma}(z).$ (110)

where, $P_{\sigma}(z)$ is defined in Eq. (107) and $\delta(u)$ is the delta
function. Performing the integral over $z$ we get,

$\displaystyle P_{\text{st}}(x)=\frac{1}{a}\sum_{\sigma=\pm
1}P_{\sigma}\Big{(}\frac{x}{a}+\frac{b}{a}\sigma\Big{)}\Theta\Big{(}\frac{\nu}{\kappa}a+b\sigma+x\Big{)}$
(111) $\displaystyle\Theta\Big{(}\frac{\nu}{\kappa}a-b\sigma-x\Big{)}.$ (112)

To order $O(\tau_{1})$, $a=1$ and we arrive at the expression quoted in Eq.
(67).

By definition, the velocity $v$ is given by the time derivative of Eq. (106).
To the leading order in $\tau$, the Langevin equation for the velocity $v$ is
then given by,

$\displaystyle v=\dot{z}+O(\tau_{1})$ (113)

Consequently, the stationary velocity distribution is given by,

$\displaystyle P_{\text{st}}(v)=\sum_{\sigma=\pm
1}\intop_{-\kappa/\nu}^{\kappa/\nu}dz\leavevmode\nobreak\ \delta(v+\kappa
z-\nu\sigma)P_{\sigma}(z).$ (114)

where, $P_{\sigma}(z)$ is defined in Eq. (107). Performing the integral over
$z$ we obtain Eq. (16).

## References

  * Cavagna and Giardina (2014) A. Cavagna and I. Giardina, Annu. Rev. Condens. Matter Phys. 5, 183 (2014).
  * Feinerman  _et al._ (2018) O. Feinerman, I. Pinkoviezky, A. Gelblum, E. Fonio, and N. S. Gov, Nat. Phys. 14, 683 (2018).
  * Pavlov and Kasumyan (2000) D. S. Pavlov and A. O. Kasumyan, J. Ichthyology 40, S163 (2000).
  * Mukundarajan  _et al._ (2016) H. Mukundarajan, T. C. Bardon, D. H. Kim, and M. Prakash, J. Exp. Biol. 219, 752 (2016).
  * Jiang and Granick (2012) S. Jiang and S. Granick, eds.,  _Janus Particle Synthesis, Self-assembly and Applications_ (Royal Society of Chemistry, 2012).
  * Bechinger  _et al._ (2016) C. Bechinger, R. Di Leonardo, H. Löwen, C. Reichhardt, G. Volpe, and G. Volpe, Rev. Mod. Phys. 88, 045006 (2016).
  * Walther and Muller (2013) A. Walther and A. H. Muller, Chem. Rev. 113, 5194 (2013).
  * Nelson  _et al._ (2010) B. J. Nelson, I. K. Kaliakatsos, and J. J. Abbott, Annu. Rev. Biomed. 12, 55 (2010).
  * Basu  _et al._ (2018) U. Basu, S. N. Majumdar, A. Rosso, and G. Schehr, Phys. Rev. E 98, 062121 (2018).
  * Howse  _et al._ (2007) J. R. Howse, R. A. Jones, A. J. Ryan, T. Gough, R. Vafabakhsh, and R. Golestanian, Phys. Rev. Lett. 99, 048102 (2007).
  * Malakar  _et al._ (2018) K. Malakar, V. Jemseena, A. Kundu, K. V. Kumar, S. Sabhapandit, S. N. Majumdar, S. Redner, and A. Dhar, J. Stat. Mech.: Theory Exp. 2018, 043215 (2018).
  * Tailleur and Cates (2008) J. Tailleur and M. E. Cates, Phys. Rev. Lett. 100, 218103 (2008).
  * Bonilla (2019) L. L. Bonilla, Phys. Rev. E 100, 022601 (2019).
  * Martin  _et al._ (2021) D. Martin, J. O’Byrne, M. E. Cates, É. Fodor, C. Nardini, J. Tailleur, and F. Van Wijland, Phys. Rev. E 103, 032607 (2021).
  * Mandal  _et al._ (2019) S. Mandal, B. Liebchen, and H. Löwen, Phys. Rev. Lett. 123, 228001 (2019).
  * Su  _et al._ (2021) J. Su, H. Jiang, and Z. Hou, New J. Phys. 23, 013005 (2021).
  * De Karmakar and Ganesh (2020) S. De Karmakar and R. Ganesh, Phys. Rev. E 101, 032121 (2020).
  * Liao  _et al._ (2021) J.-j. Liao, F.-j. Lin, and B.-q. Ai, Phys. A: Stat. Mech. Appl. 582, 126251 (2021).
  * Negro  _et al._ (2022) G. Negro, C. B. Caporusso, P. Digregorio, G. Gonnella, A. Lamura, and A. Suma, Eur. Phys. J. E 45, 75 (2022).
  * Marconi  _et al._ (2021) U. M. B. Marconi, L. Caprini, and A. Puglisi, New J. Phys. 23, 103024 (2021).
  * Caprini  _et al._ (2021) L. Caprini, C. Maggi, and U. Marini Bettolo Marconi, J. Chem. Phys. 154 (2021).
  * Caprini and Marconi (2021) L. Caprini and U. M. B. Marconi, Soft Matter 17, 4109 (2021).
  * Khali  _et al._ (2024) S. S. Khali, F. Peruani, and D. Chaudhuri, Phys. Rev. E 109, 024120 (2024).
  * Patel and Chaudhuri (2023) M. Patel and D. Chaudhuri, New J. Phys. 25, 123048 (2023).
  * (25) M. Patel and D. Chaudhuri, arXiv:2404.01107 .
  * Sprenger  _et al._ (2023) A. R. Sprenger, L. Caprini, H. Löwen, and R. Wittmann, J. Phys. Condens. Matter 35, 305101 (2023).
  * Scholz  _et al._ (2018) C. Scholz, S. Jahanshahi, A. Ldov, and H. Löwen, Nat. Commun. 9, 5156 (2018).
  * Breoni  _et al._ (2020) D. Breoni, M. Schmiedeberg, and H. Löwen, Phys. Rev. E 102, 062604 (2020).
  * Lisin  _et al._ (2022) E. Lisin, O. Vaulina, I. Lisina, and O. Petrov, Phys. Chem. Chem. Phys. 24, 14150 (2022).
  * Smith and Farago (2022) N. R. Smith and O. Farago, Phys. Rev. E 106, 054118 (2022).
  * Caprini and Marini Bettolo Marconi (2021) L. Caprini and U. Marini Bettolo Marconi, J. Chem. Phys. 154 (2021).
  * Muhsin and Sahoo (2022) M. Muhsin and M. Sahoo, Phys. Rev. E 106, 014605 (2022).
  * (33) P. M. Obreque, O. Garrido, D. Romero, H. Löwen, and F. Guzmán-Lastra, arXiv:2403.02569 .
  * (34) A. P. Antonov, L. Caprini, C. Scholz, and H. Löwen, arXiv:2404.06615 .
  * Muhsin and Sahoo (2023) M. Muhsin and M. Sahoo, Phys. Rev. E 107, 054601 (2023).
  * (36) F. Adersh, M. Muhsin, and M. Sahoo, arXiv:2402.05538 .
  * Dhar  _et al._ (2019) A. Dhar, A. Kundu, S. N. Majumdar, S. Sabhapandit, and G. Schehr, Phys. Rev. E 99, 032132 (2019).
  * (38) DLMF, “NIST Digital Library of Mathematical Functions,” https://dlmf.nist.gov/, Release 1.2.0 of 2024-03-15, f. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller, B. V. Saunders, H. S. Cohl, and M. A. McClain, eds.
  * Santra  _et al._ (2021) I. Santra, U. Basu, and S. Sabhapandit, Phys. Rev. E 104, L012601 (2021).
  * Jiménez-Aquino and Romero-Bastida (2006) J. Jiménez-Aquino and M. Romero-Bastida, Rev. Mex. Fis. 52, 182 (2006).
  * Hänggi and Jung (1994) P. Hänggi and P. Jung, Adv. Chem. Phys. 89, 239 (1994).