Beyond Itô versus Stratonovich

The stochastic differential equation, or the Langevin equation in physics, is usually a more precise description of reality than the purely deterministic one [1]–[3]. Additive noise is frequently considered for physical systems, but for biological or other complex systems a state variable dependent (multiplicative) noise is often encountered. In this paper we use the physicists' notation with a multiplicative noise

$$\begin{eqnarray} &&\dot {\mathbf{q}}=\mathbf{f}(\mathbf{q})+N(\mathbf{q})\xi (t), \end{eqnarray} \tag{ 1 }$$

where q, f are n-dimensional vectors and f is a nonlinear function of the state variable q. The noise ξ(t) is k-dimensional Gaussian white noise with zero mean, 〈ξ(t)〉 = 0, and the covariance 〈ξ(t)ξ^τ(t')〉 = δ(t − t')I_k. The superscript τ denotes the transpose of a matrix, δ(t − t') is the Dirac delta function, 〈⋯〉 indicates the average over the noise distribution, I_k is the k-dimensional identity matrix. An element of the n × k matrix N(q) can be a nonlinear function of q. This matrix is further described by

$$\begin{eqnarray} &&N(\mathbf{q}){N}^{\tau }(\mathbf{q})=2 \epsilon D(\mathbf{q}), \end{eqnarray} \tag{ 2 }$$

where is a constant quantifying the noise strength and D(q) is an n × n positive semi-definite diffusion matrix. Note that if the noise has less than n independent components, k < n, D(q) has zero eigenvalue(s).

During the study of a biological switch [5], a structured form equivalent to equation (1) was discovered,

$$\begin{eqnarray} &&[S(\mathbf{q})+A(\mathbf{q})]\dot {\mathbf{q}}=-\nabla \phi (\mathbf{q})+\hat {N}(\mathbf{q})\xi (t), \end{eqnarray} \tag{ 3 }$$

where S(q) is symmetric and positive semi-definite, A(q) is antisymmetric, ϕ(q) is the potential function, ξ(t) is identical to that in equation (1) and the matrix $\hat {N}(\mathbf{q})$ is constrained by (see also section 3.1)

$$\begin{eqnarray} &&\hat {N}(\mathbf{q}){\hat {N}}^{\tau }(\mathbf{q})=2 \epsilon S(\mathbf{q}). \end{eqnarray} \tag{ 4 }$$

We note that when S = 0, equation (3) has the same structure as the Hamiltonian equation in physics. The (n = 2l)-dimensional system (l generalized coordinates and l generalized momenta are combined in the n-dimensional system) with

$$\begin{eqnarray*} A=\left (\begin{array}{@{}cc@{}} \displaystyle 0&\displaystyle {I}_{l}\\ \displaystyle -{I}_{l}&\displaystyle 0 \end{array}\right ) \end{eqnarray*}$$

and ϕ = H (the Hamiltonian) shows the embedded symplectic structure.

Two assumptions are implied in this framework. The first one is the equivalence of equations (1) and (3), that is, they are describing the same dynamical process. It is straightforward to go from equation (3) to equation (1), once [S(q) + A(q)] is nonsingular, which holds when the components q_i (i = 1,...,n) of the state variable q = (q₁,...,q_n)^τ are independent. The reverse, however, is much more difficult, including the need to obtain S(q), A(q) and ϕ(q) for general dynamics. Nevertheless, we can still assume, without a rigorous mathematical proof but rigorous enough from a working scientist's view, that the reverse part holds based on a case by case construction, that is, considering it as a protocol and not a theorem. Hence, we may replace $\dot {\mathbf{q}}$ with the right-hand side of equation (1) in (3)

$$\begin{eqnarray} &&\left [S(\mathbf{q})+A(\mathbf{q})\right ]\left [\mathbf{f}(\mathbf{q})+N(\mathbf{q})\xi (t)\right ]=-\nabla \phi (\mathbf{q})+\hat {N}(\mathbf{q})\xi (t). \end{eqnarray} \tag{ 5 }$$

The second assumption is the separated equality of the deterministic and stochastic dynamics in equation (5)

$$\begin{eqnarray} &&\left [S(\mathbf{q})+A(\mathbf{q})\right ]\mathbf{f}(\mathbf{q})=-\nabla \phi (\mathbf{q}), \end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray} &&\left [S(\mathbf{q})+A(\mathbf{q})\right ]N(\mathbf{q})=\hat {N}(\mathbf{q}). \end{eqnarray} \tag{ 7 }$$

Intuitively, this assumption on separation is plausible for two different reasons: first, the zero mean noise function is nowhere differentiable but the deterministic forces are usually smooth (C^∞) functions, and hence two very different mathematical objects; second, the stochastic force and the deterministic forces describe different timescales of a physical system with different physical origins. By replacing equation (6) with an equivalent form (8) and plugging equations (7) and (2) into (4), we obtain the potential condition (8) and the generalized Einstein relation equation (9)

$$\begin{eqnarray} &&\nabla \times \left \{ \left [S(\mathbf{q})+A(\mathbf{q})\right ]\mathbf{f}(\mathbf{q})\right \} =0, \end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray} &&[S(\mathbf{q})+A(\mathbf{q})]D(\mathbf{q})[S(\mathbf{q})-A(\mathbf{q})]=S(\mathbf{q}). \end{eqnarray} \tag{ 9 }$$

In principle, the potential function ϕ(q) can be derived analytically by solving the n(n − 1)/2 partial differential equations (under proper boundary conditions) of equation (8), together with the n(n + 1)/2 equations given by equation (9) (n² unknowns in [S(q) + A(q)] and n² equations). It can also be calculated numerically through a gradient expansion [6]. Detailed results near fixed points can be found in [9]. Explicit constructions of potential functions for limit cycles are contained in [17, 20].

In the one-dimensional case, A = 0, let  = k_BT, if the friction γ is a constant, then S = γ/k_BT, equation (9) reduces to SD = γD/k_BT = 1, namely, the product of the friction and diffusion coefficients is a constant, discovered by Einstein [31] one hundred years ago. Equation (9) is a generalized form of the Einstein relation in two ways: the diffusion matrix can be nonlinearly dependent of the state variable (a verification experiment has been reported [32]) and the detailed balance condition can be broken (A(q) ≠ 0, see also section 2.2.2).

From the generalized Einstein relation equation (9), we know that the symmetric part of [S(q)+A(q)]⁻¹, D(q) = 1/2[(S(q)+A(q))⁻¹ + ((S(q)+A(q))⁻¹)^τ], is the diffusion matrix defined in equation (2). Hence we have [S(q) + A(q)][D(q) + Q(q)] = I (I is the identity matrix, singularity of D(q), S(q) or A(q) can usually be tolerated), where Q(q) is antisymmetric. Rewriting equation (6) as

$$\begin{eqnarray} &&\mathbf{f}(\mathbf{q})=-\left [D(\mathbf{q})+Q(\mathbf{q})\right ]\nabla \phi (\mathbf{q}), \end{eqnarray} \tag{ 10 }$$

and by multiplying ∇ϕ(q) on both sides of equation (10), we reach the Hamilton–Jacobi equation

$$\begin{eqnarray} &&\mathbf{f}(\mathbf{q})\boldsymbol{ \cdot }\nabla \phi (\mathbf{q})+\nabla \phi (\mathbf{q})\boldsymbol{ \cdot }D(\mathbf{q})\boldsymbol{ \cdot }\nabla \phi (\mathbf{q})=0. \end{eqnarray} \tag{ 11 }$$

Finally, using the knowledge from statistical physics, we can consider that the potential function leads to a Boltzmann–Gibbs distribution

$$\begin{eqnarray} &&{\rho }_{s}(\mathbf{q})=\frac{1}{Z(\epsilon )}\exp \left (-\frac{\phi (\mathbf{q})}{\epsilon }\right ), \end{eqnarray} \tag{ 12 }$$

as the steady state distribution (if it exists) for the stochastic process described by equation (3): as the 'temperature' and ϕ(q) as the 'Hamiltonian'. Experimental verification is available for some cases, e.g.  [18].

Hence, the two assumptions to implement the equivalence between equations (1) and (3) lead to a potential function for the general Langevin equation and to a corresponding Boltzmann–Gibbs distribution for the final steady state independent of both matrices A(q) and S(q). While such a distribution function seems special, it should be pointed out that it is one of the most successful equations in physics, even broader than that of the stationary Schrödinger equation. The framework can be applied for systems with an arbitrary noise strength: equations (8) and (9) do not contain . In particular, deterministic systems can be regarded as the weak noise limit ( → 0) of stochastic systems. The framework can also be extended directly to include parameter dependent potential and corresponding dynamical components [13, 20], e.g. time t and temperature dependent ϕ(q,t,). We note that equations (6) and (9) have a close connection to Lyapunov's direct method in engineering [17]. This evidence makes the assumptions more plausible and appealing.

This construction further leads to a stochastic integration (A-type) different from conventional ones such as Itô or Stratonovich. One may enquire about the mathematical consistency in this new approach. We point out that the consistency of A-type integration may be seen from two aspects. First, equation (1) or (3) by themselves are not complete descriptions of a stochastic process. A proper stochastic interpretation is needed, hence A-type is possible in principle. Second, we will show in the next subsection that there is an important situation in which one starts without the need to differentiate various known stochastic integrations, and reaches the consequence that an A-type interpretation is the natural choice for the SDE equations (3) and (1).

In this section, we summarize the main steps to demonstrate that A-type integration is mathematically consistent. The rationale is as follows. The starting point is the important stochastic differential equation in physics, Newton equation in 2n dimensions with noise and friction. Its corresponding 2n-dimensional Fokker–Planck equation (FPE) is uniquely determined, known as the generalized Klein–Kramers equation. The valuable feature here is that all approaches to stochastic integration discussed in the present paper lead to the same 2n-dimensional FPE. When taking the zero mass limit, the 2n-dimensional Newton equation reduces to an n-dimensional equation in exactly the same form as equation (3), while the corresponding n-dimensional FPE derived from the generalized Klein–Kramers equation has the ability to permit the Boltzmann–Gibbs distribution as its stationary solution. After this zero mass limit, the connection between the resulting n-dimensional FPE and the stochastic differential equation is the A-type integration.

2.2.1. The generalized Klein–Kramers equation.

We begin with a general situation via doubling the degrees of freedom, the 2n-dimensional equation with noise, proposed in [6],

$$\begin{eqnarray} &&\dot {\mathbf{q}}=\frac{\mathbf{p}}{m},\\ &&\dot {\mathbf{p}}=-[S(\mathbf{q})+A(\mathbf{q})]\frac{\mathbf{p}}{m}-{\nabla }_{\mathbf{q}}\phi (\mathbf{q})+\hat {N}(\mathbf{q})\xi (t), \end{eqnarray} \tag{ 13 }$$

where S(q), A(q), ϕ(q), $\hat {N}(\mathbf{q})$ and ξ(t) are identical to those in equation (3), and the subscript of ∇_q means that it operates on q only, the same as in equation (3). We note that equation (13) has two properties: first, by taking the zero mass limit m → 0, p → 0, we can recover equation (3) from (13); second, we wish to emphasize that starting with equation (13), we can find that different stochastic interpretations, Itô, Stratonovich, or A-type, will reach the same result, since ∇_p ⋅[S(q) + A(q)] = 0.

The probability density function in the (q,p) phase space can be calculated through the path integral as

$$\begin{eqnarray} &&\rho (\mathbf{q},\mathbf{p},t)\equiv \langle \delta \left (\mathbf{q}-\bar {\mathbf{q}}\left (t,\{ \xi (t)\} \right )\right )\delta \left (\mathbf{p}-\bar {\mathbf{p}}\left (t,\{ \xi (t)\} \right )\right )\rangle , \end{eqnarray} \tag{ 14 }$$

where $\bar {\mathbf{q}}(t,\{ \xi (t)\} )$ and $\bar {\mathbf{p}}(t,\{ \xi (t)\} )$ represent the solution of the ordinary differential equation (ODE) reduced from equation (13) under a fixed sample path {ξ(t)} of the sample space of ξ(t) (a given noise configuration). By taking a partial derivative of time on both sides of equation (14), we obtain

$$\begin{eqnarray} &&\fl {\partial }_{t}\rho =\left \langle {\nabla }_{\bar {\mathbf{q}}}\delta \left (\mathbf{q}-\bar {\mathbf{q}}\right )\boldsymbol{ \cdot }\frac{\mathrm{d}\bar {\mathbf{q}}}{\mathrm{d}t}\delta \left (\mathbf{p}-\bar {\mathbf{p}}\right )+\delta \left (\mathbf{q}-\bar {\mathbf{q}}\right ){\nabla }_{\bar {\mathbf{p}}}\delta \left (\mathbf{p}-\bar {\mathbf{p}}\right )\boldsymbol{ \cdot }\frac{\mathrm{d}\bar {\mathbf{p}}}{\mathrm{d}t}\right \rangle \nonumber\\ &&\fl \phantom {\partial _t\rho } =-\int \delta \left (\mathbf{q}-\bar {\mathbf{q}}\right )\delta \left (\mathbf{p}-\bar {\mathbf{p}}\right )\left \{ {\nabla }_{\bar {\mathbf{q}}}\boldsymbol{ \cdot }\left [\dot {\bar {\mathbf{q}}}P(\{ \xi \} )\right ]+{\nabla }_{\bar {\mathbf{p}}}\boldsymbol{ \cdot }\left [\dot {\bar {\mathbf{p}}}P(\{ \xi \} )\right ]\right \} \hspace{0.167em} \mathrm{d}\{ \xi \} \nonumber\\ &&\fl \phantom {\partial _t\rho } =-{\nabla }_{\mathbf{q}}\boldsymbol{ \cdot }\frac{\mathbf{p}}{m}\rho -{\nabla }_{\mathbf{p}}\boldsymbol{ \cdot }\Biggl \{ -\left [S(\mathbf{q})+A(\mathbf{q})\right ]\boldsymbol{ \cdot }\frac{\mathbf{p}}{m}\rho -{\nabla }_{\mathbf{q}}\phi (\mathbf{q})\rho \nonumber\\ &&\fl \phantom {\partial _t\rho xx} +~\hat {N}(\mathbf{q})\langle \xi (t)\delta \left (\mathbf{q}-\bar {\mathbf{q}}\right )\delta \left (\mathbf{p}-\bar {\mathbf{p}}\right )\rangle \Biggr \} , \end{eqnarray} \tag{ 15 }$$

in which ρ, $\bar {\mathbf{q}}$ and $\bar {\mathbf{p}}$ have the same meanings as in equation (14), and P({ξ}) is the probability density function for the given noise configuration {ξ(t)}. From step 1 to step 2 of equation (15), integration by parts is used. According to the Novikov identity [33], for a given functional f[ξ(t),t] of noise ξ(t),

$$\begin{eqnarray} \langle \xi (t)f[\xi (t),t]\rangle &=&\int \langle \xi (t){\xi }^{\tau }({t}^{\prime})\rangle \left \langle \frac{\delta f[\xi (t),t]}{\delta \xi ({t}^{\prime})}\right \rangle \hspace{0.167em} \mathrm{d}{t}^{\prime}\nonumber\\ &=&\left \langle \frac{\delta f[\xi (t),t]}{\delta \xi (t)}\right \rangle , \end{eqnarray} \tag{ 16 }$$

and using

$$\begin{eqnarray} &&\frac{\delta [\int \nolimits _{0}^{t}\xi ({t}^{\prime})\hspace{0.167em} \mathrm{d}{t}^{\prime}]}{\delta \xi (t)}=\frac{1}{2}, \end{eqnarray} \tag{ 17 }$$

noting that the solution of equation (13) can be expressed as

$$\begin{eqnarray} &&\bar {\mathbf{q}}(t)=\mathbf{q}(0)+\int \nolimits _{0}^{t}\frac{\bar {\mathbf{p}}}{m}\hspace{0.167em} \mathrm{d}{t}^{\prime}, \end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray} &&\bar {\mathbf{p}}(t)=\mathbf{p}(0)+\int \nolimits _{0}^{t}\left \{ -\left [S(\bar {\mathbf{q}})+A(\bar {\mathbf{q}})\right ]\frac{\bar {\mathbf{p}}}{m}-{\nabla }_{\bar {\mathbf{q}}}\phi (\bar {\mathbf{q}})+\hat {N}(\bar {\mathbf{q}})\xi ({t}^{\prime})\right \} \hspace{0.167em} \mathrm{d}{t}^{\prime}, \end{eqnarray} \tag{ 19 }$$

we obtain the relations

$$\begin{eqnarray} &&\frac{\delta \bar {\mathbf{q}}(t)}{\delta \xi (t)}=0, \end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray} &&\frac{\delta \bar {\mathbf{p}}(t)}{\delta \xi (t)}=\frac{1}{2}{\hat {N}}^{\tau }(\bar {\mathbf{q}}). \end{eqnarray} \tag{ 21 }$$

The last left term in equation (15) can be handled as

$$\begin{eqnarray} &&\hat {N}(\mathbf{q})\langle \xi (t)\delta \left (\mathbf{q}-\bar {\mathbf{q}}\right )\delta \left (\mathbf{p}-\bar {\mathbf{p}}\right )\rangle =-2 \epsilon S(\mathbf{q})\boldsymbol{ \cdot }\frac{1}{2}{\nabla }_{\mathbf{p}}\rho . \end{eqnarray} \tag{ 22 }$$

We finally reach the generalized version of the Klein–Kramers equation [2]

$$\begin{eqnarray} &&{\partial }_{t}\rho ={\nabla }_{\mathbf{p}}\boldsymbol{ \cdot }\left \{ \left [A(\mathbf{q})+S(\mathbf{q})\right ]\boldsymbol{ \cdot }\frac{\mathbf{p}}{m}+{\nabla }_{\mathbf{q}}\phi (\mathbf{q})+\epsilon S(\mathbf{q})\boldsymbol{ \cdot }{\nabla }_{\mathbf{p}}\right \} \rho -\frac{\mathbf{p}}{m}\boldsymbol{ \cdot }{\nabla }_{\mathbf{q}}\rho . \end{eqnarray} \tag{ 23 }$$

Note that equation (23) has a solution (if the partition function Z() is finite)

$$\begin{eqnarray} &&{\rho }_{s}(\mathbf{q},\mathbf{p})=\frac{1}{Z(\epsilon )}\exp \left (-\frac{{p}^{2}/2 m+\phi (\mathbf{q})}{\epsilon }\right ). \end{eqnarray} \tag{ 24 }$$

Equation (24) is not explicitly time dependent, ∂_tρ_s = 0, it is a steady state distribution. The state variables q and p are separated, which simplifies the process of taking zero mass limit to be discussed below.

2.2.2. Zero mass limit.

We have mentioned that after taking zero mass limit m → 0, equation (13) reduces to equation (3). Therefore, by acting the same limit process on equation (23), the corresponding FPE to equation (3) can be obtained [10],

$$\begin{eqnarray} &&{\partial }_{t}\rho (\mathbf{q},t)=\nabla \boldsymbol{\cdot }\left [D(\mathbf{q})+Q(\mathbf{q})\right ]\boldsymbol{ \cdot }[\epsilon \nabla +\nabla \phi (\mathbf{q})]\rho (\mathbf{q},t), \end{eqnarray} \tag{ 25 }$$

where D(q) is the diffusion matrix, Q(q) is defined in equation (10), ∇ ≡ ∇_q, and the derivative ∇ in ∇ϕ(q) does not operate on ρ(q,t). This FPE has a drift velocity  − [D(q) + Q(q)]∇ϕ(q) = f(q) and a generalized symmetric form ∇⋅[D(q) + Q(q)] ⋅∇ρ(q,t). The commonly defined probability current density j(q,t) = (j₁(q,t),...,j_n(q,t))^τ is

$$\begin{eqnarray} &&{j}_{i}(\mathbf{q},t)=({f}_{i}+\Delta {f}_{i})(\mathbf{q})\rho (\mathbf{q},t)-{\partial }_{j}\left [\epsilon {D}_{i j}(\mathbf{q})\rho (\mathbf{q},t)\right ], \end{eqnarray} \tag{ 26 }$$

where f(q) = (f₁(q),...,f_n(q))^τ, Δf_i(q) = ∂_j[D_ij(q) + Q_ij(q)] (see also section 3.2), with i = 1,...,n and j = 1,...,n, D_ij(q) the ith-row and jth-column element of the diffusion matrix D(q), and Q_ij(q) the element of Q(q) in equation (25). For steady state, equation (12), ∇⋅j_s(q) = 0. One can calculate that when Q = 0 then j_s = 0; if Q(q) ≠ 0 then generally j_s(q) ≠ 0, since ∂_j[Q_ij(q)ρ_s(q)] ≠ 0. Therefore the framework encompasses the cases without detailed balance. The term 'detailed balance' used here is defined for abstract phase space, that is, the net current between any two states in the phase space is zero [34], identical to that for the Markov process in mathematics.

There are various ways to take the zero mass limit (eliminate the fast degrees p). A derivation adopting the standard projection operator method [1] is given in [10], that is, to project the (q,p) space to the q space. We note that another commonly used approach of degree reduction, adiabatic elimination [1, 4, 35], is based on the postulation that the velocity is a fast variable, which is consistent with the zero mass limit. Three other derivations for one-dimensional systems are discussed in [11]. This FPE (equation (25)) defines the A-type integration for an SDE of the form of equation (3). Since generally equation (1) is equivalent to equation (3), equation (25) defines a stochastic integration for general Langevin equations. The steady state distribution in the phase space (if it exists) is the Boltzmann–Gibbs distribution, equation (12). We should clarify that this is not to say that A-type is the only way to integrate equation (3), but it is a stochastic interpretation with physical significance.

We can consider a time evolution system described by an SDE as a particle motion inside the phase space. From a physical point of view, the motion of an object should be driven by some underlying forces. We consider this particle as a charged and massless one driven by the force ${\mathbf{F}}_{\mathrm{total}}=m\ddot{\mathbf{q}}=0$ (m = 0). The driving force can be separated generally in physics into three parts, dissipative, conservative and random

$$\begin{eqnarray} &&{\mathbf{F}}_{\mathrm{total}}={\mathbf{F}}_{\mathrm{dissipative}}+{\mathbf{F}}_{\mathrm{conservative}}+{\mathbf{F}}_{\mathrm{random}}=0. \end{eqnarray} \tag{ 27 }$$

Equation (3) has exactly this structure. In the usual study of two- and three-dimensional cases, equation (3) corresponds to a known fundamental equation in physics. We can generally use a frictional force as the dissipative component: ${\mathbf{F}}_{\mathrm{dissipative}}=-S(\mathbf{q})\dot {\mathbf{q}}$, where S(q) is symmetric and positive semi-definite; and a Lorentz force together with an energy induced force (for example, the electrostatic force) as the conservative component, ${\mathbf{F}}_{\mathrm{conservative}}=e\dot {\mathbf{q}}\times \mathbf{B}(\mathbf{q})+[-\nabla \phi (\mathbf{q})]$; the random component is described by a Gaussian white noise ${\mathbf{F}}_{\mathrm{random}}=\hat {N}(\mathbf{q})\xi (t)$ with zero mean that has a common origin with the frictional force, as formulated by the fluctuation-dissipation theorem [36]–[38], $\hat {N}(\mathbf{q}){\hat {N}}^{\tau }(\mathbf{q})=2 \epsilon S(\mathbf{q})$, hence

$$\begin{eqnarray} &&-S(\mathbf{q})\dot {\mathbf{q}}+e\dot {\mathbf{q}}\times \mathbf{B}(\mathbf{q})-\nabla \phi (\mathbf{q})+\hat {N}(\mathbf{q})\xi (t)=0. \end{eqnarray} \tag{ 28 }$$

The friction matrix S(q) is positive semi-definite (guaranteed by the fluctuation-dissipation theorem), keeping the resistant property of the dissipative force (the allowed values are confined in the negative half space). It can be a non-diagonal matrix, describing an anisotropic frictional force. The potential function ϕ(q) plays a similar role to the Hamiltonian in dissipative systems, leading to a Boltzmann–Gibbs distribution on the final steady state (if it exists) of the stochastic process. For higher dimension, we generalize $\mathbf{B}(\mathbf{q})\times \dot {\mathbf{q}}$ as $A(\mathbf{q})\dot {\mathbf{q}}$, where A(q) is an antisymmetric matrix. There is a direct and simple physical realization based on this intuitive explanation of the SDE (3); for instance, practical experiments to verify the generalized Einstein relation were designed in section 5.3 of [13]. A recent experiment [18] has been implemented under a similar physical setting.

We have mentioned that the A-type integration reduces to an α-type one with α = 1 for a one-dimensional system, hence it is different from Itô or Stratonovich's integration. A more detailed demonstration of this point is provided in [11]. For general cases, when the dimension is higher than one, things become more complicated: A-type integration is not even α-type, which means that at different time intervals the choices of α are not the same, or the parameter α is position dependent, α(q). A relation between A-type integration and Itô integration is formulated in [16] by considering the correspondence of their FPEs; that is, the SDE

$$\begin{eqnarray} &&\dot {\mathbf{q}}=(\mathbf{f}+\Delta \mathbf{f})(\mathbf{q})+N(\mathbf{q})\xi (t), \end{eqnarray} \tag{ 29 }$$

where Δf(q) = (Δf₁(q),...,Δf_n(q))^τ, using Itô integration describes the same process as that of the SDE

$$\begin{eqnarray} &&\dot {\mathbf{q}}=\mathbf{f}(\mathbf{q})+N(\mathbf{q})\xi (t) \end{eqnarray} \tag{ 30 }$$

using A-type integration. Here Δf_i(q) = ∂_j[D_ij(q) + Q_ij(q)]. We note that there is an obvious distinction between Δf(q) here and the additional (to Itô's) drift term for α-type (α chosen as a constant) integrations $\Delta {f}_{i}^{(\alpha )}(\mathbf{q})=\alpha [{\partial }_{k}{N}_{i j}(\mathbf{q})]{N}_{k j}(\mathbf{q})$.

An illustrative example has been provided in [16] to show a significant difference in the final steady state distributions between Itô and A-type integrations, the dynamics $\dot {\mathbf{q}}=\mathbf{f}(\mathbf{q})+N(\mathbf{q})\xi (t)$ is

$$\begin{eqnarray} && \dot {x}=2 x-x\left ({x}^{2}+{y}^{2}\right )+\sqrt{{x}^{2}+{y}^{2}}{\xi }_{x}(t),\\ && \dot {y}=2 y-y\left ({x}^{2}+{y}^{2}\right )+\sqrt{{x}^{2}+{y}^{2}}{\xi }_{y}(t), \end{eqnarray} \tag{ 31 }$$

with q = (x,y)^τ, $\epsilon =\frac{1}{2}$, the diffusion matrix D(x,y) = (x² + y²)I₂, ξ(t) = (ξ_x,ξ_y)^τ has zero mean, and 〈ξ(t)ξ^τ(t')〉 = δ(t − t')I₂. The final steady state distribution of this system with A-type integration can be calculated based on the framework as ρ_sA(x,y) = (1/Z_A())exp[ − ϕ/] with

$$\begin{eqnarray} &&\phi (\mathbf{q})=-\ln \left ({x}^{2}+{y}^{2}\right )+\frac{{x}^{2}+{y}^{2}}{2}. \end{eqnarray} \tag{ 32 }$$

The distribution of Itô integration for this system, f = f' + Δf, is identical to the A-type distribution for the dynamics f' (in order to apply the result provided in equations (29) and (30), we decompose the original system f (equations (31)) into the form f' + Δf)

$$\begin{eqnarray} && \dot {x}=-x\left ({x}^{2}+{y}^{2}\right )+\sqrt{{x}^{2}+{y}^{2}}{\xi }_{x}(t),\\ && \dot {y}=-y\left ({x}^{2}+{y}^{2}\right )+\sqrt{{x}^{2}+{y}^{2}}{\xi }_{y}(t), \end{eqnarray} \tag{ 33 }$$

whose expression can be similarly calculated as ρ_sI(x,y) = (1/Z_I())exp[ − ψ/] with

$$\begin{eqnarray} &&\psi (\mathbf{q})=\frac{1}{2}\left ({x}^{2}+{y}^{2}\right ). \end{eqnarray} \tag{ 34 }$$

The two distributions ρ_sA and ρ_sI have obvious differences; for instance, ρ_sA(0,0) = 0 but the origin (0,0) is the most probable state for ρ_sI. These theoretical results have been verified by numerical experiments in [16] (shown in figure 1). Note that the potential ϕ(q) in equation (32) serves as a global Lyapunov function [17] for the deterministic dynamics (when  = 0) of the system equations (31). The Lie derivative [39] of ϕ(q) (∇ϕ(q) ⋅f(q)) stays non-positive. Therefore, there is an exact correspondence between the deterministic dynamics and the steady state distribution derived using A-type integration for SDEs; for example, the stable fixed points are locally most probable states, in accordance with one's intuition. This direct correspondence cannot be kept by applying Itô or Stratonovich integration: after using Itô integration, the unstable point (0,0) of equations (31) ( = 0) becomes stable in ψ(q); moreover, stable points on x² + y² = 2 in equations (31) ( = 0) disappear in equation (34). The unique advantage of using A-type integration compared to Itô or Stratonovich integration is shown here. It enables a straightforward calculation of the transition probability from one fixed point to another after taking into account the noise influence for ODE models which is essential in many applications [5, 7, 8, 12, 14]. For Itô or Stratonovich integrations, this is not direct and can even be impossible, since one cannot recognize the original stable fixed points from the long time sampling distribution, like what is shown in this example.

Zoom In Zoom Out Reset image size

A comprehensive discussion on one-dimensional systems has been given previously [11]. The corresponding FPEs and the steady state distributions for three different types of stochastic integration (contained in [11]) are listed in table 1. As has been shown in [40], there exist real processes choosing Stratonovich integration. A recent experiment [18] on a one-dimensional physical process suggested that A-type integration of the corresponding SDE is consistent with the experimental data. We show in the following that the zero mass limit is well established based on the physical setting of the experiment.

Table 1.  Fokker–Planck equations and steady state distributions with  = 1.

α	Fokker–Planck equation	Steady state distribution
0	$\frac{\partial }{\partial t}\rho (q,t)=\left [\frac{\partial }{\partial q}\frac{\partial }{\partial q}D(q)+\frac{\partial }{\partial q}D(q)\phi (q)\right ]\rho (q,t)$	$\frac{1}{{Z}_{I}}\left \{ \frac{1}{D(q)}\exp \left [-\phi (q)\right ]\right \}$
$\frac{1}{2}$	$\frac{\partial }{\partial t}\rho (q,t)=\left [\frac{\partial }{\partial q}{D}^{1/2}(q)\frac{\partial }{\partial q}{D}^{1/2}(q)+\frac{\partial }{\partial q}D(q)\phi (q)\right ]\rho (q,t)$	$\frac{1}{{Z}_{S}}\left \{ \frac{1}{{D}^{\frac{1}{2}}(q)}\exp \left [-\phi (q)\right ]\right \}$
1	$\frac{\partial }{\partial t}\rho (q,t)=\left [\frac{\partial }{\partial q}D(q)\frac{\partial }{\partial q}+\frac{\partial }{\partial q}D(q)\phi (q)\right ]\rho (q,t)$	$\frac{1}{{Z}_{A}}\exp \left [-\phi (q)\right ]$

For the colloidal particle studied in [18], ρ = 1510 kg m⁻³, R = 655 × 10⁻⁹ m, the volume V = 4/3πR³ = 4/3 × π × (655×10⁻⁹)³ m³ = 1.18 × 10⁻¹⁸ m³ and the mass of the particle is m = ρV = 1510 × 1.18 × 10⁻¹⁸ kg = 1.78 × 10⁻¹⁵ kg. The symbol dt denotes the sampling time interval in [18]; within dt ≤ 10 ms, the authors claim that 'the force acting on the particle can be treated as locally constant'. In different experiments, they choose dt around the magnitude of 1 ms.

We use the simplified formula for the friction coefficient from Mannella et al's discussion of the same experiment [4], with η = 8.5 mPa s and z₀ = 700 nm,

$$\begin{eqnarray} &&\gamma (z)=6 \pi \eta R\frac{(z+{z}_{0})}{z}. \end{eqnarray} \tag{ 35 }$$

For z → ∞, we have γ_∞ = 6πηR = 6π × 8.5 × 10⁻³ × 655 × 10⁻⁹ N s m⁻¹ = 1.1 × 10⁻⁷ N s m⁻¹, then the time needed for the particle to return to equilibrium after being acted on by an external force F(z) is t(z) = [F(z)/γ(z)]/[F(z)/m] = m/γ(z), t_∞ = m/γ_∞ = 1.78 × 10⁻¹⁵/1.1 × 10⁻⁷ s = 1.6 × 10⁻⁸ s ≪ 1 × 10⁻³ s = 1 ms ≈ dt. In figure 2 of [18], for instance, at z = 200 nm, γ(z) = 1.1 × 10⁻⁷ N s m⁻¹ × (200 + 700)/200 = 5.0 × 10⁻⁷ N s m⁻¹. Since γ(z) is always larger than γ_∞, then t(z) < t_∞ ≪ dt. We conclude that during the sampling time interval the particle has been in equilibrium, hence the mass can be considered as zero in the experiment, realizing the zero mass limit.

Some related previous attempts have been discussed in our former works [6, 9, 13]. In this subsection we provide a brief summary. No effort is made to give a complete list, but we hope to have picked some of the major ones.

To the best of our knowledge, the one-dimensional case of the present framework, α = 1 type integration, has been mathematically discussed by Wong and Zakai in the 1960s [41] (see also [9, 13]). Their work provided a method for α∈[0,1].

From a physics perspective Grabert et al [42] provided a method for nonlinear irreversible processes, which was discussed by Graham in a long review which we commented on in [6]. Their method has a restriction based on the singularity of the diffusion matrix D(q) (see equations 2.14–16 in [42]) and handles the phase variables separately according to a reversibility defined by this singularity. By contrast, our framework can be used directly for a given D(q) in its general form without a separated treatment, a global potential is then obtained in the whole phase space.

Theoretical considerations for the situations with detailed balance have been proposed by Klimontovich [43] (see also [13]), while our approach applies to situations both with and without detailed balance. In the weak noise limit, previous works have reached the same Hamilton–Jacobi equation as our equation (11), the nondifferentiable potentials reported in [44, 45] are then included within our framework; moreover, our potential function is valid for an arbitrary noise strength. A major difference between our construction of the potential and those in the literature such as the Graham–Haken construction [46] is that no assumption on the stationary distribution function is needed in our approach in the limit t → ∞. In particular, the potential in our framework can be time dependent (see also [6, 20]). Other results beyond Itô and Stratonovich have been proposed in [47]–[49].