Drift-diffusion processes from elimination of fast variables under inhomogeneous conditions

TL;DR

This paper develops a unified method to eliminate fast variables in inhomogeneous environments, deriving effective drift-diffusion processes for particles, applicable to both passive and active motion, with detailed generator calculations.

Contribution

It introduces a straightforward approach to derive effective drift-diffusion models from complex systems with position-dependent velocities and fast variables, enhancing understanding of environmental effects on particle motion.

Findings

Unified framework for eliminating fast variables

Explicit drift-diffusion generators for various models

Insights into inhomogeneity-induced drift in passive and active particles

Abstract

The problem of eliminating fast-relaxing variables to obtain an effective drift-diffusion process in position is solved in a uniform and straightforward way for models with velocity a function jointly of position and fast variables. A more unified view is thereby obtained of the effect of environmental inhomgeneity on the motion of a diffusing particle, in particular, whether a drift is induced, covering both passive and active particles. Infinitesimal generators (equivalently, drift-diffusion fields) for the contracted processes are worked out in detail for several models.

Equations

$$(P^{t}f)(p)={\mathbb{E}}_p[f(X_t)]={⟨f(t)⟩}_p,$$

$$\mathcal{Z}=u^i(x)\frac{\partial }{\partial x^i}+D^{ij}(x)\frac{{\partial }^2}{\partial x^i\partial x^j},$$

$$u^i(x_0)$$

$$2D^{ij}(x_0)$$

$$u^{\alpha }(y)=\frac{\partial y^{\alpha }}{\partial x^j}{u}^j(x)+\frac{{\partial }^2{y}^{\alpha }}{\partial x^j\partial x^k}{D}^{jk}(x).$$

$${\mathbb{E}}_{\mu }[f_2(X_{t_1+t_2})f_1(X_{t_1})]$$

$$=\int \mu (dx)P^{t_1}{f}_1{P}^{t_2}{f}_2(x).$$

$${\mathcal{Z}}_{A_0}-\gamma (x)v^i\frac{\partial }{\partial v^i}+D(x){\Delta }_v+v^i\frac{\partial }{\partial x^i},$$

$${\mathcal{Z}}_{B_0}=D(x){\Delta }_{\mathbb{S}}+v_0(x){\xi }^i\frac{\partial }{\partial x^i},$$

$$\mathcal{Z}={\mathcal{Z}}_{\mathrm{fst}}+v^i(\xi ,x)\frac{\partial }{\partial x^i},$$

$$P_{\text{fst}}^{t}f={⟨f⟩}_{\text{fst}}+\mathcal{O}\left(e^{-t/\tau }\right).$$

$${\mathcal{Z}}_{\text{fst,new}}:=\frac{1}{ϵ}{\mathcal{Z}}_{\text{fst,old}},$$

$$v_{\text{new}}(\xi ,x)-{⟨v_{\text{old}}⟩}_{\text{fst}}:={ϵ}^{1/2}(v_{\text{old}}(\xi ,x)-{⟨v_{\text{old}}⟩}_{\text{fst}})$$

$$\delta \mathcal{Z}=v\cdot {\partial }_x:=v^i(x,\xi )\frac{\partial }{\partial x^i},$$

$$P^T$$

$$=P_{\text{fst}}^T+{\int }_{0\le t_1\le T}{P}_{\text{fst}}^{t_1}(\delta \mathcal{Z})P_{\text{fst}}^{T-t_1}$$

$$+{\int }_{0\le t_2\le t_1\le T}{P}_{\text{fst}}^{t_2}(\delta \mathcal{Z})P_{\text{fst}}^{t_1-t_2}(\delta \mathcal{Z})P_{\text{fst}}^{T-t_1}+\cdots$$

$${\int }_{0\le t_n\le \cdots \le t_1\le T}$$

$$\cdots P_{\text{fst}}^{t_2-t_3}(\delta \mathcal{Z})P_{\text{fst}}^{t_1-t_2}(\delta \mathcal{Z})P_{\text{fst}}^{T-t_1}.$$

$$\overline{\mathcal{Z}}$$

$$S_1$$

$$S_{n+1}$$

$${\overline{\mathcal{Z}}}_{\text{markovian}}=$$

$$+{\int }_0^{\infty }dt{⟨\left(\tilde{v}\cdot {\partial }_x\right)P_{\text{fst}}^t\tilde{v}\cdot {\partial }_x⟩}_{\text{fst}},$$

$$\tilde{v}:=v-{⟨v⟩}_{\text{fst}}.$$

$$\overline{u}={⟨v⟩}_{\text{fst}}+{\int }_0^{\infty }dt{⟨\tilde{v}(0)\cdot ({\partial }_x\tilde{v}(t))⟩}_{\text{fst}},$$

$$\overline{D}=2{\int }_0^{\infty }dt\mathrm{Sym}{⟨\tilde{v}(0)\otimes \tilde{v}(t)⟩}_{\text{fst}},$$

$${\tilde{v}}^i(\xi ,x)={{\sigma }^i}_j(x){\tilde{v}}_0^j(\xi );$$

$${\mathcal{Z}}_{\mathrm{fst}}{\tilde{v}}^i=-{{\Lambda }^i}_j{\tilde{v}}^j,$$

$$P_{\text{fst}}^t\tilde{v}=e^{-t\Lambda }\tilde{v}.$$

$$\overline{u}-{⟨v⟩}_{\text{fst}}=$$

$$={\int }_0^{\infty }{⟨(\tilde{v}(0)\cdot {\partial }_x)\exp \left(-t\Lambda \right)\tilde{v}⟩}_{\text{fst}}$$

$$={⟨(\tilde{v}\cdot {\partial }_x){\Lambda }^{-1}\tilde{v}⟩}_{\text{fst}}.$$

$$K:={⟨\tilde{v}\otimes \tilde{v}⟩}_{\text{fst}},$$

$${\overline{u}}^i-{⟨v^i⟩}_{\text{fst}}$$

$$2{\overline{D}}^{ij}$$

Full text

1 Introduction

Thirty years ago, van Kampenvan Kampen (1988) investigated the effects of environmental inhomogeneity in models of random motion of passive particles. In general, the problem continues to attract interestChristensen and Pedersen (2003); Bringuier (2009); Yang and Ripoll (2013); Livshits (2016); Bhattacharyay (2019); Baldovin et al. (2019); Abdoli et al. (2019), and has been especially revived for active matterHong et al. (2007); Butler et al. (2015); Byun et al. (2017); Zhao et al. (2018); Popescu et al. (2018). Cates and TailleurCates and Tailleur (2013) and Vuijk et aliaVuijk et al. (2019), in particular, studied variants of a model to be taken up here. Consider, for instance, a micron-sized particle propelling itself via chemical reactions using fuel from the environment, at a speed increasing with fuel concentration. Random torques cause its direction of motion to wander (time scale $\sim 1$ sec), so that its motion is diffusive on longer time scales. Will a fuel concentration gradient cause drift? If so, the mechanism would seem to be rather different than for a passive particle.

This paper presents a unified perspective for velocity models — models for which velocity is given as a function of continuous position and fast-relaxing variables. I show how to eliminate the fast variables to obtain the infinitesimal generator for a positional drift-diffusion in the markovian limit by a uniform and natural method. Drift is just a part of this contraction problem. Although the infinitesimal generator is interconvertible with a Fokker-Planck-type equation, being a dual, or adjoint formulation, the method given here seems significantly simpler than methods which work with a probability density from the beginning. Significant generalizations of both the active particle model mentioned above and the Langevin-Kramers model of a brownian particle discussed by van Kampen are treated in detail, together with a coupling of them. All those belong to the class of quasilinear velocity models, for which practical formulas for drift-diffusion fields are given, requiring a handful of tensor fields on position space expressing simple properties of the original velocity model and its fast equilbria. Intuition suggests that, all other things being equal, a passive brownian particle drifts toward lower velocity damping while the simple active particle drifts toward higher fuel concentration. Suitably embellished, these are the only two general mechanisms of gradient-induced drift for quasilinear models, and the same is arguably roughly true for general velocity models.

Zoom out now to situate that contraction problem in a broader idea of coarse-graining. Drift-diffusions, characterized by drift ($u$) and diffusion ($D$) coefficients, and markovian random processes with continuous sample paths, are essentially the same thing. Therefore, to make a drift-diffusion approximation (of a simple random walk, for an elementary example) is to approximate by a markovian process with continuous sample paths. If the process as given lacks markovianity, it may display it approximately over longer time scales — a temporal coarse-graining is called for. The processes treated in this paper have a velocity determined jointly by position and some fast variables, that is, with quickly decaying correlations. They belong to this first general category. If the process as given has discontinuous sample paths (a so-called hopping modelvan Kampen (1988, 1992); Christensen and Pedersen (2003); Smith et al. (2017), for instance) it might be approximable as continuous after a spatial coarse-graining. Those are outside the purview of this paper. A seeming third category comprises discrete time markovian processes, such as the archetypal random walk. However, viewed as time-homogeneous continuous-time processes, as they must be to prepare for a drift-diffusion approximation, they are nonmarkovian. Therefore, they too conform to the rule of temporal/spatial coarse-graining to obtain markovianity/continuous paths.

Here is a content guide. Section 2 reviews the aforementioned equivalence of drift-diffusions with markovian with continuous paths, and the observable-centered (“Heisenberg picture”) semigroup formalism we will use. Section 3 solves the contraction problem — elimination of fast variables — for a general velocity model, leading to expressions for the infinitesimal generator of the contracted process in terms of time correlation functions in fast-equilibrium. [See Eqs. (15), (16), (17)]. Section 4 elaborates the Langevin-Kramers model, to one with arbitary tensor parameters, and a simple model of an active particle, to one with fully anisotropic rotational diffusion characteristics, as well as a coupling of the two. The elaborated models are fit into the framework of (quasilinear) models describable in terms of intrinsic drift ($\left\langle{v}\right\rangle_{\text{fst}}$), decay rate ($\Lambda$), velocity transformation ($\sigma$), and kinetic ($K$), tensor fields. A general conclusion here is that drift which is not already implied by fast equilibrium is due to derivatives of either $\Lambda$ or $\sigma$. Section 5 makes contact with the more common Fokker-Planck perspective, discussing the relationship of $(\overline{u},\overline{D})$ with drift and diffusive currents is discussed, as well as position-conditional mean instantaneous velocity.

2 Background: Feller processes and semigroups on manifolds

This Section is a quick review of some standard loreGardiner (1985); Kallenberg (1997), with two principal objectives. First, to introduce the semigroup formulation which will be used subsequently. This is a description which focuses on the evolution of observables rather than a probability density. It stands in relation to what might loosely be called a Fokker-Planck-type formulation as the Heisenberg picture in quantum mechanics stands to the Schrödinger picture. Secondly, and concurrently, to recall the equivalence between Markov processes with continuous sample paths and drift-diffusions.

For a Markov process $X$ on a manifold ${\cal M}$, and with $X_{t}$ denoting the state of the system (a point in ${\cal M}$) at time $t$, operators $P^{t}$ are defined via

[TABLE]

this being the expectation of $f(X_{t})$ assuming initial state $p$. The markovian no-memory property is expressed by the Chapman-Kolmogorov equation $P^{s}P^{t}=P^{t+s}$ ($s,t\geq 0$), which can also be read as saying that $(P^{t})_{t\in[0,\infty)}$ comprises a semigroup. A Feller semigroup is one which maps $C_{0}({\cal M})$ (continuous functions tending to zero at infinity) into itself and is strongly continuous with respect to supremum norm, and its (infinitesimal) generator ${\cal Z}$ is defined by ${\cal Z}f{\kern 1.99997pt}:{\kern-5.0pt}={\kern 1.99997pt}\lim_{t\downarrow 0}\frac{1}{t}(P^{t}f-f)$, with domain the subspace of $C_{0}({\cal M})$ on which the limit exists. The generator and semigroup are related by $\frac{d}{dt}P^{t}f={\cal Z}P^{t}f$, solution to which is formally written as $P^{t}=\exp t{\cal Z}$. Assuming as well that $C_{c}^{\infty}({\cal M})\subset\mathrm{dom}\,{\cal Z}$ (so far, all is all consistent with a jump process) and that all sample paths are continuous, ${\cal Z}$ is guaranteed to be a second-order differential operator. In a coordinate chart ${x}\colon{{\cal U}}\rightarrow{{\mathbb{R}}^{n}}$,

[TABLE]

where

[TABLE]

and $\chi\in C_{c}^{\infty}({\cal U})$ equals one on some neighborhood of $x_{0}$. Here and throughout, summation convention is followed for matching upper and lower indices. Eqs. (2) provide the interpretation of $(u,D)$. Exactness of the second-order form goes back to a 1931 paper of KolmogorovKolmogoroff (1931). Now, while $u^{i}\partial_{i}$ is covariant, $D^{ij}\partial_{i}\partial_{j}$ is not. This has the slightly confusing consequence that $D^{ij}$ transforms as a twice-contravariant tensor, but the drift has an anomalous term in its coordinate-transformation rule:

[TABLE]

This curious fact indicates that $(u,D)$ should be thought of as a unit, not two fully distinct things.

In the converse direction to the conclusion represented by (2), suppose given a differential operator as generator of a semigroup $(P^{t})$. If $P^{t}$ bears interpretation as the semigroup associated with a stochastic process on ${\cal M}$ — a probability measure or coherent system of starting-point-indexed probability measures on paths $[0,\infty)\rightarrow{\cal M}$ — then the process is perforce markovian (Chapman-Kolmogorov equation) and has continuous sample paths since ${\cal Z}f(x)=0$ whenever $f$ vanishes on any neighborhood of $x$. Further, as already seen, ${\cal Z}$ must in that case be second order.

The semigroup formulation is convenient for calculating multi-time correlation functions. For example, with initial measure $\mu$ and $t_{1},t_{2}\geq 0$,

[TABLE]

In the final expression, time ordering must be observed, and $P^{t_{1}}$ acts on everything to its right (i.e., $f_{1}P^{t_{2}}f_{2}$).

A general Markov process on ${\cal M}$ can be thought of as a continuous motion occasionally interupted by unanticipated jumps. Later, we will discuss adding jumps to our prototype models, but will do no concrete calculations, hence do not introduce any notation for them.

3 Contraction of velocity models

This section constructs the contraction of a general velocity model, with special attention to the markovian limit. The generator $\overline{{\cal Z}}$, equivalently the drift-diffusion field $(\overline{u},\overline{D})$ of the contracted process, are expressed in terms of time correlation functions in the fast equilibrium. As concrete examples to keep in mind, simple prototypes will be described first.

3.1 Prototypes

We introduce some very simple concrete models as examples. They live on product spaces ${\mathbb{R}}_{x}^{d}\times F$, where the first factor is for position and the second for some internal, fast degrees of freedom.

$A_{0}$: $F={\mathbb{R}}_{v}^{d}$ is velocity space and the generator is

[TABLE]

where $\Delta$ is the $d$-dimensional euclidean laplacian. This is nothing but the Langevin-Kramers model for a particle in a thermal bath, except that there is no assumption that anything is thermal, so the fluctuation-dissipation link is broken and $\gamma$ and $D$ are functions of $x$ (as emphasized in the notation). The drift in Model $A_{0}$ was investigated already decades agoErmak and McCammon (1978); van Kampen (1988).

$B_{0}$: $F={\mathbb{S}}^{d-1}$, the unit sphere in ${\mathbb{R}}^{d}$, so a point $\xi$ in $F$ can be thought of as a unit vector in ${\mathbb{R}}^{d}$. $B_{0}$ has infinitesimal generator

[TABLE]

where $\Delta_{\mathbb{S}}$ is the Laplace-Beltrami operator on the sphere. As motivation, consider that this can model a self-powered colloidal particle which propels itself by consuming fuel in its environment. The concentration of fuel, hence its speed $v_{0}$ varies with the position-dependent concentration. At the same time, its direction of motion $\xi$ is subject to Brownian wandering represented by the first term of ${\cal Z}$. Variants of this model were investigated relatively recentlyCates and Tailleur (2013); Vuijk et al. (2019).

In Section 4 these models will be significantly elaborated. But the simple examples are a good thing to keep in mind during the following abstract development.

3.2 Slow and fast degrees of freedom

The general and natural geometrical setting for our investigation is as follows. ${Q}$ is a manifold, the configuration, or base, manifold, representing slow degrees of freedom of the system — position in the concrete models. ${\pi_{E}}\colon{E}\rightarrow{Q}$ is a fiber bundle over $Q$ with typical fiber $F$. $F$ represents fast degrees of freedom. Conventionally, coordinates or points in the fibers will be denoted by $\xi$. In the case of models $A_{0}$ and $B_{0}$, $E$ is just a product bundle, with $Q={\mathbb{R}}_{x}^{d}$. Several motivations can be adduced for working in the manifold setting. Certainly there are interesting cases of particles moving in noneuclidean geometries, in a lipid bilayer membrane, for example. Fast variables, such as orientational variables, are perhaps even more likely to inhabit noneuclidean spaces. Model $B$ of this paper is an example. Also, since the questions at issue are local in position, the mathematical overhead involved is fairly light, while nonuniform parameters are very natural in arbitary coordinates.

What determines the rate of change of position is the velocity, of course. $\xi$ may simply be velocity, as for Model $A_{0}$. In any case, the velocity is a function $v(x,\xi)$ of the fast and slow coordinates. Thus, $v\colon E\rightarrow TQ$ is a bundle map. The basic assumption made about the fundamental Markov process is that it have a generator of the form

[TABLE]

where ${\cal Z}_{\mathrm{fst}}$ generates a fiber-preserving semigroup $P_{\text{fst}}^{t}$ (the state never moves out of the fiber it started in) in each fiber. The structure of (8) explains the name “velocity model”. In addition, the fast processes are assumed to attain equilbrium on the time scale $\tau$:

[TABLE]

That is, $\left\langle{\cdot}\right\rangle_{\text{fst}}$ represents a fiber-by-fiber equilibrium depending on base position (in $Q$).

The smaller $\tau$ compared to our observation time scale for the contracted process, the closer the latter is to markovian. In order to access the markovian limit in a formal way, we scale ${\cal Z}_{\mathrm{fst}}$ up by a large factor, scaling down the time required to reach fast equilibrium. That alone would have the effect of rendering all fluctuations negligible, so a rescaling of the velocity is also needed. With temporary “old/new” labels for clarity, we take

[TABLE]

For infinitesimal $\epsilon$, the contracted process is markovian, and therefore, according to Section 2, fully characterised by drift-diffusion fields $(\overline{u},\overline{D})$ on $Q$. An overbar on a symbol indicates that it pertains to a contracted process.

3.3 Contracted semigroup

The fast generator ${\cal Z}_{\mathrm{fst}}$ generates a fast semigroup $P_{\text{fst}}$ and drives the system toward the fast equilibrium $\left\langle{\;}\right\rangle_{\text{fst}}$. These are all black boxes to us. In such generality, there are very few options of how to proceed. Take the clutching term

[TABLE]

from ${\cal Z}$, and put it into the Duhamel-Dyson formula as a perturbation:

[TABLE]

So far, this resembles perturbation expansions for evolution operators as they occur in many fields, quantum mechanics notably. However, we wish to apply this to functions of position (in $Q$) alone, thereby constructing a contracted semigroup $\overline{P}^{t}$ with generator $\overline{{\cal Z}}$. To that end, examine the structure of the $n$th term of the expansion:

[TABLE]

The integrand is a string of operators, $P_{\text{fst}}^{s}$ and $(\delta{\cal Z})$, each of which potentially acts on everything to its right. If $\epsilon\tau\ll{t_{k}-t_{k-1}}$, that action is effectively the same as taking the fast-equilibrium value. Thus, the integrand reduces to a product of clusters of the form $\left\langle{\text{ sequence of }\delta{\cal Z}\text{'s within times }\sim\epsilon\tau}\right\rangle_{\text{fst}}$. The entire perturbation expansion thus reduces, formally, to a sum over all (ordered) products of clusters as the cluster positions range over the interval $[0,T]$. That is, $\overline{P}^{T}=\exp(T\overline{{\cal Z}})$, with

[TABLE]

Insofar as we are interested in the markovian limit, however, this is overkill. Each factor of $\delta{\cal Z}$ carries a factor $\epsilon^{-1/2}$ in the $v$, and each integration $\int_{0}^{\infty}P_{\text{fst}}^{t}\cdots dt$ brings another factor $\epsilon$. Formally, then, $\left\langle{S_{1}}\right\rangle_{\text{fst}}=\left\langle{v}\right\rangle_{\text{fst}}\cdot{\partial_{x}}={\cal O}\left({1}\right)$, and, for $n\geq 0$, $\left\langle{S_{n+2}}\right\rangle_{\text{fst}}={\cal O}\left({\epsilon^{n/2}}\right)$. Therefore, in the limit $\epsilon\downarrow 0$, we obtain the markovian generator

[TABLE]

using the abbreviation

[TABLE]

(The subtraction of the first $\left\langle{v}\right\rangle_{\text{fst}}$ on the second line is allowed since what it multiplies has zero fast expectation value.) Henceforth, we are interested only in this markovian limit and therefore drop the subscript.

Now, recalling once again that the derivatives in (15) act on everything to the right and matching against the general form (2), the drift-diffusion fields of the contracted process in the markovian limit are revealed:

[TABLE]

and

[TABLE]

where ‘Sym’ indicates symmetrization on the tensor indices. Verification that these formulas conform to the coordinate-transformation law (4) is straightforward.

4 Prototypes elaborated

This Section explores various ways to elaborate the prototype models $A_{0}$ and $B_{0}$. Rather than doing so in an unconstrained fashion, the aim is to preserve a simplicity of the relevant correlation functions. Thus, in preparation we define a restricted class of velocity models with a convenient and fairly intuitive set of parameters: intrinsic drift ($\left\langle{v}\right\rangle_{\text{fst}}$), decay ($\Lambda$), velocity transformation ($\sigma$), and kinetic ($K$) tensors. Definitions follow.

4.1 Quasilinear models

quasilinear models are defined by two conditions.

(a) The velocity function has a factorized form:

[TABLE]

(b) $\tilde{v}$ decays exponentially, on average:

[TABLE]

where $\Lambda$ is a position-dependent linear operator with spectrum in the right complex half-plane $\mathrm{Re}\,z\gtrsim\tau^{-1}$.

Condition (b) implies that

[TABLE]

The calculations for $\overline{u}$ and $\overline{D}$ under these assumptions are similar, the latter being slightly easier. For $\overline{u}$ (16), we have

[TABLE]

Now, according to the factorization condition (a), $\partial_{x}$ here is oblivious to $\tilde{v}_{0}$, so we can rearrange factors and put it before the derivative. The result will be correct in any specific coordinate system. (Vide infra for comment about transformation properties.) With the definition

[TABLE]

we obtain

[TABLE]

The expression for the generator corresponding to this $(\overline{u},\overline{D})$ pair is actually hardly less compact than either. It is

[TABLE]

Here, the alternative notation for derivatives is meant to emphasize that $\partial_{j}$ is an operator that operates on everything to its right, not just the parenthesized expression immediately following. Written in terms of fast-equilibrium correlations, this is

[TABLE]

This expression is valid with only condition (b) imposed, but that is not the reason it is written here. (26) exhibits general coordinate covariance, whereas (25) does not. The reason is that, in the latter, a factor $\tilde{v}_{0}$ was moved in front of a derivative so that we could write things in terms of $K$. $\tilde{v}_{0}$ itself is not position-dependent, but with a change of coordinates, it will generate a jacobian factor which may be, but which is invisible in (25).

Before moving on to concrete models, a few observations are in order about the significance of (23) and (24). (i) Apart from $\left\langle{v}\right\rangle_{\text{fst}}$, there are two basic mechanisms to produce drift in a quasilinear model: gradients of either the decay tensor $\Lambda$ or the velocity scaling tensor $\sigma$. These two are precisely the mechanisms whereby drift is generated in the prototypes $A_{0}$ and $B_{0}$, respectively. (ii) In contrast to the case of a generic velocity mdoel, the quasilinear parametrization allows a precise statement of that sort. Still, it seems reasonable to say that something roughly like that holds in general. (iii) Derivatives of $K$ do not enter. Therefore, we cannot say that inhomogeneity generically induces drift. This is not in conflict with the existence of thermophoresisPiazza and Parola (2008). Temperature-dependent free energy of interaction with the surrounding medium can contribute to $\left\langle{v}\right\rangle_{\text{fst}}$. (iv) The only information required directly about the fast equilibria is that encoded in $\left\langle{v}\right\rangle_{\text{fst}}$ and $K$, which is to say, not very much.

4.2 Model $A$

Now, we briefly investigate what sorts of physically plausible elaborations can be made of Model $A_{0}$ (later, $B_{0}$) from Section 3.1 without spoiling the quasilinearity conditions (18) and (19). The first, most obvious way to modify Model $A_{0}$ is simply to make all the scalar coefficients tensors. Thus, Model A is defined by the fast generator ($\gamma$, $f$, and $D$ are functions of $x$ alone)

[TABLE]

Here, $f$ has the interpretation of some kind of external force. $\gamma$ is allowed to be anisotropic and may have partially rotational character, and $D$ is symmetric and positive, but otherwise unconstrained. With regard to the anisotropies, this is not written down with any specific physical system in mind, but one may imagine motion in an anisotropic medium, e.g. nematic liquid crystal, a magnetic field, and random excitation coming from an anisotropic nonthermal source such as ultrasonic pulses.

Three of the quasilinear parameters follow immediately:

[TABLE]

Finally, demanding $\left\langle{{\cal Z}_{A}(\tilde{v}\otimes\tilde{v})}\right\rangle_{\text{fst}}=0$ gives the equation

[TABLE]

which determines $K$ implicitly, and has a unique solution given the spectral restrictions on $\Lambda$. In fact, writing $K$ and $D$ as bilinear forms, (29) says $K({\Lambda}v,w)+K(v,{\Lambda}w)=2D(v,w)$. The solution is $K(v,w)=\int_{0}^{\infty}2D(e^{-t{\Lambda}}v,e^{-t{\Lambda}}w)\,dt$, which is manifestly symmetric and positive, since $D$ is.

Multiplying (29) by $\Lambda^{-1}$ and ${{\Lambda}^{\top}}^{-1}$ from left and right, respectively, yields

[TABLE]

bypassing $K$. Unfortunately, $K$ is still needed for $\overline{u}$, and no tidy formula appears available in general. Three special cases will be considered.

Model $A_{1}$: $\gamma=\gamma_{s}I$ is proportional to the identity (subscript ‘$s$’ for ‘scalar’). Then, the foregoing quickly leads to

[TABLE]

We might try to further specialize by supposing $D$ is scalar. That makes no sense in general, but if $Q$ has a riemannian metric $g$, we may suppose $D^{ij}$ and $K^{ij}$ proportional to $g^{ij}$. Write $K^{ij}=Tg^{ij}$, so $D^{ij}=(T/\gamma_{s})g^{ij}$; this can be viewed as a definition of $T$. Essentially, this is Model $A_{0}$ with an external force. The result here is easily seen to agree with that obtained by van Kampen in Section 4 of Ref. van Kampen, 1988 after notational translation.

Model $A_{2}$: $Q$ is equipped with a riemannian metric $g$ and $K^{ij}=Tg^{ij}$. The interpretation of this is that fast-equilibrium is thermal equilibrium at position-dependent temperature $T(x)$, and relies on something outside the model itself for its justification. Thus, the formula for $K$ is given by fiat, and implies the restriction

[TABLE]

where one index of $\gamma$ is raised with $g$ in the usual way. From this, we easily obtain

[TABLE]

Model $A_{3}$: $\Lambda$ has a full complement of smoothly varying position-dependent (complex!) eigenvectors labelled by $\alpha=1,\ldots,d$:

[TABLE]

With projection operators defined as

[TABLE]

by use of the dual basis $(\theta^{\alpha})_{1}^{d}$ to the basis $(\theta_{\alpha})_{1}^{d}$,

[TABLE]

we obtain

[TABLE]

Now, we can proceed to write everything else in terms of the eigenvalues $\lambda_{\alpha}$ and projections $P_{\alpha}$. The results are not particularly transparent, but from (40) below, one at least gets some sense of how derivatives of both the eigenvalues and the eigendirections contribute to the drift. The kinetic tensor is

[TABLE]

the contracted diffusion tensor is

[TABLE]

and the contracted drift is

[TABLE]

Consider briefly some additional elaborations that would preserve quasilinearity without writing formulas for the parameters. First, $D$ can be made velocity dependent. Secondly, a variety of jumps can be included — a mean zero jump distribution would contribute to $K$, while jumps with probability density depending on a scattering angle and fractional decrease of $|v|$ would contribute to $\Lambda$.

4.3 Model $B$

The elaboration of Model $B_{0}$ is perhaps more interesting. As for Model $A_{2}$, we assume here that $Q$ is a riemannian manifold. The prototype, $B_{0}$, was defined on the sphere bundle $Q\times{\mathbb{S}}^{d-1}$; $\xi$ was just a unit vector in that case. Implicitly, the particle was supposed to have only a single distinguished axis, that along which it propels itself. The full Model $B$ will be defined on the associated orthonormal frame bundle; $\xi$ will be an entire orthonormal frame keeping track of all the orientational degrees of freedom. Consider a particle shaped like a brick. This is easy to talk about and visualize, but an ellipsoid, or even more complicated shape with appropriate reflection symmetries, will have the same type of rotational characteristics. A brick has three principal axes for its hydrodynamic properties with, as a result, distinct rotational diffusion coefficients about each of them. A dynamically inert dot painted on one of each pair of opposite faces of the brick allows unambiguous definition of an oriented axis and hence the association of a distinct orthonormal frame with each particle configuration. In the $d$-dimensional case, the body-frame $\alpha\beta$ plane has its own distinct orientational diffusion coefficient $D_{\alpha\beta}$ ($D_{\beta\alpha}=D_{\alpha\beta}$, diagonal entries zero).

This picture is intuitively simple. However, some notation needs to be set up to write down formulas. Given a reference frame $({\bm{e}}_{1},\ldots,{\bm{e}}_{d})$, each body frame can be identified with the rotation matrix $\xi$ in $SO(d)$ which brings the reference frame into coincidence with it. That is, $\xi\in SO(d)$ is identified with the frame $(\xi_{1},\ldots,\xi_{d})$, where $\xi_{n}$ is the $n$th column of $\xi$ with components ${\xi^{i}}_{n}$. Now, let $\llbracket\alpha\beta;\theta\rrbracket$ denote the rotation by angle $\theta$ in the reference $\alpha\beta$ plane, and define the operators $\tilde{L}_{\alpha\beta}$ for $\alpha\neq\beta$ via

[TABLE]

This differentiates with respect to rotation about the local, configuration-dependent, $\alpha\beta$ plane for each $\xi$. For computational purposes, it is useful to note that $\tilde{L}_{\alpha\beta}$ is a derivation: $\tilde{L}_{\alpha\beta}(fg)=(\tilde{L}_{\alpha\beta}f)g+f(\tilde{L}_{\alpha\beta}g)$; and

[TABLE]

where $E_{\alpha\beta}$ is a matrix with $\alpha\beta$ component equal to $1$, $\beta\alpha$ component equal to $-1$ and all other components zero.

Finally, we can write down the generator for fully anisotropic orientational diffusion:

[TABLE]

To finish the specification of the model, assume that the velocity is along ${\bm{e}}_{1}$ in the reference configuration, so that its direction in configuration $\xi$ is ${\xi}_{1}$, and the velocity map is

[TABLE]

It may seem that there is no possible physical interpretation for nonscalar $\sigma$, but it will be put to good use in Section 4.4.

For Model $B$, the fast process induced on the tangent space $T_{x}Q$ via the velocity map $v(\cdot,x)$ is not generally markovian. To see that, suppose that $D_{1,2}$ is much larger than all other orientational diffusion coefficients. And, suppose that we have observations of $v$ at two times with a delay of order $D_{1,2}^{-1}$. Then, we have a good idea of where the local $2$-axis is and can expect the velocity to keep moving approximately in the $1$-$2$ plane for the near future. Nevertheless, the induced process is quasilinear, as follows from a straightforward calculation yieldng

[TABLE]

with the definition

[TABLE]

In case $D_{1\beta}=D_{\text{rot}}$ for $2\leq\beta\leq d$, the motion of the orientation $\xi_{1}$ is ergodic, other rotational degrees of freedom are decoupled, and the model effectively reduces to the prototype $B_{0}$.

Along with $\sigma(x)$, given a priori, the quasilinear parameters are then

[TABLE]

Substitution into the quasilinear formula (25) yields

[TABLE]

The formula for $K$ follows from $\left\langle{{\xi^{j}}_{i}{\xi^{k}}_{i}}\right\rangle_{\text{fst}}=\delta_{jk}/d$, or more strongly from equidistribution of $\xi_{1}$ over ${\mathbb{S}}^{d-1}$ in fast equilibrium, which the reader may consider evident enough. For a more formal derivation, consider $0=\left\langle{{\cal Z}_{B,\text{fst}}{\xi^{j}}_{i}{\xi^{k}}_{i}}\right\rangle_{\text{fst}}=\sum_{n}D_{in}[\left\langle{{\xi^{j}}_{n}{\xi^{k}}_{n}}\right\rangle_{\text{fst}}-\left\langle{{\xi^{j}}_{i}{\xi^{k}}_{i}}\right\rangle_{\text{fst}}]$. This implies that $\left\langle{{\xi^{j}}_{i}{\xi^{k}}_{i}}\right\rangle_{\text{fst}}=(\sum_{n}D_{in})^{-1}\sum_{n}D_{in}\left\langle{{\xi^{j}}_{n}{\xi^{k}}_{n}}\right\rangle_{\text{fst}}$ is independent of $i$ (each of a set of expectations is a weighted average of the others). Two appeals to orthogonality of $\xi$ after summing, first over $i$ and then over $j=k$, show that $\left\langle{{\xi^{j}}_{i}{\xi^{k}}_{i}}\right\rangle_{\text{fst}}=\delta_{jk}/d$.

Due to the interest for run-and-tumble particlesCates and Tailleur (2013), addition of jumps to the model are now briefly considered. In keeping with the operative symmetries, the jump rate measure $\nu(\xi;d\eta)$ (rate of transitions $\xi\mapsto d\eta$) is required to have a well-defined form with respect to the body frame, and to be invariant under reflection $r_{i}$ through the hyperplane normal to the body axes $\xi_{i}$, $i=2,\ldots,d$:

[TABLE]

With $|\nu|=\int\nu_{0}(d\eta)$ denoting the total jump rate, the associated contribution to the generator is defined via

[TABLE]

In particular, using reflection symmetry in the second equality,

[TABLE]

where $0\leq\alpha\leq 2$. For scattering which sends ${\bm{e}}_{1}$ uniformly over the sphere or pure reversal, $\alpha=$ $1$ or $2$, respectively. In the markovian limit, therefore, jumps simply make an additive contribution to the effective diffusion coefficient (and thence to $\Lambda$). $D_{\text{rot}}$ is already an effective quantity insofar as the various $D_{1\beta}$ are not necessarily equal. With jumps, it should be replaced by

[TABLE]

4.4 A two-stage contraction

Now, consider combining two velocity models, 1 and 2. A trivial way to do this is to define a (sum) velocity model via $v(\xi_{1},\xi_{2},x)=v_{1}(\xi_{1},x)+v_{2}(\xi_{2},x)$, where the fast parts of models 1 and 2 are independent. Then the sum model is not quasilinear, but that is not a problem, because each of models 1 and 2 makes its own separate contribution to $\bar{u}$ and $\bar{D}$ (i.e., we should add $\bar{u}$ and $\bar{D}$).

If models 1 and 2 have very different relaxation times, however, a more interesting possibility — two-stage contraction — emerges. Let us consider an example constructed from Models $A_{2}$ and $B$. It will be essentially the same model as treated by Vuijk et aliaVuijk et al. (2019), except that we allow inhomogeneous $T$ and $D_{\text{rot}}$. Following that source, the model is motivated as follows. Consider an active colloidal particle ($\sim 1\mu$m) a self-generated thrust $f_{0}$ in the direction ${\xi}_{1}$ moving in a liquid. Ordinarily, the velocity would simply be some scalar multiple of that, in fact $(f_{0}/m\gamma){\xi}_{1}$, where $\gamma^{-1}$ is the viscous relaxation time for velocity, which, for a micron-sized particle in aqueous medium, is of order 1 $\mu$s. Really, we should understand that velocity as the output of a Langevin-Kramers model with the external driving force being the particle’s self-generated thrust. But, what if the damping tensor $\gamma$ is not a simple scalar? The scenario envisaged by Vuijk et alia (loc. cit.) is a charged particle in a magnetic field, but we do not commit to that and allow general $\gamma$.

Since the relaxation time of the Model $A_{2}$ velocity in this scenario is of order $10^{-6}$ times that of the thrust direction, we first eliminate the Model $A_{2}$ variables. According to Eq. (33) ($m=1$), that gives

[TABLE]

The second term does not depend on $\xi$ so we set it aside temporarily. The coefficient in the first term is now taken as the velocity function of Model $B$: $v^{i}(x,{\xi}_{1})=f_{0}(\gamma^{-1}{\xi}_{1})^{i}$, equivalently, $\sigma=f_{0}\gamma^{-1}$, where both $f_{0}$ and $\gamma$ are position-dependent. Inserting this data into (47) immediately yields

[TABLE]

Adding this to the set-aside piece of $\overline{{\cal Z}}_{A}$ gives the fully contracted generator

[TABLE]

Note that all parameters here ($f_{0}$, $\gamma$, $T$, $D_{\text{rot}}$, and $g$) are potentially position-dependent. In the case of 3-dimensional euclidean space with uniform $T$, $D_{\text{rot}}$, and $f_{0}$, this agrees with Eq. (12) of Ref. Vuijk et al., 2019, where it was obtained by very different methods. (Note: ${\gamma^{i}}_{j}$ here corresponds to $\gamma{\Gamma^{i}}_{j}$ of Vuijk et alia.)

5 Contacting the Fokker-Planck perspective

This section makes contact with the Fokker-Planck perspective. Section 5.1 takes a look at the current and asks what, if anything, the approach taken here says about splitting it into drift and diffusion parts. Section 5.2 returns to the methods of Section 3 to derive the conditional (on position) mean instantaneous velocity. Throughout, $Q$ is assumed to be equipped with a riemannian metric $g$. That allows the definition of a scalar probability density $\rho$ with respect to volume measure $d\Omega=\sqrt{|g|}dx^{1}\cdots dx^{n}$. Thus, $\rho$ is independent of coordinates.

5.1 Fokker-Planck equation and currents

This subsection considers a drift-diffusion which may not have come from contraction, so the overbars on $(u,D)$ have been dropped.

Via integration by parts, the time evolution of the expectation of test function $f(x)$ is found to be

[TABLE]

This is simply an implicit way of writing the the Fokker-Planck equation $\frac{\partial\rho}{\partial t}=-\operatorname{\mathrm{Div}}J$. Here,

[TABLE]

is the probability current density with ${\bm{e}}_{i}$ the unit vector field along the $i$th coordinate direction and $\operatorname{\mathrm{Div}}V={|g|}^{-\frac{1}{2}}\partial_{i}({|g|}^{\frac{1}{2}}V^{i})$ denotes the covariant divergence of the contravariant vector field $V$. A possibly more aesthetically pleasing way to express the current density is implicitly via

[TABLE]

where $\nabla_{j}$ denotes a covariant derivative and $\Gamma^{i}_{jk}$ is a Christoffel symbol. Since the anomalous terms in the transformations of $u$ and $\Gamma$ cancel, the object in square brackets is a genuine covariant vector field.

The expression for the current is trivially rearranged as

[TABLE]

The literature seems to show a preference for rearranging the right side of (53) as (55) and identifying the first term as the total drift, with $-{\partial{D}^{ij}}/{\partial x^{j}}$ being an “extra drift”van Kampen (1988). I argue now against that identification, recognizing that, insofar as only the total current is physically relevant, it is a relatively low-stakes issue. Let us move to euclidean space to keep matters simple, and visualize the density $\rho$ as an ensemble of noninteracting particles. “Drift” ought to refer to a propensity of any one of these particles to move in a particular direction. According to the defintion (2), $u\,dt$ provides us with an unbiased estimate of where a randomly selected member of the ensemble will go, while $D\,dt$ quantifies the uncertainty. Furthermore, there is nothing strange about a diffusive current arising from a gradient of $D$ with uniform $\rho$. Consider the familiar, elementary explanation of ordinary diffusive current: a plane $P$ separates a low-density region on the left from a high-density region on the right. Particles move randomly away from wherever they are. More particles randomly move across from right-to-left than vice versa, hence a diffusive current. If now we have equal densities, but more vigorous random motion on the right (larger $D$), a right-to-left current will again arise. Thus, it is perfectly proper to consider $(\nabla D)\rho$ a contribution to diffusive current. All this may have a whiff of “Itô versus Stratonovich”van Kampen (1981) about it. However, since the approach taken here has neither the need for an interpretational rule that a Langevin equation does, not the complication of having a general probability distribution from the beginning, it is difficult to see any interpretational wiggle room.

5.2 Conditional mean instantaneous velocity

Now suppose a velocity model with contracted probability density $\rho(x)$ which is well-aged (i.e., has been freely evolving for many fast-relaxation times $\tau$). What is the joint distribution of position and velocity? Refs. Vuijk et al., 2019; Cates and Tailleur, 2013 used a direction calculation of that joint distribution to find Fokker-Planck equations for Model $B$ variants. The method of Section 3 avoided it, but it is nevertheless a question of independent interest. One’s first thought might be that the joint distribution takes the factorized form

[TABLE]

where $\rho_{\text{fst}}^{x}$ denotes the fast-equilibrium distribution at $x$. Actually, the distribution over $v$ is slightly skewed from this, even in the markovian limit. We will not work out the full distribution, but only the first moment $\left\langle{\tilde{v}}\right\rangle_{\rho}^{x}=\left\langle{v}\right\rangle_{\rho}^{x}-\left\langle{v}\right\rangle_{\text{fst}}^{x}$. Since $\rho$ is a distribution only over $x$, this requires some explanation. In accordance with the preceding discussion, the expectation with respect to $\rho(x)$ of a function $g$ of both $x$ and $v$ should be read as

[TABLE]

while for $g$ a function on $T_{x}Q$, $\left\langle{g}\right\rangle_{\rho}^{x}$ is the corresponding distribution conditioned on position being $x$. The $P^{T}$ operator in (57) means that the velocity has had a chance to equilibrate with the “environment” described by the position distribution $\rho$. A little reflection should convince one that the current density is related to the instantaneous velocity expectation via

[TABLE]

(For one-stage contraction; two stages is a bit more complicated.) The only real question is whether the framework of Section 3 is able to deliver this conclusion. In order to calculate the first moment $\left\langle{\tilde{v}}\right\rangle_{\rho}^{x}$ of the conditional mean instantaneous velocity, it will be enough to compute

[TABLE]

for an arbitrary test function $f(x)$ of position only. Now, recall (12)

[TABLE]

$P_{\text{fst}}$ is oblivious to $f$, so the first term gives a contribution $\int\rho(x)f(x)\left\langle{\tilde{v}^{i}}\right\rangle_{\text{fst}}^{x}\,d\Omega=0$ to (59). The next term is more interesting. Since $T$ is infinitely large compared to the fast relaxation time, we obtain (with an integration-by-parts in the final step)

[TABLE]

Since $f$ is an arbitrary test function, we conclude that

[TABLE]

Comparing this to Eq. (53) for the probability current density, one see that (58) is indeed verified.

All this prompts one to ask which kind of measurement would be expected to yield the drift velocity $\overline{u}$, which the instantaneous velocity $\left\langle{v}\right\rangle_{\rho}$. The calculations used here to find them both involved evolving the system for a time $T\gg\tau$, however, that evolution served very different purposes in the two cases. Suppose the system is started (or found) at $x$. Then, $\overline{u}$ is the expected average velocity of the system over the subsequent time $T\gg\tau$. If now, the system has a probability density $\rho$ (perhaps we have calculated this for a time $T^{\prime}\gg T$ via the Fokker-Planck equation from an initial point mass at $x$), then our expectation of the instantaneous velocity of the system, regardless of its actual position, is $\left\langle{v}\right\rangle_{\rho}$. “Instantaneous” here means that the measurement is done over a time small (or at worst, comparable to) $\tau$.

6 Conclusion

By means of a “Heisenberg picture” approach, elimination of fast variables was carried out in general for velocity models, resulting in a general expression for the infinitesimal generator $\overline{{\cal Z}}$, equivalently, the drift-diffusion field $(\overline{u},\overline{D})$, of the contracted drift-diffusion process on position space in terms of time correlation functions in fast-equilibrium. Compared to Chapman-Enskog type expansions, which work in the “Schrödinger picture”, the method seems to be simpler and more automatically general. Information about the distribution of fast variables, which comes naturally with Chapman-Enskog expansions, is also available in the current approach, if desired (Sec. 5.2). Examination of the subclass of quasilinear models showed that drift, other than that already present in fast-equilbrium, traces to derivatives of either velocity decay ($\Lambda$) or the dependence of velocity on fast variables ($\sigma$). The first of these is exemplified by the humble brownian particle in a viscosity gradient, and the second by the active particle in a fuel concentration gradient.

Acknowledgements.

I thank Prof. Vincent H. Crespi for critical reading and suggestions. This work was funded by the Penn State MRSEC, Center for Nanoscale Science, under award National Science Foundation DMR-1420620.