Three-dimensional diffusion with helical persistence

Persistent random walks have been proposed as models for a wealth of physical and biological systems [1–6] . In the simplest scenario, this process is characterized by the tendency to continue moving in the direction of the last step. In certain circumstances, this tendency can be viewed as a kind of inertial effect that keeps the particle moving in its direction of motion [7]. On the other hand, in the context of biology, the tendency may arise from the fact that animals move in the direction in which they are facing [8].

For one-dimensional systems the probability distribution function for the position of the walker independent of the direction it is moving, satisfies the telegrapher's equation in the continuous limit [7]. In two-dimensions, the probability distribution function for the position of the walker can be shown to satisfy a telegrapher's equation only approximately, and little is known for three and higher dimensions [9, 10]. However, in all cases, the net effect of persistence is that it induces a short term memory in the motion in the sense that the direction of the path is closely related to previous directions over a correlation time, but in the long run, the transport process remains diffusive.

An interesting generalization of this process is the case in which the tendency to move is not in the same direction as the previous step, but rather, at a fixed deviation from the previous direction. In two dimensions, this modification induces a tendency to 'loop' if, say, the direction of each step is narrowly centered at a fixed, non zero angle from the previous step. This modification was shown to give rise to non-monotonic behavior of the walker's effective diffusion constant as a function of the width of the distribution of angle changes [10]. For very small widths, the walker would perform many loops that brings it back close to its initial position, giving rise to small diffusions. If, on the other hand, the distribution is very broad, then the persistence effects are quickly lost and the process is essentially a normal random walk. However, between both limits, large loop segments can become uncorrelated before returning near the starting point, giving rise to a large effective diffusion coefficient [10].

In this work we address the three-dimensional case, namely we study the transport properties of persistent diffusion with curvature and torsion. This three-dimensional motion induces a tendency for the trajectories to spiral, a behavior that has been observed in many animals, indeed: '...helical trajectories are nearly universal for organisms that are less than 0.5 mm long and live in water.' [1].

Further, in spite of its limitations, persistent random walks are useful models for the shape of linear polymers with some degree of bond rigidity; thus persistent diffusion with curvature and torsion would be appropriate models for polymers that tend to form helices [4, 6].

In [5] a similar problem was addressed with the additional complication of chemotaxis, which we shall not consider here. In that work a set of stochastic equations was proposed for the description of stochastic helical motion and the persistence time for the effective persistent random walk was determined in the weak noise limit.

In this work we approach the problem from a different perspective: we derive a Fokker–Planck equation for the evolution of the probability distribution of the Frenet–Serret reference frame—for a thorough discussion of this approach, see e.g. [11]—and we evaluate the transport properties of the system from that equation. The model is shown to depend essentially on four parameters: the mean curvature, the mean torsion and two parameters characterizing the fluctuations of both these quantities. We analyze in detail the parameter dependence of the mean position and orientation of the walker after a time t, as well as the diffusion constant of the walker and various detailed characteristics of the diffusion, such as the rate at which the orientation of the helical axis—the so-called Darboux vector—relaxes to an isotropic distribution. Finally, we point out a peculiar transition in the behavior of the time-autocorrelation functions of the velocity as the parameters change.

As mentioned above, we formulate the problem in terms of the Frenet–Serret equations with noise terms added to the curvature and torsion, as in [5]. Thus, we are actually dealing with a set of nine coupled stochastic differential equations, however, there are only two independent noises in the system, thus, the fluctuations induced by the curvature noise on the evolution of the tangent vector, say, are not independent of the fluctuations of the normal vector, since these are partly due to the same noise; a similar thing happens with the binormal vector. Further, the evolution of the Frenet frame, including the noises in the curvature and torsion, must always always consist of pure rotations. To keep track of these details we find it useful to consider a discrete version of the process that tends to a unique continuum limit. Actually, this discrete process, with appropriately small values of the time step, is what is used to perform numerical simulations of the system.

Having a Markov chain that describes a helical random walk, we can write directly the master equation. From this discrete master equation we then proceed to take the continuum limit in a way that is analogous to the derivation of the diffusion equation from the random walk [7]. We thus construct a Fokker–Planck equation [12] for the full helical persistent diffusion process, from which we can calculate the correlation functions of the Frenet frame. This should be distinguished from the Fokker–Planck equation reported in [5], which describes the probability distribution of the helical axis. In particular, using our equation, we can evaluate the velocity autocorrelation function from which the effective diffusion constant of the process is obtained.

To formulate the persistent diffusive process, we consider the Frenet description of a trajectory [1–3, 5], which consists in the evolution of the right-handed frame of unit vectors given by $\hat{T}$, the vector tangent to the curve, the normal $\hat{N}$, in the direction of the derivative of $\hat{T}$, and the binormal $\hat{B}=\hat{T}\times \hat{N}$. The Frenet–Serret equations then state that

$$\begin{eqnarray}&&\frac{{\rm d}}{{\rm d}s}\left( \begin{array}{ccccccccccccccc} {\hat{T}} \\ {\hat{N}} \\ {\hat{B}} \\ \end{array} \right)=\left( \begin{array}{ccccccccccccccc} 0 & K & 0 \\ -K & 0 & T \\ 0 & -T & 0 \\ \end{array} \right)\left( \begin{array}{ccccccccccccccc} {\hat{T}} \\ {\hat{N}} \\ {\hat{B}} \\ \end{array} \right),\end{eqnarray} \tag{ 1 }$$

where s is the arc length, K and T are the curvature and torsion of the curve respectively, which can, in principle, depend on s too. The time evolution of the particle following the trajectory is given by

$$\begin{eqnarray}&&\frac{{\rm d}\vec{X}}{{\rm d}t}=v\hat{T},\end{eqnarray} \tag{ 2 }$$

where v is the speed of the particle, which we will take as a constant throughout this work.

To simplify the construction of the Fokker–Planck equation, we begin by considering a discrete version of the Frenet–Serret equation as follows. We write the arc length as $s=n\Delta s$ and take the following discretized version of the evolution equation:

$$\begin{eqnarray}&&{{\left( \begin{array}{ccccccccccccccc} {\hat{T}} \\ {\hat{N}} \\ {\hat{B}} \\ \end{array} \right)}_{n+1}}={{R}_{n}}{{\left( \begin{array}{ccccccccccccccc} {\hat{T}} \\ {\hat{N}} \\ {\hat{B}} \\ \end{array} \right)}_{n}},\end{eqnarray} \tag{ 3 }$$

where

$$\begin{eqnarray}&&{{R}_{n}}=\frac{1+\frac{1}{2}{{\Omega }_{n}}}{1-\frac{1}{2}{{\Omega }_{n}}},\end{eqnarray} \tag{ 4a }$$

$$\begin{eqnarray}{{\Omega }_{n}}=\left( \begin{array}{ccccccccccccccc} 0 & {{K}_{n}} & 0 \\ -{{K}_{n}} & 0 & {{T}_{n}} \\ 0 & -{{T}_{n}} & 0 \\ \end{array} \right).\end{eqnarray} \tag{ 4b }$$

The rationale for this choice is that at each step, the Frenet frame undergoes a pure rotation. Equivalently, inverting (3) we have

$$\begin{eqnarray}&&{{\left( \begin{array}{ccccccccccccccc} {\hat{T}} \\ {\hat{N}} \\ {\hat{B}} \\ \end{array} \right)}_{n}}=R_{n}^{-1}{{\left( \begin{array}{ccccccccccccccc} {\hat{T}} \\ {\hat{N}} \\ {\hat{B}} \\ \end{array} \right)}_{n+1}}.\end{eqnarray} \tag{ 5 }$$

Finally, we introduce noise by taking

$$\begin{eqnarray}&&{{K}_{n}}\to k\Delta s+{{\sigma }_{1}}{{\eta }_{1}}(n){{(\Delta s)}^{1/2}},\qquad {{T}_{n}}\to \tau \Delta s+{{\sigma }_{2}}{{\eta }_{2}}(n){{(\Delta s)}^{1/2}},\end{eqnarray} \tag{ 6 }$$

where k and τ will be assumed to be constant, while ${{\eta }_{1}}(n)$ and ${{\eta }_{2}}(n)$ are random variables satisfying

$$\begin{eqnarray}&&\langle {{\eta }_{j}}(n)\rangle =0,\qquad \langle {{\eta }_{j}}(n){{\eta }_{i}}(m)\rangle ={{\delta }_{i,j}}{{\delta }_{n,m}}.\end{eqnarray} \tag{ 7 }$$

To order $\Delta s$ we can evaluate the matrix R_n as follows

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} R_{n}^{-1} & \approx & 1-{{\Omega }_{n}}+\frac{1}{2}\Omega _{n}^{2} \\ \qquad & = & \left( \begin{array}{ccccccccccccccc} 1-\frac{1}{2}K_{n}^{2} & -{{K}_{n}} & {{K}_{n}}{{T}_{n}} \\ {{K}_{n}} & 1-\frac{1}{2}(K_{n}^{2}+T_{n}^{2}), & -{{T}_{n}} \\ {{K}_{n}}{{T}_{n}} & {{T}_{n}} & 1-\frac{1}{2}T_{n}^{2} \\ \end{array} \right) \\ \qquad & \approx & \left( \begin{array}{ccccccccccccccc} 1-\frac{\Delta s}{2}\sigma _{1}^{2}\eta _{1}^{2} & -k\Delta s-{{\sigma }_{1}}{{\eta }_{1}}{{(\Delta s)}^{1/2}} & \Delta s\;{{\sigma }_{1}}{{\sigma }_{2}}{{\eta }_{1}}{{\eta }_{2}} \\ k\Delta s+{{\sigma }_{1}}{{\eta }_{1}}{{(\Delta s)}^{1/2}} & 1-\frac{\Delta s}{2}(\sigma _{1}^{2}\eta _{1}^{2}+\sigma _{2}^{2}\eta _{2}^{2}) & -\tau \Delta s-{{\sigma }_{2}}{{\eta }_{2}}{{(\Delta s)}^{1/2}} \\ \Delta s\;{{\sigma }_{1}}{{\sigma }_{2}}{{\eta }_{1}}{{\eta }_{2}} & \tau \Delta s+{{\sigma }_{2}}{{\eta }_{2}}{{(\Delta s)}^{1/2}} & 1-\frac{\Delta s}{2}\sigma _{2}^{2}\eta _{2}^{2} \\ \end{array} \right), \\ \end{array}\end{eqnarray} \tag{ 8 }$$

where the n-dependencies have been suppressed from the last expression in the interest of typographical clarity.

Thus, if we denote ${{P}_{n}}(\hat{T},\hat{N},\hat{B})$ the probability density function for the Frenet frame being given by the vectors $\hat{T}$, $\hat{N}$ and $\hat{B}$ at step n, then, correct to order $\Delta s$, we derive the Fokker–Planck equation from the relation

$$\begin{eqnarray}&&{{P}_{n+1}}(\hat{T},\hat{N},\hat{B})=\left\langle {{P}_{n}}\left( R_{n}^{-1}(\hat{T},\hat{N},\hat{B}) \right) \right\rangle .\end{eqnarray} \tag{ 9 }$$

Here the angular brackets denote the average over the random noise ${{\eta }_{1}}(n)$ and ${{\eta }_{2}}(n)$ and $R_{n}^{-1}$ is as in (8).

To lighten the expressions to come, let us define the drifts

$$\begin{eqnarray}&&{{\vec{v}}_{T}}=\frac{\sigma _{1}^{2}}{2}\hat{T}+k\hat{N},\end{eqnarray} \tag{ 10a }$$

$$\begin{eqnarray}&&{{\vec{v}}_{N}}=\frac{\sigma _{1}^{2}+\sigma _{2}^{2}}{2}\hat{N}-k\hat{T}+\tau \hat{B},\end{eqnarray} \tag{ 10b }$$

$$\begin{eqnarray}&&{{\vec{v}}_{B}}=\frac{\sigma _{2}^{2}}{2}\hat{B}-\tau \hat{N}.\end{eqnarray} \tag{ 10c }$$

Putting in (9), ${{P}_{n}}(\hat{T},\hat{N},\hat{B})=P(\hat{T},\hat{N},\hat{B};s)$, with $s=n\Delta s$, expanding to linear order in $\Delta s$ and finally dividing by $\Delta s$, we obtain

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} \frac{\partial P}{\partial s} & = & -{{{\vec{v}}}_{T}}\cdot {{\nabla }_{T}}P-{{{\vec{v}}}_{N}}\cdot {{\nabla }_{N}}P-{{{\vec{v}}}_{N}}\cdot {{\nabla }_{B}}P \\ {} & {} & +\frac{1}{2}\mathop{\sum }\limits_{i}\mathop{\sum }\limits_{j}\left[ {{N}_{i}}{{N}_{j}}\left( \sigma _{1}^{2}\frac{{{\partial }^{2}}P}{\partial {{T}_{i}}\partial {{T}_{j}}}+\sigma _{2}^{2}\frac{{{\partial }^{2}}P}{\partial {{B}_{i}}\partial {{B}_{j}}} \right)+\left( \sigma _{1}^{2}{{T}_{i}}{{T}_{j}}+\sigma _{2}^{2}{{B}_{i}}{{B}_{j}} \right)\frac{{{\partial }^{2}}P}{\partial {{N}_{i}}\partial {{N}_{j}}} \right] \\ {} & {} & -\mathop{\sum }\limits_{i}\mathop{\sum }\limits_{j}\left[ \sigma _{1}^{2}{{N}_{i}}{{T}_{j}}\frac{{{\partial }^{2}}P}{\partial {{T}_{i}}\partial {{N}_{j}}}+\sigma _{2}^{2}{{N}_{i}}{{B}_{j}}\frac{{{\partial }^{2}}P}{\partial {{B}_{i}}\partial {{N}_{j}}} \right], \\ \end{array}\end{eqnarray} \tag{ 11 }$$

where ${{\nabla }_{T}}$, ${{\nabla }_{N}}$ and ${{\nabla }_{B}}$ denote the gradient in the T, N and B coordinates respectively. The equation above describes the evolution of the probability distribution for the components of the Frenet frame along the arc length of the particle's trajectory and becomes exact in the limit $\Delta s\to 0$. Equation (11) can be rewritten in a slightly more symmetric form as

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} \frac{\partial P}{\partial s} & = & -{{{\vec{v}}}_{T}}\cdot {{\nabla }_{T}}P-{{{\vec{v}}}_{N}}\cdot {{\nabla }_{N}}P-{{{\vec{v}}}_{N}}\cdot {{\nabla }_{B}}P \\ {} & {} & +\frac{\sigma _{1}^{2}}{2}\mathop{\sum }\limits_{i}\mathop{\sum }\limits_{j}\left( {{N}_{i}}\frac{\partial }{\partial {{T}_{i}}}-{{T}_{i}}\frac{\partial }{\partial {{N}_{i}}} \right)\left( {{N}_{j}}\frac{\partial }{\partial {{T}_{j}}}-{{T}_{j}}\frac{\partial }{\partial {{N}_{j}}} \right)P \\ {} & {} & +\frac{\sigma _{2}^{2}}{2}\mathop{\sum }\limits_{i}\mathop{\sum }\limits_{j}\left( {{N}_{i}}\frac{\partial }{\partial {{B}_{i}}}-{{B}_{i}}\frac{\partial }{\partial {{N}_{i}}} \right)\left( {{N}_{j}}\frac{\partial }{\partial {{B}_{j}}}-{{B}_{j}}\frac{\partial }{\partial {{N}_{j}}} \right)P; \\ \end{array}\end{eqnarray} \tag{ 12 }$$

either way, the equation seems to be rather intractable. However, it can be used to study the evolution of the moments of $\hat{T}$, $\hat{N}$ and $\hat{B}$; from these, the transport properties of $\vec{x}(t)$ can be inferred. In particular, by carrying out the integrations, one finds that

$$\begin{eqnarray}&&\frac{{\rm d}}{{\rm d}s}\left( \begin{array}{ccccccccccccccc} \langle {{{\hat{T}}}_{s}}|{{{\hat{T}}}_{0}},{{{\hat{N}}}_{0}},{{{\hat{B}}}_{0}}\rangle \\ \langle {{{\hat{N}}}_{s}}|{{{\hat{T}}}_{0}},{{{\hat{N}}}_{0}},{{{\hat{B}}}_{0}}\rangle \\ \langle {{{\hat{B}}}_{s}}|{{{\hat{T}}}_{0}},{{{\hat{N}}}_{0}},{{{\hat{B}}}_{0}}\rangle \\ \end{array} \right)=\left( \begin{array}{ccccccccccccccc} -\sigma _{1}^{2}/2 & k & 0 \\ -k & -[\sigma _{1}^{2}+\sigma _{2}^{2}]/2 & \tau \\ 0 & -\tau & -\sigma _{2}^{2}/2 \\ \end{array} \right)\left( \begin{array}{ccccccccccccccc} \langle {{{\hat{T}}}_{s}}|{{{\hat{T}}}_{0}},{{{\hat{N}}}_{0}},{{{\hat{B}}}_{0}}\rangle \\ \langle {{{\hat{N}}}_{s}}|{{{\hat{T}}}_{0}},{{{\hat{N}}}_{0}},{{{\hat{B}}}_{0}}\rangle \\ \langle {{{\hat{B}}}_{s}}|{{{\hat{T}}}_{0}},{{{\hat{N}}}_{0}},{{{\hat{B}}}_{0}}\rangle \\ \end{array} \right).\end{eqnarray} \tag{ 13 }$$

Here the dependence on the initial condition has been written explicitly to emphasize the fact that these statistics give us access to the time correlation functions. The solution can be written as

$$\begin{eqnarray}&&\left( \begin{array}{ccccccccccccccc} \langle {{{\hat{T}}}_{s}}|{{{\hat{T}}}_{0}},{{{\hat{N}}}_{0}},{{{\hat{B}}}_{0}}\rangle \\ \langle {{{\hat{N}}}_{s}}|{{{\hat{T}}}_{0}},{{{\hat{N}}}_{0}},{{{\hat{B}}}_{0}}\rangle \\ \langle {{{\hat{B}}}_{s}}|{{{\hat{T}}}_{0}},{{{\hat{N}}}_{0}},{{{\hat{B}}}_{0}}\rangle \\ \end{array} \right)={{{\rm e}}^{\Gamma s}}\left( \begin{array}{ccccccccccccccc} {{{\hat{T}}}_{0}} \\ {{{\hat{N}}}_{0}} \\ {{{\hat{B}}}_{0}} \\ \end{array} \right),\end{eqnarray} \tag{ 14 }$$

where Γ is the matrix appearing in (13). Thus, the matrix of correlation functions can be written explicitly as

$$\begin{eqnarray}\mathcal{M}(s)=\left( \begin{array}{ccccccccccccccc} \langle {{{\hat{T}}}_{s}}\cdot {{{\hat{T}}}_{0}}\rangle & \langle {{{\hat{T}}}_{s}}\cdot {{{\hat{N}}}_{0}}\rangle & \langle {{{\hat{T}}}_{s}}\cdot {{{\hat{B}}}_{0}}\rangle \\ \langle {{{\hat{N}}}_{s}}\cdot {{{\hat{T}}}_{0}}\rangle & \langle {{{\hat{N}}}_{s}}\cdot {{{\hat{N}}}_{0}}\rangle & \langle {{{\hat{N}}}_{s}}\cdot {{{\hat{B}}}_{0}}\rangle \\ \langle {{{\hat{B}}}_{s}}\cdot {{{\hat{T}}}_{0}}\rangle & \langle {{{\hat{B}}}_{s}}\cdot {{{\hat{N}}}_{0}}\rangle & \langle {{{\hat{B}}}_{s}}\cdot {{{\hat{B}}}_{0}}\rangle \\ \end{array} \right)={{{\rm e}}^{\Gamma s}}.\end{eqnarray} \tag{ 15 }$$

The Fourier–Laplace transform of the correlation matrix is thus given by

$$\begin{eqnarray}&&\mathcal{M}(\omega )={{({\rm i}\omega -\Gamma )}^{-1}}.\end{eqnarray} \tag{ 16 }$$

If we now define $\Delta (\omega )$ as follows

$$\begin{eqnarray}&&\Delta (\omega )=(\sigma _{1}^{2}+2{\rm i}\omega )(\sigma _{2}^{2}+2{\rm i}\omega )(\sigma _{1}^{2}+\sigma _{2}^{2}+2{\rm i}\omega )+4{{k}^{2}}(\sigma _{2}^{2}+2{\rm i}\omega )+4{{\tau }^{2}}(\sigma _{1}^{2}+2{\rm i}\omega ),\end{eqnarray} \tag{ 17 }$$

we find the matrix $\Delta (\omega )\mathcal{M}(\omega )/2$ to be given by

$$\begin{eqnarray}\left( \begin{array}{ccccccccccccccc} (\sigma _{2}^{2}+2{\rm i}\omega )(\sigma _{1}^{2}+\sigma _{2}^{2}+2{\rm i}\omega )+4{{\tau }^{2}} & 2k(\sigma _{2}^{2}+2{\rm i}\omega ) & 4k\tau \\ -2k(\sigma _{2}^{2}+2{\rm i}\omega ) & (\sigma _{1}^{2}+2{\rm i}\omega )(\sigma _{2}^{2}+2{\rm i}\omega ) & 2\tau (\sigma _{1}^{2}+2{\rm i}\omega ) \\ 4k\tau & -2\tau (\sigma _{1}^{2}+2{\rm i}\omega ) & (\sigma _{1}^{2}+2{\rm i}\omega )(\sigma _{1}^{2}+\sigma _{2}^{2}+2{\rm i}\omega )+4{{k}^{2}} \\ \end{array} \right).\end{eqnarray} \tag{ 18 }$$

Note that the matrix is not symmetric as the process is not invariant under simple time inversion. Rather, $\hat{T}$ and $\hat{B}$ are odd, whereas $\hat{N}$ is even, leading to the observed alternation between symmetric and antisymmetric elements in the matrix.

The explicit inversion of the above results to obtain the time-dependent correlation functions can, of course, be carried out in principle, as it only requires the solution of a third order polynomial. Nevertheless, the results are messy and unenlightening. Instead, we focus on particular cases which illustrate different aspects of the process. One simple situation is reached in the planar case, in which both the torsion τ and its associated noise vanish (${{\sigma }_{2}}=0$). In this case, the binormal vector is constant and the motion is restricted to a plane. In this situation, equation (18) simplifies to

$$\begin{eqnarray}\quad {{\mathcal{M}}_{2d}}(\omega )=\frac{2}{\Delta (\omega )}\left( \begin{array}{ccccccccccccccc} 2{\rm i}\omega (\sigma _{1}^{2}+2{\rm i}\omega ) & 4{\rm i}\omega k & 0 \\ -4{\rm i}\omega k & 2{\rm i}\omega (\sigma _{1}^{2}+2{\rm i}\omega ) & 0 \\ 0 & 0 & {{(\sigma _{1}^{2}+2{\rm i}\omega )}^{2}}+4{{k}^{2}} \\ \end{array} \right),\end{eqnarray} \tag{ 19 }$$

where $\Delta (\omega )$ here reduces to

$$\begin{eqnarray}&&{{\Delta }_{2d}}(\omega )=2{\rm i}\omega [{{(\sigma _{1}^{2}+2{\rm i}\omega )}^{2}}+4{{k}^{2}}].\end{eqnarray} \tag{ 20 }$$

The s-dependent correlations can also be determined explicitly in this case. Indeed one readily finds

$$\begin{eqnarray}{{{\rm e}}^{\Gamma s}}=\left( \begin{array}{ccccccccccccccc} {\rm cos} (ks){{{\rm e}}^{-\sigma _{1}^{2}\,s/2}} & {\rm sin} (ks){{{\rm e}}^{-\sigma _{1}^{2}\,s/2}} & 0 \\ -{\rm sin} (ks){{{\rm e}}^{-\sigma _{1}^{2}\,s/2}} & {\rm cos} (ks){{{\rm e}}^{-\sigma _{1}^{2}\,s/2}} & 0 \\ 0 & 0 & 1 \\ \end{array} \right),\end{eqnarray} \tag{ 21 }$$

which corresponds to circular motion which is damped due to the noise.

Further, we can use $\mathcal{M}$ to write the Fourier–Laplace transforms for the average of the Frenet vectors; for example, for the mean tangent vector we have:

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} \langle {{{\hat{T}}}_{\omega }}|{{{\hat{T}}}_{0}},{{{\hat{N}}}_{0}},{{{\hat{B}}}_{0}}\rangle & = & \frac{2}{\Delta (\omega )}\left\{ \left[ (\sigma _{2}^{2}+2{\rm i}\omega )(\sigma _{1}^{2}+\sigma _{2}^{2}+2{\rm i}\omega )+4{{\tau }^{2}} \right]{{{\hat{T}}}_{0}} \right. \\ {} & {} & +\left. 2k(\sigma _{2}^{2}+2{\rm i}\omega ){{{\hat{N}}}_{0}}+4k\tau {{{\hat{B}}}_{0}} \right\}. \\ \end{array}\end{eqnarray} \tag{ 22 }$$

This quantity connects to the behavior of the walker's average position by

$$\begin{eqnarray}&&\langle {{\vec{x}}_{\omega }}|{{\hat{T}}_{0}},{{\hat{N}}_{0}},{{\hat{B}}_{0}}\rangle =\frac{1}{{\rm i}\omega }\langle {{\hat{T}}_{\omega }}|{{\hat{T}}_{0}},{{\hat{N}}_{0}},{{\hat{B}}_{0}}\rangle .\end{eqnarray} \tag{ 23 }$$

Thus, in the long time limit, the average position of the persistent walker saturates at

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} \langle {{{\vec{x}}}_{t}}|{{{\hat{T}}}_{0}},{{{\hat{N}}}_{0}},{{{\hat{B}}}_{0}}\rangle & \to & \langle {{{\vec{x}}}_{\infty }}|{{{\hat{T}}}_{0}},{{{\hat{N}}}_{0}},{{{\hat{B}}}_{0}}\rangle \\ {} & = & 2\frac{\left\{ \sigma _{2}^{2}(\sigma _{1}^{2}+\sigma _{2}^{2})+4{{\tau }^{2}} \right\}{{{\hat{T}}}_{0}}+2k\sigma _{2}^{2}{{{\hat{N}}}_{0}}+4k\tau {{{\hat{B}}}_{0}}}{\sigma _{1}^{2}\sigma _{2}^{2}(\sigma _{1}^{2}+\sigma _{2}^{2})+4{{k}^{2}}\sigma _{2}^{2}+4{{\tau }^{2}}\sigma _{1}^{2}}. \\ \end{array}\end{eqnarray} \tag{ 24 }$$

Another interesting quantity is the Darboux vector [1], which in our case, will be written as $\vec{A}=\tau \hat{T}+k\hat{B}$. This vector is akin to the angular momentum of the walker. In the absence of noise, this vector is constant, as can be seen directly from (13), and lies on the axis of the helical motion. In the presence of noise, however, the Darboux vector evolves according to

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} \quad \langle {{{\vec{A}}}_{\omega }}|{{{\hat{T}}}_{0}},{{{\hat{N}}}_{0}},{{{\hat{B}}}_{0}}\rangle & = & \frac{2}{\Delta (\omega )}\left\{ (\sigma _{1}^{2}+\sigma _{2}^{2}+2{\rm i}\omega )\left( \sigma _{2}^{2}\tau {{{\hat{T}}}_{0}}+\sigma _{1}^{2}k{{{\hat{B}}}_{0}} \right) \right. \\ {} & {} & +2\left[ {\rm i}\omega \left( \sigma _{1}^{2}+\sigma _{2}^{2}+2{\rm i}\omega \right)+2({{k}^{2}}+{{\tau }^{2}}) \right]{{{\vec{A}}}_{0}} \\ {} & {} & +\left. 2k\tau \left( \sigma _{2}^{2}-\sigma _{1}^{2} \right){{{\hat{N}}}_{0}} \right\}. \\ \end{array}\end{eqnarray} \tag{ 25 }$$

Here ${{\vec{A}}_{0}}=\tau {{\hat{T}}_{0}}+k{{\hat{B}}_{0}}$ is the initial value of $\vec{A}$. In figure 1 through 3 we show the evolution of the components of the Darboux vectors in relation to the initial frame of reference. It is also worth noting that in the small noise limit, the saturation value of the average position, given in (24) is given by

$$\begin{eqnarray}&&\langle {{\vec{x}}_{\infty }}|{{\hat{T}}_{0}},{{\hat{N}}_{0}},{{\hat{B}}_{0}}\rangle \approx \frac{2\tau {{{\vec{A}}}_{0}}}{\sigma _{2}^{2}{{k}^{2}}+\sigma _{1}^{2}{{\tau }^{2}}},\end{eqnarray} \tag{ 26 }$$

which reflects the intuitively appealing fact that in the low-noise limit, the net displacement is in the direction of the axis of the helix.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Furthermore, from the correlation functions we can compute the effective diffusion constant. Assuming, again, that the instantaneous speed v is constant, we can write

$$\begin{eqnarray}&&\left\langle \vec{x}(t)\cdot \vec{x}(t) \right\rangle ={{v}^{2}}\int _{0}^{t}\int _{0}^{t}\left\langle {{{\hat{T}}}_{{{s}_{1}}}}\cdot {{{\hat{T}}}_{{{s}_{2}}}} \right\rangle \;{\rm d}{{t}_{1}}\;{\rm d}{{t}_{2}},\end{eqnarray} \tag{ 27 }$$

recalling that s_i = vt_i. The effective diffusion constant D will be given by

$$\begin{eqnarray}&&D=\mathop{{\rm lim} }\limits_{t\to \infty }\frac{\left\langle \vec{x}(t)\cdot \vec{x}(t) \right\rangle }{6t}=\frac{v}{3}\int _{0}^{\infty }\left\langle {{{\hat{T}}}_{s}}\cdot {{{\hat{T}}}_{0}} \right\rangle {\rm d}s=\frac{2v}{3}\;\frac{\sigma _{2}^{2}(\sigma _{1}^{2}+\sigma _{2}^{2})+4{{\tau }^{2}}}{\sigma _{1}^{2}\sigma _{2}^{2}(\sigma _{1}^{2}+\sigma _{2}^{2})+4{{k}^{2}}\sigma _{2}^{2}+4{{\tau }^{2}}\sigma _{1}^{2}}.\end{eqnarray} \tag{ 28 }$$

This expression's behavior becomes easier to visualize after introducing scaled versions of the curvature, the torsion and the diffusion constant as well as the adimensional parameter κ, which measures the relative importance of the curvature and the torsion:

$$\begin{eqnarray}&&\bar{k}=\frac{k}{{{\sigma }_{1}}\sqrt{\sigma _{1}^{2}+\sigma _{2}^{2}}},\end{eqnarray} \tag{ 29a }$$

$$\begin{eqnarray}&&\bar{\tau }=\frac{\tau }{{{\sigma }_{2}}\sqrt{\sigma _{1}^{2}+\sigma _{2}^{2}}},\end{eqnarray} \tag{ 29b }$$

$$\begin{eqnarray}&&\bar{D}=\frac{3\sigma _{1}^{2}}{2v}\;D,\end{eqnarray} \tag{ 29c }$$

$$\begin{eqnarray}&&\kappa =\frac{4{{{\bar{k}}}^{2}}}{1+4{{{\bar{\tau }}}^{2}}}.\end{eqnarray} \tag{ 29d }$$

In terms of these parameters, equation (28) implies

$$\begin{eqnarray}&&\bar{D}=\frac{1+4{{{\bar{\tau }}}^{2}}}{1+4{{{\bar{k}}}^{2}}+4{{{\bar{\tau }}}^{2}}}=\frac{1}{1+\kappa }.\end{eqnarray} \tag{ 30 }$$

Thus $\bar{D}$ grows with $\bar{\tau }$ and decreases with $\bar{k}$. Intuitively, this corresponds to the fact that torsion favors a rectilinear motion in the direction of the Darboux vector, thus favoring diffusion by enhancing transport, whereas curvature leads to looping around in circles, which confines the walker's motion. For the true diffusion constant, we can write $\delta =2v/(3\sigma _{1}^{2})$ so that together with the dimensionless parameter κ, we have

$$\begin{eqnarray}&&D=\frac{\delta }{1+\kappa },\end{eqnarray} \tag{ 31 }$$

which is increasing in δ and decreasing in κ, but may show complex non-monotonic behavior in other variables.

Let us turn to the limiting behaviors of the diffusion constant. First note that if we set $\tau =0$ in equation (28), the result is a continuous function of ${{\sigma }_{2}}$, and that as ${{\sigma }_{2}}\to 0$, the diffusion constant approaches

$$\begin{eqnarray}&&{{D}_{2d}}=\frac{2\sigma _{1}^{2}}{\sigma _{1}^{4}+4{{k}^{2}}}\;\frac{v}{3},\end{eqnarray} \tag{ 32 }$$

which is (up to a factor of $1/3$ due to a difference in the definition of the diffusion constant in two and three-dimensions) the two-dimensional result previously derived in [10]. We note that, as mentioned above, D_2d is a non-monotonic function of ${{\sigma }_{1}}$.

On the other hand, a quite different, intrinsically three dimensional result, is obtained if we take the limits in opposite order. Indeed, if we let ${{\sigma }_{2}}=0$ first, we see that $\kappa =0$ and therefore

$$\begin{eqnarray}&&{{D}_{H}}=\frac{2}{\sigma _{1}^{2}}\;\frac{v}{3}.\end{eqnarray} \tag{ 33 }$$

Which is independent of all the other parameters of the problem, in particular of τ. From the scaling form (30), it follows quite straightforwardly that the same result also holds if k = 0. Intuitively, however, it is far from obvious why, in the former case, the diffusion constant only depends on ${{\sigma }_{1}}$ and v. In this case, it follows from (33) that the diffusion constant differs by a finite amount from the one that obtains when $\tau =0$; in other words, the diffusion constant is discontinuous in τ. This suggests that there exists a crossover time ${{t}_{c}}(\tau )$ going to infinity as $\tau \to 0$ such that an effective value of the diffusion coefficient close to D_2d is observed for times less than ${{t}_{c}}(\tau )$, whereas the asymptotic diffusion coefficient D_H is observed for $t\gg {{t}_{c}}(\tau )$. The geometric meaning of ${{t}_{c}}(\tau )$ is the time required for the orientation of the plane of the motion to change significantly.

Indeed, let us consider ${{\mathcal{M}}_{11}}(\omega )$ for the case in which ${{\sigma }_{2}}=0$ and both τ and ω are small compared to k and ${{\sigma }_{1}}$. This means we are considering large arc lengths (corresponding to long times) and small torsions. Then, introducing a new variable $w=\omega \sigma _{1}^{2}/2{{\tau }^{2}}$, we get

$$\begin{eqnarray}&&{{\mathcal{M}}_{11}}(\omega )\approx \frac{2\sigma _{1}^{2}(1+{\rm i}w)}{\sigma _{1}^{4}+{\rm i}w(4{{k}^{2}}+\sigma _{1}^{4})},\end{eqnarray} \tag{ 34 }$$

which is quite compatible with the picture presented above: one value (related to D_H as defined by (33)) arises for w = 0, but as w takes on values of order one, that is, on length scales of order $\sigma _{1}^{2}/2{{\tau }^{2}}$ and below, the Fourier transform of the velocity autocorrelation function changes rapidly and takes on values related to D_2d as defined in (32).

As a final remark, we observe a peculiar transition in the large-time behavior of the average trajectory. Indeed, the average position of the trajectory at arc length s is determined by the matrix ${\rm exp} (\Gamma s)$. For large s it must therefore be dominated by that eigenvalue of Γ having the smallest real part in absolute value. As is readily seen, the three eigenvalues of Γ in the limit of small noise are given by

$$\begin{eqnarray}&&{{\lambda }_{0}}=-\frac{{{\tau }^{2}}\sigma _{1}^{2}+{{k}^{2}}\sigma _{2}^{2}}{2({{k}^{2}}+{{\tau }^{2}})},\end{eqnarray} \tag{ 35a }$$

$$\begin{eqnarray}&&{{\lambda }_{\pm }}=-\frac{\sigma _{1}^{2}(2{{k}^{2}}+{{\tau }^{2}})+\sigma _{2}^{2}({{k}^{2}}+2{{\tau }^{2}})}{4({{k}^{2}}+{{\tau }^{2}})}\pm {\rm i}\sqrt{{{k}^{2}}+{{\tau }^{2}}},\end{eqnarray} \tag{ 35b }$$

up to terms of higher order in the noise. Thus, in the low-noise limit, we have two groups of eigenvalues, one tending to 0 and corresponding essentially to the Darboux vector as an eigenvector. This thus describes the motion in the direction of the helix axis. The other two eigenvalues correspond to eigenvectors determining the rotational motion. As we see, their real parts are not equal, and an easy calculation shows that ${{\lambda }_{0}}$ dominates if

$$\begin{eqnarray}&&\Phi =(2{{k}^{2}}-{{\tau }^{2}})\sigma _{1}^{2}+(2{{\tau }^{2}}-{{k}^{2}})\sigma _{2}^{2}\gt 0\end{eqnarray} \tag{ 36 }$$

and ${{\lambda }_{\pm }}$ otherwise. Thus, according to the sign of Φ, we have two qualitatively different motions: if $\Phi \gt 0$, the circular motion decays more rapidly than the progressive motion, and the asymptotic shape of the particle's trajectory is a very loosely wound helical motion, in which the particle in one rotation advances by a distance far larger than the helix' radius. If $\Phi \lt 0$, on the other hand, the opposite situation prevails, and the asymptotic motion is on a circle of a radius much larger than the distance by which the particle advances for each rotation. In other words, we might say that the asymptotic pitch of the average motion diverges at large times for $\Phi \gt 0$, whereas it vanishes for $\Phi \lt 0$. We show in figure 4 typical walks for two cases in which $\Phi \gt 0$ and $\Phi \lt 0$ respectively; in figure 5 we show the corresponding autocorrelation functions of the tangent vector $\hat{T}$ for large times. The difference is quite striking.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Finally, let us point out that the above-mentioned transition does not always occur. In particular, if k and τ are equal, or of similar size, then we always have $\Phi \gt 0$. For both signs of Φ to be possible, we require

$$\begin{eqnarray}&&{\rm max} \left( \frac{k}{\tau },\frac{\tau }{k} \right)\geqslant \sqrt{2}.\end{eqnarray} \tag{ 37 }$$

When $\Phi =0$, the model maintains a truly helicoidal dynamics for all times.

Let us briefly summarize what we have done: we began by defining a discretization of the Frenet–Serret equations describing a curve in three-dimensional space. Assuming that the curvature and torsion are both constant (we take them to be proportional to the discretization parameter $\Delta s$) these equations define a discretized helical motion. We add uncorrelated Gaussian noise of the order of ${{(\Delta s)}^{1/2}}$ to the curvature and torsion. The resulting process is, in fact, a random walk in which the probability of every step depends on the step's relative position with respect to the preceding two steps. Indeed, the curvature dictates the average angle to the preceding step, whereas the torsion determines the average angle with respect to the plane generated by the two earlier steps. The process is fully determined by the following parameters: the average curvature k, the average torsion τ, the noise intensity for the curvature $\sigma _{1}^{2}$, the noise intensity for the torsion $\sigma _{2}^{2}$ and finally the speed v at which the particle moves. While this last parameter only fixes an otherwise arbitrary time scale and can thus in principle be disregarded, the other four parameters are essential.

Next we derive the Fokker–Planck equation for the process. Using this equation, we then proceeded to determine the Fourier–Laplace transform of the various correlation functions of the orthogonal triad $\hat{T}$, $\hat{N}$ and $\hat{B}$. From these it is possible to determine the effective diffusion constant. An explicit formula for it is derived in terms of the model's parameters, see (28). Interestingly, there is a scaling parameter κ, see (3), such that the appropriately rescaled diffusion constant depends only on κ, see (30).

Other quantities which can be evaluated explicitly are the asymptotic average position (as a vector) of the walk, as well as the correlation function of the Darboux vector, which corresponds to the axis of the helix which would be generated in the noiseless limit. We show in particular that the asymptotic average position of the walk indeed coincides with the Darboux vector in the low-noise limit.

Various special cases can also be treated more or less explicitly. Among these are the two-dimensional case (${{\sigma }_{2}}=0$ and $\tau =0$). Finally, we found an interesting change in conformation as the parameters are varied. See section 4 for details. The connection of this phenomenon to other related transitions reported in the literature—see for example the transition reported by Marenduzzo et al in [13]—is still not clear and should be looked into in future work.

A still open question concerns the response of the process when an external field is applied; the difficulties for advancing in this question arise from the fact that a spatially fixed field rotates in the Frenet frame. A related problem would be to understand the effects in the transport properties due to periodic or strategic variations in the parameters, attempting to steer the particle in a given direction [1–3, 5].

FL acknowledges the support of UNAM-PAPIIT-DGAPA grant IN114014 as well as CONACyT grant 154586.