Exact moments and re-entrant transitions in the inertial dynamics of active Brownian particles

In this study, we investigate the behavior of free inertial active Brownian particles in the presence of thermal noise. While finding a closed-form solution for the joint distribution of positions, orientations, and velocities using the Fokker–Planck equation is generally challenging, we utilize a Laplace transform method to obtain the exact temporal evolution of all dynamical moments in arbitrary dimensions. Our expressions in d dimensions reveal that inertia significantly impacts steady-state kinetic temperature and swim pressure while leaving the late-time diffusivity unchanged. Notably, as a function of activity and inertia, the steady-state velocity distribution exhibits a remarkable re-entrant crossover from ‘passive’ Gaussian to ‘active’ non-Gaussian behaviors. We construct a corresponding ‘phase diagram’ using the exact expression of the d-dimensional kurtosis. Our analytic expressions describe steady states and offer insights into time-dependent crossovers observed in moments of velocity and displacement. Our calculations can be extended to predict up to second-order moments for run-and-tumble particles and the active Ornstein–Uhlenbeck process (AOUP). Additionally, the kurtosis shows differences from AOUP.

In this paper, we consider the motion of non-interacting ABPs in the presence of translational inertia.Using a Laplace transform approach, we show that the Fokker-Planck equation governing the dynamics can be utilized to evaluate the exact time dependence of all moments of dynamical variables in arbitrary d-dimensions.This is the first main contribution of this paper.Offering a unified approach for computing the precise time evolution of all conceivable dynamical variables in arbitrary dimensions, our method represents a noteworthy analytical advancement in investigating inertial effects in active matter.The technique was originally proposed back in 1952 to study equilibrium properties of worm-like chains [53] and was recently put in use to study the dynamics of overdamped ABPs in Ref. [31,[54][55][56].
We present exact expressions of several dynamical moments, which agree with the direct numerical integration results.The system eventually reaches a steady state in velocity but not displacement, with mean-squared displacement undergoing dynamic crossovers before reaching asymptotic diffusion with inertia-independent diffusivity.In contrast, the kinetic temperature and swim pressure depend on inertia; with inertia, while the kinetic temperature increases to saturate eventually, the swim pressure decreases monotonically to vanish in the infinite mass limit.
The most striking result is the occurrence of a re-entrant crossover between active and passive states as inertia varies, quantified through deviations from the equilibriumlike Gaussian velocity distribution.Analytically, this is captured by the d-dimensional kurtosis of velocity, which shows a non-monotonic variation with inertia and a decrease to negative values followed by a saturation with increasing activity.We use the results to obtain a "phase diagram" in the activity-inertia plane.This is the second main contribution of this study.Note that in this single-particle system, there is no genuine phase transition; instead, it represents a crossover between active and passive behaviors.Further, we propose approximate analytic forms of probability distributions of velocity that show reasonable agreement with numerical results at small and large inertia limits and can qualitatively capture the re-entrant active-passive crossovers.Our calculations for ABP can be extended to predict up to second-order moments for RTPs and the AOUPs.Additionally, the non-zero kurtosis for intermediate inertia shows differences from AOUP.
The paper is structured as follows: In Section 2, we introduce the model and Langevin dynamics.The subsequent section outlines a Laplace-transform-based approach to compute dynamic moments in any dimension.Sections 4 and 5 present detailed results related to velocity and displacement moments.Section 4 delves into discussions about kinetic temperature, the "phase diagram", and velocity distributions, while Section 5 explores topics like active diffusivity and swim pressure.Using appropriate mapping, we extend the results to the related RTP and AOUP systems in Section 6.Finally, in Section 7, we provide a conclusion summarizing our key findings.

Model
The underdamped motion of ABPs in d-dimensional space moving with a constant active speed v a in the orientation û = (u 1 , u 2 , ..., u d ) is described by its position evolving with time t ′ .We consider the dynamics in the presence of translational and orientational Brownian noise with diffusivities D and D r , respectively.The time and length scales are set by τ r = D −1 r and ℓ = D/D r , to give the unit of velocity v = √ DD r .We use dimensionless position r = r ′ /ℓ, time t = t ′ /τ r , and velocity v = v ′ /v to express the Langevin equation of motion in d dimensions in the following Ito form [57][58][59] where we used the dimensionless mass M = τ /τ r with inertial relaxation time τ = m/γ, m denoting the translational inertia and γ the viscous damping, and Péclet number Pe = v a / √ DD r .The Gaussian noise in translation and rotation are uncorrelated, and their components obey ⟨dB i ⟩ = 0 and ⟨dB i dB j ⟩ = δ ij dt.The first term in Eq. ( 3) denotes a projection operator for the noise onto a (d − 1)-dimensional surface perpendicular to û.The second term ensures that the unit vector remains normalized after displacement.The equations can be directly integrated numerically using the Euler-Maruyama scheme.

Calculation of moments from the Fokker-Planck equation
Noting that the heading direction û undergoes an orientational diffusion on a (d − 1)dimensional hypersurface and v performs a drift-diffusion, the probability distribution P (r, v, û, t) evolves following the Fokker-Planck equation where ∇ and ∇ v denote the gradient operator in d-dimensional position and velocity space, respectively.The Laplacian ∇ 2 u in (d − 1)-dimensional orientation space can be expressed in terms of Cartesian coordinates y as y∂ y ] by defining u i = y i /y with y = |y|.Using a Laplace transform P (r, v, û, s) = ∞ 0 dt e −st P (r, v, û, t), the Fokker-Planck equation can be expressed as Therefore, the mean of any observable ⟨ψ⟩ s = dr dv d û ψ(r, v, û) P (r, v, û, s) satisfies the moment equation where ⟨ψ⟩ 0 = dr dv d û ψ(r, v, û)P (r, v, û, 0) is given by the initial condition, which without any loss of generality, can be expressed as P (r, v, û, 0) = δ(r)δ(v−v 0 )δ( û− û0 ).Equation ( 5) can be used to determine any dynamical moment of interest in arbitrary d -dimensions as a function of time by performing an inverse Laplace transform.In the following, we denote the steady-state values lim t→∞ ⟨ψ⟩(t) ≡ ⟨ψ⟩ st = lim s→0 s⟨ψ⟩ s .Our theoretical approach, offering a unified method for computing the precise time evolution of all dynamical variables in arbitrary dimensions, differs from the recent works on inertial ABPs [19,35,51,52].Moreover, we present approximate analytic expressions for velocity distributions parallel and perpendicular to the heading direction in Sec.4.9.

Velocity moments, steady states, and "phase transitions"
We first employ the above equation to obtain the time evolution of several velocity moments, including mean, variance, quartic moment, and lag function towards the steady state.The dynamics of inertial passive particles lead to the equilibrium Maxwell-Boltzmann distribution of velocity.Activity changes this behavior [19,41].Here, we calculate up to fourth order cumulant of velocity to identify such deviations and hence determine a "phase diagram" describing departures from the equilibrium Gaussian characteristic.
In the presence of inertia, the orientation of the velocity vector v can be different from the heading direction û.Orientation fluctuation of the heading direction over the inertial relaxation time leads to this difference.In the limit of vanishing inertia, one expects ⟨v ∥ ⟩ = Pe.Let us define instantaneous velocity components parallel and perpendicular to the heading direction, Its inverse Laplace transform yields, Thus ⟨v ∥ ⟩ starts from û0 • v 0 at t = 0 to asymptotically reach the steady-state value Pe (d−1)M +1 (see Fig. 1(a) ).As expected, in the limit of vanishing inertia, ⟨v ∥ ⟩ = Pe, the velocity vector always remains oriented along the heading direction.The larger the inertia, the bigger the deviation from the heading direction (also see Fig. C1 in the Appendix).By symmetry ⟨v ⊥ ⟩ = 0.

Lag function
While analyzing the motion of inertial vibrobots, a lag function was used in Ref. [19], • v(0)⟩, which quantifies the average difference between the projection of velocity to initial heading direction and the projection of heading direction to the initial velocity.In equilibrium, time-reversal symmetry ensures C(v(t), û(t)) = 0. We obtain the following closed-form expression for inertial ABPs, using v(0) = Pe û0 .In the absence of inertia, C(v(t), û(t)) = 0.In the inertial system, C(v(t), û(t)) starts from zero at t = 0 to returns to zero as t → ∞ after reaching a maximum at an intermediate time  9), ( 10), (12), and (A. in d = 2, 3 are shown in Fig. 1(b), with a comparison against numerical simulations performed in two dimensions.A similar non-monotonic variation was observed in vibrobots [19].

Mean squared velocity
Using ψ = v 2 in equation ( 5) and ⟨ψ⟩(0 Using equation( 8), we get The inverse Laplace transform gives We plot this expression in Fig. 1(c) for d = 2, 3 and show a comparison with numerical simulations in 2d.To analyze the observed dynamical crossovers from ∼ t to ∼ t 2 growth before reaching the steady state, we expand ⟨v 2 ⟩(t) around t = 0 using v 0 = 0, The expansion shows a ⟨v 2 ⟩ ∼ t scaling at a short time before a crossover to see Appendix C. In the asymptotic limit of long time t ≫ M 2 and t ≫ M M (d−1)+1 , the mean-squared velocity reaches the steady state value Given that ⟨v⟩ st = 0, the velocity fluctuation ⟨δv 2 ⟩ st = ⟨v 2 ⟩ st .

Effective temperature and violation of equilibrium fluctuation-dissipation relation
The expression of mean squared velocity can be used to obtain the steady-state kinetic temperature of the system Remarkably, unlike in equilibrium, the kinetic temperature is a function of inertia.In the over-damped limit k B T kin = 1.The kinetic temperature increases with M linearly at a small M to saturate to 1 + Pe 2 /[d(d − 1)] at large M .The kinetic temperature in the coexistence of phase separating ABPs [40] and in active Ornstein-Uhlenbeck particles in two dimensions [60] were calculated before.With a careful mapping of the active force and persistence time of ABPs to that of Ornstein-Uhlenbeck particles, the above prediction of kinetic temperature agrees with that of Ref. [60] in d = 2 (see Sec.
Clearly, the violation is maximum in the overdamped limit of M = 0.It decreases linearly with M at small M , to vanish asymptotically as M → ∞.Thus, with increasing M , ABPs can return to equilibrium-like behavior.

Fluctuation of velocity component in the heading direction
As we have shown before, the steady state value of ⟨v ∥ ⟩ st decreases from Pe to 0 monotonically, with increasing mass (Fig. C2(a) in Appendix C).Here we calculate the variance While it is straightforward to calculate the full time-evolution of ⟨v 2 ∥ ⟩(t) performing an inverse Laplace transform of the above expression, we show the steady-state expressions In the limit of vanishing inertia, the above expression gives M ⟨δv 2 ∥ ⟩ st = 1 and matches with the overdamped limit of k B T kin = 1.In the other limit of and ⊥ ⟩ in the small and large M limits.

Fourth moment of velocity
The calculation of ⟨v 4 ⟩ s proceeds as before using equation( 5) and involves the following steps, where ⟨ û • v⟩ s , ⟨v 2 ⟩ s and ⟨( û • v) 2 ⟩ s are already calculated in equations ( 8), ( 11) and ( 18) respectively.Using them, we get the exact form of ⟨v 4 ⟩ s as shown in Appendix A. The inverse Laplace transform gives the time-dependence ⟨v 4 ⟩(t) shown in equation(A.2) of Appendix.Here, we present the steady-state result, The short-time behavior with initial condition v 0 = 0 can be analyzed in terms of a series expansion of ⟨v 4 ⟩(t) around t = 0 yielding, As shown in Fig. 1(d), ⟨v 4 ⟩(t) undergoes a crossover from ∼ t 2 scaling at short times to ∼ t 3 scaling near t I = 4d(d+2) Further, it crosses over to a ∼ t 4 scaling at t II = A 1 /B 1 before saturating to a steady state.Time-dependence of velocity moments before reaching a steady state depends on the initial condition, a point illustrated further in Appendix C.

Kurtosis: deviations from the Gaussian process
The fourth moment ⟨v 4 ⟩ for a general Gaussian process with ⟨v⟩ = 0 in terms of lower order moments can be expressed as [31] µ The departure from Gaussian behavior can be evaluated using the kurtosis At steady-state, it has the form Thus, kurtosis shows a non-monotonic variation with M (Fig. 2(a) ).It vanishes in the two limits of small and large inertia -in the limit of M → 0 as and in the limit of M → ∞ as K ∼ M −1 , indicating the two passive states with asymptotically Gaussian distributions of velocity.For intermediate M values, K becomes negative and reaches a minimum.This suggests a re-entrance transition from passive to active to passive behavior with the increase in inertia.
On the other hand, with activity Pe, kurtosis decreases monotonically to saturate (Fig. 2(b) ).For Pe → 0, it vanishes as K ∼ Pe 4 .In the limit of large activity Pe → ∞, it saturates to The large amplitude of kurtosis characterizes a strongly non-Gaussian velocity distribution denoting an active state.As shown in Fig. 2(d), the velocity distribution in two dimensions (2d) corresponding to the active state obtained from numerical simulations displays a distinct ring shape.

"Phase diagram" describing active-passive crossovers
A heat map of equation (26) in Fig. 2(c) produces a "phase diagram" identifying crossovers between active (dark, large amplitude of K) and passive (light, K ≈ 0) regions in 2d.It clearly shows a re-entrant crossover from passive (Gaussian) to active (non-Gaussian) to passive (Gaussian) with increasing inertia at a constant activity, e.g., near Pe ≈ 100.The "phase diagram" in d = 3 also shows similar behavior (figure not shown).
The re-entrant crossover can be understood using the following heuristic argument.The velocity distribution becomes Gaussian if the persistence of the heading direction is lost over the inertial relaxation time, i.e., τ ≫ τ r , which is equivalent to M ≫ 1.As a result, the active fluctuations can be approximated as additional 'thermal noise.'The "phase diagram" shows the actual boundary is near M = 1.At M < 1, i.e., τ < τ r , persistence dominates; therefore, the velocity distribution shows maximum on a (d − 1)dimensional hypersphere of radius v = Pe, with the spherical shape emerging due to isotropy in the initial orientation.Also, the thermal noise thickens the hypersurface shell by ⟨δv 2 ⟩ ≈ d/M .The v = Pe peak of the distribution gets unrecognizable when the velocity fluctuation gets larger than the ring size, ⟨δv 2 ⟩ ≫ Pe, leading to the criterion M < Pe −2 for a re-entrance to the passive Gaussian regime.This criterion is consistent with equation ( 27).

Probability distributions of velocity in 2d
While finding a closed-form solution for the joint distribution of positions, orientations, and velocities using the Fokker-Planck equation is generally challenging, it can be analyzed in the two limits of small and large M , particularly for the marginal distributions of velocities.The velocity fluctuations are governed by the kinetic temperature k B T kin in Eq. (15).For M ≪ 1, i.e., τ r ≫ τ , the heading direction can be assumed to be fixed at some angle θ.This regime is equivalent to the dynamics of a particle under external force.Moreover, k B T kin ≈ 1, as M ≪ 1.Thus, the steady-state solution is The velocity distribution P (v x , v y ) can be obtained by integrating over all possible θ, P (v x , v y ) = (1/2π) 2π 0 dθP (v x , v y , θ).The marginal distribution, by symmetry has the same form as As shown in Fig. 3(a), this distribution already describes the small inertia crossover between the passive Gaussian to active non-Gaussian (bimodal symmetric) nature, capturing the numerical simulation results at M = 10 −4 , 10 −2 for Pe = 100.Thus, the previously proposed perturbative calculation [35] is not necessary to qualitatively describe the non-Gaussian feature.Moreover, as the perturbative calculation is valid only for small M limits, it cannot capture the large M re-entrance to Gaussian behavior.
At large M ≫ 1, orientational relaxation leads to effective equilibrium-like Gaussian distributions with fluctuations governed by the kinetic temperature k B T kin ≈ 1 + Pe 2 /[d(d − 1)] (see Eq.( 15)) In the last expression, v i denotes both v x and v y .This expression captures the simulation results at M = 10 2 and Pe = 100, thereby describing the crossover to passive Gaussian behavior at large enough M .Thus, equations ( 30) and ( 31) together capture the reentrance observed as a function of M as evidenced from the numerical simulations shown in Fig. 3. Similar arguments as above can be used to describe the probability distribution of v ∥ and v ⊥ in the two limits of small and high inertia.At small M their distributions take the form These expressions agree with the numerical result at M = 10 −4 (Fig. 3).Equation (32) also suggests ⟨v ∥ ⟩ = Pe, ⟨v ⊥ ⟩ = 0, ⟨δv 2 ∥ ⟩ ∼ 1/M and ⟨δv 2 ⊥ ⟩ ∼ 1/M at small M , agreeing with equations ( 17) and ( 19).In the large M limit, distributions are described by the Pe-dependent kinetic temperature,   32) and that at M = 100 (i) are plots of equation (33).The lines at intermediate M = 0.2 (ii) are guides to the eye.
The above expressions are tested against the numerical result at M = 10 2 and show good agreements (see Fig. 3).Equation (33) suggests ⟨v ∥ ⟩ = 0 = ⟨v ⊥ ⟩ and the variance ⟨δv 2 ∥ ⟩ = ⟨δv 2 ⊥ ⟩ ∼ 1/M at the large M limit, agreeing with equations ( 17) and (19).The velocity evolution under active drive Eq.( 2) can formally be mapped to the position evolution of an overdamped ABP in a harmonic trap dr i = (v a u i − βr i )dt + √ 2 dB i (t) [54], with dimensionless activity v a and trap stiffness β.The mapping gives M = β −1 , Pe = v a /β and the strength of the thermal noise β −2 .Thus, although the algebra is equivalent, extracting the physical meaning from such mapping requires care.
A further subtlety involves the coupling of positional evolution to velocity for inertial ABPs, which we explore next.

Displacement moments, active diffusivity, and swim pressure
In this section, we calculate moments of displacement vector to explore their time evolution.We begin by using ψ = r and the initial condition ⟨ψ⟩ 0 = 0 in equation( 5) to get ⟨r⟩ s = ⟨v⟩ s /s .Further using equation( 6), we get Its inverse Laplace transform gives In the limit of M → 0 the above expression gives ⟨r(t)⟩ = Pe û0 (d−1) [1−e −(d−1)t ] in agreement with the evolution of overdamped ABPs [31].

Mean squared displacement (MSD) and asymptotic diffusivity
Let us now consider ψ = r 2 and the initial condition r(t = 0) = 0 in equation (5).It is easy to see ⟨v v ψ = 0, and ∇ 2 u ψ = 0. Thus, equation( 5) leads to ⟨r 2 ⟩ s = 2⟨v • r⟩ s /s.The calculation of the position-velocity cross-correlation ⟨v • r⟩ s proceeds in the same way using ψ = v • r in equation( 5), using initial condition To complete the calculation, one needs to evaluate the position-orientation crosscorrelation ⟨ û • r⟩ s , using ψ = û • r in equation (5) with Note that ⟨v 2 ⟩ s and ⟨ û • v⟩ s were already calculated in equations ( 11) and ( 8).Thus, we find The inverse Laplace transform of equation (38) gives the full time dependence of MSD, As shown in Fig. 4(a), the above equation captures the numerical simulation results of MSD.Moreover, the MSD shows a curious ∼ t 3 scaling at short times before reaching the asymptotic diffusive scaling.Note that the long-time diffusive behavior is independent of inertia, and the diffusion constant has the same form as overdamped ABPs [19,31].To analyze the short-time behavior, we expand ⟨r 2 ⟩(t) around t = 0  where 41) agrees with the observation of shorttime ballistic behavior of MSD [19] if the initial velocity is drawn from the equilibrium distribution v 2 0 = d/M .However, a choice of v 0 = 0 leads to ⟨r 2 ⟩ ∼ t 3 at short times in agreement with Fig. 4, and the series expansion simplifies to The expansion explains the short-time crossover from ⟨r 2 ⟩ ∼ t 3 to ⟨r 2 ⟩ ∼ t 4 observed from numerical simulations near t I = 8dM/[3(Pe 2 M − 2d)] provided Pe 2 M > 2d, a criterion fulfilled by the lines (i) (M = 100, Pe = 100) and (ii) (M = 0.2, Pe = 100) in Fig. 4(a).Note that for lines (iii) (M = 10 −4 , Pe = 100) and (iv) (M = 0.2, Pe = 1), the criterion Pe 2 M > 2d is not obeyed and as a result beyond the initial ∼ t 3 scaling, ⟨r 2 ⟩ directly approaches the asymptotic ∼ t behavior.

Position orientation cross-correlation and swim pressure
While calculating ⟨r 2 ⟩, we computed the cross-correlation ⟨ û • r⟩.Here, we consider its role in determining the swim pressure of ABPs.Note that the virial term for an ideal gas of N ABPs due to activity π a = (1/d V ) N i=1 ⟨r ′ i • γPe ′ ûi ⟩ with V denoting an effective volume [61][62][63].The dimensionless pressure has the form P = ρ(Pe/d)⟨ û • r⟩ st , where ρ denotes the dimensionless particle density and ⟨ û • r⟩ st = Pe (d−1)(1+(d−1)M ) .Thus, the swim pressure Unlike in passive particles or overdamped ABPs, the swim pressure turns out to be a function of inertia M .Note that P can be re-expressed in terms of the excess diffusivity eff .This result in the overdamped limit agrees with earlier calculations for ABPs in d = 3 [62].However, for inertial ABPs, the swim pressure P vanishes as M −1 in the limit of large inertia.The decrease of swim stress with translational inertia was observed before in Ref. [32,34].

Fourth moment of displacement
The calculation of the fourth moment in the Laplace space is long but straightforward.We show the essential steps of this calculation in Appendix B. The inverse Laplace transform can be used to obtain ⟨r 4 ⟩(t).As the expression is too lengthy to show here, we plot the expression ⟨r 4 ⟩(t) as a function of time and compare it with the result of simulation for various M and Pe values in Fig. 4(b) to find good agreement.Here, we analyze the expression in the short and long time limits.The asymptotic behavior of ⟨r 4 ⟩ in the long time limit gives with the coefficient of the t 2 scaling independent of inertia.The t 6 scaling in the t → 0 limit observed in Fig. 4(b) can be understood performing a series expansion of the expression around t = 0 with initial velocity v 0 = 0 where Note that the shortest time ∼ t 6 scaling is due to equilibrium fluctuations, independent of activity.This behavior with initial velocity v 0 = 0 starkly contrasts the shortest-time ⟨r 4 ⟩(t) ∼ t 2 scaling observed for overdamped ABPs [31].The short-time ∼ t 2 scaling for overdamped ABPs is recovered if one keeps a non-zero v 0 while taking the limit M → 0. The first non-equilibrium influence comes in ⟨r 4 ⟩(t) determining the crossover to ∼ t 7 scaling at t I = 4dM/[3(Pe 2 M − 2d)].This crossover will be observable if Pe 2 M > 2d.
As a result, we see this crossover in lines (i) and (ii) of Fig. 4(b), but not in lines (iii) and (iv).

Discussion: Comparison with RTP and AOUP
As the introduction mentions, the ABP model is closely related to the RTP and AOUP models.They differ only by the statistical properties of the active forces.However, up to the second moment, these properties are equivalent.The averages over the active forces over a time scale longer than the correlation time vanish.The two-time correlations for the active forces show exponential decay with characteristic correlation times.As a result, our predictions for ABPs up to second-order moments remain valid for these other two models and, with suitable mapping, can be used to predict exact results for them.
To illustrate this, let us denote the position and velocity vectors (r, v) to express the dynamics of an active particle in d-dimensions, where the only difference between the three models is in the dynamics of the active force γv a .ABPs differ from RTPs in their rotational relaxation.For the orientational diffusion of ABPs, the heading direction auto-correlation ⟨ û(t) • û(t ′ )⟩ = e −(d−1)Dr|t−t ′ | using equation (5).RTPs undergo discrete tumbling events with rate γ r [28,64,65], such that the orientation propagator P ( û, t| û0 , 0) = e −γrt δ( û − û0 ) + (1 − e −γrt ) 1 Ω with Ω the area of a d-dimensional unit sphere.Thus, the orientational correlation controlling the active forces follows ⟨ û(t)• û(t ′ )⟩ = e −γr|t−t ′ | .In contrast, within the AOUP model, using τ a va = −v a + v ao 2τ a /d η v (t) in d dimensions with Gaussian white noise η v (t) [31,60], the steady-state auto-correlation follows ⟨v a (t) • v a (t ′ )⟩ = v 2 ao e −|t−t ′ |/τa .Thus, the active force in the three models is characterized by ⟨v a ⟩ = 0 and In the above relations, we used v r and v a to denote the amplitude of active velocity in the RTP and ABP models, respectively.Up to the second moment, all the models are equivalent with the mapping v a = v r = v ao and (d − 1)D r = γ r = 1/τ a .The connection between ABP and RTP with γ r ↔ (d − 1)D r was noted before in Ref. [28].
The units of time and length are set in the three models by (i) Using these mappings and the definitions M = m/γτ r , Pe= v a / D/τ r the expressions derived for ⟨v⟩, ⟨v 2 ⟩, ⟨r⟩, ⟨r 2 ⟩ within the ABP model remains valid for RTP and AOUP models.The same is true for the expressions of kinetic temperature, diffusivity, and pressure.The ensuing results agree with, e.g., the earlier AOUP calculations of Ref. [60] in d = 2. Finally, the velocity distribution in RTP can show Gaussian-non-Gaussian crossovers, as in ABP.In contrast, given the inbuilt Gaussian nature, AOUP can not give rise to any non-Gaussian distribution, and all its higher-order dynamical moments can be calculated using the first two moments.Whether experimental observations demonstrate a non-Gaussian nature or not can be used to determine the potential application of the AOUP model in their analysis.

Conclusion
In this paper, we conducted a thorough analysis of the dynamics of free inertial Active Brownian Particles (ABP) in the presence of thermal noise.To achieve this, we employed a method based on the Laplace transform of the Fokker-Planck equation, allowing us to obtain expressions of time-dependent dynamical moments of arbitrary observables in any d-dimensional space.We obtained explicit expressions for several moments of velocity and displacement and equal time cross-correlations.The analytical predictions were rigorously validated through direct numerical simulations in two dimensions.
Our study has uncovered several intriguing characteristics of inertial active particles, setting them apart from their overdamped and passive counterparts.Notably, the inertia-dependent characteristics manifested in the steady-state kinetic temperature and swim pressure while leaving the late-time diffusivity unchanged.Specifically, the kinetic temperature initially increased before saturation, while the pressure consistently demonstrated a monotonic decrease.
Furthermore, our study revealed fascinating re-entrant crossovers in the steadystate velocity distribution, transitioning from "passive" Gaussian to "active" non-Gaussian behaviors.To illustrate this phenomenon, we constructed a detailed "phase diagram" in the activity-inertia plane, utilizing the exact expression of d-dimensional kurtosis.Additionally, we proposed approximate velocity distribution formulas that are closely aligned with numerical simulations, especially in small and large inertia limits.These formulations also captured the transition from low inertia passive to active states.Moreover, our findings can be expanded to encompass related models such as RTP and AOUP, enabling predictions up to second moments of dynamic variables through suitable mappings.
Our exact expressions, beyond describing steady states, could be used to explain time-dependent crossovers observed in various moments of velocity and displacement.Overall, this study provides significant insight into the intricate interplay between activity, inertia, and thermal noise governing the behavior of inertial active particles.
to understand the differences in short-time scaling behaviors from the v 0 = 0 case illustrated in the main text.Here, in the limit t → 0, ⟨v 2 ⟩ = v 2 0 and ⟨v 4 ⟩ = v 4 0 , as one would expect.The first t-dependence in ⟨v 2 ⟩ and ⟨v 4 ⟩ both show ∼ t scaling.Note the difference from v 0 = 0 case particularly for ⟨v 4 ⟩, as the linear in t term vanishes for v 0 = 0 in equation (C.2).Further, as shown in Fig. C1, ⟨v ∥ ⟩ can even decrease with time before reaching the steady state.
Fig. C2 illustrates the steady-state variation of ⟨v ∥ ⟩ st , ⟨δv 2 ∥ ⟩ st and ⟨δv 2 ⊥ ⟩ st with M in 2d using Pe = 100.The lines are analytic functions, and the points denote simulation results.
6).At equilibrium I = D − µk B T = 0, due to the equilibrium fluctuation-dissipation relation.Using the expression of effective diffusivity D eff = 1 + Pe 2 /[d(d − 1)] derived later in equation (40) in Sec. 5, we obtain the dimensionless form of the violation of equilibrium fluctuation-dissipation scalings in the small and large M limits are observable in Fig. C2(b) of Appendix C. Further, using the definition ⟨v 2

Figure 2 .
Figure 2. The steady-state kurtosis K as a function of M at Pe = 10 (a) and as a function of activity Pe at M = 1 (b).The solid (dashed) line shows the plot in d = 2 (d = 3).(c) "Phase diagram": A heat map of steady-state kurtosis as a function of activity Pe and mass M in d = 2. Passive and active states are represented by light (yellow) and dark (black) regions, respectively.(d) Heat maps of P (v x , v y ) in the passive state identified by point (iv) in the "phase diagram" (c) and (e) P (v x , v y ) in the active state identified by point (ii) in (c).

Figure 3 .
Figure 3. Velocity distributions in 2d at Pe = 100: (a) Probability distributions of components of velocity v i with i = x, y.Lines denote analytic expressions, and points denote simulation results.The small M Gaussian to non-Gaussian crossover is captured by equation(30) showing good agreement with simulation results at M = 10 −4 and M = 10 −2 .The large M = 10 2 result from simulations agrees with equation(31).Probability distributions of velocity components parallel P (v ∥ ) and perpendicular P (v ⊥ ) to the heading direction are shown in (b) and (c).The distributions denoted by (i), (ii), and (iii) in the legends correspond to the points marked in Fig.2(b) at fixed Pe = 100 and M = 100, 0.2, 10 −4 , respectively.The lines at M = 10 −4 (iii) are plots of equation(32) and that at M = 100 (i) are plots of equation(33).The lines at intermediate M = 0.2 (ii) are guides to the eye.
d−1)M ] .In the limit of M → 0, swim pressure is the maximum, P = ρD (a)