# Optimal navigation of microswimmers in complex and noisy environments

Lorenzo Piro  Max Planck Institute for Dynamics and Self-Organization, 37077
Göttingen, Germany  Benoît Mahault  Max Planck Institute for Dynamics and
Self-Organization, 37077 Göttingen, Germany  Ramin Golestanian
<EMAIL_ADDRESS>Max Planck Institute for Dynamics and Self-
Organization, 37077 Göttingen, Germany  Rudolf Peierls Centre for Theoretical
Physics, University of Oxford, Oxford OX1 3PU, United Kingdom

###### Abstract

We design new navigation strategies for travel time optimization of
microscopic self-propelled particles in complex and noisy environments. In
contrast to strategies relying on the results of optimal control theory, these
protocols allow for semi-autonomous navigation as they do not require control
over the microswimmer motion via external feedback loops. Although the
strategies we propose rely on simple principles, they show arrival time
statistics strikingly similar to those obtained from stochastic optimal
control theory, as well as performances that are robust to environmental
changes and strong fluctuations. These features, as well as their
applicability to more general optimization problems, make these strategies
promising candidates for the realization of optimized semi-autonomous
navigation.

##  I Introduction

The problem of finding the most efficient route towards a desired target has
important technological and medical applications, such targeted drug delivery
at the microscale [1, 2], environmental monitoring [3], and optimization of
plane routes [4]. With the advent of theoretical and experimental
prescriptions for how to make microswimmers [5, 6, 7, 8], the fictional idea
of making microscopic devices capable to delivering molecular cargo at the
desired location in human body has moved closer to reality.

A solution of this classical problem can be traced back to the work of E.
Zermelo [9], who derived the steering policy which minimizes the travel time
of a vessel moving at constant speed in presence of wind. In the context of
micoswimmers, Zermelo’s work has been extended to more general contexts,
including fuel consumption and dissipation minimization [10], time-varying
flows [11], motion on non-Euclidean spaces [12], or the role of hydrodynamic
interactions [13]. The approach of Zermelo, however, does not account for
thermal fluctuations, which play a prominent role at the micro-scale.

Optimal navigation in the presence of noise falls into the class of problems
addressed by stochastic optimal control theory [14, 15, 16, 17, 18].
Considering a cost function defined on the configuration space (e.g. the mean
arrival time to a given target starting from a specific position), a
stochastic optimization principle can be used to show that it obeys a so-
called Hamilton-Jacobi-Bellman equation, from which the optimal control map
can be obtained [16, 19]. In parallel, machine learning algorithms such as
reinforcement learning provide convenient and increasingly popular routes to
determining the optimal control landscapes [20, 21, 22, 23, 24, 25]. In
practice, implementations of such optimal policies can then be achieved via
external feedback loops for the actuation of a microswimmer motion [26, 27,
28, 29, 30, 31, 32].

On the other hand, a number of natural and artificial microswimmers exhibit
tactic behaviour [33], i.e. are able to adapt their motility in response to
external stimuli such as light [34, 35, 36], chemical concentration [37, 38,
39, 40] or viscosity [41, 42] gradients, as well as magnetic [43, 6, 44, 45]
or gravitational [46] fields, in an autonomous fashion. Harnessing guidance
provided by taxis allows microswimmers to perform complex tasks [38, 47, 48,
49] in a semi-autonomous way as it does not rely on real time external
feedback mechanisms.

Here, we show how these ideas can be applied to the problem of optimal
navigation in complex and noisy environments. Considering a minimal but non-
trivial optimization problem in two dimensions, we propose novel navigation
policies inspired by the control maps provided by the exact global
optimization problem, which can be implemented in a semi-autonomous fashion.
Using extensive Brownian dynamics simulations, we show that the new policies
show performances comparable to that obtained from stochastic optimal control
theory, and demonstrate their robustness upon changes in the environment, as
well as to positional and rotational fluctuations. Lastly, we illustrate how
the semi-autonomous policies can be conveniently adapted to a broader class of
problems such as navigation on curved manifolds.

Figure 1: Optimal navigation in a Taylor-Green flow. (a) The Taylor-Green flow
map, where the colour codes for the intensity and the black arrows indicate
the direction. The green curve indicates Zermelo’s optimal trajectory at
vanishing noise connecting the initial point
$\bm{r}_{0}=-(\ell/2)\hat{\bm{e}}_{x}$ and the target $\bm{r}_{T}=\bm{0}$,
marked with a green and magenta circle, respectively. The red and cyan lines
show two representative stochastic trajectories following respectively the
Optimal and Straight policies at $\gamma=0.7$ and ${\rm Pe}=400$. (b) The
optimal control obtained by solving Eqs. (2) and (3) for same parameters as in
panel (a). (c) Schematic defining the quantities used for the implementation
of the semi-autonomous navigation strategies. The red dot indicates the active
particle position $\bm{r}$ and the red arrow its heading direction
$\bm{\hat{u}}$. $\bm{r}_{c}$ (red cross) marks the closest point to $\bm{r}$
on the optimal trajectory, while the two vectors $\bm{\hat{n}}$ and
$\bm{\hat{t}}$ are respectively normal and tangent to the optimal trajectory
in $\bm{r}_{c}$.

##  II Results

###  II.1 Supervised optimal navigation

We consider an overdamped self-propelled particle moving at a fixed speed
$v_{0}$ in presence of a stationary force field $\bm{f}(\bm{r})$, which may in
general include a contribution due to advection by the solvent flow velocity,
and translational diffusion with diffusivity $D$. For simplicity, we set the
friction coefficient (and mobility) to unity. The position $\bm{r}$ of the
self-propelled particle obeys the following stochastic differential equation

$\dot{\bm{r}}=v_{0}\bm{\hat{u}}+\bm{f}(\bm{r})+\sqrt{2D}\,\bm{\xi}\ ,$ (1)

where $\bm{\hat{u}}$ is the unit vector setting the direction of self-
propulsion, and $\bm{\xi}$ is a Gaussian white noise vector whose components
have zero mean and unit variance. Within this setting, the only degree of
freedom accessible to the self-propelled particle for navigation is its
orientation $\bm{\hat{u}}$. For the sake of presentation, we will first assume
a full control over $\bm{\hat{u}}$, either from external sources or by the
particle itself, while this constraint will be relaxed later.

We now want to determine the optimal navigation protocol for the particle to
reach a target position $\bm{r}_{T}$ given a force profile and an initial
position $\bm{r}_{0}$. Following standard techniques that invoke the backward
Fokker-Planck equation [50], the mean first-passage time (MFPT) ${\cal
T}(\bm{r})$ to reach $\bm{r}_{T}$ while starting initially at $\bm{r}$ can be
found as the solution of the following equation (see Methods for a
derivation):

$-v_{0}|\bm{\nabla}{\cal T}|+\bm{f}(\bm{r})\cdot\bm{\nabla}{\cal
T}+D\bm{\nabla}^{2}{\cal T}=-1\ ,$ (2)

which is to be solved subject to the boundary condition ${\cal
T}(\bm{r}_{T})=0$. Here, the prescription for optimal control comprises tuning
the orientation according to the the following rule

$\bm{\hat{u}}_{\rm opt}(\bm{r})=-\bm{\nabla}{\cal T}/|\bm{\nabla}{\cal T}|\ ,$
(3)

at every point in space. We note that Eqs. (2) and (3) provide strategies
which, by design, lead to the fastest possible trajectories on average as can
be seen in the vanishing diffusivity limit where Zermelo’s solution [9] is
recovered. Hereafter, we will refer to the strategy that corresponds to the
solution of Eqs. (2) and (3) as the _optimal policy_ (OP), and the $D=0$
optimal path as the _Zermelo path_.

We now illustrate OP by considering a simple but nontrivial setup in which a
self-propelled particle navigates in the two dimensional plane spanned by the
unit vectors $\\{\hat{\bm{e}}_{x},\hat{\bm{e}}_{y}\\}$ between neighbouring
stationary points of a Taylor-Green vortex flow (see the colour map and black
arrows in Fig. 1a). This configuration corresponds to
$\bm{f}(\bm{r})=v_{f}[\cos(ky)\sin(kx)\,\hat{\bm{e}}_{x}-\cos(kx)\sin(ky)\,\hat{\bm{e}}_{y}]$
with $k=2\pi/\ell$ and $\ell$ being the characteristic length scale of the
flow. Rescaling space and time as $\bm{r}\to\ell\bm{r}$ and $t\to\ell
t/v_{0}$, the dynamics (1) is characterized by only two nondimensional
parameters: the ratio between flow intensity and self propulsion $\gamma\equiv
v_{f}/v_{0}$, and the Péclet number ${\rm Pe}=\ell v_{0}/D$. Here, we will
focus only on cases where the self-propulsion is always stronger than the
flow, namely, $0\leq\gamma\leq 1$.

With the setup shown in Fig. 1a, the most direct route between the departure
and arrival points requires travelling counter-flow all the way. Consequently,
the straight path becomes increasingly disadvantageous as the flow amplitude
grows, such that for $\gamma\gtrsim 0.4$ the Zermelo path makes use of the
flow profile and takes a bell-shaped curve (see the green line in Fig. 1a).
For finite Péclet number ($D\neq 0$), such a feature is moreover consistent
with the control map provided by OP, since the latter generally orients the
self-propulsion away from the straight path (Fig. 1b). Simulations of the
Brownian dynamics (1) with the control map solving Eqs. (2) and (3) indeed
reveal that the OP trajectories tend to remain close to the Zermelo path for a
broad range of Péclet number values, provided that the flow strength and
particle self-propulsion dominate over fluctuations (see e.g. the density map
in Fig. 2d and Methods for more details on the stochastic dynamics
simulations).

Figure 2: The semi-autonomous navigation policies. (a)-(c) The control
$\bm{\hat{u}}$ obtained for the AAP from Eqs. (4) and (5) in the Taylor-Green
flow at $\gamma=0.7$ with $\varepsilon=0.05$(a), 0.2(b) and 1(c). (d,e) Heat
maps of $10^{3}$ stochastic trajectories obtained from numerical simulations
of the OP(d) and AAP(e). The solid green lines in (a)-(e) represent the $D=0$
Zermelo path. (f) Example of curves leading to the determination of the
optimal sensitivities for the AP and AAP ($\gamma=0.7$ and ${\rm Pe}=400$)
with ${\langle t\rangle}$ and $t_{\rm opt}$ respectively denoting the mean
arrival time and the optimal travel time in absence of noise. The optimal
values of $\alpha$(AP) and $\varepsilon$(APP) are indicated by the vertical
dashed lines. Note the small range of values for the normalized mean arrival
time, showing the robustness of the two policies. Here, both curves are
normalized by their minimum value.

###  II.2 Semi-autonomous optimal navigation

The above observations suggest that optimized navigation in the finite Péclet
number regime may be achievable using only the local information of the
relative positions of the stochastic swimmer and the Zermelo path, as opposed
to the OP which requires a global control map $\sim\bm{\nabla}{\cal
T}(\bm{r})$.

We now explore a number of such local strategies and probe their efficiencies
in comparison with OP. For a given particle position $\bm{r}$, we define
$\bm{r}_{c}\equiv\min_{\bm{r}^{\prime}}|\bm{r}-\bm{r}^{\prime}|$ as the
corresponding closest point on the Zermelo path. Moreover, we suppose that the
latter is smooth and can be parametrized by the moving frame
$\\{\hat{\bm{t}},\hat{\bm{n}}\\}$, of tangent and normal vectors, with
$\hat{\bm{t}}$ heading towards the target as shown in Fig. 1c. Assuming that
the swimmer is able to measure its relative position to the Zermelo path, it
can regulate it by steering its self propulsion direction $\bm{\hat{u}}$ via
the following rule

$\bm{\hat{u}}\cdot\bm{\hat{n}}={\cal G}(\Delta r,\bm{f}(\bm{r}))\ ,$ (4)

where $\Delta r\equiv(\bm{r}-\bm{r_{c}})\cdot\bm{\hat{n}}$ and the function
${\cal G}\in[-1,1]$ depends on the amount of information available to the
swimmer. As the rhs of Eq. (4) depends on $\bm{r}$ solely through $\Delta r$
and the external flow $\bm{f}$, it defines a class of optimal navigation
policies relying only on the swimmer’s local knowledge of its environment.

In the simplest case where the swimmer can only determine the direction
$\hat{\bm{n}}$ to the Zermelo path (from its current position), it can choose
to keep a constant angle $\alpha$ between its self-propulsion direction and
$\hat{\bm{n}}$. Such _aligning policy_ (AP) corresponds to a protocol ${\cal
G}_{\rm AP}=\pm\cos(\alpha)$, where the $\pm$ sign ensures that
$\bm{\hat{u}}\cdot\hat{\bm{t}}\geq 0$.

Although AP allows the swimmer to remain in the vicinity of the Zermelo path,
it also slows it down as it imposes a finite angle between $\bm{\hat{u}}$ and
$\hat{\bm{t}}$ even for (arbitrarily) small separations. For swimmers able to
evaluate their relative distances to the Zermelo path, AP can thus be refined
by allowing ${\cal G}$ to depend on $\Delta r$. This defines the _adaptive
aligning policy_ (AAP). Here, we choose, for simplicity, a piecewise linear
form for ${\cal G}(\Delta r)$, namely

${\cal G}_{\rm AAP}(\Delta r)=\begin{cases}+1&\Delta r<-\varepsilon\\\ -\Delta
r/\varepsilon&|\Delta r|<\varepsilon\\\ -1&\Delta r>\varepsilon\end{cases},$
(5)

where the parameter $\varepsilon$ sets a cut-off scale above which the
stochastic particle points normally to the Zermelo path.

Figure 3: Comparison of the performances of the different policies. (a)
Arrival time distribution as a function of the normalized time $t/{\langle
t\rangle}$ for the four navigation strategies: straight policy (blue circles),
optimal policy (yellow triangles), adaptive aligning policy (green diamonds)
and aligning policy (red squares). (b) Kullback-Leibler divergence $D_{\rm
KL}$ between the numerical and theoretical (Eq. (6)) distributions as a
function of the relative flow strength $\gamma$. (c) Mean arrival time
(normalized by the noiseless optimal time $t_{\rm opt}$) as a function of
$\gamma$. Inset: Zoom-out showing the larger values taken by SP. In (a-c) the
symbols show the data obtained from Langevin simulations, while in (a) the
solid lines indicate the theoretical prediction (6) calculated with ${\langle
t\rangle}$ and $\sigma$ obtained from the numerical data. All data in (a)-(c)
are averaged over $10^{5}$ independent trajectories.

The parameters $\alpha$ and $\varepsilon$ introduced above essentially play
the role of sensitivities for AP and AAP, respectively. Their optimal values
(that minimize the mean travel time in this case) cannot be selected a priori,
and need to be determined empirically. However, the existence of such optimal
values can be intuitively understood from the control maps obtained from Eq.
(5) for various $\varepsilon$ values (with the generalization to AP being
straightforward). As shown in Figs. 2a-c, exceedingly small $\varepsilon$
values force the swimmer to mostly point normally to the Zermelo path, while
for excessively large $\varepsilon$ stochastic trajectories are less
efficiently confined and can visit less favourable flow regions. Therefore, it
is natural to expect an intermediate value of $\varepsilon$ providing the
optimal trade-off between efficient confinement and tangential motion along
the Zermelo path.

This heuristic picture is confirmed by numerical simulations showing that the
mean arrival time $\langle t\rangle$ indeed exhibits a minimum at a value
$\varepsilon=\varepsilon_{\rm opt}$ (see Fig. 2f). We moreover note that
$\langle t\rangle$ varies relatively little with $\varepsilon$, such that in
practice the policy implementation does not require a fine tuning of this
parameter. The heat map of trajectories obtained from simulations of AAP at
optimal $\varepsilon$ shows that they globally follow the Zermelo path (Fig.
2e), similarly to the OP case. Contrary to OP, however, they are not
distributed symmetrically with respect to the desired path due to a non-zero
transverse component of the flow (see Fig. 1a). Better conformity to the OP
results thus requires additional features such as a function $\cal G$ in Eq.
(4) that depends explicitly on the local flow field $\bm{f}(\bm{r})$. Here, we
restrict to the case where the swimmer is unaware of the local flow structure
around it and will address such more elaborate policies in a forthcoming
publication [51].

###  II.3 Performance assessment of the navigation policies

We now compare the performances of the two policies introduced above (AP and
AAP) with that of OP by simulating Eq. (1) with the controls defined by (3)
and (4) in the Taylor-Green flow setup (Fig. 1a). To illustrate the relevance
of nontrivial policies, we moreover consider the _straight policy_ (SP) for
which the swimmer always points towards the direction of the target regardless
of its current position (see the cyan curve in Fig. 1a for a representative
trajectory).

Figure 4: Robustness of the navigation policies. (a) Colour maps showing the
performances of OP (left) and AAP (right) as functions of the relative flow
amplitude $\gamma$ and the Péclet number ${\rm Pe}$. (b) Normalized mean
arrival time as a function of the Péclet number at $\gamma=0.7$ for the four
policies (horizontal axis in log scale). (c) Inverse of the mean arrival time
as a function of the rotational noise strength $\eta^{2}$ at ${\rm Pe}=400$
and $\gamma=0.7$ for OP, AAP and SP (caption is the same as in (b)). Here the
data is normalized by the value ${\langle t\rangle}$ measured at $\eta=0$.
Solid lines show linear fits in the small $\eta$ regime. (d), (e) Heat maps of
$10^{3}$ stochastic trajectories obtained from numerical simulations at
$\gamma=0.7$ and ${\rm Pe}=400$ of OP and AAP, respectively. In both panels
the green curve shows the Zermelo path connecting the initial point
$\bm{r}_{0}=-(3\ell/2)\hat{\bm{e}}_{x}$ (green circle) and the target
$\bm{r}_{T}=\bm{0}$ (magenta circle). (f) Arrival time probability
distributions as functions of $t/{\langle t\rangle}$ corresponding to (d,e) as
well as for the SP case under similar conditions; the caption is the same as
(b). Solid lines show the theoretical curves obtained similarly to Fig. 3(a).
All data in (a)-(c) and (f) are averaged over $10^{5}$ independent
trajectories.

We first work at fixed ${\rm Pe}=400$ and vary the relative flow amplitude
$\gamma\in[0,1]$. Figure 3a shows the arrival time distributions $P(\tau)$
with $\tau\equiv t/\langle t\rangle$ for each of the four policies (OP, AP,
AAP, and SP) at $\gamma=0.7$. Remarkable overlap between the OP and AAP
distributions can be observed. We moreover find that they are both well
described by a so-called inverse Gaussian distribution of the form

$P(\tau)=\sqrt{\frac{{\langle
t\rangle}^{2}}{2\pi\sigma^{2}\tau^{3}}}\exp\biggr{[}-\frac{{\langle
t\rangle}^{2}(\tau-1)^{2}}{2\sigma^{2}\tau}\biggr{]}\ ,$ (6)

with variance $\sigma^{2}$. To quantify this correspondence we furthermore
calculate the Kullback-Leibler divergence $D_{\rm KL}=\langle\ln[P_{\rm
num}(\tau)/P(\tau)]\rangle_{P_{\rm num}}$ between the numerically obtained
distribution $P_{\rm num}$ and the prediction of Eq. (6), with ${\langle
t\rangle}$ and $\sigma$ determined from the data. As shown in Fig. 3b, $D_{\rm
KL}$ remains almost zero for both OP and AAP over a wide range of $\gamma$
values, highlighting the robustness of (6). As the inverse Gaussian
corresponds to the First Passage Time distribution of a driven Brownian
particle in one dimension [52], Eq. (6) is closely related to the confinement
of the OP and AAP trajectories along the Zermelo path as shown in Figs. 2d-e.
In fact, for both OP and AAP the loss of correspondence with the inverse
Gaussian coincides with the regime of strong fluctuations that prevent the
swimmers from being efficiently guided along the Zermelo path (details in SI).

In contrast, for sufficiently large flow amplitudes the simpler aligning
policy shows arrival time distributions that do not follow the inverse
Gaussian law (red squares in Figs. 3a,b). These distributions indeed exhibit a
crossover from inverse Gaussian-like behaviour at $\tau<1$ to an exponential
decay at $\tau>1$ with a characteristic time $\tau_{\rm AP}$ that is
systematically larger than the value $2\sigma^{2}/{\langle t\rangle}^{2}$
predicted by Eq. (6). As detailed in the SI, these deviations from the inverse
Gaussian law correspond to asymmetric trajectory distributions around the
Zermelo path, similarly to the case observed for OP and AAP in the large noise
regime. In particular, for the AP case most trajectories land on the left of
the target such that in most cases the swimmer has to navigate counter-flow in
the last part of its journey.

Finally, SP shows arrival time distributions globally compatible with the
inverse Gaussian (6), with only slight deviations at large flow amplitudes
(see blue circles in Figs. 3a,b), indicating that in this case too
trajectories are nearly one dimensional. As SP trajectories are mostly
oriented against the flow, they are characterized by a lower effective drift
on average, resulting in a larger ratio $\sigma/{\langle t\rangle}$.

Figure 5: Optimal navigation on a sphere. (a) Illustration of the chosen
setup. The black arrows indicate the direction of the flow while the solid
green curve connecting the points
$\bm{r}_{T}=(2\pi/3)\hat{\bm{e}}_{\theta}+(3\pi/2)\hat{\bm{e}}_{\phi}$
(magenta circle) and $\bm{r}_{0}=(\pi/6)\hat{\bm{e}}_{\theta}$ (green circle)
corresponds to the Randers geodesic at $\gamma=0.5$. Here the colour map shows
the trajectory distribution for AAP. (b) The arrival time distribution as a
function of $t/{\langle t\rangle}$ for straight (blue circles), adaptive
aligning (green diamonds) and aligning (red squares) policies at $\gamma=0.5$.
(c) Mean arrival time as a function of $\gamma$ for the three policies; the
caption is the same as (b). In (a-c) the value of the Péclet number is set to
${\rm Pe}=10^{3}$. All data in (b,c) are averaged over $10^{5}$ independent
trajectories.

As customary in this context, we compare the performances of the navigation
policies by measuring the mean arrival time to reach the target $\langle
t\rangle$ 111Other measures of performances which take into account additional
features of the arrival time distribution can be defined. As detailed in the
SI, we find that while those allow for more refined evaluation of the
efficiency, they lead to similar results to those obtained from the mean
arrival time ${\langle t\rangle}$., here normalized using the optimal value
$t_{\rm opt}$ at $D=0$ obtained from Zermelo’s solution. The corresponding
results shown in Fig. 3c reveal that, naturally, OP performs the best with the
mean arrival times always remaining higher than $t_{\rm opt}$ by only a few
percent. Conversely, the trivial SP performances strongly deteriorate as
$\gamma$ increases. For small flow amplitudes where the Zermelo path is almost
straight ($\gamma\lesssim 0.4$) SP performs similarly to OP with $\langle
t\rangle\gtrsim t_{\rm opt}$, whereas for sufficiently large $\gamma$ values
it exhibits mean arrival times reaching five to six times $t_{\rm opt}$ (see
the inset of Fig. 3c). Despite the presence of fluctuations, the performances
of the policies are thus primarily set by their ability to make efficient use
of the stationary flow profile. This feature is moreover illustrated by both
AP and AAP, which show mean arrival times no more than $10\%$ higher than that
of OP, regardless of the relative flow amplitude. As expected, in the non-
trivial cases ($\gamma>0.4$) the performances of the different strategies
reflect the amount of information they require for navigation, such that, in
the order of increasing efficiency, one finds SP, AP, AAP and OP.

###  II.4 Robustness of the new protocols

The above analysis shows that AAP displays arrival time statistics similar to
that of OP. Both AP and AAP moreover exhibit performances comparable to OP,
despite them relying only on local information. We now assess the generality
of these results, focusing on AAP, by discussing more general situations with
different model parameters and evaluation setups.

The role of ${\rm Pe}$. The two colour maps of Fig. 4a show how the ratio
$\langle t\rangle/t_{\rm opt}$ varies with the Péclet number and relative flow
amplitude $\gamma$ for OP and AAP. In agreement with previous results,
$\langle t\rangle/t_{\rm opt}$ for OP does not significantly depend on
$\gamma$ while we observe a slight increase with decreasing ${\rm Pe}$. The
AAP case, on the other hand, exhibits two distinct regimes. At small flow
strengths ($\gamma\lesssim 0.4$), $\langle t\rangle/t_{\rm opt}$ remains
nearly constant upon varying ${\rm Pe}$ such that the AAP performances are not
significantly altered by the amplitude of noise. Conversely, at larger
$\gamma$ values where the Zermelo path is more curved, the mean arrival time
is more affected by translational noise. As shown in Fig 4b for $\gamma=0.7$,
all non-trivial strategies show a slight decrease in performance as ${\rm Pe}$
is lowered, whereas SP becomes slightly more favourable upon increasing the
noise, since in this case stronger fluctuations lead the swimmer to visit less
unfavourable flow regions.

Misalignment of self-propulsion. We have so far assumed full control over the
self-propulsion orientation $\bm{\hat{u}}$. In reality, however,
$\bm{\hat{u}}$ is subject to fluctuations—e.g. due to rotational Brownian
motion or inaccuracies in the evaluation of the desired direction—which affect
the performances of the policies using it as a control. To model the effect of
rotational noise, we applied random rotations $\bm{\hat{u}}\to{\cal
R}(\beta)\bm{\hat{u}}$ to the controls (Eqs. (3) and (4)), where the angle
$\beta$ was sampled from a uniform distribution in $(-\eta\pi,\eta\pi]$; we
present the corresponding results in Fig. 4c. For small $\eta$ values, we find
that the inverse of the mean arrival time normalized by its value at $\eta=0$,
${\langle t\rangle}/{\langle t\rangle}_{\eta}$, decays linearly with
$\eta^{2}$ with a policy-dependent slope. In particular, OP and AAP show
similar trends and appear to be much more robust to the effect of rotational
noise than the SP case, as for the latter the mean arrival time has increased
by a factor $10$ at $\eta=\tfrac{1}{2}$, while the corresponding drop in
performance for OP and AAP is about $40\%$. We moreover show in Methods how
the scaling ${\langle t\rangle}/{\langle t\rangle}_{\eta}-1\sim\eta^{2}$ can
be derived from a simple argument in the limit of small noises, and how it is
related to the quasi-one dimensional nature of the problem. This scaling is
thus not expected to hold for large noises, as suggested by the deviations
from the linear decay observed for the largest $\eta$ values in Fig. 4c.

Complex navigation tasks. Increasing the distance between the initial and
target points in the Taylor-Green flow allows us to design more complex paths.
Upon translating the initial swimmer position along $\hat{\bm{e}}_{x}$, as
shown in Figs. 4d and 4e, both OP and AAP lead to trajectories focused around
the Zermelo path on average. Consequently, the corresponding arrival time
distributions remain well characterized by the inverse Gaussian law (Figs.
4f). The performances of AAP and OP are moreover stable upon increasing the
total travel distance, with mean arrival times ${\langle t\rangle}$ not higher
than $t_{\rm opt}$ by more than a few percent.

###  II.5 Optimal navigation on a manifold

We next show how AP and AAP navigation protocols are applicable to motion on
curved landscapes. Self-propelled motion on curved surfaces has recently
earned growing attention both at individual [54, 55, 56, 12] and collective
[57, 58, 59, 60, 61] levels. As stochastic motion taking place on a generic
Riemannian manifold involves multiplicative noise, solving the corresponding
MFPT equation (2) requires advanced computational techniques [62], which will
introduce additional challenges for determining the stochastic optimal control
(3). On the other hand, Zermelo’s approach was recently generalized to self-
propelled motion on curved surfaces using a mapping to Finsler geometry [12].
The corresponding noiseless optimal path, also known as the _Randers geodesic_
, can then straightforwardly be used to extend AP and AAP policies to non-
Euclidean spaces.

For the sake of illustration, let us consider the case of active motion on a
sphere in the presence of a unidirectional flow
${\bm{f}}(\theta,\phi)=v_{f}\sin\theta\,\hat{\bm{e}}_{\phi}$, where $\theta$
and $\phi$ respectively denote the polar and azimuthal angles in the spherical
coordinate system. As shown in Fig. 5a (see the black arrows), this flow—which
is characterized by a pair of vortices at the poles and is maximum at the
equator—generally leads to non-trivial Randers geodesics (solid green line)
between two arbitrary points on the sphere.

Simulating the counterpart of the Langevin equation (1) on the sphere (details
in Methods) at fixed ${\rm Pe}=10^{3}$, we are able to compare the
performances of AP, AAP and SP. The corresponding arrival time distributions
are shown in Fig. 5b. As for the Taylor-Green flow in flat space, we find that
they are all in good agreement with the inverse Gaussian law (6). In the AP
case, this result is probably due to the rather large value of ${\rm Pe}$
chosen for convenience, which allows all policies to exhibit trajectories well
distributed around a one dimensional path. Figure 5c moreover shows that at
small flow strengths all policies perform similarly with $\langle
t\rangle\gtrsim t_{\rm opt}$, while for larger $\gamma$ values leading to more
complex Randers geodesics SP becomes increasingly disadvantageous. On the
contrary, both AP and AAP always remain close to optimality, as they exploit
the information of the noiseless optimal path.

##  III Discussion

We have introduced a class of policies that allow for semi-autonomous optimal
navigation of microswimmers in complex and noisy environments. These policies
rely on the swimmer knowledge of its local environment, such as its relative
position and distance to a desired path, and their performances were found to
improve with the amount of information available to the swimmer. In
particular, extensive numerical simulations reveal that the adaptive aligning
policy (for which the self-propulsion orientation varies with the distance to
the Zermelo path) shows performances comparable to that obtained from the
optimal MFPT control. The two strategies lead to statistically similar
trajectories and nearly identical inverse Gaussian arrival time distributions.
Our analysis moreover shows that the best performing strategies are also the
most robust to environmental changes, such as stronger translational diffusion
and the introduction of rotational noise. Finally, it was shown that the newly
introduced navigation strategies have the additional advantage of being easily
applicable to the problem of optimal navigation on curved surfaces. In an
illustrative example on spherical geometry, the semi-autonomous navigation
strategies were once again found to perform significantly better than the
trivial one consisting of pointing straight to the target.

Although the analysis carried out here focused on the problem of travel time
optimization, the policies we proposed are fully determined by the noiseless
optimal path. Therefore, they are straightforwardly generalizable to a broader
class of optimization problems like energy dissipation or fuel consumption
minimization. Furthermore, even though the new policies show remarkable
performances as compared to OP, some differences persist. In particular, while
trajectories following OP are symmetrically distributed with respect to the
Zermelo path, this is not the case for AP and AAP (Figs. 2(d,e)). Our ongoing
work investigates how further improvement could be reached by designing
policies based on the ability of some swimmers to adapt their swimming
direction according to the local flow field [39, 41].

Lastly, we note that the constraint of imposing a constant sensitivity
(represented by the parameters $\alpha$ and $\varepsilon$) throughout the
swimmer motion might restrict the performances of the policies. While making
these parameters explicitly space dependent would break the semi-autonomous
nature of the policies, determining the functions $\alpha(\bm{r})$ or
$\varepsilon(\bm{r})$ is certainly a simpler problem than calculating the full
optimal control map of the swimmer orientation. Therefore, the framework
presented here could serve as a basis for reinforcement learning based
approaches applied to complex navigation problems [20, 21, 25].

##  IV Methods

The optimal control for an active particle in a stationary flow. The scope of
this section is to provide a derivation of the mean first passage time
equation and the corresponding optimal control (Eqs. (2) and (3) in the main
text). The derivation relies on standard techniques presented in Ref. [50].

First of all, we aim to derive a policy that minimizes the mean arrival time
at the target given an initial position $\bm{x}$. Let us therefore denote
${\cal T}(\bm{x})$ as the mean travel time. Now, we consider the joint
probability $p(\bm{y},t|\bm{x},0)$ that at time $t$ the swimmer is at position
$\bm{y}$, given that it started from $\bm{x}$ at $t=0$. This probability obeys
the backward Fokker-Planck equation (from (1)):

$-\partial_{t}p+(v_{0}{\bm{\hat{u}}}(\theta)+{\bm{f}}({\bm{x}}))\cdot\bm{\nabla}p+D\nabla^{2}p=0\
,$ (7)

where the gradients are taken with respect to $\bm{x}$. Moreover, $p$ is
related to ${\cal T}$ via

${\cal T}(\bm{x})=\int{\rm d}t{\rm d}\bm{y}\,p(\bm{y},t|\bm{x},0),$

such that we obtain

$(v_{0}{\bm{\hat{u}}}(\theta)+{\bm{f}}({\bm{x}}))\cdot\bm{\nabla}{\cal
T}+D\nabla^{2}{\cal T}=-1\ .$ (8)

The optimal choice for the heading direction $\theta$ (our control parameter)
can be obtained by taking the variational derivative of both sides of (8) with
respect to $\theta$ itself, leading to

$\tan\theta=\frac{\partial_{y}{\cal T}}{\partial_{x}{\cal
T}}\iff\bm{\hat{u}}(\theta)=-\frac{\bm{\nabla}{\cal T}}{|\bm{\nabla}{\cal
T}|}\ ,$ (9)

which corresponds to Eq. (3) in the main text. Putting together Eqs. (8) and
(9), we get the MFPT equation of our problem (Eq.(2) in the main text).

Numerical methods. The Brownian dynamics simulations of Eq. (1) have been
performed using an Euler-Mayurama scheme with a time step $dt=10^{-3}$. We
have verified that the selected time step is sufficiently small, such that the
results presented here do not depend on its value. In all our simulations a
given run ends when the active particle is within a distance $\delta
r=0.025\ell$ from the target. We have also checked that the choice of the disk
radius $\delta r$ does not significantly influence the results as long as the
thermal fluctuations length scale is kept relatively small, i.e.
$\sqrt{2Ddt}\ll\delta r\ll\ell$.

Details on the implementation of the protocols. Here we provide details
regarding the implementation of the navigation strategies presented in the
main text.

The MFPT equation (2) was numerically solved using the Finite Elements Method
implemented in the NDSolve routine of Wolfram Mathematica 12.3.1 [63]. The
Optimal policy was then implemented from the corresponding solution (3) by
discretizing the simulation domain on a square grid of step $l=0.01$, and
assigning to each box the optimal control orientation $\bm{\hat{u}}_{\rm
opt}(\bm{r}_{b})$, with $\bm{r}_{b}$ being the position of the centre of the
box. In the stochastic simulations, the swimmer following OP was then aligning
its direction of motion with the orientation associated with its current
position on the grid.

Both AP and AAP rely on the evaluation of the point ${\bm{r}}_{c}$ on the
Zermelo path closest to the particle position $\bm{r}$. For numerical
efficiency, the Zermelo path was thus discretized and the distance between the
particle and the curve was calculated from the positions of the mid-point of
each segments. In all simulations the initial particle orientation was taken
to be equal to the one prescribed by the Zermelo solution.

The scaling of the mean arrival time with rotational noise. Here, we show how
the linear scaling of the mean arrival time ${\langle t\rangle}_{\eta}$ with
the rotational noise strength $\eta^{2}$ shown in Fig. 4c can be understood
from an effective one dimensional model of driven Brownian motion.

Assuming that the particle remains in the vicinity of the mean path and
neglecting effects due to the curvature of the latter, we consider the
following dynamics

$\dot{r}_{\|}=v_{\|}(r_{\|})\cos\theta+\sqrt{2D}\xi_{\|},$ (10)

where the subscript $\|$ stands for quantities projected along the mean path.
The first term on the rhs of (10) thus accounts for the total mean velocity of
the particle along the mean path, which includes the combined effects of flow
and self-propulsion. In general, the angle $\theta$ obeys a nontrivial and
policy dependent dynamics. However, in the limit of small $\eta$ and $D$ where
the particle remains close to the mean path, we approximate $\theta$ as a
Gaussian noise with zero mean and variance $\propto\eta^{2}$. Expanding the
cosine and performing an average over the noises we thus obtain

$\langle\dot{r}_{\|}\rangle\simeq v_{\|}(r_{\|})(1-\kappa\eta^{2}),$

with $\kappa>0$ being a constant that depends on the navigation details, e.g.
the protocol used. In the one-dimensional approximation and assuming that
$v_{\|}$ varies little with $r_{\|}$, the mean travelling time to reach an
absorbing barrier at distance $L$ scales as ${\langle t\rangle}_{\eta}\sim
L/\langle\dot{r}_{\|}\rangle$, such that we obtain

${\langle t\rangle}/{\langle t\rangle}_{\eta}\simeq 1-\kappa\eta^{2}\ ,$

which corresponds to the scaling observed in Fig. 4c.

Langevin simulations on the sphere. To describe the motion of an overdamped
particle on the sphere, the Langevin equation (1) has to be adjusted to take
into account the multiplicative noise induced by the space curvature. Namely,
it is given by

$\dot{\hat{\bm{r}}}=v_{0}\bm{\hat{u}}+\bm{f}(\bm{r})+\sqrt{2D}\hat{\bm{r}}\times\bm{\xi},$
(11)

where $\hat{\bm{r}}\equiv{\bm{r}}/|{\bm{r}}|$, while the noise $\bm{\xi}$
shares the same statistics as in Eq. (1) and is interpreted in the
Stratonovich sense. Furthermore, in contrast with the Taylor-Green flow case
studied in the main text, for this setup the characteristic length scale of
the flow is comparable with the sphere radius: $\ell\sim R$, such that the
Péclet number is here defined as ${\rm Pe}=Rv_{0}/D$.

Lastly, the extension of AP and AAP to the case where motion takes place on a
sphere straightforwardly follows from the presentation in the text, as it only
requires us to generalize the definitions of the distance and relative
direction between two points of interest. On a sphere of radius $R$, the
shortest distance between the points $\bm{r}$ and $\bm{r}_{c}$—also known as
great-circle distance—is defined as

$|\Delta r|\equiv R\arccos\left(\frac{\bm{r}\cdot\bm{r}_{c}}{R^{2}}\right).$
(12)

The direction $\hat{\bm{n}}$ from $\bm{r}$ to $\bm{r}_{c}$ at the point
$\bm{r}$ is likewise defined from the shortest arc linking the two points,
namely

$\hat{\bm{n}}\equiv\frac{(\bm{r}\times\bm{r}_{c})\times\bm{r}}{|(\bm{r}\times\bm{r}_{c})\times\bm{r}|}.$
(13)

The desired heading direction $\bm{\hat{u}}$ for AP and AAP was then obtained
from rotations of $\hat{\bm{n}}$ around the axis set by $\bm{r}$ using
Rodrigues rotation formula [64]:

$\bm{\hat{u}}={\cal G}(\Delta
r)\hat{\bm{n}}\pm(\bm{r}\times\hat{\bm{n}})\sqrt{1-{\cal G}^{2}(\Delta
r)}+\bm{r}(\bm{r}\cdot\hat{\bm{n}})(1-{\cal G}(\Delta r))\ ,$

where the protocol function ${\cal G}(\Delta r)$ is defined by Eq. (4) and the
$\pm$ sign ensures that $\bm{\hat{u}}\cdot\hat{\bm{t}}\geq 0$.

## Acknowledgements

This work has received support from the Max Planck School Matter to Life and
the MaxSynBio Consortium, which are jointly funded by the Federal Ministry of
Education and Research (BMBF) of Germany, and the Max Planck Society.

Figure 6: Examples of the curves leading to the determination of the optimal
sensitivities $\varepsilon$ (AAP, left) and $\alpha$ (AP, right) at ${\rm
Pe}=400$ and $\gamma=0.7$. Blue circles: mean arrival time $\langle t\rangle$
in units of the optimal time in absence of noise $t_{\rm opt}$. Orange
triangles: probability of getting to the target in a time $t<t_{\rm opt}$.
Here, both observables are normalized by their maximum/minimum value. In both
plots the corresponding optimal sensitivity values are indicated by the
vertical dashed lines.

##  Appendix A An alternative policy performance indicator

As customary in the context of optimal navigation, the primary performance
indicator is the mean arrival time at the target ${\langle t\rangle}$. The
corresponding results shown in the main text, reveal that the new simple
protocols introduced exhibit performances actually very close to that of OP
regardless of the relative flow strength. However, as our numerical
simulations give us access to the full arrival time distribution we may get
further insights into the performance of a navigation protocol by considering
additional observables which also take explicitly into account the effect of
fluctuations.

The specific choice we have made here stems from the following remark: owing
to the presence of thermal fluctuations, the active particle may reach the
target in less time than at $D=0$. The frequency of these events can be
quantified by looking at the probability of arriving before the optimal time
$t_{\rm opt}$ in the absence of noise (defined by the Zermelo solution, see
the main text): ${\rm Prob}(t<t_{\rm opt})\equiv\pi_{<}$. The latter is a
measure of how the policies manage to optimize the effect of fluctuations by
maximizing the frequency of small arrival time events. Note that since in
general ${\langle t\rangle}\geq t_{\rm opt}$, $\pi_{<}$ is bounded by
$\tfrac{1}{2}$ from above.

The new navigation policies introduced in the main text both depend on a free
parameter representing the protocol sensitivity. Figure 6 shows the policies
performances measured from ${\langle t\rangle}$ and $\pi_{<}$ as function of
sensitivity for both the adaptive aligning policy (AAP, left) and the aligning
policy (AP, right). Similarly to the mean arrival time, the performance
indicator $\pi_{<}$ allows to obtain a clear optimal sensitivity value for
both AP and AAP. Moreover, the optimal values obtained independently from
${\langle t\rangle}$ and $\pi_{<}$ generally coincide for AAP while some small
differences are observed for AP. In the latter case, choosing either of the
two estimates does not lead to significant variations of the values of
${\langle t\rangle}$ and $\pi_{<}$, and the corresponding arrival time
distributions do not differ significantly.

Figure 7: Probability of reaching the target before the optimal time in
absence of noise ($t_{\rm opt}$) as function of (a) the relative flow strength
$\gamma$ at ${\rm Pe}=400$ and (b) the inverse of the Péclet number at
$\gamma=0.7$ (horizontal axis in log scale). (c) Colour maps showing the
performance of OP (left) and AAP (right) as function of both the relative flow
strength $\gamma$ and the Péclet number ${\rm Pe}$. The colour gradient
represents the probability of reaching the target faster than in absence of
noise from low (blue) to high (red). All data in (a)-(c) were averaged over
$10^{5}$ independent trajectories. Figure 8: (a) Arrival time variance to
square mean ratio as function of the relative flow strength $\gamma$ at ${\rm
Pe}=400$ for the four policies: straight policy (blue circles), optimal policy
(yellow triangles), adaptive aligning policy (green diamonds) and aligning
policy (red squares). The inverted triangles are the values obtained from the
fit of the large-time tails of AP distribution (vertical axis in log scale).
(b) Arrival time probability distributions as function of the normalized time
$t/{\langle t\rangle}$ corresponding to AP for three different values of the
relative flow strength. Solid lines show the exponential fit at large times,
while the dashed lines correspond to the tails predicted by the inverse
Gaussian law (14) using the values ${\langle t\rangle}$ and $\sigma$ from the
data. (c) Heat map of $10^{3}$ stochastic trajectories obtained from numerical
simulations at $\gamma=0.7$ and ${\rm Pe}=400$ of AP. The green curve shows
the Zermelo path connecting the initial point
$\bm{r}_{0}=-\ell/2\hat{\bm{e}}_{x}$ (green circle) and the target
$\bm{r}_{T}=\bm{0}$ (magenta circle). All data in (a) and (b) are averaged
over $10^{5}$ independent trajectories.

Performance and robustness assessment. In Fig. 7a we show the probability
$\pi_{<}$ as a function of the relative flow strength $\gamma$. On the one
hand, the policies performances show trends analogous to those reported in the
main text for the analysis of the mean arrival time, with a preserved
hierarchy at strong flows. As expected, in order of increasing performance we
also find here SP, AP, AAP and OP. On the other hand, the differences between
the various protocols appear more striking. For example, the probability that
a swimmer following AAP reaches the target in a shorter time than $t_{\rm
opt}$ is around $35\%$ lower than that of a swimmer following OP at
$\gamma=1$, while the corresponding mean arrival times deviate by only a few
percent (see main text). As shown in Fig. 7b, similar conclusions can be
reached examining the behaviour of $\pi_{<}$ as function of the Péclet number.
Namely, all OP, AAP and AP show a slight decrease in performance upon
increasing the strength of fluctuations, while SP becomes more advantageous at
large noises. For completeness, we also show in Fig. 7c the heat maps
comparing the values of $\pi_{<}$ for OP and AAP as function of $\gamma$ and
${\rm Pe}$. Here again, the behaviour of the performance indicator is similar
to that of the mean arrival time.

Overall, the new indicator $\pi_{<}$ is largely dominated by the mean arrival
time, although it accounts explicitly for the effect of fluctuations and
provides a more refined evaluation of the protocols performances. $\pi_{<}$
therefore leads to qualitatively analogous conclusions regarding the policies
efficiency, which highlights the robustness of the results presented in the
main text.

##  Appendix B Crossover from Inverse Gaussian to exponential decay for AP

In the main text, we showed that for most flow and noise strengths, the
arrival time distributions of AAP and OP essentially follow inverse Gaussian
laws, whose expression we report here for convenience:

$P(\tau)=\sqrt{\frac{{\langle
t\rangle}^{2}}{2\pi\sigma^{2}\tau^{3}}}\exp\biggr{[}-\frac{{\langle
t\rangle}^{2}(\tau-1)^{2}}{2\sigma^{2}\tau}\biggr{]}\ ,$ (14)

where $\sigma^{2}$ denotes the corresponding variance and $\tau=t/{\langle
t\rangle}$.

However, we observed that while the the arrival time distributions for AP are
well described by the inverse Gaussian law at small $\gamma$, significant
deviations were found when increasing the intensity of the flow. In
particular, these deviations mainly lie in the large-time tail of the
distribution that shows an exponential decay with characteristic time
$\tau_{\rm AP}$ generally larger than the value $2\sigma^{2}/{\langle
t\rangle}^{2}$ predicted by the inverse Gaussian law (14) (see Fig. 8b).
Defining $\tau^{*}\equiv 2\sigma^{2}/{\langle t\rangle}^{2}$ as the variance
to square mean ratio of arrival time for each protocol, its scaling with the
flow strength is shown in Fig. 8a for all policies. As a sign that AAP and OP
on average lead the swimmer to travel faster as the flow strength is
increased, $\tau^{*}$ decays with $\gamma$ for these two policies. In
contrast, swimmers following SP always travel counter-flow and are thus get
slower on average as $\gamma$ increases. They are thus more subject to
fluctuations, such that for SP $\tau^{*}$ grows with $\gamma$. Lastly, for AP
$\tau^{*}$ undergoes a crossover from a decay with flow strength at small
$\gamma$, to a growth with $\gamma$ at large flows (red squares in Fig. 8a).
On the other hand, the value $\tau_{\rm AP}$ obtained from the large-time
tails of the distribution always grows with $\gamma$, similarly to SP (purple
inverted triangles in Fig. 8a). These observations can be rationalized from
the heat map of AP trajectories shown in Fig. 8c. Namely, it shows that most
of the trajectories end at the left of the target, such that the swimmers have
to travel counter-flow and are thus generally slower in the final stretch.
This suggests that, for AP, the events characterized by $\tau>1$ are dominated
by this last part of the swimmer’s journey, ultimately caused by the overall
pronounced asymmetry with respect to the Zermelo path.

Figure 9: (a) Kullback-Leibler divergence $D_{\rm KL}$ between the numerical
and theoretical arrival time distributions as function of the inverse of the
Péclet number at $\gamma=0.7$ (horizontal axis in log scale). (b) Heat maps of
$10^{3}$ stochastic trajectories obtained from numerical simulations of OP
(left), AAP (center) and SP (right) at ${\rm Pe}=40$ and $\gamma=0.7$. The
solid green lines here represent the $D=0$ Zermelo path connecting the initial
point $\bm{r}_{0}=-\ell/2\hat{\bm{e}}_{x}$ and the target $\bm{r}_{T}=\bm{0}$,
marked with a green and magenta circle, respectively. The solid red line is
instead the straight path linking the two points.

##  Appendix C Large fluctuations and the inverse Gaussian law

The correspondence observed between the inverse Gaussian law and the arrival
time distributions of OP, AAP and SP is due to the fact that, as long as the
thermal fluctuations are reasonably small, these policies guide the particle
along a fictional path, thus reducing the system to an active Brownian
particle navigating in a quasi-1D environment.

However, this analogy may break down in the regime of strong fluctuations.
This can be quantified as in the main text computing the Kullback-Leibler
divergence $D_{\rm KL}$ between the distributions obtained from the direct
numerical simulations $P_{\rm num}$ and those predicted by the inverse
Gaussian law (14) with parameters ${\langle t\rangle}$ and $\sigma$ extracted
from the data. More specifically, given a discretization ${t_{i}}$ of the
relevant time interval, the Kullback-Leibler divergence has been calculated as
$D_{\rm KL}=\sum_{i}P_{\rm num}(t_{i})\ln[P_{\rm num}(t_{i})/P(t_{i})]$.
Figure 9a shows the corresponding results obtained for all the policies.
Remarkably, the arrival time distributions of both OP and AAP turn out to be
inverse Gaussian for a wide range of noise amplitudes, with significant
deviations arising only from ${\rm Pe}^{-1}\gtrsim 0.01$. As shown in the heat
maps in Fig. 9b (left and central panels), when fluctuations are strong the
stochastic trajectories are in fact less focused around the optimal path and
also more asymmetrically distributed around it. As already discussed in the
previous section, this asymmetry leads to larger tails in the probability
distributions (data not shown) and therefore to the observed deviations from
the inverse Gaussian law.

In this scenario, the shape of SP distribution turns out to be the most robust
to fluctuations. This can be better understood by looking at the corresponding
heat map in Fig. 9b (rightmost panel). Despite being more dispersed, the
stochastic trajectories still look symmetrically distributed around the
straight path. This is strictly related to the flow symmetry with respect to
the line connecting the starting point with the target (to this end, please
refer to Fig. 1a in the main text).

## References

## References

  * Park  _et al._ [2017] B.-W. Park, J. Zhuang, O. Yasa, and M. Sitti, ACS Nano 11, 8910 (2017).
  * Manzari  _et al._ [2021] M. T. Manzari, Y. Shamay, H. Kiguchi, N. Rosen, M. Scaltriti, and D. A. Heller, Nature Review Materials 6, 351 (2021).
  * Trincavelli  _et al._ [2008] M. Trincavelli, M. Reggente, S. Coradeschi, A. Loutfi, H. Ishida, and A. J. Lilienthal, in  _2008 IEEE/RSJ International Conference on Intelligent Robots and Systems_ (2008) pp. 2210–2215.
  * Szczerba  _et al._ [2000] R. Szczerba, P. Galkowski, I. Glicktein, and N. Ternullo, IEEE Transactions on Aerospace and Electronic Systems 36, 869 (2000).
  * Najafi and Golestanian [2004] A. Najafi and R. Golestanian, Phys. Rev. E 69, 062901 (2004).
  * Dreyfus  _et al._ [2005] R. Dreyfus, J. Baudry, M. L. Roper, M. Fermigier, H. A. Stone, and J. Bibette, Nature 437, 862 (2005).
  * Golestanian  _et al._ [2007] R. Golestanian, T. B. Liverpool, and A. Ajdari, New Journal of Physics 9, 126 (2007).
  * Illien  _et al._ [2017] P. Illien, R. Golestanian, and A. Sen, Chemical Society Reviews 46, 5508 (2017).
  * Zermelo [1931] E. Zermelo, ZAMM - Journal of Applied Mathematics and Mechanics / Zeitschrift für Angewandte Mathematik und Mechanik 11, 114 (1931).
  * Liebchen and Löwen [2019] B. Liebchen and H. Löwen, EPL (Europhysics Letters) 127, 34003 (2019).
  * Techy [2011] L. Techy, Intelligent Service Robotics 4, 271 (2011).
  * Piro  _et al._ [2021] L. Piro, E. Tang, and R. Golestanian, Phys. Rev. Research 3, 023125 (2021).
  * Daddi-Moussa-Ider  _et al._ [2021] A. Daddi-Moussa-Ider, H. Löwen, and B. Liebchen, Commun Phys 4, 10.1038/s42005-021-00522-6 (2021).
  * Pontryagin [1987] L. S. Pontryagin,  _Mathematical Theory of Optimal Processes_ (Routledge, 1987).
  * Bellman [1954] R. Bellman, Bulletin of the American Mathematical Society 60, 503 (1954).
  * Yong and Zhou [1999] J. Yong and X. Y. Zhou,  _Stochastic controls: Hamiltonian systems and HJB equations_, Vol. 43 (Springer Science & Business Media, 1999).
  * Kappen [2005] H. J. Kappen, IOP Publishing 2005, 11011 (2005).
  * Pinti  _et al._ [2020] J. Pinti, A. Celani, U. H. Thygesen, and P. Mariani, Theoretical Ecology 13, 583 (2020).
  * Bertsekas [2012] D. P. Bertsekas,  _Dynamic programming and optimal control - 4th Edition_, Vol. 2 (Athena scientific Belmont, MA, 2012).
  * Colabrese  _et al._ [2017] S. Colabrese, K. Gustavsson, A. Celani, and L. Biferale, Phys. Rev. Lett. 118, 158004 (2017).
  * Schneider and Stark [2019] E. Schneider and H. Stark, EPL (Europhysics Letters) 127, 64003 (2019).
  * Biferale  _et al._ [2019] L. Biferale, F. Bonaccorso, M. Buzzicotti, P. C. D. Leoni, and K. Gustavsson, Chaos 29, 103138 (2019).
  * Yang  _et al._ [2020] Y. Yang, M. A. Bevan, and B. Li, Advanced Intelligent Systems 2, 1900106 (2020).
  * Cichos  _et al._ [2020] K. Cichos, F. nnd Gustavsson, B. Mehlig, and G. Volpe, Nat. Mach. Intell. 2, 94 (2020).
  * Muiños-Landin  _et al._ [2021] S. Muiños-Landin, A. Fischer, V. Holubec, and F. Cichos, Science Robotics 6, eabd9285 (2021).
  * Tierno  _et al._ [2008] P. Tierno, R. Golestanian, I. Pagonabarraga, and F. Sagues, The Journal of Physical Chemistry B 112, 16525 (2008).
  * Bregulla  _et al._ [2014] A. P. Bregulla, H. Yang, and F. Cichos, ACS Nano 8, 6542 (2014).
  * Mano  _et al._ [2017] T. Mano, J.-B. Delfau, J. Iwasawa, and M. Sano, Proc. Natl. Acad. Sci. 114, 2580 (2017).
  * Das  _et al._ [2015] S. Das, A. Garg, A. I. Campbell, J. Howse, A. Sen, D. Velegol, R. Golestanian, and S. J. Ebbens, Nature Communications 6, 10.1038/ncomms9999 (2015).
  * Gomez-Solano  _et al._ [2017] J. R. Gomez-Solano, S. Samin, C. Lozano, P. Ruedas-Batuecas, R. van Roij, and C. Bechinger, Scientific Reports 7, 14891 (2017).
  * Chen  _et al._ [2018] X.-Z. Chen, B. Jang, D. Ahmed, C. Hu, C. De Marco, M. Hoop, F. Mushtaq, B. J. Nelson, and S. Pané, Advanced Materials 30, 1705061 (2018).
  * Fernandez-Rodriguez  _et al._ [2020] M. A. Fernandez-Rodriguez, F. Grillo, L. Alvarez, M. Rathlef, I. Buttinoni, G. Volpe, and L. Isa, Nature Communications 11, 4223 (2020).
  * Tsang  _et al._ [2020] A. C. Tsang, E. Demir, Y. Ding, and O. S. Pak, Advanced Intelligent Systems 2, 1900137 (2020).
  * Bennett and Golestanian [2015] R. R. Bennett and R. Golestanian, Journal of The Royal Society Interface 12, 20141164 (2015).
  * Lozano  _et al._ [2016] C. Lozano, B. ten Hagen, H. Löwen, and C. Bechinger, Nature Communications 7, 12828 (2016).
  * Dai  _et al._ [2016] B. Dai, J. Wang, Z. Xiong, X. Zhan, W. Dai, C.-C. Li, S.-P. Feng, and J. Tang, Nature Nanotechnology 11, 1087 (2016).
  * Hong  _et al._ [2007] Y. Hong, N. M. K. Blackman, N. D. Kopp, A. Sen, and D. Velegol, Phys. Rev. Lett. 99, 178103 (2007).
  * Lagzi  _et al._ [2010] I. Lagzi, S. Soh, P. J. Wesson, K. P. Browne, and B. A. Grzybowski, Journal of the American Chemical Society 132, 1198 (2010).
  * Saha  _et al._ [2014] S. Saha, R. Golestanian, and S. Ramaswamy, Phys. Rev. E 89, 062316 (2014).
  * Kranz  _et al._ [2016] W. T. Kranz, A. Gelimson, K. Zhao, G. C. L. Wong, and R. Golestanian, Phys. Rev. Lett. 117, 038101 (2016).
  * Palacci  _et al._ [2015] J. Palacci, S. Sacanna, A. Abramian, J. Barral, K. Hanson, A. Y. Grosberg, D. J. Pine, and P. M. Chaikin, Science Advances 1, e1400214 (2015), https://www.science.org/doi/pdf/10.1126/sciadv.1400214 .
  * Kantsler  _et al._ [2014] V. Kantsler, J. Dunkel, M. Blayney, and R. E. Goldstein, eLife 3, e02403 (2014).
  * Klumpp and Faivre [2016] S. Klumpp and D. Faivre, The European Physical Journal Special Topics 225, 2173 (2016).
  * Matsunaga  _et al._ [2017] D. Matsunaga, F. Meng, A. Zöttl, R. Golestanian, and J. M. Yeomans, Phys. Rev. Lett. 119, 198002 (2017).
  * Meng  _et al._ [2018] F. Meng, D. Matsunaga, and R. Golestanian, Phys. Rev. Lett. 120, 188101 (2018).
  * ten Hagen  _et al._ [2014] B. ten Hagen, F. Kümmel, R. Wittkowski, D. Takagi, H. Löwen, and C. Bechinger, Nature Communications 5, 4829 (2014).
  * Zheng  _et al._ [2017] J. Zheng, B. Dai, J. Wang, Z. Xiong, Y. Yang, J. Liu, X. Zhan, Z. Wan, and J. Tang, Nature Communications 8, 1438 (2017).
  * Frangipane  _et al._ [2018] G. Frangipane, D. Dell’Arciprete, S. Petracchini, C. Maggi, F. Saglimbeni, S. Bianchi, G. Vizsnyiczai, M. L. Bernardini, and R. Di Leonardo, eLife 7, e36608 (2018).
  * Huang  _et al._ [2019] H.-W. Huang, F. E. Uslu, P. Katsamba, E. Lauga, M. S. Sakar, and B. J. Nelson, Science Advances 5, eaau1532 (2019), https://www.science.org/doi/pdf/10.1126/sciadv.aau1532 .
  * Risken [1996] H. Risken,  _The Fokker-Planck Equation_ (Springer Berlin Heidelberg, Berlin, Heidelberg, 1996).
  * Piro  _et al._ [2022] L. Piro, B. Mahault, and R. Golestanian, in preparation (2022).
  * Folks and Chhikara [1978] J. L. Folks and R. S. Chhikara, Journal of the Royal Statistical Society. Series B (Methodological) 40, 263 (1978).
  * Note [1] Other measures of performances which take into account additional features of the arrival time distribution can be defined. As detailed in the SI, we find that while those allow for more refined evaluation of the efficiency, they lead to similar results to those obtained from the mean arrival time ${\langle t\rangle}$.
  * Fily  _et al._ [2016] Y. Fily, A. Baskaran, and M. F. Hagan, Active particles on curved surfaces (2016), arXiv:1601.00324 [cond-mat.soft] .
  * Apaza and Sandoval [2018] L. Apaza and M. Sandoval, Soft Matter 14, 9928 (2018).
  * Castro-Villarreal and Sevilla [2018] P. Castro-Villarreal and F. J. Sevilla, Phys. Rev. E 97, 052605 (2018).
  * Sknepnek and Henkes [2015] R. Sknepnek and S. Henkes, Phys. Rev. E 91, 022306 (2015).
  * Henkes  _et al._ [2018] S. Henkes, M. C. Marchetti, and R. Sknepnek, Phys. Rev. E 97, 042605 (2018).
  * Lin  _et al._ [2020] S.-Z. Lin, Y. Li, J. Ji, B. Li, and X.-Q. Feng, Soft Matter 16, 2941 (2020).
  * Napoli and Turzi [2020] G. Napoli and S. Turzi, Phys. Rev. E 101, 022701 (2020).
  * Shankar  _et al._ [2021] S. Shankar, A. Souslov, M. J. Bowick, M. C. Marchetti, and V. Vitelli, Topological active matter (2021), arXiv:2010.00364 [cond-mat.soft] .
  * Dziuk and Elliott [2013] G. Dziuk and C. M. Elliott, Acta Numerica 22, 289 (2013).
  * Inc. [2021] W. R. Inc., Mathematica, Version 12.3.1 (2021).
  * Rodrigues [1840] O. Rodrigues, Journal de Mathématiques Pures et Appliquées , 380 (1840).