{'fraction_non_alphanumeric': 0.047822784417678156, 'fraction_numerical': 0.01550171293927979, 'mean_word_length': 4.339347790100977, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Swimmer dynamics in externally-driven fluid flows: The role of noise

We theoretically investigate the effect of random fluctuations on the motion of elongated microswimmers near hydrodynamic transport barriers in externally-driven fluid flows. Focusing on the two-dimensional hyperbolic flow, we consider the effects of translational and rotational diffusion as well as tumbling, i.e. sudden jumps in the swimmer orientation. Regardless of whether diffusion or tumbling are the primary source of fluctuations, we find that noise significantly increases the probability that a swimmer crosses one-way barriers in the flow, which block the swimmer from returning to its initial position. We employ an asymptotic method for calculating the probability density of noisy swimmer trajectories in a given fluid flow, which produces solutions to the time-dependent Fokker-Planck equation in the weak-noise limit. This procedure mirrors the semiclassical approximation in quantum mechanics and similarly involves calculating the least-action paths of a Hamiltonian system derived from the swimmer's Fokker-Planck equation. Using the semiclassical technique, we compute (i) the steady-state orientation distribution of swimmers with rotational diffusion and tumbling and (ii) the probability that a diffusive swimmer crosses a one-way barrier. The semiclassical results compare favorably with Monte Carlo calculations.


I. INTRODUCTION
The advection of self-propelled particles in externally-driven fluid flows presents many surprises when compared to passive advection. Perhaps the biggest surprise is that the transport efficiency of swimmers does not simply increase as swimming speed increases.
For example, when swimmers are placed in a two-dimensional (2D) oscillating vortex array exhibiting chaotic mixing, faster swimming does not always lead to faster transport [1,2].
Even in a steady 2D vortex array, swimmer trapping inside vortices may be enhanced when the particles swim faster, depending on the shape of the particle [3]. Similarly, transport efficiency does not simply increase as a swimmer's rotational diffusivity increases, either. In fact, the opposite occurs in the 2D oscillating vortex array [2,4]. It is reasonable to expect that the more a swimmer's propulsion direction fluctuates, the smaller its net displacement in a fixed amount of time, and hence the lower the transport efficiency. Unexpectedly however, the presence of both rotational noise and shear flow can effectively trap swimmers in certain regions, as has been experimentally observed for swimming bacteria [5] and swimming phytoplankton [6] in a channel flow. While numerous studies have investigated the spatial distributions of noisy swimmers in a variety of flows [7][8][9][10][11], a basic understanding of how rotational noise alters swimmer trajectories is lacking. Our objective in this paper is to develop a theory that quantifies the effect of noise on swimmer dynamics in externally driven fluid flows, especially near transport barriers.
Recently, transport barriers analogous to separatrices-and the related invariant manifoldsof passive advection were identified for self-propelled particles in fluid flows [12]. Perfectly smooth-swimming particles are blocked by so-called swimming invariant manifolds (SwIMs) in position-orientation space. The SwIMs project to one-way barriers, called SwIM edges, to swimmer motion in position space. Swimmers with orientational noise, on the other hand, can cross SwIM edges, but they are still blocked by one-way barriers known as burning invariant manifolds (BIMs), which were originally introduced as barriers for propagating chemical reaction fronts in fluid flows [13,14]. By one-way barriers, we mean that swimmers can pass through a BIM or SwIM edge in one direction, but not the other.
This theory was applied to analyze the experimental trajectories of smooth-swimming and run-and-tumble Bacillus subtilis bacteria in a microfluidic cross-channel featuring a hyperbolic fluid flow, illustrated in Fig. 1a. Whereas the run-and-tumble bacteria exhibit strong rotational noise in the form of sporadic, abrupt changes in swimming direction (Fig. 1c), the smooth-swimming bacteria tend to swim straight in the absence of a flow, with minimal rotational noise (Fig. 1b). The vertical lines in Figs. 1b and 1c are the SwIM edges blocking inward swimming particles, while the horizontal lines are the SwIM edges blocking outward swimming particles. In the hyperbolic flow, the BIMs coincide exactly with the SwIM edges, and hence these curves are one-way barriers for both the smooth swimming and run-andtumble bacteria. Here, all experimentally recorded trajectories are rectified so they appear to enter the flow from above and exit to the right. Therefore, there can be no trajectories to the left of the SwIM edge x = −1, beyond which any trajectory would be swept to the left, as is evident from Figs. 1b and 1c. However, we observe a gap between the left SwIM edge and the measured trajectories of run-and-tumble bacteria in Fig. 1c, compared to the smooth-swimmers in Fig. 1b that can just graze the left SwIM edge before swimming off to the right. Because the gap represents a depletion of the density of trajectories near the SwIM edge relative to the smooth swimmer case, we refer to it as the depletion effect. The depletion effect is caused by the orientation fluctuations of the run-and-tumble bacteria. A run-and-tumble swimmer near this SwIM edge and initially swimming to the right is very likely to tumble and end up crossing the left SwIM edge, forcing it to escape to the left. At the same time, we observe that the run-and-tumble bacteria swim much closer to the lower SwIM edge than their smooth-swimming counterparts. This again is due to tumbling. While smooth swimmers get aligned with the extensional x direction of the flow and thus cannot swim very far below the line y = 0, run-and-tumble swimmers can tumble out of alignment and swim towards the lower SwIM edge. These stark differences between the trajectories of smooth versus run-and-tumble swimmers have motivated the present work.
In this paper, we show how to calculate the probability of particular noisy swimmer trajectories in a given fluid flow, taking the hyperbolic flow as a case study. Our approach focuses on computing solutions to the time-dependent Fokker-Planck equation of a swimmer (alternately known as a master equation or Smoluchowski equation) in the weak-noise limit [15][16][17][18][19]. This differs from traditional approaches to the swimmer Fokker-Planck equation, which are focused on the stationary (time-independent) solution and are in the Eulerian frame-of-reference [5,[8][9][10]. In contrast, we construct a time-dependent swimmer probability density function by following the Lagrangian paths of a swimmer. This procedure is derived from the weak-noise limit in a manner that is nearly identical to the semiclassical approximation in quantum mechanics [20], so we refer to it as the semiclassical approximation to the Fokker-Planck equation. We use this approach to quantify the depletion of noisy swimmers near a BIM, and compare the results of our semiclassical calculations with Monte Carlo calculations, i.e. direct numerical simulations of the swimmer equations of motion.
The paper is organized as follows. In Sec. II, we provide background information on the model for swimmer motion employed here and the semiclassical approximation to the Fokker-Planck equation. In Sec. III, we review the dynamics of a deterministic smooth swimmer in the hyperbolic flow, in particular the role of the SwIMs and BIMs. In Sec. IV, we compute the position-independent orientation distributions of a swimmer in the hyperbolic flow, obtaining results analogous to the orientation distributions of magnetotactic and viscotactic swimmers in external fields [21][22][23], and we apply the semiclassical approximation to calculate the orientation distribution of swimmers with both rotational diffusion and run-and-tumble dynamics. In Sec. V, we compare Monte Carlo and semiclassical calculations of the depletion effect. Concluding remarks are in Sec. VI. In the appendix, we present a complete derivation of the semiclassical approximation to the Fokker-Planck equation.

II. NOISY SWIMMER DYNAMICS
We consider the motion of an ellipsoidal swimmer in two dimensions, with position r = (x, y) and orientationn = (cos θ, sin θ). The stochastic differential equations describing a noisy swimmer in a fluid flow u(r) are [2,4] where ω =ẑ · (∇ × u) is the vorticity,n ⊥ = (− sin θ, cos θ), and E = (∇u + ∇u T )/2 is the symmetric rate-of-strain tensor. The shape parameter α equals (a 2 − 1)/(a 2 + 1), where a is the aspect ratio of the ellipse; α varies from −1 to 1, where α = 0 is a circle, and |α| = 1 is an infinitely thin rod. Positive (negative) values of α correspond to swimming parallel (perpendicular) to the major axis. Each equation of (1) contains a deterministic drift term, proportional to dt, and noise terms. In Eq. (1a), the noise terms are the independent Wiener processes w r = (w x , w y ) and account for translational diffusion with diffusivity D T . Note that for certain swimmers, the strength of the translational diffusivity along the particle's major axis may differ from the translational diffusivity along the minor axis, and their may be additional correlations between translational and rotational noise [24]. We ignore these issues here for simplicity.
Equation (1b), on the other hand, contains two stochastic terms describing fluctuations in the swimming direction. We distinguish between two types of rotational noise: rotational diffusion and tumbling. Rotational diffusion refers to continual random perturbations in the swimmer orientation, such that in the absence of a flow, the orientation θ would exhibit free diffusion (given by the Wiener process w θ ) with a rotational diffusivity D R . This can arise due to random fluctuations in the propulsion force of the swimmer [24,25] or from the thermal fluctuations of the surrounding fluid. In the former case, the noise mechanism would lead to correlations between translational and rotational diffusion, but here we neglect those for simplicity. Tumbling refers to the sudden resetting of θ to a random orientation, independent of its present value, which occurs sporadically as a Poisson process L(ν) with frequency ν. This kind of sudden, random reorientation is seen in swimming bacteria in the "run-and-tumble" mode of swimming. In practice, the distribution of new orientations may depend on the previous value, as is the case for wild-type strains of the swimming bacteria E. coli [26], but we neglect this here for simplicity. Hence, the new angle after a tumble is uniformly and randomly distributed between 0 and 2π.
Our goal in this paper is to estimate the probability of various swimmer trajectories of Eq. (1). This can certainly be accomplished by Monte Carlo simulations, i.e. direct numerical simulations of Eq. (1), but we also develop an analytical approach to calculating such probabilities, which is less computationally costly and provides deeper theoretical insight into the swimmer dynamics. To study the probability of swimmer trajectories, we investigate the Fokker-Planck equation for the probability density of the particle P (r, θ, t) [27], where we have abbreviated the deterministic drift terms from Eq. (1) as F, with Here, ∇ = (∂/∂x, ∂/∂y, ∂/∂θ). We have non-dimensionalized the coordinates using a length scale L and velocity scale U, so that is the strength of rotational diffusion, γ = D T /(L 2 D R ) is the ratio of translational diffusion to rotational diffusion (usually γ < 1), and λ = νL/U is the non-dimensional tumbling rate.
Note that the rotational Péclet number is Pe = 2/ε [9]. The first two terms on the righthand side of Eq. (2) are the usual drift and diffusion terms. The third term proportional to λ accounts for tumbling, with the first term in brackets describing the loss of probability due to tumbling out of the present angle, and the second term describing the gain of probability from the swimmers at all other angles that have tumbled into the present angle. Equation (2) is difficult to attack in general, so we focus on special cases where exact or approximate analytical (or semi-analytical) solutions may be found.
Of particular interest is the λ = 0 case, describing non-tumbling swimmers or, alternatively, the evolution of the probability density in between tumble events. In this case we seek asymptotic solutions in the weak diffusion (ε 1) limit, which have the WKB form quantum mechanics, and hence we refer to it as the semiclassical approximation. Substituting this approximation into Eq. (2) with λ = 0 leads to a Hamilton-Jacobi equation for W and a related transport equation for A.
Here, we briefly describe the mathematical theory of approximation (5), while a detailed derivation and discussion are contained in Appendix A. The solution W to the Hamilton-Jacobi equation (A11) is related to the classical action accumulated along the trajectories of a particular Hamiltonian system associated with Eq. (2). The Hamiltonian is given by where q = (x, y, θ) and p = (p x , p y , p θ ). Equation (6) follows from the Hamiltonian (A12) for a general Fokker-Planck equation (A2), of which Eq. (2) (with λ = 0) is a special case.
The function W is equivalent to the Onsager-Machlup-Freidlin-Wentzell action function [17,28,29] which arises in nonequilibrium statistical mechanics [15,16] and rare event modeling [30]. At each point (r, θ), the action W (r, θ, t) can be expressed as an integral along the trajectory of the Hamiltonian system with Hamiltonian (6) that arrives at that point at time t. This makes Eq. (5) a Lagrangian, as opposed to Eulerian, description of the probability density. The trajectories associated with the minima of W , i.e. the minimumaction paths, correspond to the most likely trajectories of a noisy swimmer, because at lowest order in ε, the probability density (5) is peaked at these points. Hence, finding the minimum-action paths is the main focus of most works involving the Onsager-Machlup-Freidlin-Wentzell action function, including recent work on escape paths of active particles in potential wells [19]. In contrast, we consider all possible paths, in order to get a global approximation to the probability density (5).
Throughout the paper, we focus on the hyperbolic flow u = (Bx, −By). Therefore, We have taken the velocity scale U = v 0 and the length scale L = v 0 /B in the nondimensional equation (7). The typical values of ε, γ, λ, and α depend on the system being considered. In the hyperbolic flow experiments leading to Fig. 1, B = 0.44 s −1 [12].
The measured rotational diffusivity for several species of swimming phytoplankton is D R = 0.15-0.27 rad 2 /s [6], which in the hyperbolic flow would yield ε = 0.68-1.2. For wild-type E. coli, which exhibit run-and-tumble behavior, the rotational diffusivity during runs is D R = 0.06 rad 2 /s [31] and the tumbling frequency is approximately ν = 1 s −1 [26], which in the hyperbolic flow would yield ε = 0.27 and λ = 2.3. Assuming that the translational diffusion for wild-type E. coli is due only to thermal fluctuations so that D T = 0.2 µm 2 /s Before investigating the dynamics of noisy swimmers in the hyperbolic flow, we study the deterministic dynamics of Eq. (7) with ε = 0 and λ = 0.

III. DETERMINISTIC DYNAMICS IN THE HYPERBOLIC FLOW
The deterministic dynamics of Eq. (7) is best understood through the system's fixed points and invariant manifolds, previously studied in Ref. [12]. The system possesses four fixed points, which we refer to as swimming fixed points (SFPs) to distinguish them from the passive fixed points of the fluid flow. Denoting the phase-space coordinate q = (x, y, θ), the fixed points are illustrated in Fig. 2. Each of the SFPs is a saddle. When α > 0, and in particular when α = 1, the q out ± fixed points have stable-unstable-unstable (SUU) linear stability, and the q in ± fixed points have stable-stable-unstable (SSU) linear stability. Hence, the q out ± SFPs possess 2D unstable manifolds, while the q in ± fixed points possess 2D stable manifolds. We refer to these 2D manifolds as swimming invariant manifolds (SwIMs), to distinguish them from the invariant manifolds of passive advection [12].
Taken as a whole, the stable and unstable SwIMs consist of two interlocking S-shaped sheets, plotted in Fig. 2a. The stable SwIMs (the blue surface) attached to q in + and q in − share common boundaries along the lines {(x, y, θ) | x = 0, θ = ±π} which are the 1D stable manifolds of q out ± (dark blue lines). Hence, the union of the 2D stable SwIMs with the 1D stable manifolds of q out ± is a surface (blue surface in Fig. 2a) which separates phase space into two pieces. By symmetry, the same geometric shape can be constructed by taking the union of the unstable SwIMs of q out ± with the 1D unstable manifolds of q in ± , leading to the red surface in Fig. 2a. The shape of the stable SwIMs is independent of the y coordinate and, similarly, the shape of the unstable SwIMs is independent of the x coordinate. This occurs because in the hyperbolic flow, the xθ equations Eq. (7a) and (7c) are decoupled from y and similarly the yθ equations Eq. (7b) and (7c) are decoupled from x. The stable and unstable SwIMs intersect along heteroclinic orbits going from one fixed point to another, indicated by the yellow curves in Fig. 2a.
Cross-sections of the SwIMs are shown in Fig. 3, along with the phase portraits of the xθ dynamics and the yθ dynamics. Figure 3a shows that swimmers on the left of the stable SwIM ultimately exit the flow to the left, while swimmers on the right ultimately exit right.
Similarly, the unstable SwIM (Fig. 3b) separates swimmers that entered the flow from the top from those which entered from the bottom. The SwIMs are therefore transport barriers to swimmers in the hyperbolic flow, because they carve out the xyθ phase space into distinct, qualitatively different families of trajectories. Importantly, these barriers are nonporous in phase space, meaning no swimmer trajectory can cross them (in the absence of noise).
On the other hand, in position space, the SwIMs project to one-way barriers, allowing swimmers to cross in one direction but not the other. Figure 2b shows the singularities of the projection of the SwIMs into the xy plane-that is, the folds of the S-shapes which bound the projection of the 2D surfaces into the plane. We refer to these curves as SwIM edges [12]. The stable SwIM edges at x = ±1 (blue curves in Fig. 2b) block inward swimming particles, while allowing outward swimming particles through. To see this, note that for x = −1,ẋ ≤ 0 for all θ, and for x = 1,ẋ ≥ 0 for all θ, as shown in Fig. 3a. Along the stable SwIM edges, the outward fluid flow overpowers the swimmers and they are swept away from the center of the flow. Similarly, the unstable SwIM edges (red curves in Fig. 2b) block outward swimming particles, while inward swimming particles can pass through them.
Here, for y = 1,ẏ ≤ 0, and for y = −1,ẏ ≥ 0. On the unstable SwIM edges, it is the inward flow which overpowers the swimmers and pushes them towards the center of the flow.
The SwIM edges in the hyperbolic flow coincide exactly with the BIMs-the 1D invariant manifolds of the SFPs when α = −1 [12,32]. This is important because SwIM edges are only guaranteed to be one-way barriers for purely deterministic swimmers. BIMs, on the other hand, have stronger barrier properties, in that they are also one-way barriers for swimmers with rotational diffusion or run-and-tumble dynamics in the limit of negligible translational diffusion [12]. Thus, in the hyperbolic flow, the SwIM edges also act as one-way barriers for swimmers with rotational noise. This explains why the run-and-tumble bacteria in the hyperbolic flow experiment remain bounded by the unstable SwIM edge at y = −1 in This is the origin of the depletion effect we observe when comparing the smooth swimming bacteria data (Fig. 1a) to the run-and-tumble data (Fig. 1b). The orientation fluctuations of tumbling bacteria make it very likely that a bacterium near the SwIM edge at x = −1, for example, swims across it, precluding the possibility that it subsequently exits the flow to for the probability density P (θ, t) restricted to the orientation degree-of-freedom is We focus on the stationary distributions of the θ variable, which are the stationary solutions (∂P/∂t = 0) of Eq. (9). We first treat the two limiting cases (i) no tumbling (λ = 0), and (ii) no rotational diffusion (ε = 0), before proceeding to the case where there is both tumbling and rotational diffusion. Note that the orientation dynamics of a noisy swimmer in the hyperbolic flow is very similar to the orientation dynamics of swimmers in other types external fields, such as magnetotactic swimmers in external magnetic fields [21,22] and swimmers in viscocity gradients [23]. The main difference here compared to the preceding examples, aside from the source of the torque on the swimmer, is that Eq. (9) is invariant under the symmetry θ → θ + π.
A. Orientation dynamics with rotational diffusion only (λ = 0) Without tumbling, the θ dynamics is governed by The deterministic part of the equation has the form of the gradient of a potential V (θ), dynamics is equivalent to that of an overdamped particle in the potential V (θ) with noisy driving. In this case, the long-time probability distribution of θ evolves towards a stationary state which is peaked at the potential wells at θ = 0 and θ = π (for α > 0). This probability distribution P ε (θ) can be found by solving for the stationary state of Eq. (9) with λ = 0. For gradient systems, the solution is simply P ε (θ) ∝ exp[−2V (θ)/ε], which is simple to verify, and hence we have Clearly, the distribution depends on a single dimensionless parameter, which is the ratio of the rate of alignment with the extensional direction of the flow, Aα, to the intensity of the noise. Normalizing the probability distribution, we obtain where I 0 (x) is a modified Bessel function of the first kind. The stationary distribution (13) is invariant under the shift symmetry θ → θ + π, as is the underlying stochastic process (10). Equation (13) is plotted in Figs. 4a and 4b, along with histograms from Monte Carlo simulations of Eq. (10). In Fig. 4, we map the Monte Carlo data onto the interval θ ∈ (−π/2, π/2) using symmetry and only plot Eq. (13) in this range. Similar results were previously obtained for magnetotactic swimmers in an external magnetic field with rotational diffusion [22]. (c) ε = 0, λ = 1.6. (d) ε = 0, λ = 5.

B. Orientation dynamics with tumbling only (ε = 0)
Here we consider the case of the stationary θ distribution under tumbling only. Every time a swimmer tumbles, its orientation is drawn from the uniform distribution. If it tumbles at time τ = 0, then until the next tumble, its probability density P (θ, τ ) evolves according to the Liouville equation with the initial condition P (θ, 0) = P 0 (θ) = 1/2π. Intuitively, we thus expect that the steady-state distribution under tumbling only, P λ (θ), should consist of the superposition of probability distributions P (θ, τ ) describing the relaxation of θ in between tumbles, weighted by the probability λe −λτ dτ that the last tumble occurred a time τ in the past. In other words, the stationary probability distribution must be [33] It is straightforward to verify that Eq. (15) is a stationary solution of Eq. (9) with ε = 0.
An explicit solution to Eq. (14) can be obtained using the method of characteristics, because the θ equation of motion (10) in the absence of noise (ε = 0) can be solved analytically.
The expression for the deterministic trajectory θ * (t) is where θ 0 = θ(0) is the initial condition, and here it is assumed θ 0 ∈ [−π/2, π/2]. This condition arises due to the use of the tan −1 function; when θ 0 is outside this range, this solution needs to be shifted either up or down by π, depending on θ 0 . Then, the solution to Eq. (14) is where P 0 (θ) = P (θ, 0) is an arbitrary initial orientation distribution (see Appendix B for the derivation). Again, this form of the solution is valid for θ ∈ [−π/2, π/2], and a shift by π in the argument of P 0 in Eq. (17) adapts the solution to the excluded range of θ.
Next, we obtain the stationary θ distribution P λ (θ) under tumbling with rate λ by substituting Eq. (17) into Eq. (15), with P 0 = 1/2π. Rescaling the time in Eq. (15) by the tumbling rate λ, we obtain a complicated integral that depends on a single dimensionless parameter that we call the tumbling number Tu, This is essentially the ratio of the tumbling rate to the relaxation rate of a swimmer's orientation to its equilibrium (parallel to the extensional x-direction) in the hyperbolic flow. Note, the latter relaxation rate is distinct from the relaxation rate of the orientation distribution of a swimmer with rotational diffusion to the stationary state given by Eq. (13).
It is conceivable that this is the more relevant time scale for defining Tu in the case where we have both tumbling and rotational diffusion. The stationary distribution P λ can be shown to be equal to where 2 F 1 (a, b; c; z) is the ordinary hypergeometric function. Like P ε (θ) (Eq. (13)), Eq. (19) is invariant under the shift symmetry θ → θ + π. To be sure, the hypergeometric function makes the expression (19) of the tumbling swimmer's stationary distribution more opaque than its counterpart for rotational diffusion.
In Figs. 4c and 4d, we plot Eq. (19), superimposed over histograms from numerical simulations of tumbling swimmers in the hyperbolic flow, to obtain a basic intuition for how the distribution depends on the parameters. We see an excellent agreement between the theory and simulations. For sufficiently small tumbling rates such that Tu ≤ 1 in Eq. (18) ( Fig. 4c), it can be shown that P λ is singular at the orientation equilibrium θ = 0 and relatively flat for all other orientations. On the other hand, for Tu > 1 (Fig. 4d), the peak at the equilibrium becomes finite, and the difference between the probabilities near the stable and unstable equilbria becomes more modest. Again, our results mirror those obtained for magnetotactic run-and-tumble bacteria in Ref. [22].

C. Orientation distribution with rotational diffusion and tumbling
Having treated the limiting cases of one type of noise versus another, we now seek the stationary θ distribution of a run-and-tumble swimmer with rotational diffusion, which we denote P ε λ (θ). This requires slightly modifying the approach of Sec. IV B, where the stationary distribution is the weighted time average of the probability distributions describing relaxation to equilibrium, P (θ, τ ) [see Eq. (15)]. Namely, when we have both tumbling and rotational diffusion, we need to obtain P (θ, τ ) by solving the time-dependent Fokker-Planck with initial condition P 0 (θ) = 1/2π, instead of the Liouville equation (14). An exact analytical solution of Eq. (20) is unavailable, so we resort to a short-time, small ε asymptotic approximation based on the semiclassical techniques detailed in Sec. A 5.
1. Approximation of the Fokker-Planck propagator for small ε P (θ, τ ) is the probability density of reaching θ at time τ under Eq. (10) with a random initial condition drawn from a uniform distribution. We approximate P (θ, τ ) by making use of the semiclassical approximation to the Fokker-Planck equation (20). We note that this distribution can be expressed as an integral over the propagator K(θ, θ 0 , τ ), which is the probability distribution of reaching θ in time τ from a fixed initial condition θ 0 , as Our strategy is to approximate P (θ, τ ) by approximating K analytically and performing the integral (21) numerically.
The propagator K can be approximated in the small ε limit using the semiclassical approach described in Appendix A. In particular, we are satisfied with an approximation valid for short times, because the long-time behavior of P (θ, τ ) is suppressed in the integral (15) for the steady-state distribution. Therefore, we make use of the Gaussian approximation of K about a deterministic trajectory, derived for a general 1D Fokker-Planck equation of the form (20) in Sec. A 5. In this case, K is peaked around the trajectory θ * (τ ) initiated at θ 0 , given by Eq. (16). The final expression for K follows from Eqs. (A64) and (A65), which after lengthy but straightforward calculations yields where We illustrate the validity of the approximate probability distribution (22) by comparing the prediction to numerical simulations of Eq. (10). One comparison is shown in Fig. 5, where at t = 0 we initialized the swimmers with the orientation θ 0 = 1.07, not terribly far from the turning point θ = π/2 ≈ 1.57, with a modest noise strength of ε = 0.25. These parameters were selected to push the limits of our approximations; not only do we assume ε is small, but our calculation of K(θ, θ 0 , τ ) neglects entirely the contributions of paths that cross the turning point at θ = π/2 ≈ 1.57 and relax to the equilibrium at θ = π instead of θ = 0. Despite these limitations, we see that the approximate K(θ, θ 0 , τ ) given by Eq. (22) does a good job both of tracking the center of the distribution of trajectories and accounting for their spread as a function of time. After some time, the variance of the approximate distribution (∂ 2 R/∂θ 2 ) −1 saturates and the centroid converges onto θ = 0, yielding a steadystate. This distribution is reasonably close to the numerical one at t = 4. However, we know that if we continue the numerical simulations for very long times, then eventually the distribution should approach the exact steady-state given by Eq. (13) (Figs. 4a and 4b). In contrast to our approximate distribution, which only has one peak that eventually converges to θ = 0, the true steady-state distribution is symmetrically peaked about θ = 0 and θ = π.
Over long times, this is achieved as the noise drives some swimmers' orientations over the potential barriers at θ = ±π/2, causing them to settle down around θ = π for long times.
This process is reflected by the growing peak in the density of simulated trajectories at θ = π in Fig. 5c and 5d. Our approximate distribution manifestly neglects this process, because the action associated with such trajectories is larger than for trajectories near the deterministic path, which are the only trajectories accounted for in our approximation.

Approximation of the stationary state
We can now compute the stationary orientation distributions of swimmers with both tumbling and rotational diffusion. To recap, we have an explicit approximation (22) of K(θ, θ 0 , τ ), the time-dependent probability distribution of θ for a rotationally-diffusing swimmer with orientation θ 0 at time t = 0 (which also requires Eqs. (23) and (16) to evaluate).
Thus, we are able to numerically evaluate our expression (21) for the time-dependent probability distribution P (θ, τ ) of θ for a rotationally-diffusing swimmer with an initially uniform orientation distribution. This initial state corresponds to the swimmer's orientation distribution after a tumble. Therefore, we can finally evaluate Eq. (15) for P * λ (θ), the stationary θ distribution of a tumbling and rotationally-diffusing swimmer.
We proceed by evaluating Eqs. (21) and (15) numerically, and we compare the results  Fig. 6, we conclude rotational diffusion smooths out these peaks. We observe good agreement between the stochastic simulations and the semiclassical theory in all cases. Thus, our semiclassical method for evaluating P ε λ (θ) can in principle be used to fit experimental data, allowing the determination of the effective rotational diffusivity and tumbling rate of swimmers in the hyperbolic flow.

V. DEPLETION EFFECT
Here, we present Monte Carlo and semiclassical calculations quantifying the depletion effect. We quantify the depletion effect by calculating the probability Pr(x 0 ) that a swimmer ultimately exits right with a given initial position −1 < x 0 < 1 and a given intensity of the noise. For an x 0 near the BIM x = −1, the signature of the depletion effect is a decreasing Pr(x 0 ) for increasing noise intensity. A low probability of right-exiting swimmer trajectories initialized near x = −1 would be consistent with the absence of such trajectories for runand-tumble bacteria in the experimental data shown in Fig. 1. Conversely, for an x 0 near the BIM x = 1, Pr(x 0 ) should increase with increasing noise intensity. This is simply due to the symmetry of the hyperbolic flow, which requires that We focus solely on the dynamics in the xθ plane, because it is independent of the y variable, as discussed in Sec. III. Hence, Eq. (2) becomes where is the drift restricted to the xθ plane. In Eq. (25) and throughout this section, we also take ∇ = (∂/∂x, ∂/∂θ).
We restrict our attention to the case where rotational diffusion dominates translational diffusion, i.e. γ 1, and we fix γ = 0.1. When γ = 0, all swimmers which cross the line x = 1 must ultimately exit right, due to the BIM at x = 1 blocking inward swimming particles. Therefore, the probability to exit right may be calculated by integrating the probability current through x = 1. We assume this remains approximately true for small γ.
Defining Pr(x 0 , t) as the probability that a swimmer has exited right by time t, we have Differentiating Eq. (27) with respect to time and using Eq. (25), we obtain the probability where is the probability current density excluding tumbling, with The tumbling contribution in Eq. (28) vanishes upon integration over θ. Using the divergence theorem, the probability current (28) becomes and thus the probability to exit right is given by  Figure 7 shows the results for Pr(x 0 ) for swimmers with diffusion only (λ = 0, Fig. 7a) and swimmers with tumbling only (ε = 0, Fig. 7b). For the λ = 0 swimmers, θ 0 was drawn from the stationary distribution P ε (θ 0 ) given by Eq. (13). For the ε = 0 swimmers, θ 0 was drawn from the stationary distribution P λ (θ 0 ) given by Eq. (19). We also show Pr(x 0 ) for deterministic swimmers (ε = λ = 0) initialized with a uniform distribution of θ 0 in Fig. 7c. Here, Pr(x 0 ) is obtained by calculating the fraction of trajectories on the right side of the SwIM at a given x 0 (see Fig. 3a). Figure 7 shows that as the intensity of noise increases, Pr(x 0 ) increases for x 0 > 0 and decreases for x 0 < 0. This occurs both for swimmers with diffusion only (Fig. 7a) and for swimmers with tumbling only (Fig. 7b), where the intensity of noise effectively increases when the tumbling frequency λ increases. The reduction of Pr(x 0 ) for x 0 < 0 for noisier swimmers is consistent with the depletion effect observed in the experimental data shown in Fig. 1. For smooth swimming bacteria (Fig. 1a), which behave like swimmers with weak diffusion, the exit-right probability Pr(x 0 ) is substantial for most values of x 0 , even those relatively close to the BIM at x = −1 (Fig. 7a). Therefore, it is not unlikely to observe bacteria trajectories which graze the BIM at x = −1, as we indeed see in Fig. 1a. For run-and-tumble bacteria on the other hand (Fig. 1b), Pr(x 0 ) is very small near x 0 = −1 for sufficiently large λ (Fig. 7b). Therefore, it is very unlikely to observe bacteria trajectories

B. Semiclassical approximation for diffusion
We use the semiclassical approximation to compute Pr(x 0 ) when λ = 0 and investigate how accurately it matches the Monte Carlo calculations. For the xθ dynamics in the hyperbolic flow, Hamiltonian (6) becomes We evaluate Eq. (32) for Pr(x 0 ) using our semiclassical approximation for P (x, θ, t). This essentially requires integrating over a subset of trajectories of Eq. (33), which begin at x 0 at t = 0 and hit x = 1 at a later time. One advantage of the semiclassical approximation is that this set of trajectories is independent of ε, so once these trajectories are computed, Eq. (32) can be evaluated for any arbitrary value of ε. Another advantage is that this set of trajectories provides insight into the actual paths in xθ space that noisy swimmers take on their way to exiting right.
To illustrate the relationship between the trajectories of Hamiltonian (33) and the noisy trajectories, we first consider the semiclassical evolution of a probability density initially concentrated at a single point (x 0 , θ 0 ). This corresponds to an initial condition which is the initial condition for the propagator of the Fokker-Planck equation (Eq. (A6)).
The semiclassical solution for such an initial condition (Eq. (A39)) requires one to integrate all trajectories of Hamiltonian (33) beginning at (x 0 , θ 0 ), which means considering all possible initial momenta (p x0 , p θ0 ) at that point (Appendix A 1). This 2D surface of initial conditions of the Hamiltonian system is called a Lagrangian manifold [20]. Along the way, one keeps track of the accumulated action R(x, θ, x 0 , θ 0 , t) along each trajectory (Eqs. (A16) and(A19)).
For Hamiltonian (33), the accumulated action is where the integral is along the trajectory connecting (x 0 , θ 0 ) to (x, θ) in time t. In the case of the propagator initial condition, the function W of the semiclassical probability density (5) is simply equal to the accumulated action, W (x, θ, t) = R(x, θ, x 0 , θ 0 , t). The exponential dependence of the semiclassical probability density on W makes the probability density peaked around the local minima and valleys of W . The Hamiltonian trajectories reaching these minima or valleys can be thought of as prototypical noisy trajectories.
For example, we consider swimmers exiting right from (x 0 , θ 0 ) = (0.5, 2), as shown in  Fig. 8. We selected one trajectory that hits x = 1 at t ≈ 1, and a second trajectory that hits at t ≈ 2; aside from these prescribed hitting times, the trajectories were selected at random. We can calculate the trajectories of the system with Hamiltonian (33) which hit x = 1 at those same times. There are infinitely many, each hitting with a different final θ. Out of this set of trajectories, we find the ones which minimize the action at x = 1, equivalent to the condition ∂W (1, θ, t) where the last equality follows from Eq. (A13). In other words, for a given t, the minimumaction trajectory is the one which hits x = 1 with p θ (t) = 0. Equation (36) is the condition for a valley of W (x, θ, t) because it is a local minimum of W with x (and t) held fixed. The minimum-action trajectories corresponding to the two noisy trajectories are plotted as the dotted curves in Fig. 8.
The resemblance between the noisy paths and the minimum-action paths in Fig. 8  at t = 0 to x = 1 at specified times t. In the asymptotic ε → 0 limit, the probability density becomes increasingly concentrated along the minimum-action paths. However, for any finite ε, the probability density has a finite width around these minimum-action paths, so any actual noisy trajectory will deviate from the minimum-action path, as seen in Fig. 8.  Table I).
Next, we turn to the semiclassical calculation of Pr(x 0 ), given that the swimmer's initial orientation θ 0 is distributed according to Eq. (13) as in Sec. V A. This requires the solution of the Fokker-Planck equation (A2) for P (x, θ, t) with initial condition This is a hybrid propagator-WKB initial condition of the form (A9), where, in the notation of Appendix A, A 0 = 2πI 0 α ε −1 and U (θ) = −α cos 2θ. We use the semiclassical probability density (A52) to evaluate Eq. (32). This means that for each x 0 , we must integrate over all  Hamiltonian paths which hit x = 1, with initial conditions on the Lagrangian manifold Therefore, the initial Lagrangian manifold may be parametrized by the coordinates z = (p x0 , θ 0 ). The probability current integral Eq. (32) becomes where A is given by Eq. (A50) and R is given by Eq. (35). Equation (39) Table I.
This converts Eq. (39) into an initial value representation [34,35]. We also neglect the ∂(ln A)/∂x term at the end of the first line of Eq. (39), because it is of order ε relative to the other terms. This means it is of higher order in ε than we account for in our asymptotic expression Eq. (5) (see also Eq. (A3)), and thus it may be neglected within the framework of the semiclassical approximation. We must also truncate the range of p x0 for numerical evaluation, so we take |p x0 | < p max . Lastly, Eq. (39) is even in θ 0 by symmetry, so we can restrict the domain θ 0 ∈ [0, π] and double the result. Hence, the final expression that we evaluate numerically is The domain for the integral (41) is the set of initial conditions (p x0 , θ 0 ) that eventually exit right, i.e. those initial conditions that reach x = 1 at some time t. Figure 9 shows an example integration domain for x 0 = 0.5.
We evaluate Eq. (41) numerically using the trapezoidal rule. For each x 0 , we discretize the set of initial conditions (p x0 , θ 0 ) on the Lagrangian manifold (38) with a uniform grid.
We numerically integrate each trajectory until it hits x = 1, up to a maximum integration time of t = 6, consistent with the corresponding Monte Carlo calculations in Fig. 7a. We x 0 that exit right, including those which begin to the right side of the SwIM (Fig. 3a) and would have exited right even without noise, as well as those which begin on the left side of the SwIM and cross it due to fluctuations (Fig. 8).
As |x 0 | gets closer to 0, however, the semiclassical predictions begin to deviate from the Monte Carlo calculations with increasing ε. In particular, Eq. (24) requires that Pr(0) = 0.5, that is, a swimmer starting in the middle of the flow has an equal probability of going left or right. While the semiclassical prediction appears consistent with this property for ε = 0.1, as ε increases further, we see the semiclassical Pr(0) increase above 0.5 in the inset of Fig. 10.
This fact, combined with our use of Eq. (24) to obtain the semiclassical Pr(x 0 ) for x 0 < 0, causes the apparent kink at x 0 = 0 in our semiclassical predictions plotted in Fig. 10.
We believe that at least part of the discrepancy between the semiclassical Pr(x 0 ) and the Monte Carlo calculations for x 0 close to 0 is due to the presence of caustics, a technical issue that we have ignored until now. Our semiclassical approximation assumes the uniqueness of Hamiltonian trajectories that originate on the initial Lagrangian manifold and go from (x 0 , θ 0 ) to (x, θ) in time t. This means, for example, that when evaluating the accumulated action R(x, θ, x 0 , θ 0 , t) in Eq. (35), there is a unique such trajectory. The uniqueness is the critical property that makes R a well-defined function. However, uniqueness is guaranteed only for sufficiently short times, meaning that once t is sufficiently large, there will be multiple Hamiltonian trajectories connecting the two points. In this case, R becomes multi-valued [36]. Geometrically, the uniqueness breaks down when the evolving Lagrangian manifold develops fold singularities, such that when it is projected into q space, parts of the projection overlap with each other. These overlap regions are the regions where multiple Hamiltonian trajectories can arrive at a single point. In the context of the semiclassical approximation in quantum mechanics, these fold singularities are called caustics, and they have been extensively studied (see [20] and references therein). The formulas for the semiclassical wave function need to be corrected to account for the occurrence of caustics through the inclusion of a Maslov index, a phase factor that essentially counts the number of caustics encountered by each trajectory.
We know of no general prescription for dealing with caustics in the semiclassical approximation to the Fokker-Planck equation, even though they commonly occur. Previously, caustics have been investigated in the semiclassical formulation of the noise-driven dynamics of a nonlinear oscillator in two dimensions [37]. When calculating the steady-state probability density of this system, one must account for switching lines in phase space, i.e. curves on either side of which the least-action path changes discontinuously due to the presence of multiple Hamiltonian trajectories arriving at those locations. Near the switching line, where the distinct paths with coincident endpoints have nearly the same action, the semiclassical probability density needs to account for each of the distinct paths; in fact, the signature of multiple paths leading to the same point in phase space has been observed experimentally in a noisy electronic oscillator [16]. Caustics also arise in the theoretical description of noise-induced transitions in non-gradient dynamical systems with metastable fixed points [38,39]. Specifically, the quasi-stationary probability densities underlying escape from the metastable fixed points may be approximated semiclassically, though one must go beyond the standard WKB approximation in the vicinity of the caustics. It is not obvious how to generalize these previously obtained results to the noisy dynamics of the swimmer in the hyperbolic flow, because the swimmer phase space does not possess stable fixed points and is fundamentally transient, i.e. time-dependent.
Such a generalization is needed because caustics do indeed occur in the dynamics of noisy swimmers in the hyperbolic flow. To demonstrate this, we track the number of caustics encountered by the trajectories underlying our semiclassical calculation of Pr(x 0 ) up until the point that they hit x = 1. At points where caustics occur, the projection from the Largrangian manifold into q space is not invertible, meaning the Jacobian determinant det ∂q/∂z must be zero at that point. Thus, we count the number of caustics encountered along a trajectory by tracking zero-crossings of det ∂q/∂z . In Fig. 11, we plot the fraction of trajectories used in our numerical semiclassical approximation that encounter at least one caustic on the way to x = 1. While the fraction of caustic-crossing trajectories is around 6% or smaller for x 0 > 0.5, it rapidly rises to nearly 15% as x 0 decreases from 0.5 to 0. This trend is reasonable, because caustics only begin to occur after a sufficiently long time. Trajectories beginning closer to x = 1 will tend to reach it sooner, potentially before many caustics have occurred. At the same time, in the x 0 > 0.5 range, we see good agreement between the semiclassical and Monte Carlo calculations in Fig. 10, while in the x 0 < 0.5 range we observe a discrepancy with increasing ε. This correlation is evidence that the semiclassical approximation works well when few trajectories have encountered caustics, while a discrepancy can be caused by improper treatment of the caustics, which is important when considering sufficiently long-time processes.

VI. CONCLUSION
To summarize, we have quantified the effect of noise on swimmer dynamics in a steady, two-dimensional hyperbolic fluid flow. In such a flow, swimmers are ultimately forced to escape to the left or the right, with their transient dynamics near the passive unstable fixed point determining which way they go. Without noise, a swimmer's fate is sealed based on its position relative to the SwIM in the xθ phase space. With noise, the swimmer's motion is a stochastic process. We calculated the steady-state orientation distributions of diffusive, run-and-tumble, or mixed swimmers in the hyperbolic flow. The fluctuations give some swimmers greater opportunity to cross the SwIM and exit on the opposite side than they would have without noise. There is however a maximal distance that swimmers can get on either side of the passive fixed point and still be able to swim back to the other side-this is where the stable BIMs block inward swimming particles.
Fluctuations make it increasingly likely that a swimmer close to one of these BIMs does indeed end up crossing it, causing irreversible changes to the fluctuating swimmers' trajectories (assuming negligible translational diffusion). We quantified this probability using Monte Carlo calculations and a semiclassical approximation to the swimmer Fokker-Planck equation. The semiclassical approximation accurately predicts the probability Pr(x 0 ) a swimmer exits right given that it began at a position x 0 relative to the passive fixed point, especially for x 0 close to the BIM. It also predicts the probability distribution of paths that fluctuating swimmers take in phase space given specified boundary conditions. SwIMs and BIMs are present in nonlinear flows as well, such as alternating vortex flows [12]. Thus, we expect the depletion effect to occur in the vicinity of the BIMs of such flows as well.
This study demonstrates the utility of the semiclassical approximation for understanding the noisy dynamics of a non-trivial active matter system. However, it also reveals a key shortcoming of the existing semiclassical theory for Fokker-Planck dynamics. In particular, a general approach is needed for taking into account the occurence of caustics, i.e. multiple branches of Hamiltonian paths connecting points in configuration space. While this issue has been examined in a few specific cases [16,[37][38][39], no general theory is currently available to the best of our knowledge. A procedure for coherently summing the contributions of multiple paths, similar to the Maslov theory in quantum mechanics [20,40], would be highly desirable, both for accurate numerical computations and for the theoretical analysis of most-likely noisy paths of a dynamical system.
Finally, the semiclassical approximation may be a valuable tool for analyzing experimental data of noisy swimmers in fluid flows. For example, with a sufficiently large number of experimentally-recorded trajectories of the type shown in Fig. 1, it would be possible to test the semiclassical predictions of the exit-right probability Pr(x 0 ). It should also be possible to investigate the distribution of experimentally-measured trajectories satisfying specific boundary conditions [37,41]. The semiclassically-predicted distributions may be used to fit the experimental data in order to extract physical parameters, such as rotational diffusivity and swimmer shape [42].

ACKNOWLEDGMENTS
We thank Tom Solomon for stimulating discussions. This material is based upon work supported by the National Science Foundation under Grant No. CMMI-1825379.
Appendix A: Semiclassical approximation for the Fokker-Planck equation We consider the stochastic process in a d-dimensional phase space where w = (w 1 , w 2 , . . . , w n ) is a set of uncorrelated Wiener processes and C is a d × n matrix, assumed to be constant for simplicity. The G term is included in Eq. (A1) to account for potential noise-induced drift [17]. The corresponding Fokker-Planck equation for the probability density P (q, t) is where D = CC T is the diffusion tensor (up to a factor of 2). The diffusion tensor is required to be positive-definite.
Our goal is to find an approximate solution to Eq. (A2) in the weak-noise (ε 1) limit.
We use an approach closely related to the semiclassical approximation of quantum mechanics [20], and our derivation closely follows Ref. [17]. Similar techniques have been applied to stochastic dynamics in a variety of settings [15,16,18,29,[37][38][39]43]. We consider an asymptotic expansion of the solution where N is the maximum number of terms in the expansion and the S n are functions to be determined. We restrict our attention to N = 2 and rewrite the solution as where W = S 0 and A = e −S 1 . Equation (A4) constitutes a WKB approximation to Eq. (A2).
It is in the same spirit as the semiclassical approximation to the Schrödinger equation; with the substitution ε → i , the semiclassical wave function is expressed as ψ = Ae iW/ . In that case, W is the action of the classical system associated to the quantum Hamiltonian, and A 2 = |ψ| 2 is the probability density of the system. In the Fokker-Planck case, we shall see that W also corresponds to the action of a particular classical Hamiltonian derived from the Fokker-Planck equation, while A is essentially a normalization function.
The functions A and W (equivalently the S n ) are determined by substituting Eq. (A3) into Eq. (A2). This yields the equation The goal is to equate terms of equal order in ε in Eq. (A5), which leads to equations for the S n . We are only interested in the ε 0 and the ε 1 terms.
The solutions to Eq. (A5) depend on the initial condition P 0 (q) ≡ P (q, 0). We derive the solutions for three types of initial conditions. The first is a Dirac δ function centered on an arbitrary phase-space point q 0 , where we use the "0" subscript to refer to an initial condition fixed at some particular value. In this case, the solution is denoted P = K(q, q 0 , t), where K is called the propagator, with initial condition The propagator encodes the probability for the system to make a transition from q 0 to q in time t. This fact, combined with the linearity of the Fokker-Planck equation, allows its solution with any initial probability distribution P 0 to be written as However, evaluating Eq. (A7) in practice is very computationally costly, which leads us to consider initial conditions which vary smoothly over some or all of the q coordinates. In particular, we consider an initial condition already in WKB form (A4) [43], i.e.
Last, we consider initial conditions that are a hybrid between the WKB form (A8) and the propagator form (A6), of the type Here, the variables are split into two groups, q = (q f , q k ), where the q f variables are fixed at a specific value q f 0 at t = 0 by the δ function in Eq. (A9), while the q k variables are distributed according to the WKB part of Eq. (A9). We are able to solve for W and A for this hybrid initial condition in the special case that D is block-diagonal, with one block corresponding to the q f variables and the other to the q k variables, i.e.
where d f is the size of q f and d k is the size of q k . The equations satisfied by W and A in Eq. (A4) are the same in each of the three cases, but the solutions must be selected such that the initial condition is satisfied. In Secs. A 1-A 3, we derive the solutions for the propagator initial condition (A6), and in Sec. A 4, we derive the solutions for the WKB and hybrid initial conditions.

Hamilton-Jacobi equation
From the zeroth order of Eq. (A5), we find W satisfies H(q, p) = 1 2 p · Dp + p · F(q). Here, q is the coordinate of the original stochastic dynamical system (A1), and p is the momentum canonically conjugate to q. The solution of Eq. (A11) is obtained by using the method of characteristics, which aims to find solutions of the form W (q(t), t) along particular paths (the characteristics) q(t). For the case of Eq. (A11), the characteristics turn out to be the projections of trajectories (q(t), p(t)) of the Hamiltonian system into configuration space. The relationship between the canonical momentum and W is The characteristics obey the equationṡ The overdots signify the total time-derivative d/dt along the characteristics q(t).
Equation (A18) is alternately known as the Onsager-Machlup action functional [28] or the Freidlin-Wentzell action functional [29]. i.e. a d-dimensional surface in phase space containing the initial conditions (q , p ), must be selected so that the initial condition is satisfied, i.e. lim t→0 A(q, t) exp[−W (q, t)/ε] = P 0 (q).
We use primed variables to refer to the space of initial conditions. Though Eq. (A13) explicitly gives the relationship between W and the Lagrangian manifold for times t = 0, it is ill-defined for any P 0 (q) for which W is singular in the t → 0 limit. This occurs when the projection of the initial Lagrangian manifold into configuration space is singular, which is the case for the propagator and hybrid initial conditions.
For the propagator initial condition (A6), the initial Lagrangian manifold is the phasespace surface q = q 0 . This means the Lagrangian manifold includes all possible initial momenta p ∈ R d . The projection of this Lagrangian manifold into q space is singular because all points of the Lagrangian manifold project to the same configuration space point, q 0 . The solution of Eq. (A11) in this case is W (q, t) = R(q, q 0 , t), where R is Hamilton's principal function, given by where (q(τ ), p(τ )) is the Hamiltonian path [i.e. solution of Eqs. (A14) and (A15)] going from q 0 to q in time t. Equation (A19) clearly follows from Eq. (A16). Because this solution is obtained by evolving the surface q = q 0 forward in time, each path has a distinct initial momentum p 0 . In fact, p 0 may be expressed in terms of q 0 , q, and t as To prove it, we consider the change in R as we make an infinitesimal change to the initial coordinate q 0 , while keeping the final coordinate q and transit time t fixed. We get In Sec. A 3, we show that Eq. (A19) indeed gives a probability density that satisfies the initial condition (A6).
All of our arguments make the assumption that there is only one Hamiltonian path going from q 0 to q in time t, so that R is single-valued (for all times t > 0). In general, however, at longer times there are multiple paths connecting q 0 and q in the same time t, with distinct initial momenta p 0 . Thus, the action becomes multi-valued. This situation also arises in semiclassical quantum mechanics. In that case, the quantum propagator consists of a sum of terms of the form of Eq. (A4), one for each branch of solutions, that are stitched together in such a way that the propagator is continuous [20,40]. We are not aware of a similar procedure for the Fokker-Planck equation. Thus, we continue to assume that there is a unique characteristic for any (q 0 , q, t) and hence, R is single-valued.

Transport equation
Next, we look at the equation arising from the first order terms of Eq. (A5). These lead to a transport equation for S 1 , which can be rearranged into a transport equation for A.
This results inȦ where we recallȦ = dA/dt = ∂A/∂t+q·∇A, withq given by Eq. (A14). We solve Eq. (A21) by integrating along the characteristics q(t) defined be Eqs. (A14)-(A16). Rearranging, we Equation (A22) gives the change of A as one moves along a characteristic path from a point on the initial Lagrangian manifold, with configuration space coordinate q , to the terminal coordinate q in time t. By integration, we obtain where the integration on the right-hand side of Eq. (A23) is carried out along the characteristic path. Equation (A23) applies to any initial Lagrangian manifold, but it must be handled carefully for initial Lagrangian manifolds with a singular projection, such as the propagator case for which q = q 0 . The propagator initial condition (A6) is itself singular at q = q 0 , and it turns out that A is also singular at this point. We introduce the quantity A 0 ≡ A(q 0 , 0) as a placeholder for now, and it is properly accounted for in Sec. A 3.
To clarify the first term on the right-hand side of Eq. (A23), we introduce the phase space functions (q(z , t), p(z , t)). The functions (q(z , t), p(z , t)) express the positions and momenta (q, p) at time t of Eqs. (A14) and (A15) as a function of their initial coordinate z on the Lagrangian manifold. We allow for an arbitrary parametrization z of the initial Lagrangian manifold, and we express the initial conditions as (q (z ), p (z )). For the propagator initial condition, the initial Lagrangian manifold can be simply parametrized as (q (z ), p (z )) = (q 0 , z ).
The functions (q, p) can be expressed using the flow functions (Q(q , p , t), P(q , p , t)), which map initial conditions (q , p ) to their values Q and P at time t and satisfy Hamilton's equations (A14) and (A15). It is obvious that where we have used Eq. (A27). From this it follows The right-hand side of Eq. (A30) is the first term that appears in Eq. (A23) up to a factor of −1/2. In Eq. (A30), For the case of the propagator initial condition, q → q 0 as t → 0, independent of z , and therefore ∂q/∂z = 0 in the limit, implying D 0 = 0. This is related to the divergence of A as t → 0, so we keep D 0 as a placeholder here and return to this point in Sec. A 3. Using Eq. (A30), we may now solve Eq. (A23) for A, giving where we have replaced z by p by virtue of Eq. (A24).

Normalization of the propagator
To fix the value of the constant A 0 √ D 0 , we must consider the limit t → 0 and impose the initial condition (A6). In this limit we have the estimates Solving Eq. (A33) for q, we find Hence, the semiclassical propagator in this limit becomes Equation (A37) is a Gaussian approximation to the δ function initial condition (A6), which approaches the δ function in the t → 0 limit. Therefore, the solutions for W and A satisfy the initial condition, provided that so that the Gaussian is properly normalized. The final semiclassical expression for the propagator is (A39)

a. WKB initial condition
For the WKB initial condition (A8), the initial Lagrangian manifold is the surface defined by p = ∂W 0 ∂q (q ).
We parametrize the Lagrangian manifold by (q (z ), p (z )) = z , ∂W 0 ∂q (z ) , and let the phase space functions be (q(z , t), p(z , t)). The position coordinate part of this function is assumed to be invertible, so that the initial configuration space coordinate can be expressed as q 0 = q 0 (q, t). Then, the solution to Eq. (A11) is [20,36] W (q, t) = W 0 (q 0 (q, t)) + R(q, q 0 (q, t), t), where R is evaluated along the Hamiltonian trajectory with initial coordinate q 0 (q, t) and initial momentum given by Eq. (A40) evaluated at q 0 (q, t).
The solution to the transport equation (A21) is almost identical to the propagator initial condition case. The main differences are the parametrization of the Lagrangian manifold (A41) and the specific initial condition A(q, 0) = A 0 (q). Taking these into account, we obtain A(q, t) = A 0 (q 0 (q, t)) det ∂q ∂q The final expression for the semiclassical probability density is P (q, t) ≈ A 0 (q 0 ) det ∂q ∂q For the hybrid initial condition (A9), the initial Lagrangian manifold is the surface defined by We parametrize the Lagrangian manifold by the coordinates z ≡ (p f , q k ), and let the phase space functions be (q(z , t), p(z , t)). We assume an inverse to q(z , t) exists, which we denote (p f 0 , q k0 ) = z 0 = z 0 (q, q f 0 , t). Then, the solution to the Hamilton-Jacobi equation is where q 0 = (q f 0 , q k0 (q, q f 0 , t)), and the initial momentum of the trajectory terminating at q at time t is p k0 (q, q f 0 , t) = ∂U ∂q k (q k0 ).
Solving the transport equation is again similar to the propagator case. The quantity det ∂q/∂z still satisfies Eq. (A29), while in Eq. (A23) we write A(q , 0) = A 0 (q k )A * , where A * is a placeholder constant. We therefore obtain A(q, t) = A * D 0 A 0 (q k0 (q, q f 0 , t)) det ∂q ∂z Taking the t → 0 limit of the full solution and forcing it to satisfy the initial condition (A9) leads to Thus, the full semiclassical probability density is P (q, t) = 1 (2πε) d f /2 A 0 (q k0 ) det ∂q ∂z For sufficiently small noise and times, it is useful to approximate the semiclassical propagator as a Gaussian centered on the deterministic trajectory q * (t). Because of the exponential form of Eq. (A39), for short times and small noise, the most important part to capture is the part near the absolute minimum of the action R(q, q 0 , t) = 0, which due to Eq. (A18), occurs at the deterministic solution q = q * (t). We illustrate the approximation for the 1D case for simplicity, so that q → q and F → F , and we let the diffusion tensor D → 1 and G → 0. Thus, Hamilton's equations becomė K(q, q 0 , t) ≈ A(q * ) exp − 1 2ε where In Eq. (A56), the Jacobian matrix and integral are evaluated along the deterministic trajectory q * (t), so that A(q * ) is a function of time. Equation (A55) thus constitutes a Gaussian approximation to the propagator, which would be properly normalized provided that Next, we show that this is indeed the case.
Recalling that p = ∂R/∂q, we have that We now rewrite Eq. (A58) in terms of partial derivatives with respect to the initial conditions (q , p ), which we can then compute using the tangent flow equations of Eqs. (A53) and (A54).