Inertial self-propelled particles in anisotropic environments

Self-propelled particles in anisotropic environments can exhibit a motility that depends on their orientation. This dependence is relevant for a plethora of living organisms but dif ﬁ cult to study in controlled environments. Here, we present a macroscopic system of self-propelled vibrated granular particles on a striated substrate that displays orientation-dependent motility. An extension of the active Brownian motion model involving orientation-dependent motility and inertial effects reproduces and explains our experimental observations. The model can be applied to general n -fold symmetric anisotropy and can be helpful for predictive optimization of the dynamics of active matter in complex environments


INTRODUCTION
The survival of organisms in complex environments essentially depends on their fitness and strategy to react and adapt to external conditions.In particular, a realistic environment is never isotropic but typically anisotropic, i.e., its traversability depends on the direction of motion [1].Anisotropy can be caused on various scales by many different means: by an external force arising from gravity [2,3], viscosity [4], light [5], and chemical gradients [6], electromagnetic fields [7], through steric confinement by channels, veins, and anisotropic porous media [8,9], or by motion in a liquid-crystalline [10][11][12][13][14] or crystalline [15][16][17] medium.Anisotropic environments can have a pronounced impact on the motion of self-propelled particles.These "active" particles convert energy from their environment into directed motion and comprise both living organisms and artificial inanimate objects, like activated colloids [18][19][20], granules [21][22][23][24][25], and robots [26][27][28].Standard models of self-propelled particles [29] assume that the propulsion force is isotropic in the sense that it always points into the direction of the particle orientation with a constant self-propulsion speed even in an inhomogeneous environment [30][31][32][33][34][35][36].In anisotropic environments, a dependence of the selfpropulsion speed of the particle on its orientation is frequently observed, i.e., some biological organisms react to their environment in a sense that the propulsion force depends on their orientation relative to the environment.For instance, microorganisms can move faster towards light sources [37] or in the direction of food sources [38].Additionally, flying animals such as bees and birds control their flying speed by relative changes of their environment, which in turn leads to anisotropic flying velocities within structured environments [39][40][41].Similarly, anisotropic movement is also observed for smaller insects like ants in guiding structures [42,43].Those macroscopic self-propelled particles in low-friction environments (e.g., such as flying insects) where the effect of anisotropy is most prominent, are also governed by inertial effects [44].This poses a challenging problem because inertia introduces correlations that can persist for longer times [45][46][47][48][49][50][51][52][53][54].
In this communication, we present an experimental realization of a self-propelled granular particle on an anisotropically structured substrate, which exhibits orientation-dependent motility.We observe pronounced anisotropy in the motion of the particle, which is well explained analytically by an extension of the active Brownian motion model with inertia and orientation-dependent motility.The orientation-dependence can be written in terms of a Fourier series which allows a general solution for anisotropic motility that can be applied to our experiments.Our findings establish a class of active matter models useful for anisotropic environments and shed light on the potential self-propulsion strategies of organisms in such anisotropies.The analytical results of our model can be particularly useful for predictive optimization of control parameters of artificial active agents, such as robots [26][27][28], to better explore anisotropic environments [55].

Experimental observation of anisotropic self-propulsion
Macroscopic active matter with orientation-dependent motility can be realized from self-propelled 3D-printed agents called vibrobots (see Fig. 1a) on structured substrates.These particles are excited by vertical vibrations generated by a rectangular acrylic baseplate attached to an electromagnetic shaker.The particles stand When orientated diagonally, the particle moves with average velocity v ∥ along its orientation while simultaneously experiencing active propulsion v ⊥ perpendicular to it (g).The particle moves with decreased velocity v ∥ when perpendicularly aligned to the grooves (h).i-k Three representative trajectories with an excitation amplitude A = 1.60 g.The persistence length is noticeably shorter for perpendicularly aligned particles than for parallel aligned particles.Length ratios and velocity contributions are not to scale.
on slightly tilted legs, which causes the particles to hop forward.These legs are all tilted equally along the orientation (or symmetry axis) of the particle.The baseplate is covered with a lenticular plastic sheet on top, which is the source for the anisotropic motility.The experimental setup is depicted in Fig. 1b.An illustration with a side-view of the particle resting on such a grooved surface are shown in Fig. 1c.The vibration frequency is set to f = 80 Hz.In this frequency range, the plate vibration is sufficiently uniform [47].Three different peak acceleration amplitudes A = 1.28 g, 1.44 g, and 1.60 g are investigated, which varies the motility and motion properties of the vibrobot.For this choice of f and A, the vibration is strong enough to ensure stable vibrobot motion, but not too strong to prevent particles from falling over.
We find pronounced anisotropy in the motion of the particle and observe a modulation of the velocity parallel but also perpendicular to the orientation of the grooves, as well as an increased activity with increasing excitation amplitude.The motion of the particles is illustrated in Supplementary videos 1 -6, where we show a montage of all measured trajectories for each excitation amplitude as well as for parallel and perpendicular initial orientation, respectively.From the trajectories, the anisotropy is already visible by the naked eye, in particular when comparing parallel and perpendicular starting orientations.
This anisotropy is best illustrated when displaying all recorded trajectories (integrated and smoothed) and distinguishing parallel and perpendicular initial orientations, as shown in Fig. 1d, e.For particles starting parallel to the grooves, we observe that the peak of the density (which is linked to the starting position of the particles) is broad along and narrow perpendicular to the starting orientation since the particles tend to move faster parallel to the grooves and therefore propagate further before they reorient.In the case of perpendicular starting orientation, the density spreads more around the peak, since particles reorient near to the starting position.Hence the persistence length depends on the orientation of the particle.Surprisingly, from individual particle trajectories, we also identify a driving-force component perpendicular to the orientation, whenever a particle is not moving exactly parallel or perpendicular to the grooves.
The anisotropic self-propulsion is caused by the grooved surface of the vibrating plate.Our conjecture is that this is due to the strong dependence of the particle speed on the relative inclination angle between legs and surface [56].When resting on the vibrating plate, the legs are bent along the orientation of the particle.This deformation stores elastic energy.Then, after detaching from the base, the energy is released and the vibrobot jumps forward.When the particle is oriented perpendicular to the grooves, the legs face an elliptical half-cylinder and the relative inclination angle of the legs is decreased (see Fig. 1c).As a result, the legs will bend less compared to the case where the particle is oriented along the grooves.If the particle is diagonally aligned with the grooves, the legs will not bend along the orientation and the particle experiences a force perpendicular to its orientation.This in fact results in propulsion perpendicular to the orientation of the particle.In Fig. 1f-h, we illustrate the two velocity contributions for three different orientations of the particle.
As described in the literature, we also observe orientational fluctuations, caused by an instability of the driving mechanism to the microscopic surface roughness, and inertial delay effects due to the mass of the particles [47,57].When vibrobots are excited above a certain amplitude threshold, they begin to tumble [58].As a result, they randomly reorient while moving and eventually change the direction of their path.Figure 1i-k shows three representative trajectories with different initial orientations.Clearly, the particle does not show a deterministic motion, apart from short-time correlations due to initial orientation and inertia.The particle rather undergoes an anisotropic two dimensional random walk with a certain persistence length.
Due to the simplicity of our particles, compared to living active matter, our experiment allows us to investigate kinetic properties of particles with orientation-dependent motility, which can be useful for optimization of motion and search strategies of active matter in general.This requires an analytical description of the motion that captures the essential properties of the particle and must be applicable to general cases of anisotropic motility.

Langevin dynamics model
Finding an analytical description for macroscopic selfpropelled systems can be challenging due to the complex interaction of particles and environment.Here, we model those interactions with an effective driving force and thereby introduce a minimal model, where the interplay of orientation-dependent motility, inertia, and fluctuations, is treated in terms of a generalized active Langevin dynamics model.Our model reproduces the experimental observations quantitatively despite its complex anisotropic nature.
We assume that the particle has non-negligible mass M and moment of inertia J.The motion of such an underdamped particle is in general characterized by the translational center-of-mass velocity ṙ(t) with the center-ofmass position r(t) and the time variable t as well as by the angular velocity φ(t) and the angle of orientation ϕ(t), which denotes the angle between the orientation vector n = (cos ϕ, sin ϕ) and the positive x-axis.By taking the above considerations into account, the translational and rotational motion of the particle is governed by the force balance between inertial, frictional, self-propulsive driving, and random forces and torques M r(t) + γ t ṙ(t) = γ t v ϕ(t) + 2D t γ t ξ(t), (1) Here, γ t and γ r denote the translational and rotational friction coefficients, respectively.To take translational and rotational diffusion into account, the Langevin equations contain independent Gaussian white noise terms ξ(t) and η(t), with zero means ⟨ξ(t)⟩ = 0 and ⟨η(t)⟩ = 0 and delta-correlated variances , where i, j ∈ {x, y}.
Therein, D t and D r are the translational and rotational short-time diffusion coefficients of the particle, respectively.The brackets ⟨. ..⟩ denote the noise average in the stationary state (meaning after losing correlation with initial conditions [52]) and δ ij is the Kronecker delta.Most importantly, v(ϕ) denotes an arbitrary orientation-dependent motility which accounts for the interaction between the particle and environment.For mathematical convenience, we represent v(ϕ) as a Fourier series where c k is the Fourier-coefficient vector of the mode k, and i denotes the imaginary unit.This representation lets us solve the model for any type of orientationdependence and then apply the results to our experimental system.In particular, this description can be used for different experimental realizations ranging from anisotropic illuminated Janus particles, triangular microparticles in traveling ultrasound waves, and the motion of living insects in guiding structures to the specific setup studied in this communication [59,60].In general, for a given propulsion velocity v(ϕ), these Fourier coefficients can be calculated as The seminal case of isotropic propulsion is recovered for the two non-zero coefficients c ±1 = v(1, ∓i)/2.Note that we exclude the mode k = 0 in Eq. ( 3), which would correspond to a drift velocity induced by a constant external force (e.g., gravity) not measured in the experiment.
Moreover, as typical 3D-printed particles are not perfectly symmetrical, they tend to perform circular motions on long time scales.To capture this behaviour, we assume a systematic torque which acts on the particle and leads to an angular speed ω.In contrast to v(ϕ), we measured no orientational dependency in the angular speed which could in principle be caused by the anisotropic substrate.
Concluding, our theoretical model depends on a number of parameters: the angular velocity ω, the rotational diffusion coefficient D r , the rotational friction time τ r = J/γ r , the set of Fourier coefficients {c k } describing the anisotropic motility, the translational diffusion coefficient D t and the translational friction time τ t = M/γ t .In the context of the experimental observations, we assume that the vibrobot is moving with an orientationdependent velocity (4) where n(ϕ) = (cos ϕ, sin ϕ) is pointing parallel and n⊥ (ϕ) = (− sin ϕ, cos ϕ) is pointing perpendicular to the particle's orientation.The sine and cosine terms in Eq. ( 4) reflect the orientation dependence of the particle velocity and the symmetry of the system.This adds the parallel speed v ∥ , the parallel speed anisotropy δv ∥ , and the perpendicular speed anisotropy δv ⊥ , leading to a total of 8 independent parameters.The four non-zero Fourier coefficients of Eq. ( 4) read These parameters are determined from analytic fits to the experimental results.We use temporal correlation functions, like the orientational correlation function C(t) = ⟨n(t) • n(0)⟩ and the velocity correlation function Z(t) = ⟨ṙ(t) • ṙ(0)⟩, to determine the relevant timescales and diffusion coefficients.Further stationary observables, like the mean translational velocity v 0 = ⟨ṙ(0)⟩ and the mean angular velocity ⟨ φ(0)⟩, are used to estimate all motility parameters.More information on the parameter estimation can be found in the Methods section and the parameter values are listed in Tab.I.In the following, we compare the experimental data with analytic predictions derived from the theoretical model and discuss the anisotropy found in several observables.

Comparison between analytical results and experiment
As described above, the mean self-propulsion strongly depends on the relative orientation of the particle with respect to the groove direction.The model describes this via two orthogonal velocity components.In Fig. 2, we separately show the mean velocity along the body-axis v ∥ = v 0 • n and perpendicular to it v ⊥ = v 0 • n⊥ as functions of the orientation ϕ.The parallel contribution v ∥ in Fig. 2a shows considerably greater propulsion along the grooves than perpendicular to them.For the perpendicular contribution (see Fig. 2b) we find the assumed sin(2ϕ)-modulation (see Eq. ( 4)), which has an alignment effect on the overall velocity direction in favor of the groove direction.Overall, we measure increased activity for larger excitation amplitudes while the degree of anisotropy remains almost the same for all three measurements.From the theoretical side, the mean instantaneous velocity v 0 = ⟨ṙ(0)⟩ at a specific orientation ϕ 0 can be computed in general as instatanteous with the dimensionless coefficients and the generalized incomplete gamma function Γ(s, x 1 , x 2 ) = x2 x1 t s−1 e −t dt.The analytic result is plotted in Fig. 2 and yields good agreement with the experimental data.In contrast to overdamped motion, where the particle's mean velocity is simply equal to the internal self-propulsion velocity, here the particle moves on average with a smaller velocity due to inertial delay effects, i.e., |v 0 (ϕ)| ≤ |v(ϕ)|.Further, the faster varying contributions (i.e., the higher Fourier modes) of the propulsion are more affected by these inertial delay effects, resulting in a more isotropic mean velocity for in- creasing mass M .Conversely, the anisotropy is restored for increasing moment of inertia: lim J→∞ v 0 (ϕ) = v(ϕ).
A suitable quantifier for the presence of inertial effects is the delay function 48,54].This function quantifies the average difference between the projection of the initial velocity on the orientation and the projection of the initial orientation on the velocity.In overdamped systems, this function is zero at all times.Here, we find that this function is significantly different from zero in particular for large excitation amplitudes A (see the Methods section).The standard delay function can be generalized to resolve anisotropy in the system by conditioning the average with a specific initial orientation ϕ 0 at time t = 0.In Fig. 3 we plot the anisotropic delay function d ϕ0 (t) both as a function of ϕ 0 for given t and as a function of t for given ϕ 0 .We compare the experimental data with simulations which follow Eqs.( 1) and ( 2) and are initialized similar to the experiments.The delay function is a highly fluctuating quantity making the experimental data difficult to interpret.The simulated data suggests an isotropic delay for short times and a larger delay along the grooves as time proceeds mimicking the modulation of the self-propulsion velocity.The simulated data always fits within the standard error of the experimental data.
For stochastic processes, it is common to analyze the first and seconds moments of the motion, i.e., the mean and mean square displacement.In anisotropic systems, these quantities will strongly depend on the initial orientation of a particle.In Fig. 4, we compare the experimental mean displacement ⟨∆r(t)⟩ conditioned at different initial orientations ϕ 0 with that resulting from our theoretical model.To demonstrate the effect of the orientation-dependent motility, we show the mean displacement as a function of the initial orientation ϕ 0 after fixed times t forming elliptic-like shapes in the xy−plane (see Fig. 4a).In Fig. 4b, we plot the absolute mean displacement |⟨∆r(t)⟩| as a function of time t for particles which are initially orientated along the grooves (blue) and for those starting perpendicular to the grooves (red).The experimental data fit within theoretical results for short time, where the particle moves linearly in time with ⟨∆r(t)⟩ = v 0 t + O(t 2 ).For longer time, confinement effects play an increasing role.Since recordings are stopped once a particle hits the boundary, events where the particle reorients beforehand dominate the statistic.As a consequence, the measured mean displacement decreases for times larger than the mean first-passage time of hitting the boundary.We perform simulations with absorbing boundaries and find an excellent agreement for all experimental accessible time scales (indicated by the black dashed curves in Fig. 4b).Without confinement, the theoretical mean displacement saturates to an anisotropic persistence length L p = lim t→∞ ⟨∆r(t)⟩ for long times with the persistence time of mode k and Ω k = D r τ r k 2 + iωτ r k.The persistence length L p consists of two contributions: the first term is given by the mean stationary velocity v 0 which is damped over the translational friction time τ t .The second term in Eq. ( 6) describes the active propulsion getting decorrelated due to the rotational noise D r .Again, the degree of anisotropy increases as a function of the moment of inertia J.For vanishing angular speed ω = 0, we find the following asymptotic behavior for small and large J, respectively: Note that for large J the contribution of higher modes decays only linearly instead of quadratically, demonstrating the relevance of the moment of inertia as an important control parameter.Last, we address the mean-square displacement, which is most commonly investigated for passive and active Brownian motion.In Fig. 5, we compare the experimentally determined mean-square displacement with the corresponding theoretical result.For short times, the particle is moving ballistically, as ⟨∆r 2 (t)⟩ = ⟨ṙ 2 (0)⟩ t 2 +O(t 3 ) (see Fig. 5a).For larger times, the particle transitions towards a diffusive regime ⟨∆r 2 (t)⟩ ∼ 4D L t, which is characterized by the long-time diffusion coefficient Similar to the mean displacement, the mean-square displacement is affected by the confinement for long times which hinders the particle to reach a diffusive state.In Fig. 5b, we show the mean-square displacement parallel and perpendicular to the grooves comparing experiment, theory, and simulation.The mean-square displacement is non-monotonic in time due to the confinement.At longer times, the particle needs to reorient before hitting the wall.The non-monotonic behavior results from the persistency of the particle and therefore is not observed for passive particles.The particle makes larger displacements along the grooves than perpendicular to them.In the absence of confinement, this anisotropy can persist even in the long-time limit characterized by the long-time diffusion matrix for i, j ∈ {x, y}.The eigenvalues of this matrix are given as D ± = D L ± ∆D L , with the long-time anisotropy which describes the long-time diffusion along the principal axes of maximal and minimal diffusion, respectively.The existence of a long-time anisotropy ∆D L ̸ = 0 will depend in general on the specific form of v(ϕ).

DISCUSSION
Anisotropic motility has a strong impact on the motion of active particles both on short and long time scales.Our experiments demonstrate this explicitly for short and intermediate times and implicitly for long time-scales through simulations.Anisotropy persists for long times in the mean and mean-square displacement.We derived an analytical description that explains this behavior in terms of the Fourier series of the anisotropic driving term.The Fourier modes of the motility are linked to different time scales that add up and have an effect on the stationary mean velocity, persistence length and long-time diffusion.Specifically, these quantities are mostly affected by the low-order Fourier coefficients.
Our theoretical results predict that the degree of anisotropy is not only set by the orientation-dependent motility itself but depends non-trivially on all time scales 1/D r , 1/|ω|, τ t , and τ r of the model.In Fig. 6, we depict the anisotropy of the stationary mean velocity, persistence length, and long-time diffusion for different values of the moment of inertia J and two exemplary orientation-dependent motilities v(ϕ) = v(1 + cos(nϕ))n(ϕ) with 2-fold symmetry (n = 2) and 3-fold symmetry (n = 3).In general, the mass and the moment of inertia have contrary effects on the anisotropy  for short and intermediate times.For increasing mass, the dynamics of the particle involves stronger delay effects, smoothing the trajectories of the particle and effectively decreasing the anisotropy.On the other hand, increasing the moment of inertia leads to more resistance to reorientation and subsequently to higher persistence.The stationary parallel velocity in Fig. 6a,b shows an increasing degree of anisotropy (being the ratio of outermost points to the innermost points on these curves) for increasing moment of inertia J.For the persistence length (see Fig. 6c,d), the degree of anisotropy remains fairly invariant with increasing J but overall we find a large persistence length (recalling Eq. ( 8)).Note that the mean displacement and thus the persistence length inherit the symmetry of the driving velocity v(ϕ).This symmetry is in general lost for long times, since the longtime diffusion can either follow a 2-fold symmetric modulation or behaves fully isotropic in every direction (see Fig. 6e,f).In fact, for motilities with higher rotational symmetry than two-fold, the long-time diffusion is always isotropic.Thus, we like to stress that even a system showing isotropic diffusion can hide anisotropic dynamics on shorter time scales.
As an outlook, we want to highlight the intriguing possibilities that arise from combining position-and orientation-dependent motility.This opens up avenues to explore migration in gradients of anisotropy.Our experimental system relies on pre-molded lenticular sheets to create the anisotropic substrate.However, by employing a larger 3D printer or an engraving tool, more complex substrates could be generated, for instance, to introduce gradients in anisotropy.For proof of principle, we assume that the particle exhibits 2-fold symmetric motility v(r, ϕ) = (v + δv(r) cos(2ϕ)) n(ϕ), where the anisotropy δv(r) increases along the direction of ŝ.Specifically, we employ a logistic function to describe the spatial dependency, δv(r) = v/(1 + e −κr •ŝ ), with κ representing the growth rate.This function yields maximum anisotropy (δv(r) = v) for r • ŝ → ∞, and isotropic motility (δv(r) = 0) for r • ŝ → −∞, with a symmetric slope around the origin (see Fig. 7b).For simplicity, we consider only the overdamped case (m = J = 0) and examine the mean position along the gradient ⟨r(t) • ŝ⟩ for particles initially starting in the origin (which corresponds to the inflection point of δv(r)).In Fig. 7a, we present the mean position for different growth rates κ and gradient directions ŝ.We observe opposing behavior depending on whether the gradient is aligned with x or ŷ, i.e., positive displacement when ŝ = x and negative dis-placement for ŝ = ŷ.Thus, the particle exhibits motion parallel or antiparallel to the gradient towards regions where the motility increases.Specifically, for ŝ = x the particle moves towards more anisotropic regions, whereas for ŝ = ŷ it moves towards more isotropic regions.Since the rotational dynamics is independent of the particle position, there is always an equal probability of moving upward and downward the gradient.However, when particles move towards regions of higher motility, they experience simultaneous acceleration within the persistence time.Consequently, the persistence length is always greater in the direction of increasing motility compared to the opposite direction, indicating a displacement for long times towards regions of higher motility.We like to stress that this result is not contradictory to previous studies in spatial motility fields, where the stationary positional probability shows accumulation in regions of low motility [61].Here, we are considering a gradient in an unbounded space.Thus, we are not reaching stationarity within finite simulation time.
Our model could be useful to predictively optimize driving parameters for the navigation of active matter in anisotropic environments [62][63][64][65], for instance robotic systems.In particular, the persistence length is an important control parameter that strongly impacts collective phenomena, like motility-induced phase separation [66][67][68].Swarms of self-propelled particles moving with an orientation-dependent motility would be an interesting topic for future research, for which our model provides a baseline [69][70][71][72].

Particle fabrication
The particle used in this work has been manufactured by 3D-printing using a stereolithographic acrylic based photopolymer 3D printer (Formlabs Form 2, using Grey V3 material, identical to Ref. [47]). Figure 1a shows an image of the particle.It consists of a cylindrical core (diameter 9 mm, height 4 mm) and a cap (diameter 15 mm, height 2 mm).Seven tilted cylindrical legs (diameter 0.8 mm, inclination angle 4 degrees) are attached to the cap in a regular heptagon around the bottom cylinder.The legs are tilted parallel to each other defining the orientation of the particle.The length of the legs is chosen such that the bottom of the particle is lifted by 1 mm above the surface.The particle is marked with a sticker from which the orientation can be determined using computational image processing.The particle's mass is about m = 0.76 g.From the particle's mass and shape, its moment of inertia is computed to be J = 1.64 × 10 −8 kg m 2 , assuming homogeneous density.

Experimental setup and analysis
Particle motion is excited by vertical vibrations of a rectangular acrylic baseplate (side length 300 mm, thickness 15 mm) with a lenticular plastic sheet on top, attached to an electromagnetic shaker (Tira TV 51140).The sheet's surface consists of equally spaced elliptical half-cylinders with a density of 0.787 mm −1 (20 lines per inch) and a groove depth of 0.315 mm.An illustration and a cross-section of the particle resting on such a grooved surface are shown in Fig. 1c, respectively.Lenticular sheets of this kind are typically used in digital printing or displays to create images with the illusion of depth.
Here, we use it to induce an anisotropic driving of the particle parallel and perpendicular to the lines, since the speed of the particle is very sensitive to the contact angle of the legs to the surface.Note that the width and height of the grooves are chosen such that the particle legs cannot be significantly trapped (see Fig. 1c), in order to prevent the particle simply from sliding along grooves.
The tilt of the plate is adjusted with an accuracy of 0.01 • to minimize gravitational drift.The vibration frequency is set to f = 80 Hz and three different peak acceleration amplitudes A = 1.28 g, 1.44 g and 1.60 g are studied.
A mid-to-high-speed camera system (Allied Vision Mako-U130B) operating at 150 frames per second is used to record the experiment with a spatial resolution of 1024 × 1024 pixels.The particle location and orientation are determined and tracked using standard image recognition methods (Hough transform and morphological image region analysis) to a spatial accuracy of about ±3 × 10 −5 m and a orientational accuracy of ±0.74 • [47].Multiple single trajectories are recorded for each amplitude, until 20 min of data are acquired per recording.Half of the recorded time the particle starts parallel and the other half of the time it starts perpendicular to the grooves.Events involving particle-border collisions mark a trajectory's termination and are subsequently discarded, resulting in trajectories of various lengths.
The velocity was calculated from the displacement of successive positions of the particle as v(t) = (r(t + ∆t) − r(t)) /∆t, where ∆t = 1/150 s is the time between two frames.The time steps are not fully equidistant between recorded frames, therefore the experimental data were linearly interpolated to obtain equidistant points.Experimental means with respect to a specific initial orientation ϕ 0 were calculated by averaging in the interval [ϕ 0 − δϕ, ϕ 0 + δϕ].We chose δϕ = 10 • and modified the theoretical results accordingly by exp(ikϕ) → exp(ikϕ) sin(kδϕ)/(kδϕ).We took advantage of the rotational and inflection symmetries of the experiment (by rotating some trajectories by 180 degrees) to increase the angular statistics for the mean displacement.

Analytic results
Both the translational velocity ṙ(t) and the angular velocity φ(t) undergo a simple stochastic process for which a general solution is easily obtained (see Eqs. ( 1) and ( 2)).Several dynamical correlation function as well as low-order moments can be consequently calculated using standard methods of stochastic calculus [73].The orientational correlation function C(t) = ⟨n(t) • n(0)⟩ displays a double exponential decay C(t) = cos(ωt)e −Dr(t−τr(1−e −t/τr )) , (as previously discussed in Ref. [45][46][47]).The velocity correlation function Z(t) = ⟨ṙ(t) • ṙ(0)⟩ is given as where the Fourier-coefficient vectors are determined by the orientation-dependent motility, as c k = π −π v(ϕ) exp(−ikϕ)/(2π) dϕ (see Eq. ( 3)), and with Ω ± k = D r τ r k 2 ±(iωτ r k+τ r /τ t ) and S k = D r τ r k 2 .The real part is denoted by Re{. . .} and the generalized incomplete gamma function is Γ(s, x 1 , x 2 ) = x2 x1 t s−1 e −t dt.The delay function measuring the difference between the direction of the velocity and the current orientation, which coincides with the result for isotropic selfpropulsion [47] (due to the projection onto the orientation).Next, we give the mean displacement ⟨∆r(t)⟩ = ⟨r(t) − r 0 ⟩ under the condition that initially the position r 0 and the orientation ϕ 0 are prescribed, with the stationary velocity v 0 (see Eq. ( 5)), and Ω k = D r τ r k 2 + iωτ r k .Lastly, we provide the result for the mean-square displacement ⟨∆r 2 (t)⟩ = ⟨(r(t) − r 0 ) 2 ⟩ which can be expressed as with the long-time diffusion coefficient D L (see Eq. ( 9)), the velocity correlation function Z(t) (see Eq. ( 13)) and where 2 F 2 denotes the generalized hypergeometric function.Last we remark that in the overdamped limit, i.e. m → 0 and J → 0, we recover the results of orientation-dependent motility in underdamped systems [7] and similarly for an isotropic self-propulsion v(ϕ) = v 0 n(ϕ), we obtain the expressions of Ref. [52].

Parameter estimation
The underdamped active Brownian motion model depends on eight independent parameters.All parameters were obtained using the MATLAB standard optimizer fminsearch (Nelder-Mead optimization of a function of several variables on an unbounded domain).Our cost function consists of five terms covering different parameters.Each term is constructed as follows: The absolute deviation between the experimental mean and the analytical expectation is weighted with the standard error of the mean and then averaged over time or orientation.This procedure takes into account the experimental uncertainty [74].At the same time, the value of our cost function quantifies the fit itself.We call a fit sufficiently representative of the experimental mean if the mean deviation between experimental mean and analytical expectation is no greater than one standard error.We use this definition to determine an error interval for our optimal parameters.The orientational correlation function C(t) (see Eq. ( 12)) is used to determine the rotational diffusion constant D r and the rotational friction time τ r .In addition, we use the mean stationary angular velocity ⟨ φ(0)⟩ = ω to determine the angular speed ω.Further, we use the velocity correlation function Z(t) (see Eq. ( 13)) to extract values for the translational friction time τ t and the translational short-time diffusion coefficient D t .Lastly, we use the mean stationary velocity v 0 (see Eq. ( 5)), which is projected parallel (v ∥ = v 0 • n) and perpendicular (v ⊥ = v 0 • n⊥ ) to the body axis, to determine all the motility parameters v ⊥ , δv ∥ , and δv ⊥ .In Fig. 8a-d, the analytic fitting curves to the experimental data are shown and the resulting set of parameter is listed in Tab.I.For vibrobots, the delay function d(t) (see Eq. ( 15)) proved to be a sensitive measure for the quality of the determined parameter-set [47].Figure 8e shows good agreement between theory and experiment for all three measurements.Note that ω and D t are not significantly different from zero for our particles.However, they are included in the model, since they can be relevant for different experimental realizations in the literature [23,47,75].For the inertial time scales τ r and τ t we see an increase for increasing A, and τ t only becomes significant for A = 1.44 g and A = 1.60 g.This is likely caused by the reduction of friction at larger A.

Simulation
Numerical data for a self-propelled particle with orientation-dependent motility enclosed by absorbing boundaries are included in Figs. 3, 4b, and 5b.Equations (1) and (2) were discretized to perform Langevin dynamics simulations using first-order finite difference discretization.For these simulations, we chose the time step size ∆t = 10 −2 s and we performed 10 5 realizations in Fig. 4b, and 5b and 2000 realizations in Fig. 3  Figure 7 presents simulation data for a self-propelled particle with additional positional dependency in its motility v(r, ϕ).In this case, we used a time-step size of ∆t = 10 −1 s and performed 10 5 realizations to compute the ensemble averages.Each trajectory started from the origin (x 0 y 0 ) = (0, 0) with a random orientation. .Initial increase of the velocity and orientationdependence of the angular velocity.a Normalized initial velocity ⟨| ṙ(t)|⟩/⟨| ṙ(0.5 s)|⟩ of a particle starting from rest (solid curves) for initial orientations ϕ0 = 0 (cyan) and ϕ0 = π/2 (yellow).The difference between horizontal and vertical starting orientation is small and within the variation observed in simulations (dashed lines).This suggests that friction in the parallel and perpendicular direction is not significantly different.b Mean angular velocity ⟨ φ⟩ as a function of particle orientation ϕ (solid curves).The slight periodicity is an artifact of the initial conditions and can be reproduced by simulating with equal conditions (dashed curves).Both are shown for excitation amplitudes A = 1.28 g (upper row), A = 1.44 g (middle row), and A = 1.60 g (lower row).

Orientation-dependent friction and torque
In our underdamped active Brownian particle model, the influence of an anisotropic environment is effectively described by an orientation-dependent motility.However, in more general cases, self-propelled particles may undergo more complicated dynamics in anisotropic environments, resulting in orientational dependencies in various model parameters beyond motility.Here, we provide supplementary information on why we chose not to model anisotropic motion using an anisotropic friction matrix or orientation-dependent torque.
While using an anisotropic friction matrix may seem like an intuitive approach to describe anisotropic motion, it is not suitable for our experimental particles.Most of the time, our particles do not have direct contact with the substrate, and energy dissipation only occurs during collisions.Only the effective angle between legs and plate is different in perpendicular and parallel directions causing anisotropic self-propulsion.Experimental evidence supporting this conjecture is presented in Fig. 9a, showing the initial increase of the velocity, which reaches an intermediate state after approximately 0.5 s regardless of the initial orientation.From this we conclude an isotropic translational damping time and thus isotropic friction.
Furthermore, particles may experience anisotropic torques.Surprisingly, measuring the angular velocity for given orientations suggests no such torques in our experiment (see Fig. 9b).The slight modulation observed in the angular velocity, which does not align with the overall two-fold symmetry of the environment, can be attributed to initialization bias.Simulation results assuming isotropic torque, indicated by the black curves in Figure 8b, are consistent with the experiment.

Figure 1 .
Figure1.Description of experimental system and observations.a Vibrationally driven self-propelled particle (vibrobot) manufactured by 3D printing.The white cross indicates the particle orientation.Scale bar represents 1 cm.b Experimental setup: Rectangular acrylic baseplate attached to an electromagnetic shaker.The size (width × length) of the top-mounted plate equals 30 cm × 30 cm. c Cross-section of the anisotropic substrate (lenticular foil) with particle to scale.d, e Trajectory density for vibrobots starting parallel (d) and perpendicular (e) to grooves with an excitation amplitude A = 1.28 g. f -h Sketch of the two velocity contributions.The particle moves with increased velocity v ∥ when aligned along the grooves (f ).When orientated diagonally, the particle moves with average velocity v ∥ along its orientation while simultaneously experiencing active propulsion v ⊥ perpendicular to it (g).The particle moves with decreased velocity v ∥ when perpendicularly aligned to the grooves (h).i-k Three representative trajectories with an excitation amplitude A = 1.60 g.The persistence length is noticeably shorter for perpendicularly aligned particles than for parallel aligned particles.Length ratios and velocity contributions are not to scale.

Figure 2 .
Figure 2. Orientation dependence of stationary velocity.a Stationary parallel velocity v ∥ and b stationary perpendicular velocity v ⊥ plotted as a function of the orientation angle ϕ for three different excitation amplitudes A = 1.28 g (upper row), A = 1.44 g (middle row), and 1.60 g (lower row).Solid dark blue and dashed red curves show the experimental data and analytical results, respectively.Blue experimental error intervals represent the standard error of the mean.

Figure 3 .
Figure 3.Anisotropic delay function.a The anisotropic delay function d ϕ 0 (t) plotted as function of the initial orientation ϕ0 after fixed times t = 0.1 s, t = 0.4 s.Solid blue and dashed red curves show the experimental and simulated data, respectively.b The anisotropic delay function d ϕ 0 (t) plotted as a function time t for parallel ϕ0 = 0 (cyan), diagonal ϕ0 = π/4 (green), and perpendicular ϕ0 = π/2 (yellow) orientations, each.Both for excitation amplitude A = 1.28 g (upper row), A = 1.44 g (middle row) and A = 1.60 g (lower row).Solid and dashed curves correspond to the experimental and simulated data (using the parameter values given in Tab.I), respectively.

Figure 4 .
Figure 4. Mean displacement.Comparison between model and measurement with excitation amplitude A = 1.28 g (upper row), A = 1.44 g (middle row), and A = 1.60 g (lower row).a The anisotropic motion of the particle is visualized by plotting the mean displacement ⟨∆r(ϕ0)⟩ for ϕ0 ∈ [0, 2π) and fixed times t = 0.2 s, t = 0.6 s and t = 1.0 s.Solid blue and dashed red curves show the experimental data and analytical results, respectively.Light blue area expresses the standard error of the mean.b The absolute mean displacement |⟨∆r(t)⟩| is plotted as a function of time t for initial orientations ϕ0 = 0 (cyan) and ϕ0 = π/2 (yellow).Solid colored curves represent the experimental data and dashed colored curves the analytic results.In addition, dashed black curves depict simulation data for a particle in confinement.Black dots correspond to the experimental values for the fixed times of Fig.4a.Theoretical predictions and simulations use the parameters given in Tab.I.

Figure 5 .
Figure 5. Mean square displacement.Comparison between model and measurement with excitation amplitude A = 1.28 g (upper row), A = 1.44 g (middle row), and A = 1.60 g (lower row).a The total mean-square displacement as a function of time t (double logarithmic scaling).Open blue circles and dashed red curves show the experimental data and analytical results, respectively.b The mean-square displacement along the x-axis (cyan) and y-axis (yellow) as functions of time t.Solid colored curves and dashed colored curves show the experimental data and analytical results, respectively.Light colored areas represent the standard error of the mean.Dashed black curves show simulation data for a particle in confinement.Theoretical predictions correspond to the parameters given in Tab.I.

Figure 6 .
Figure 6.Anisotropy of the stationary mean velocity v0, persistence length Lp, and long-time diffusion DL for various values of the moment of inertia J evaluated for a 2-fold symmetric motility (left column) and a 3-fold symmetric motility (right column).a, b Stationary mean velocity as a function of the current orientation v0(ϕ) • n(ϕ).c, d Persistence length as a function of the initial orientation Lp(ϕ) • n(ϕ).e, f Long-time diffusion projected along different directions nT (ϕ) DL n(ϕ).The moment of inertia is set to J = 0.1 γr/Dr (orange), J = γr/Dr (red), and J = 10 γr/Dr (purple).The mass is fixed at M = γt/Dr.

Figure 7 .
Figure 7. Mean displacement of a particle in an anisotropy gradient.a Mean position in the gradient direction, ⟨r(t) • ŝ⟩ for different reduced growth rates κv/Dr.The upper plane displays simulated data for the gradient direction ŝ = x, while the lower plane shows the data for ŝ = ŷ.b The anisotropy δv(r) of the orientation-dependent motility is plotted as a function of position r • ŝ, using the same reduced growth rates as in a.The plot is sideways to align with a.

Figure 8 .
Figure 8. Determination of model parameters.a Orientational correlation function C(t), b velocity correlation function Z(t), c stationary parallel velocity v ∥ d stationary perpendicular velocity v ⊥ .Solid dark blue and dashed red curves show the experimental data and analytical results, respectively.Experimental error intervals represent the standard error of the mean.The parameter values are listed in Tab.I. e Time-dependence of the delay function d(t) validating the parameters on an independent quantity.The different vibration amplitudes are A = 1.28 g (upper row), A = 1.44 g (middle row), and A = 1.60 g (lower row).

Figure 9
Figure 9.Initial increase of the velocity and orientationdependence of the angular velocity.a Normalized initial velocity ⟨| ṙ(t)|⟩/⟨| ṙ(0.5 s)|⟩ of a particle starting from rest (solid curves) for initial orientations ϕ0 = 0 (cyan) and ϕ0 = π/2 (yellow).The difference between horizontal and vertical starting orientation is small and within the variation observed in simulations (dashed lines).This suggests that friction in the parallel and perpendicular direction is not significantly different.b Mean angular velocity ⟨ φ⟩ as a function of particle orientation ϕ (solid curves).The slight periodicity is an artifact of the initial conditions and can be reproduced by simulating with equal conditions (dashed curves).Both are shown for excitation amplitudes A = 1.28 g (upper row), A = 1.44 g (middle row), and A = 1.60 g (lower row).