Anomalous fluxes in overdamped Brownian dynamics with Lorentz force

TL;DR

This paper investigates the effects of Lorentz force on overdamped Brownian particles, revealing that the standard overdamped approximation produces unphysical fluxes despite accurately capturing position statistics.

Contribution

It demonstrates that the overdamped limit introduces an additional drift term and unphysical fluxes, highlighting limitations of the common approximation in systems with Lorentz forces.

Findings

Overdamped equation captures position statistics accurately.

Unphysical fluxes persist in the long-time limit.

Analytical flux calculations match simulations.

Abstract

We study the stochastic motion of a particle subject to spatially varying Lorentz force in the small-mass limit. The limiting procedure yields an additional drift term in the overdamped equation that cannot be obtained by simply setting mass to zero in the velocity Langevin equation. We show that whereas the overdamped equation of motion accurately captures the position statistics of the particle, it leads to unphysical fluxes in the system that persist in the long time limit; an anomalous result inconsistent with thermal equilibrium. These fluxes are calculated analytically from the overdamped equation of motion and found to be in quantitative agreement with Brownian dynamics simulations. Our study suggests that the overdamped approximation, though perfectly suited for position statistics, can yield unphysical values for velocity-dependent variables such as flux and entropy production.

Equations

$$\dot{\mathbf{r}}(t)$$

$$m\dot{\mathbf{v}}(t)$$

$$\gamma \dot{\mathbf{r}}(t)=\mathbf{F}(\mathbf{r}(t))+\sqrt{2\gamma k_{B}T}\mathbf{\xi }(t).$$

$$\dot{\mathbf{r}}(t)$$

$$m\dot{\mathbf{v}}(t)$$

$$\mathbf{S}(\mathbf{r})=\frac{k_{B}T}{\gamma }\left(\mathbf{G}(\mathbf{r})\nabla \mathbf{G}(\mathbf{r})-\nabla {\mathbf{G}}_a(\mathbf{r})\right),$$

$$\mathbf{G}=1-\frac{\gamma qB(\mathbf{r})}{{\gamma }^2+q^2{B}^2(\mathbf{r})}\mathbf{M}+\frac{q^2{B}^2(\mathbf{r})}{{\gamma }^2+q^2{B}^2(\mathbf{r})}{\mathbf{M}}^2.$$

$$\dot{\mathbf{r}}(t)=\frac{k_{B}T}{\gamma }\left(\mathbf{G}(\mathbf{r})\nabla \mathbf{G}(\mathbf{r})-\nabla {\mathbf{G}}_a(\mathbf{r})\right)+\sqrt{\frac{2k_{B}T}{\gamma }}\mathbf{G}(\mathbf{r})\mathbf{\xi }(t).$$

$$\dot{\mathbf{r}}(t)=\sqrt{\frac{2k_{B}T}{\gamma }}\mathbf{G}\mathbf{\xi }(t).$$

$$\mathbf{J}(\mathbf{r},t)=-\frac{k_{B}T}{\gamma }\mathbf{G}\nabla Q(\mathbf{r},t),$$

$$\bar{\mathbf{J}}(\mathbf{r},t)=⟨\dot{\mathbf{r}}(t){\delta }^{(3)}(\mathbf{r}(t)-\mathbf{r})⟩,$$

$$⟨{\xi }_{\alpha }(t)\mathbf{R}[\mathbf{\xi }]⟩=\int ds⟨{\xi }_{\alpha }(t){\xi }_{\beta }(s)⟩⟨\frac{\delta \mathbf{R}[\mathbf{\xi }]}{\delta {\xi }_{\beta }(s)}⟩,$$

$$\bar{\mathbf{J}}(\mathbf{r},t)=-\frac{k_{B}T}{\gamma }\left(\nabla {\mathbf{G}}_a(\mathbf{r})\bar{Q}(\mathbf{r},t)+{\mathbf{G}}_s(\mathbf{r})\nabla \bar{Q}(\mathbf{r},t)\right),$$

$${\bar{J}}_x(y)={\rho }_b{k}_{B}T\frac{\partial B}{\partial y}\frac{{\gamma }^2-(qB{)}^2}{({\gamma }^2+(qB{)}^2{)}^2},$$

$$\frac{\partial }{\partial t}Q(\mathbf{r},t)=\frac{k_{B}T}{\gamma }\nabla \cdot \left(\mathbf{G}(\mathbf{r})\nabla Q(\mathbf{r},t)\right),$$

$$\mathbf{J}(\mathbf{r},t)$$

$$\dot{\mathbf{r}}(t)=\frac{k_{B}T}{\gamma }\left(\mathbf{G}\nabla \mathbf{G}-\nabla {\mathbf{G}}_a\right)+\sqrt{\frac{2k_{B}T}{\gamma }}\mathbf{G}(\mathbf{r}(t))\mathbf{\xi }(t).$$

$${\mathbf{J}}^{(1)}(\mathbf{r})$$

$$=\frac{k_{B}T}{\gamma }\left(\mathbf{G}(\mathbf{r})\nabla \mathbf{G}(\mathbf{r})-\nabla {\mathbf{G}}_a(\mathbf{r})\right)\bar{Q}(\mathbf{r},t),$$

$$\frac{\delta }{\delta {\xi }_{\alpha }(t)}{\int }_0^{t}ds{G}_{\nu \beta }(\mathbf{r}(s)){\xi }_{\beta }(s)=\frac{1}{2}{G}_{\nu \alpha }(\mathbf{r}(t)).$$

$${\bar{J}}_{\alpha }^{(2)}(\mathbf{r},t)$$

$$=\sqrt{\frac{2k_{B}T}{\gamma }}⟨\frac{\delta }{\delta {\xi }_{\beta }(t)}\left(G_{\alpha \beta }(\mathbf{r}(t)){\delta }^{(3)}(\mathbf{r}(t)-\mathbf{r})\right)⟩$$

$$=\sqrt{\frac{2k_{B}T}{\gamma }}⟨\frac{\delta r_{\nu }(t)}{\delta {\xi }_{\beta }(t)}\frac{\partial }{\partial r_{\nu }(t)}\left(G_{\alpha \beta }(\mathbf{r}(t)){\delta }^{(3)}(\mathbf{r}(t)-\mathbf{r})\right)⟩$$

$$=\frac{k_{B}T}{\gamma }⟨G_{\nu \beta }(\mathbf{r}(t))\left({\delta }^{(3)}(\mathbf{r}(t)-\mathbf{r})\frac{\partial }{\partial r_{\nu }(t)}{G}_{\alpha \beta }(\mathbf{r}(t))+G_{\alpha \beta }(\mathbf{r}(t))\frac{\partial }{\partial r_{\nu }(t)}{\delta }^{(3)}(\mathbf{r}(t)-\mathbf{r})\right)⟩$$

$$=\frac{k_{B}T}{\gamma }\left[\left(G_{\nu \beta }(\mathbf{r}){\partial }_{\nu }{G}_{\alpha \beta }(\mathbf{r})\right)\bar{Q}(\mathbf{r},t)-{\partial }_{\nu }\left(G_{\nu \beta }(\mathbf{r})G_{\alpha \beta }(\mathbf{r})\bar{Q}(\mathbf{r},t)\right)\right]$$

$$=-\frac{k_{B}T}{\gamma }\left[(G_{\alpha \beta }(\mathbf{r}){\partial }_{\nu }{G}_{\nu \beta }(\mathbf{r}))\bar{Q}(\mathbf{r},t)+G_{\alpha \beta }(\mathbf{r})G_{\nu \beta }(\mathbf{r}){\partial }_{\nu }\bar{Q}(\mathbf{r},t\right]$$

$${\bar{\mathbf{J}}}^{(2)}(\mathbf{r},t)=-\frac{k_{B}T}{\gamma }\left[\left(\mathbf{G}(\mathbf{r})\nabla \mathbf{G}(\mathbf{r})\right)\bar{Q}(\mathbf{r},t)+{\mathbf{G}}_s(\mathbf{r})\nabla \bar{Q}(\mathbf{r},t)\right]$$

$$\bar{\mathbf{J}}(\mathbf{r},t)=-\frac{k_{B}T}{\gamma }\left(\nabla {\mathbf{G}}_a(\mathbf{r})\bar{Q}(\mathbf{r},t)+{\mathbf{G}}_s(\mathbf{r})\nabla \bar{Q}(\mathbf{r},t)\right),$$

$$\frac{\partial }{\partial t}P(t)=(L_{rev}+L_{irr})P(t),$$

$$L_{rev}P(t)=-\mathbf{v}\cdot \nabla P(t)+\frac{q}{m}B(\mathbf{r}){\nabla }_{\mathbf{v}}\cdot \left[\mathbf{M}\mathbf{v}P(t)\right]$$

$$L_{irr}P(t)=\frac{\gamma }{m}{\nabla }_{\mathbf{v}}\cdot \left[\mathbf{v}P(t)+\frac{k_{B}T}{m}{\nabla }_{\mathbf{v}}P(t)\right].$$

$$\frac{\partial }{\partial t}\bar{P}(t)=\left({\bar{L}}_{rev}+{\bar{L}}_{irr}\right)\bar{P}(t),$$

$$\bar{P}(t)=P(t)R(\mathbf{v}{)}^{-1/2}$$

Full text

active colloids, phase separation, wetting

I Introduction

The motion of a particle suspended in a solvent can be modeled using the Langevin equation approach Langevin (1908). In this approach one writes an equation of motion for a particle, in which its interaction with the other degrees of freedom of the system (solvent) is modeled in terms of a stochastic force with suitable statistical properties. For instance, in absence of hydrodynamics the dynamics of a Brownian particle of mass $m$ can be described by the Langevin equation for its position $\boldsymbol{r}$ and velocity $\boldsymbol{v}$:

[TABLE]

where $\boldsymbol{F}(\boldsymbol{r})$ is an external force, $\gamma$ is a friction coefficient, $k_{B}$ is the Boltzmann constant and $T$ is the temperature. The noise $\boldsymbol{\xi}(t)$ is Gaussian with zero mean and time correlation $\langle\boldsymbol{\xi}(t)\boldsymbol{\xi}^{T}(t^{\prime})\rangle=\boldsymbol{1}\delta(t-t^{\prime})$. The importance of the Langevin equation approach lies in its applicability to a wide class of nonequilibrium problems Balakrishnan (2008).

The velocity correlations decay on a time scale $\tau=m/\gamma$, which implies that for times $t>>\tau$, the inertia term $m\dot{\boldsymbol{v}}(t)$ is negligible and can be set to zero to obtain an effective equation of motion for $\boldsymbol{r}$ as

[TABLE]

This equation, referred to as the overdamped equation of motion, is extensively used in theoretical and computational studies of nonequilibrium problems in which the correlation time $\tau$ is much smaller than the time scale of diffusion the particle Balakrishnan (2008); Gardiner (2009). This decoupling of velocity and position on time scales larger than $\tau$ makes it easier to find analytical solutions and has the advantage of significantly faster numerical computation. In fact, it has become a common practise to start with the overdamped equation of motion of the particle (Eq. (2)) as the model of the nonequilibrium system under study Brader and Krüger (2011); Cates and Tailleur (2013); Farage et al. (2015); Sharma and Brader (2016); Stenhammar et al. (2016); Sharma et al. (2017); Sharma and Brader (2017); Vuijk et al. (2018a).

The overdamped equation of motion is generally obtained in a simple way: set the inertia term to zero in the velocity Langevin equation and rearrange to describe the dynamics of the slow position variable. However, this procedure does not always yield the correct overdamped equation of motion, for instance when the noise is position dependent. In this case one must follow a systematic limiting procedure ($m\rightarrow 0$) of Eq. (1) to obtain the appropriate overdamped equation Hottovy et al. (2015); Volpe and Wehr (2016). This procedure yields an an additional drift term often referred to as the noise-induced drift in the literatureLau and Lubensky (2007); Hottovy et al. (2015); Volpe and Wehr (2016); Vuijk et al. (2018a). This additional drift term is absent if one simply sets $m=0$ in the velocity Langevin equation.

Additional drift also appears in the overdamped Langevin equation when the friction coefficient is position dependent, or more precisely, when the coefficient multiplying the velocity is position dependent. Hottovy et al. (2015); Volpe and Wehr (2016). One particularly interesting case, which is the main focus of this paper, is that of a Brownian particle subject to Lorentz force due to spatially varying magnetic field. The Lorentz force acting on a particle can be written as an antisymmetric matrix acting on $\boldsymbol{v}$, which, when added to the friction term $-\gamma\boldsymbol{v}$, results in an equation with position dependent coefficient in front of $\boldsymbol{v}$. Lorentz force is distinct from other nonconservative forces (e.g. shear) which input energy to the system. Shear forces can drive a system out of equilibrium resulting in nonequilibrium steady states. This stands in contrast to Lorentz force. Although Lorentz force generates particle currents, these are purely rotational and do no work on the system, which is thus not driven out of equilibrium. Being in equilibrium, such a system has (a) a stationary density profile given by the Boltzmann distribution and (b) no fluxes.

In this paper, we show that whereas the overdamped equation of motion for a Brownian particle in a spatially varying magnetic field accurately captures the position statistics, it leads to unphysical fluxes in the system. We first obtain the overdamped equation from the velocity Langevin equation using existing methods and show that the trajectory from the velocity Langevin equation converges on the trajectory from the overdamped equation with decreasing mass. We then show that for a particular choice of the spatially varying magnetic field, the overdamped equation fails to satisfy the no flux condition in equilibrium. This anomalous behaviour of the overdamped equation is the main result of this paper. We calculate these unphysical fluxes analytically from the overdamped equation, perform Brownian dynamics simulations of the overdamped equation of motion, measure the fluxes and show that they agree with the analytical predictions.

The anomalous behaviour does not invalidate the overdamped approximation. Rather, it is a manifestation of the subtle nature of the limiting procedure that yields the overdamped equation. The overdamped equation accurately captures the statistics of the position of the particle over finite time intervals Hottovy et al. (2015). This is seen clearly when one considers the Fokker-Planck equation for the position variable. The Fokker-Planck equation obtained from the overdamped equation is the same as that obtained from an independent alternative route. The computing of flux, however, involves taking the small-mass limit of velocity dependent terms which may result in additional terms. For instance, it has been shown that the overdamped equation does not yield the correct entropy production in the presence of a temperature gradient Celani et al. (2012); Marino et al. (2016); Birrell (2018).

The paper is organized as follows. In Sec. II, we briefly describe how the overdamped Langevin equation is obtained for a charged particle in spatially varying magnetic field. In Sec. III, we consider the special case of uniform magnetic field and demonstrate the existence of unusual curl-like fluxes. In Sec. IV, we show analytically and numerically that the overdamped equation obtained in Sec. II leads to unphysical fluxes in the system. The Fokker-Planck equation for the position variable is derived in Sec. V. Finally we present conclusions and brief discussion in Sec. VI.

II Langevin equation

We consider a single charged Brownian particle in a magnetic field $\boldsymbol{B}(\boldsymbol{r})$. The state of the particle is determined by the position vector $\boldsymbol{r}$ and velocity $\boldsymbol{v}$. Omitting hydrodynamic interactions, the dynamics of the particle are described by the following Langevin equation:

[TABLE]

where $m$ is the mass of the particle, $q$ is the charge, $k_{B}$ is the Boltzmann constant, $T$ is the temperature and $\boldsymbol{\xi}(t)$ is Gaussian white noise with zero mean and time correlation $\langle\boldsymbol{\xi}(t)\boldsymbol{\xi}^{T}(t^{\prime})\rangle=\boldsymbol{1}\delta(t-t^{\prime})$. Let $\boldsymbol{n}$ be the unit vector in the direction of the magnetic field, and $B(\boldsymbol{r})$ be the magnitude (i.e., $\boldsymbol{B}(\boldsymbol{r})=B(\boldsymbol{r})\boldsymbol{n}$). We define a matrix $\boldsymbol{M}$ with elements given by $M_{\alpha\beta}=-\epsilon_{\alpha\beta\nu}n_{\nu}$, where $\epsilon_{\alpha\beta\nu}$ is the totally antisymmetric Levi-Civita symbol in three dimensions and $n_{\nu}$ is $\nu$-component of $\boldsymbol{n}$ for the Cartesian index $\nu$. The Lorentz force can be written as $qB(\boldsymbol{r})\boldsymbol{M}\boldsymbol{v}$. One can rewrite the equation in terms of the position dependent matrix $\boldsymbol{\Gamma}(\boldsymbol{r})=(\gamma\boldsymbol{1}+qB(\boldsymbol{r})\boldsymbol{M})$ as

[TABLE]

When one is only interested in the slow degree of freedom (i.e. the position of the particle), the simulations are generally performed using an overdamped equation of motion. This equation of motion is obtained by taking the small-mass limit of Eq. (4). The limiting procedure is mathematically involved and is described in detail in Hottovy et al. (2015). The additional drift term of the overdamped equation corresponding to Eq. (4) is

[TABLE]

where $\boldsymbol{G}_{a}=(\boldsymbol{G}-\boldsymbol{G}^{T})/2$ is the antisymmetric part of $\boldsymbol{G}$, and $\boldsymbol{G}\equiv\gamma\boldsymbol{\Gamma}^{-1}$ which is

[TABLE]

It is important to note that the drift term (Eq. (5)) depends on whether the overdamped equation is interpreted in It$\rm\hat{o}$ or Stratonovich sense Lau and Lubensky (2007); Gardiner (2009). Equation (5) gives the additional drift in the Stratonovich interpretation of the overdamped equation which is given as

[TABLE]

The small-mass limit of the Langevin Eq. (4) involves a subtle limiting procedure which may be appreciated by noting that (a) the integration of Langevin Eq. (4) is independent of the interpretation (It$\rm\hat{o}$ or Stratonovich) whereas that of the overdamped equation (7) is not and (b) the term $\nabla\boldsymbol{G}_{a}(\boldsymbol{r})$ cannot be eliminated by choosing a different integration calculus; that is, this term is independent of the sense in which the overdamped equation is interpreted.

The following conventions are followed throughout the article: $\partial_{\alpha}$ stands for $\partial/\partial r_{\alpha}$, where $r_{\alpha}$ is the $\alpha$-component of $\boldsymbol{r}$. The $\alpha$-component of $\nabla\boldsymbol{G}$ is given as $(\nabla\boldsymbol{G})_{\alpha}=\partial_{\beta}G_{\beta\alpha}$, where repeated index is summed over. Similarly $(\boldsymbol{G}\nabla\boldsymbol{G})_{\alpha}=G_{\alpha\nu}\partial_{\beta}G_{\beta\nu}$.

III Uniform magnetic field

We first consider the case of uniform magnetic field. The overdamped Langevin equation can be obtained from Eq. (7) by setting $\nabla\boldsymbol{G}=0$ as

[TABLE]

In a recent study Chun et al. (2018), Chun et. al studied the small-mass limit of the Langevin equation (3) and showed explicitly that the resulting overdamped equation is not given by Eq. (8). They first calculated the noise correlation matrix from the Langevin equation (3) for a finite mass and then took the small-mass limit to obtain the correlation matrix of the noise in the overdamped equation. This procedure yielded the surprising result that the noise appearing in the overdamped equation of motion is a nonwhite Gaussian noise. This is in sharp contrast with the overdamped equation (8) which has a white Gaussian noise.

It was also shown in Ref. Chun et al. (2018) that the flux $\boldsymbol{J}(\boldsymbol{r},t)$, obtained from the correct overdamped equation, is

[TABLE]

where $Q(\boldsymbol{r},t)$ is the probability density. This flux is unusual because the matrix $\boldsymbol{G}$ cannot be interpreted as the diffusion matrix: it is not symmetric whereas a diffusion matrix is always symmetric. The flux can be written as sum of two terms: $\boldsymbol{G}_{s}\nabla Q(t)$, which we call the diffusive flux determined by the symmetric part $\boldsymbol{G}_{s}=(\boldsymbol{G}+\boldsymbol{G}^{T})/2$ of $\boldsymbol{G}$ and $\boldsymbol{G}_{a}\nabla Q(t)$, which we refer to as the curl flux determined by the antisymmetric part of $\boldsymbol{G}$. We note that Eq. (8) cannot give rise to curl flux. From Eq. (8), one only obtains the diffusive flux.

The unusual flux is a consequence of the nonwhite noise that appears in the overdamped equation Chun et al. (2018). Chun et. al studied how the nonwhite noise impacts dissipation in a system subject to nonconservative force that couples to the position. However, a direct demonstration of these fluxes was not presented. We show below that the curl flux can be measured in numerical simulations by starting with a nonequilibrium density distribution and measuring the fluxes, which arise from the density gradient. Ideally, this would be done using the overdamped equation with the nonwhite noise reported in Ref. Chun et al. (2018); however, at present it is not known how to generate the nonwhite noise that appears in this equation. Therefore, we demonstrate the presence of curl flux by numerically integrating the Langevin equation (3) with a small mass. We consider noninteracting particles that are initially uniformly distributed in the region $1\leq x,y\leq 3$ with $z=0$. We then numerically integrate the Langevin equation (3) with mass $m=2\times 10^{-3}$, $\boldsymbol{B}=5\hat{\boldsymbol{z}}$, and integration step $dt=5\times 10^{-6}$. Throughout the article we have used $k_{B}T=1$, $\gamma=1$, and $q=1$. The density distribution and flux are shown in Fig. 1 at time $t=0.05$. The velocity autocorrelation time, $m/\gamma=2\times 10^{-3}$, is much shorter than this time. As can be seen in Fig. 1 (a), density distribution becomes nonzero in the neighbourhood of the square region. The change in the distribution is due to the diffusive flux of the particles which is perpendicular to the edges of the square region. In addition to the diffusive flux there is also curl flux, which is shown in Figs. 1 (b) and (c). This flux, which is along the edges of the square region, is divergence free and therefore does not influence the time evolution of the density.

That the flux has a curl like component has also been reported in Refs. Kwon et al. (2005); Wang and Yuan (2016). However, the flux was obtained following the Fokker-Planck approach (shown below in Sec. V) which does not require the overdamped Langevin equation.

IV Inhomogeneous magnetic field

The overdamped motion of a charged particle in an inhomogeneous magnetic field has been studied in the past Cerrai and Freidlin (2011); Freidlin and Weber (2012); Hottovy et al. (2015). In this case, the overdamped Langevin equation (7) of a particle has an additional drift term. We show below that whereas this equation accurately describes the position of the particle and therefore the correct density distribution in the long-time limit, there are fluxes in steady state. We take the following approach: we compare the trajectory of the particle obtained from Eq. (3) with a small mass to the trajectory obtained from integrating Eq. (7). With decreasing mass the trajectories should converge. Figure 2(a,b) shows trajectories obtained from Eq. (3) with different masses and a trajectory obtained from Eq. (7), with magnetic field $\boldsymbol{B}=8\sin(2\pi y/L)\hat{\boldsymbol{z}}$, where $L$ is the size of the simulation box. We apply periodic boundary conditions in all directions. Figure 2(c,d) shows the comparison of trajectories in presence of a harmonic potential $(y-5)^{2}$. The magnetic field is $\boldsymbol{B}=1.6(y-5)\hat{\boldsymbol{z}}$. As can be seen in the Fig. 2, the trajectory from Eq. (3) seems to converge on the trajectory from Eq. (7) with decreasing mass.

Past studies have also relied on the comparison of trajectories to establish the accuracy of the overdamped Langevin equation of motion Hottovy et al. (2012); Wang and Yuan (2016); Volpe and Wehr (2016). This would seem to be a perfectly reasonable approach to establish the validity of the overdamped equation. If the two trajectories are matching, the overdamped equation of motion is accurately capturing the dynamics of the position of the particle. However, despite the matching trajectories, the two equations yield different particle fluxes in steady state. The flux obtained from the Langevin equation (3) with a small mass is identically zero at every spatial location; however, the flux obtained from the overdamped Langevin equation is nonzero in steady state; see Fig. 3. From an equilibrium thermodynamics standpoint, the steady state should be characterized by a Boltzmann probability density with no net fluxes. The Langevin equation (3) is consistent with thermodynamic equilibrium whereas the overdamped equation (7) is not. Nevertheless, as we discuss below, this inconsistency with equilibrium does not invalidate the overdamped stochastic differential equation.

The steady state flux can be obtained analytically by evaluating

[TABLE]

where $\boldsymbol{r}(t)$ is the position of the particle at time $t$ and $\dot{\boldsymbol{r}}(t)$ is the velocity, which is given by Eq. (7). In order to avoid any confusion, we clarify that $\boldsymbol{r}(t)$ is denoting the position of the particle and $\boldsymbol{r}$ is the position in space at which the fux is calculated. The flux can be calculated by substituting Eq. (7) for $\dot{\boldsymbol{r}}(t)$ in Eq. (10). The term containing $\langle\boldsymbol{\xi}(t)\delta^{(3)}(\boldsymbol{r}(t)-\boldsymbol{r})\rangle$ can be evaluated using the Novikov identity Novikov (1965)

[TABLE]

where $\boldsymbol{\xi}$ is Gaussian noise, $\alpha,\beta$ denote the $x,y$ or $z$ component, and $\boldsymbol{R}[\boldsymbol{\xi}]$ is a functional of the noise. The details of the calculation are shown in the appendix A. The final expression for the flux is

[TABLE]

where $\bar{Q}(\boldsymbol{r},t)$ is the probability density of the particle corresponding to the equation (7). The second term in the expression for the flux is a diffusive flux with the position-dependent diffusion coefficient $k_{B}T\boldsymbol{G}_{s}$, where we have used $\boldsymbol{G}_{s}=\boldsymbol{G}\boldsymbol{G}^{T}$. We consider the long time limit in which the probability density is homogeneously distributed implying that the diffusive flux is identically zero. It follows from Eq. (12) that for the particular choice of the magnetic field there should be a flux in the $x$-direction and no fluxes in the other directions. The $x$-component of the flux, obtained from Eq. (12), is

[TABLE]

where $\rho_{b}$ is the bulk probability density. Clearly, the numerically obtained flux is in excellent agreement with the analytical prediction (see Fig. 3).

That there are fluxes in the system is clearly inconsistent with thermal equilibrium. However, it does not invalidate the overdamped approximation. Rather, it is a manifestation of the subtle nature of the limiting procedure that yields the overdamped equation. When dealing with singular limits of equations, one can not speak broadly about the correct limiting equation in an absolute sense because the correct limit depends on the observable that one wishes to study with the limiting equation. In the case of the overdamped equation (7), the observable is the position of the particle Hottovy et al. (2015). It is clear from Fig. 2 that indeed the overdamped equation accurately captures the position statistics of the particle. The computing of flux, however, involves taking the small mass-limit of velocity dependent terms which may result in additional contributions in the limiting procedure which are not present in the overdamped equation. Another observable which involves taking the small-mass limit of velocity dependent terms is entropy production. Previous studies have shown that the overdamped equation does not yield the correct entropy production in the presence of a temperature gradient Celani et al. (2012); Marino et al. (2016); Birrell (2018).

Since the uniform magnetic field in Sec. III is a special limit of the inhomogeneous magnetic field, we believe that nonwhite noise would also emerge in the case of spatially varying magnetic field. It would be ideal to obtain the correlation matrix of the noise for a spatially varying magnetic field following the same approach as in Chun et al. (2018). However, at present our efforts have not been successful.

V Fokker Planck Equation

The Fokker-Planck equation corresponding to the overdamped equation (7) can be derived using standard methods Gardiner (2009) and is given as

[TABLE]

where $\boldsymbol{Q}(\boldsymbol{r},t)$ is the probability density of finding the particle at location $\boldsymbol{r}$ at time $t$. Alternatively, the Fokker-Planck equation can be obtained by an independent route which does not require the knowledge of the overdamped equation. The derivation, using the method described in Ref. Risken (1989), is presented in the appendix B and yields the same equation as in Eq. (14).

That one obtains the same Fokker-Planck equation establishes the validity of the overdamped equation (7) as the correct description of the position of the particle. The advantage of this alternative method of deriving the Fokker-Planck equation is that it also yields an expression for the flux $\boldsymbol{J}(\boldsymbol{r},t)$ in the small-mass limit:

[TABLE]

Though the two fluxes, Eqs. (15) and (12) have different dependence on $\boldsymbol{G}$ and $\boldsymbol{Q}$, their divergence is the same which is why both yield the same Fokker-Planck equation. The steady state density distribution can be obtained from Eq. (14) as a uniform distribution. Consistent with thermal equilibrium, the flux in Eq. (15) is identically zero for a uniformly distributed density. This is in contrast to the predictions of Eq. (12) which yields finite fluxes in a uniformly distributed system.

Note that the flux has exactly the same form as in Eq. (9) but with position dependent $\boldsymbol{G}$. It may seem that one can read off the expression for flux from Eq. (14) by casting the Fokker-Planck equation in the form of a continuity equation $\partial Q/\partial t+\nabla\cdot\boldsymbol{J}=0$. Though this approach yields the correct flux in most of the cases, there can be exceptions where it would not work. For instance if the flux has a constant divergence-free part, which would leave the Fokker-Planck equation unchanged, one cannot uniquely determine the flux from the Fokker-Planck equation alone. This is clearly seen in the case of uniform magnetic field: the Fokker-Planck equation, which is given as $\partial Q/\partial t=k_{B}T/\gamma\nabla\cdot(\boldsymbol{G}_{s}\nabla Q)$, remains unchanged due to the divergence-free flux $\boldsymbol{G}_{a}\nabla Q(\boldsymbol{r},t)$ (see Fig. 1).

VI Discussion and conclusions

In this paper we studied the motion of a Brownian particle subject to Lorentz force in the small-mass limit. We specifically considered the case in which the Lorentz force is position dependent; that is, the applied magnetic field is spatially varying. Spatially varying Lorentz force manifests itself as a position dependent coefficient in the Langevin equation for the velocity variable. One cannot then simply set the mass of the particle to zero to obtain the overdamped equation of motion Hottovy et al. (2015). When the coefficient multiplying the velocity is position dependent, the small-mass limit of the Langevin equation yields an overdamped equation of motion which has an additional drift term that depends on the gradient of the coefficient. Using existing techniques, we obtained the overdamped Langevin equation of motion of the particle with the additional drift term. We compared the trajectory obtained from the overdamped equation of motion with the trajectory from the velocity Langevin equation in the limit of small mass. We found that whereas the overdamped equation of motion accurately captures the position statistics of the particle, it leads to unphysical fluxes in the system.

That there are unphysical fluxes in the system is clearly inconsistent with thermal equilibrium. However, it does not invalidate the overdamped equation (7). The subtle limiting procedure used to obtain the overdamped equation ensures that the statistics of the position observable are accurately captured Hottovy et al. (2015). However, it is not suitable for studying velocity-dependent observables such as flux. Previously, It has been shown that the overdamped equation does not yield the correct entropy production in the presence of a temperature gradient Celani et al. (2012); Marino et al. (2016); Birrell (2018)

The Fokker-Planck equation for the position variable can be obtained by an independent route which does not require the overdamped equation. We find that the Fokker-Planck equation obtained from this method is the same as that from the overdamped equation of motion. This establishes the validity of the overdamped equation as the corect description of the position of the particle. The flux entering the Fokker-Planck equation (Eq. (14)) is unusual in the sense that a density gradient gives rise not only to a flux parallel to it (diffusive) but also perpendicular to it (curl like). The unusual form of the flux in Eq. (15) was most recently reported in Ref. Chun et al. (2018) in which the authors obtained the overdamped Langevin equation of motion for a Brownian particle in a uniform magnetic field. The authors elegantly demonstrated that this equation has nonwhite noise whose correlation matrix has antisymmetric components. Unfortunately, it is presently not known how to generate such a noise process. We have not been successful to obtain the correlation matrix of the noise in the case of spatially varying magnetic field.

Although the unusual form of flux has been previously reported, an unambiguous demonstration of such a flux using numerical simulations has been lacking. By numerically integrating the Langevin Eq. (3) with a small mass, we measured the flux directly and confirmed the theoretical predictions. When the magnetic field is uniform, the curl flux is divergence free and does not affect the time evolution of the probability distribution. By only retaining the diffusive flux, the resulting Fokker-Planck equation has a (diffusion) tensor which is symmetric. However, in a spatially varying magnetic field, the curl flux is not divergence free. Therefore, one has to retain the full tensor in the Fokker-Planck equation. This tensor can not be regarded as a diffusion tensor due to the antisymmetric components.

The Fokker-Planck equation often serves as the starting point for theoretical description of nonequilibrium problems such as spinodal decomposition Archer and Evans (2004); Vuijk et al. (2018b), linear response Sharma and Brader (2016, 2017); Merlitz et al. (2018), and first passage time problems Vuijk et al. (2019). It will be very interesting to investigate how the presence of these unusual curl like fluxes affects the dynamics of these phenomena.

Appendix A: Calculation of flux

Here we calculate the flux resulting from the overdamped equation of motion (7) using the Novikov relation Novikov (1965). We denote the position of the particle at time $t$ as $\boldsymbol{r}(t)$ to distinguish it from the spatial position $\boldsymbol{r}$ at which we calculate the flux. The Stratonovich stochastic differential equation for the position $\boldsymbol{r}(t)$ is given as (Eq. (7))

[TABLE]

The flux is calculated using $\bar{\boldsymbol{J}}(\boldsymbol{r},t)=\left\langle\dot{\boldsymbol{r}}(t)\delta^{(3)}(\boldsymbol{r}(t)-\boldsymbol{r})\right\rangle=\bar{\boldsymbol{J}}^{(1)}(\boldsymbol{r})+\bar{\boldsymbol{J}}^{(2)}(\boldsymbol{r})$, where $\bar{\boldsymbol{J}}^{(1)}$ is the contribution to the flux from the deterministic part of the equation for $\dot{\boldsymbol{r}}(t)$ and $\bar{\boldsymbol{J}}^{(2)}$ from the stochastic part. $\bar{\boldsymbol{J}}^{(1)}$ can be calculated in a straightforward fashion as

[TABLE]

where $\bar{Q}(\boldsymbol{r},t)=\left\langle\delta^{(3)}(\boldsymbol{r}(t)-\boldsymbol{r})\right\rangle$ is the probability density at $\boldsymbol{r}$.

The calculation of $\bar{\boldsymbol{J}}^{(2)}$ uses the Novikov relation and is presented below. In the derivation below we have used the following Fox (1986):

[TABLE]

We calculate the flux component wise. The $\alpha$-component of the flux $\bar{\boldsymbol{J}}^{(2)}$ can be written as

[TABLE]

where Eq. (18) is used in the fourth step of the derivation. Equation (19) can be cast in vector notation as

[TABLE]

Adding Eqs. (17) and (19), we get

[TABLE]

Appendix B: Fokker-Planck derivation

It follows exactly from the Langevin equation (3) that the probability distribution $P(t)\equiv P(\boldsymbol{r},\boldsymbol{v},t)$ evolves in time according to Gardiner (2009)

[TABLE]

where the time-evolution operator has been split up in a reversible part

[TABLE]

and an irreversible part

[TABLE]

To derive the Fokker-Planck equation for the position of the particle, we follow the method described in Ref. Chun et al. (2018); Risken (1989). We first recast Fokker-Planck equation equation (22) as

[TABLE]

where

[TABLE]

and

[TABLE]

where $L$ can be either of the operators in Eq. (22), and

[TABLE]

is the solution to $L_{irr}R(\boldsymbol{v})=0$, normalized such that the integral over $\boldsymbol{v}$ is one. The transformed operators are

[TABLE]

[TABLE]

where

[TABLE]

The eigenfunctions of the operator $b^{\dagger}_{\alpha}b_{\alpha}$, where $\alpha$ is either $x$,$y$ or $z$, are

[TABLE]

and

[TABLE]

where $H_{n}$ are Hermite polynomials. The operators $b_{\alpha}^{\dagger}$ and $b_{\alpha}$ are the raising and lowering operators of the eigenfunctions: $b_{\alpha}^{\dagger}\psi_{n}(v_{\alpha})=\sqrt{n+1}\psi_{n+1}(v_{\alpha})$ and $b_{\alpha}\psi_{n}(v_{\alpha})=\sqrt{n}\psi_{n-1}(v_{\alpha})$. The eigenfunctions are orthonormal,

[TABLE]

and can be used to expand $\bar{P}(t)$:

[TABLE]

where $c_{n_{x},n_{y},n_{z}}=c_{n_{x},n_{y},n_{z}}(\boldsymbol{r},\boldsymbol{p},t)$.

Without loss of generality, the magnetic field is oriented along the $z$ direction and $\boldsymbol{B}(\boldsymbol{r})=B(\boldsymbol{r})\hat{\boldsymbol{z}}$. Equation (25) together with the orthonormality of the eigenfunctions yields an hierarchy of equations for the functions $c_{n_{x},n_{y},n_{z}}$ called a Brinkman hierarchy Brinkman (1956):

[TABLE]

where $\boldsymbol{D}=\sqrt{\frac{k_{B}T}{m}}\nabla$.

The probability density for the position and orientation, $Q(t)\equiv Q(\boldsymbol{r},t)$, is given by the first expansion coefficient:

[TABLE]

The order of the coefficient functions $c_{n_{x},n_{y},n_{z}}=\mathcal{O}(m^{\frac{1}{2}(n_{x}+n_{y}+n_{z})})$ and up to leading order in $m$, $\partial_{t}c_{n_{x},n_{y},n_{z}}=0$ for $n_{x}+n_{y}+n_{z}>0$. Up to leading order in $m$ Eq. (Appendix B: Fokker-Planck derivation) is closed and can now be written as

[TABLE]

and

[TABLE]

where

[TABLE]

The matrix $\Gamma$ is the sum of $\gamma\boldsymbol{1}$ and the cross product with $\boldsymbol{B}(\boldsymbol{r})$, where in this case $\boldsymbol{B}(\boldsymbol{r})=B(\boldsymbol{r})\hat{\boldsymbol{z}}$. In the general case of magnetic field as $\boldsymbol{B}(\boldsymbol{r})=B(\boldsymbol{r})\boldsymbol{n}$ where $\boldsymbol{n}$ is a unit vector, the friction matrix is given as $\boldsymbol{\Gamma}(\boldsymbol{r})=(\gamma\boldsymbol{1}+qB(\boldsymbol{r})\boldsymbol{M})$. The elements of the matrix $\boldsymbol{M}$ are given as $M_{\alpha\beta}=-\epsilon_{\alpha\beta\nu}n_{\nu}$, where $\epsilon_{\alpha\beta\nu}$ is the totally antisymmetric Levi-Civita symbol in three dimensions and $n_{\nu}$ is $\nu$-component of $\boldsymbol{n}$ for the Cartesian index $\nu$.

The flux in position space is

[TABLE]

which can be calculated by using Eq. (36) and $v_{\alpha}\psi_{0}(v_{\alpha})=\sqrt{\frac{T}{m}}\psi_{1}(v_{\alpha})$:

[TABLE]

So the equation for the probability density $Q(t)\equiv Q(\boldsymbol{r},t)$ is

[TABLE]

where $\Gamma^{-1}(\boldsymbol{r})$ is given by Eq. (6).