# A perturbative theory for Brownian vortexes

Henrique W. Moyses, Ross O. Bauer, Alexander Y. Grosberg and David G. Grier
Department of Physics and Center for Soft Matter Research, New York University, New York, NY 10003
April 18, 2015
###### Abstract

Brownian vortexes are stochastic machines that use static non-conservative force fields to bias random thermal fluctuations into steadily circulating currents ((1); (2)). The archetype for this class of systems is a colloidal sphere in an optical tweezer ((3); (1); (4)). Trapped near the focus of a strongly converging beam of light, the particle is displaced by random thermal kicks into the nonconservative part of the optical force field arising from radiation pressure ((5)), which then biases its diffusion ((3); (1)). Assuming the particle remains localized within the trap, its time-averaged trajectory traces out a toroidal vortex. Unlike trivial Brownian vortexes, such as the biased Brownian pendulum, which circulate preferentially in the direction of the bias, the general Brownian vortex can change direction and even topology in response to temperature changes. Here we introduce a theory based on a perturbative expansion of the Fokker-Planck equation for weak non-conservative driving. The first-order solution takes the form of a modified Boltzmann relation and accounts for the rich phenomenology observed in experiments on micrometer-scale colloidal spheres in optical tweezers.

# I Introduction

Stochastic systems can be driven out of equilibrium in various ways. Thermal ratchets and Brownian motors use time-dependent forces to rectify Brownian motion ((6); (7); (8); (9); (10)). Other systems evolve in response to spatial and temporal variations in the temperature ((11); (10)). Even if none of these systems can reach thermodynamic equilibrium, some can achieve non-equilibrium steady states. Recently, a distinct class of stochastic machines was discovered that use noise to extract work from a static non-conservative force field ((12); (1); (13); (14)). These systems, which have been dubbed Brownian vortexes ((1); (2)), perform no work at all in the absence of stochastic forces. When activated by noise, they enter into steady-state motion characterized by toroidal vortexes in the time-averaged or ensemble-averaged probability current. These cyclic processes in principle can be coupled to external systems to extract useful work. Since their initial observation in the fluctuations of optically trapped colloidal beads ((12); (1)), Brownian vortexes have been reported in trapped colloidal rods ((13)), in trapped spheres subjected to shear flows ((14); (15)), in the response of bacterial populations to antibiotics ((16)), and in models for population dynamics ((17)).

Brownian vortexes arise in time-independent force fields that include at least one point of stable equilibrium. In the absence of thermal forces, such systems remain stationary at their fixed points and so perform no work. They differ in this respect from conventional machines that move deterministically under the influence of non-conservative forces. Random thermal forces allow a Brownian vortex to explore its force landscape. Were the force field purely conservative, the system would reach thermal equilibrium in the Boltzmann distribution, and so would have no means to perform work. In force fields with solenoidal components, however, the probability distribution can be advected by the non-conservative force. Probability currents then flow through the system under the competing influences of advection and diffusion. These currents must form closed cycles if the system’s overall probability is conserved, raising the possibility that the system can reach steady state.

Initial studies have identified two broad categories of Brownian vortex behavior: trivial Brownian vortexes that circulate in the direction dictated by the non-conservative component of the force landscape, and general Brownian vortexes that select their own topology and circulation. Both types of behavior have been explained with master equations on discrete networks ((2)). Solutions for continuous systems have been obtained for special cases ((2); (17); (16); (18)), with most examples representing trivial Brownian vortexes. Here, we introduce a perturbative theory that accounts for both trivial and general cases in continuous systems.

# II Theory

Our approach is based on the observation that systems governed by linear forces can be mapped to a generalization of the Ornstein-Uhlenbeck process, for which a general solution is known ((19)). For nonlinear forces, a common approach is to solve the system using perturbation theory ((20)). Here we develop a perturbation theory for the specific case where the non-conservative force is weak compared to the conservative one. We show that this approach captures the important features of general Brownian vortex circulation, including topological transitions and flux reversal.

Following earlier studies ((1); (2); (19); (20)), we describe a system capable of undergoing steady-state stochastic circulation as a Brownian particle of mobility $\mu$ moving through a static force landscape $\mathbf{F}(\mathbf{r})$ at absolute temperature $T$. Such a description applies naturally to the motion of a colloidal particle in an optical force field, and also may be applied to more general systems such as an ensemble of identical particles interacting through conservative forces and confined by a force landscape. Assuming that the system reaches steady state, we seek the steady-state probability distribution $\rho(\mathbf{r})$ describing the likelihood of finding the particle near position $\mathbf{r}$, and the associated steady-state flux of probability $\mathbf{j}(\mathbf{r})$ flowing through that point. The flux is generated both by advection of the probability distribution and by diffusion:

 $\mathbf{j}(\mathbf{r})=\mu\rho(\mathbf{r})\mathbf{F}(\mathbf{r})-D\nabla\rho(% \mathbf{r}).$ (1)

Because the probability density is non-negative, we may express it as the exponential of an effective potential $\phi(\mathbf{r})$:

 $\rho(\mathbf{r})=e^{-\beta\phi(\mathbf{r})},$ (2)

so that

 $\mathbf{j}(\mathbf{r})=\mu\rho(\mathbf{r})\left[\mathbf{F}(\mathbf{r})+\nabla% \phi(\mathbf{r})\right].$ (3)

The Helmholtz decomposition theorem guarantees that the force field can be separated into a conservative gradient term and a non-conservative solenoidal component,

 $\mathbf{F}(\mathbf{r})=-\nabla U(\mathbf{r})+\nabla\times\mathbf{A}(\mathbf{r}),$ (4)

where $U(\mathbf{r})$ is the potential energy and $\mathbf{A}(\mathbf{r})$ is the vector potential. The probability current may be written in terms of these potentials as

 $\mathbf{j}(\mathbf{r})=\mu\rho(\mathbf{r})\,\left\{\nabla\times\mathbf{A}(% \mathbf{r})+\nabla\left[\phi(\mathbf{r})-U(\mathbf{r})\right]\right\}.$ (5)

Equation (5) makes clear that the system can only reach thermodynamic equilibrium with $\mathbf{j}_{0}(\mathbf{r})=0$ if the non-conservative force vanishes, $\nabla\times\mathbf{A}(\mathbf{r})=0$. The remaining term in Eq. (5) then yields Boltzmann’s distribution for the probability density,

 $\rho_{0}(\mathbf{r})=e^{-\beta U(\mathbf{r})},$ (6)

with $\phi_{0}(\mathbf{r})=\beta U(\mathbf{r})$, where $\beta^{-1}=k_{B}T$ is the thermal energy scale.

Systems subject to nonconservative forces may not reach thermodynamic equilibrium, but still must satisfy the continuity equation,

 $\nabla\cdot\mathbf{j}(\mathbf{r})=-\frac{\partial\rho(\mathbf{r},t)}{\partial t},$ (7)

where $\rho(\mathbf{r},t)$ is the time-dependent probability density. Equation (7) is the Fokker-Planck equation for Brownian vortexes. Any steady-state solution of Eq. (7) satisfies

 $\nabla\cdot\mathbf{j}(\mathbf{r})=0,$ (8)

which also ensures conservation of probability.

Quite remarkably, Eq. (8) implies that the steady state distribution $\rho(\mathbf{r})$ does not depend on the particle’s mobility, $\mu$, or on any other transport property. Indeed, substituting Eq. (5) for the current density into Eq. (8) yields an equation that is independent of $\mu$. This happens because of the assumption, implicit in Eq. (1), that the particle responds in the same way to an applied force whether or not it is conservative. For a potential force, the Einstein relation $D=k_{B}T\mu$ between the diffusion coefficient $D$ and mobility $\mu$ follows from the condition that there exists an equilibrium state which obeys the detailed balance condition $\mathbf{j}(\mathbf{r})=0$. In the presence of a non-conservative force, the system does not come to equilibrium, it does not satisfy detailed balance. Nevertheless, its probability distribution remains independent of the particle’s transport properties, and depends only on the form of the force field. In this sense, the steady-state probability distribution resembles Boltzmann’s distribution despite the system’s departure from equilibrium.

To find the specific steady-state probability distribution, we introduce the projection

 $p(\mathbf{r})=\nabla U(\mathbf{r})\cdot\nabla\times\mathbf{A}(\mathbf{r})$ (9)

of the non-conservative force onto the direction of the conservative force. An exact solution to Eqs. (3) and (8) is known ((2)) only for the special case $p(\mathbf{r})=0$. In this case, the steady-state probability distribution retains the form of the Boltzmann distribution, Eq. (6), and is simply advected by the nonconservative force:

 $\mathbf{j}(\mathbf{r})=\mu e^{-\beta U(\mathbf{r})}\nabla\times\mathbf{A}(% \mathbf{r}).$ (10)

The direction of $\mathbf{j}(\mathbf{r})$ in this case is fixed by $\nabla\times\mathbf{A}(\mathbf{r})$ regardless of the temperature. Equation (10) therefore describes a trivial Brownian vortex. Examples of trivial Brownian vortexes include the biased Brownian pendulum and colloidal spheres circulating in circularly polarized optical tweezers ((21)).

To move beyond this special case, we consider systems in which $\nabla\times\mathbf{A}(\mathbf{r})$ is not necessarily aligned with $\nabla U(\mathbf{r})$, but may be treated as a perturbation whose scale is characterized by a small parameter $\epsilon$. We therefore replace $\nabla\times\mathbf{A}(\mathbf{r})$ in Eq. (4) with $\epsilon\nabla\times\mathbf{A}(\mathbf{r})$. Assuming the perturbation not to be singular, the effective potential may be expanded in orders of $\epsilon$ as

 $\phi(\mathbf{r})=\phi_{0}(\mathbf{r})+\epsilon\phi_{1}(\mathbf{r})+\mathcal{O}% \!\left\{\epsilon^{2}\right\},$ (11)

with $\phi_{0}(\mathbf{r})=U(\mathbf{r})$. The associated expansion of the probability current,

 $\mathbf{j}(\mathbf{r})=\epsilon\mathbf{j}_{1}(\mathbf{r})+\mathcal{O}\!\left\{% \epsilon^{2}\right\},$ (12)

has no contribution at zero-th order in $\epsilon$. The first-order correction,

 $\mathbf{j}_{1}(\mathbf{r})=\mu e^{-\beta U(\mathbf{r})}\left[\nabla\times% \mathbf{A}(\mathbf{r})+\nabla\phi_{1}(\mathbf{r})\right]$ (13a) retains the advective term from Eq. (10), although this now distorts the probability distribution as well as transporting it. The second term in Eq. (13a) describes diffusive relaxation of that distortion. Because these two terms depend on temperature in different ways the first-order expansion admits the possibility of temperature-dependent transitions such as those characterizing general Brownian vortexes. These transitions, moreover, may be understood to arise from competition between advection and diffusion. The balance depends on the distortion of the probability distribution away from $\rho_{0}(\mathbf{r})$. As explained in Appendix A, substituting Eq. (13a) into Eq. (8) yields a differential equation for $\phi_{1}(\mathbf{r})$, $\nabla^{2}\phi_{1}(\mathbf{r})-\beta\nabla U(\mathbf{r})\cdot\nabla\phi_{1}(% \mathbf{r})=\beta p(\mathbf{r}).$ (13b) Solutions to Eq. (13b) are difficult to find for arbitrary force fields. Provided that $U(\mathbf{r})$ has at least one minimum, however, the associated operator $\hat{H}_{U}=\nabla^{2}-\beta\nabla U(\mathbf{r})\cdot\nabla$ has eigenfunctions $\psi_{n}(\mathbf{r})$ with eigenvalues $\lambda_{n}$ that satisfy $\hat{H}_{U}\psi_{n}(\mathbf{r})=\lambda_{n}\psi_{n}(\mathbf{r}).$ (13c) Symmetrized functions of the form $e^{-\frac{1}{2}\beta U(\mathbf{r})}\psi_{n}(\mathbf{r})$ constitute a complete orthogonal basis with the orthonormalization condition $\int e^{-\beta U(\mathbf{r})}\psi_{n}(\mathbf{r})\psi_{m}(\mathbf{r})\,d% \mathbf{r}=\delta_{nm}.$ (13d) We therefore can expand $\phi_{1}(\mathbf{r})$ in this basis, $\phi_{1}(\mathbf{r})=\sum_{n}c_{n}\psi_{n}(\mathbf{r})$ (13e) with expansion coefficients $c_{n}=\frac{\beta}{\lambda_{n}}\int e^{-\beta U(\mathbf{r})}\psi_{n}(\mathbf{r% })\,p(\mathbf{r})\,d\mathbf{r}.$ (13f)

Equations (13) are the principal result of this study. They describe the lowest-order extension of a trivial Brownian vortex into a general Brownian vortex under the influence of a weakly non-conservative force. The nature of the resulting topological transformations and flux reversals depends on details of $\mathbf{F}(\mathbf{r})$, and particularly on the projection of its solenoidal component onto its conservative part.

To illustrate practical applications of Eq. (13), we next apply this formalism to the particular case of colloidal spheres moving in optical tweezers, the context in which Brownian vortexes were first observed ((3); (1)). This analysis casts new light on the nature of this system’s behavior.

# III Colloids in an optical tweezer

To facilitate comparisons with experimental studies of colloidal spheres circulating in optical tweezers, we numerically compute the forces acting on a trapped sphere with the Lorenz-Mie theory of light scattering using methods described in Appendix B. Typical results are presented in Fig. 1. Colors in Fig. 1(a) represent the relative intensity of the light in an optical tweezer, as viewed in the $(r,z)$ plane in cylindrical coordinates. The light propagates in the $+\hat{z}$ direction (upward) and comes to a focus along the axis defined by $r=0$. Arrows indicate the direction and strength of the resulting force $\mathbf{F}(\mathbf{r})$ experienced by a colloidal sphere of radius $a_{p}=0.5~{}\mathrm{\upmu}\mathrm{m}$ at position $\mathbf{r}$ within that light field. We have shifted the origin of the coordinate system to coincide with the position of the trap so that $\mathbf{F}(0)=0$. Away from the origin, $\mathbf{F}(\mathbf{r})$ directs the particle back to the stable fixed point. All forces reported in Fig. 1 are computed for a power of 1 $\mathrm{m}$$\mathrm{W}$, which is a reasonable scale for typical optical trapping experiments.

Both the axial component of the force, plotted in Fig. 1(b), and the radial component plotted in Fig. 1(c) resemble a linear restoring force over a reasonably wide range of axial displacements. The axial component of the optical force shows a non-trivial dependence on radial position, as shown in Figs. 1(d) and (e), that is reasonably modeled as a quartic polynomial. This solenoidal dependence is responsible for Brownian vortex circulation in optically trapped colloidal spheres.

The three-dimensional force field is very nearly symmetric with respect to rotations about the $\hat{z}$ axis. Small distortions along the axis of the light’s polarization have little influence on the trapped particle’s motions, and will not be considered here.

# III.1 Simulated circulation

The data in Fig. 2(a) through 2(c) show results of Brownian dynamics simulations ((22)) of a colloidal silica sphere diffusing through water at $T=300\mathrm{K}$ in the computed force field from Fig. 1(a). The sphere’s trajectory $\mathbf{r}_{p}(t)$ is computed at 100 $\mathrm{\upmu}$$\mathrm{s}$ intervals for each of three values of the light’s intensity: 1.38 $\mathrm{m}$$\mathrm{W}$ in Fig. 2(a), 0.35 $\mathrm{m}$$\mathrm{W}$ in Fig. 2(b) and 0.17 $\mathrm{m}$$\mathrm{W}$ in Fig. 2(c). These correspond to values of the radial stiffness, $k$, of 1 $\mathrm{p}$$\mathrm{N}$$\mathrm{/}$$\mathrm{\upmu}$$\mathrm{m}$, 0.25 $\mathrm{p}$$\mathrm{N}$$\mathrm{/}$$\mathrm{\upmu}$$\mathrm{m}$ and 0.12 $\mathrm{p}$$\mathrm{N}$$\mathrm{/}$$\mathrm{\upmu}$$\mathrm{m}$, respectively. These values are chosen for consistency with previous experimental studies ((3); (1)). Each trajectory then is compiled into estimates for the probability density $\rho(\mathbf{r})$ and the probability current density $\mathbf{j}(\mathbf{r})$ using an adaptive kernel density estimator ((23); (3); (1)). Streamlines of $\mathbf{j}(\mathbf{r})$ are denoted by barbs whose sharp ends point in the direction of motion and whose length is proportional to the local current density. Each barb is colored by the estimate for $\rho(\mathbf{r})$.

At the highest intensity, Fig. 2(a), the trap is stiffest and the particle is most strongly confined. A single toroidal roll is evident in the probability current, and circulates downstream along the optical axis. Reducing the intensity weakens the trap and frees the particle to explore more of the force landscape. Under these conditions, shown in Fig. 2(b), streamlines of $\mathbf{j}(\mathbf{r})$ reveal a second counter-rotating roll. Reference ((1)) identifies the appearance of this second roll as a topological transition. Reducing the intensity still further allows the second, outer roll to dominate the system’s dynamics. It subsumes the inner roll so that only a single toroidal vortex remains. The single remaining roll circulates with the flux directed upstream along the optical axis, which signals a flux reversal relative to Fig. 2(a) in addition to the topological transition relative to Fig. 2(b). Precisely this behavior was reported in ((1)), which lends credence both to the experimental results and also to the present simulations.

The particle’s residence time in the trap is effectively indefinite at the highest laser intensity. This corresponds to the force field in Fig. 1 and the streamlines in Fig. 2(a). Reducing the laser intensity reduces the residence time, which falls to 10 s under the conditions in Fig. 2(c). This is consistent with results of experimental studies. Simulations under these conditions are restarted every time the particle escapes, and the results averaged to until the computed current density converges. The streamlines in Fig. 2(c) therefore should be viewed as an ensemble average.

# III.2 Analytical model

The force field $\mathbf{F}(\mathbf{r})$ resembles a cylindrically symmetric harmonic well for small excursions away from the equilibrium point. This can be seen in Figs. 1(b) and (c). We therefore model the trapping potential as

 $U(\mathbf{r})=\frac{1}{2}k\,(r^{2}+\eta\,z^{2}),$ (14)

where $\eta$ characterizes the trap’s anisotropy. In this case, Eq. (13c) has a complete set of eigenfunctions, given in cylindrical coordinates by

 $\psi_{\mathbf{n}}(\mathbf{r})=N_{\mathbf{n}}\,L_{n}\!\left(\frac{1}{2}\beta kr% ^{2}\right)H_{n_{z}}\!\left(\sqrt{\frac{1}{2}\eta\,\beta k}\,z\right),$ (15a) where $L_{n}(\cdot)$ is the Laguerre polynomial of index $n$, $H_{n_{z}}(\cdot)$ is the Hermite polynomial of index $n_{z}$ and where $\mathbf{n}=\{n,n_{z}\}$ is a set of whole-number indexes that are related to the associated eigenvalues by $\lambda_{\mathbf{n}}=-\beta k(2n+\eta\,n_{z}).$ (15b) The harmonic well’s basis functions are normalized by $N_{\mathbf{n}}=\eta^{\frac{1}{4}}\left(\frac{\beta k}{2\pi}\right)^{\frac{3}{4% }}\frac{1}{\sqrt{2^{n_{z}}n_{z}!}}.$ (15c)

We have chosen basis functions that are independent of $\theta$ to reflect the azimuthal symmetry of $\mathbf{F}(\mathbf{r})$.

# III.2.1 Quadratic perturbation in an isotropic trap

To illustrate applications of Eq. (13), we first consider the symmetric case, $\eta=1$, subjected to a quadratic perturbation,

 $F_{z}(r)=\epsilon ka_{p}\left(1-\frac{r^{2}}{a_{p}^{2}}\right),$ (16)

with $F_{x}(r)=F_{y}(r)=0$, where $\epsilon$ is the small parameter characterizing the strength of the non-conservative force and $a_{p}$ is the radius of the sphere. This perturbation is divergence-free and thus has the solenoidal form assumed in Eq. (4). Choosing $F_{z}(r)$ to be proportional to $k$ ensures that it scales with the light’s intensity in the same way as the linear restoring force. Scaling distances by the particle’s radius is reasonable because typical realizations of Brownian vortex circulation involve motions substantially smaller than $a_{p}$ ((3); (1)). In general, both $k$ and $\epsilon$ depend on the particle’s radius relative to the wavelength of light. Parameterized in this way, the perturbation may be considered weak if

 $\epsilon<\sqrt{\frac{1}{2}\beta ka_{p}^{2}}.$ (17)

Choosing a quadratic form for $F_{z}(r)$ suggests that the non-conservative force will dominate the linear restoring force at large distances, and that the particle eventually will escape from the trap. The rate at which probability leaks from the system is limited by the smallness of $\epsilon$. This is the case, for example, in the system described by Fig. 1. Under these conditions, the particle is likely to remain trapped over the period of observation, and probability may be considered to be conserved. This also is the case for the experimental studies of optically trapped colloidal particles ((3); (1); (5)) that Eqs. (14) and (16) are intended to model.

The projection factor for this system, defined in Eq. (9), has the form

 $p(\mathbf{r})=-k^{2}a_{p}z\left(1-\frac{r^{2}}{a_{p}^{2}}\right).$ (18)

It does not explicitly include $\epsilon$ because this parameter is introduced in Eq. (11) for the effective potential and in Eq. (12) for the current density. Substituting this form for $p(\mathbf{r})$ into Eq. (13) along with the basis functions from Eq. (15) yields the first-order correction to the effective potential,

 $\beta\phi_{1}(\mathbf{r})=\frac{1}{3}\frac{z}{a_{p}}\left[4-\beta ka_{p}^{2}% \left(3-\frac{r^{2}}{a_{p}^{2}}\right)\right].$ (19)

This correction’s linear dependence on $z$ displaces the probability distribution down the optical axis, as would be expected of radiation pressure. The distribution’s width

 $\sigma=\left.\sqrt{\frac{\int_{0}^{\infty}r^{3}\rho(\mathbf{r})dr}{\int_{0}^{% \infty}r\rho(\mathbf{r})dr}}\right|_{z=0}=\sqrt{\frac{2}{\beta k}}$ (20)

is the same as in the unperturbed case.

The associated current density,

 $\mathbf{j}(\mathbf{r})=\frac{2}{3}\epsilon De^{-\beta U(\mathbf{r})}\left[% \beta krz\hat{r}+(2-\beta kr^{2})\hat{z}\right],$ (21)

has the form of a toroidal roll, centered on the axis $r=0$ and circulating around a circular core in the plane $z=0$ at radius

 $r_{0}=\sqrt{\frac{2}{\beta k}},$ (22)

which coincides with the distribution’s width, $r_{0}=\sigma$. The circulation rate is controlled by the strength $\epsilon$ of the driving force and the rate at which the particle diffuses in the harmonic potential energy well given its diffusion coefficient $D=\mu k_{B}T$.

This Brownian vortex undergoes no topological transitions or flux reversals. Applying a quadratic perturbation to a harmonic well therefore constitutes a model for a trivial Brownian vortex.

# III.2.2 Quartic perturbation in an anisotropic trap

More interesting behavior arises in a more highly structured force field of the type actually observed in optical trapping experiments ((3); (1); (5)). Figure 1(e) shows that the computed axial force is well approximated near the plane $z=0$ by a quartic polynomial,

 $F_{z}(r)=-\epsilon\frac{2+\eta}{3}k\frac{r^{2}}{a_{p}}\left(1-p_{4}r^{2}\right),$ (23)

with $F_{x}(r)=F_{y}(r)=0$, where $p_{4}$ characterizes the particle’s interaction with the focused beam of light. Like $\epsilon$ and $k$, this additional parameter also depends on the sphere’s radius, $a_{p}$. As for the quadratic case, Eq. (23) describes a purely solenoidal perturbation. We will show that $p_{4}$ governs topological transitions and flux reversals in the resulting Brownian vortex circulation. Incorporating a constant offset into $F_{z}(r)$ would move the equilibrium position along the optical axis as in the quadratic case, but does not otherwise influence the system’s behavior. We have omitted such an offset from Eq. (23) for clarity. The weak-perturbation condition for this model is

 $\epsilon<\frac{3}{2+\eta}\sqrt{\frac{1}{2}\beta ka_{p}^{2}}.$ (24)

Accounting for the trap’s anisotropy $\eta$ allows for comparison with optical trapping experiments. The force field presented in Fig. 1, for example, is intended to model the influence of Gaussian beam brought to a diffraction-limited focus and is characterized by $\eta=0.47$. Results for isotropic harmonic wells can be retrieved by setting $\eta=1$. No qualitative features of Brownian vortex circulation depend on the value of $\eta$.

The projection factor, defined in Eq. (9), associated with this force field

 $p(\mathbf{r})=\frac{2+\eta}{3}\eta k^{2}\frac{z}{a_{p}}(r^{2}-p_{4}r^{4}),$ (25)

gives rise to a first-order correction to the effective potential that also has the form of a quartic polynomial:

 $\beta\phi_{1}(\mathbf{r})=\frac{2}{3}\frac{z}{a_{p}}\,\left[(1-2\tilde{p}_{4})% (2+\tilde{r}^{2})-\frac{2+\eta}{4\eta}\,\tilde{p}_{4}\tilde{r}^{4}\right].$ (26a) For conciseness we have introduced the dimensionless parameters $\displaystyle\tilde{r}^{2}=\frac{1}{2}\eta\beta k\,r^{2}\quad\text{and}$ (26b) $\displaystyle\tilde{p}_{4}=\frac{8}{(4+\eta)\beta k}\,p_{4}.$ (26c)

As in the quadratic case, the quartic perturbation displaces probability along $\hat{z}$. The quartic term also tends to broaden the probability distribution relative to $\rho_{0}(\mathbf{r})$ by introducing a super-Gaussian tail at large $r$. This broadening may be of interest for precision measurements of optical forces because it can influence ((3)) calibration protocols based on analysis of thermal fluctuations ((24); (25); (26)). The availability of an analytical form for $\phi_{1}(\mathbf{r})$ will help to assess when nonequilibrium redistribution of the probability density may be ignored ((27)).

The first-order correction to the current density,

 $\mathbf{j}_{1}(\mathbf{r})=j_{r}(\mathbf{r})\hat{r}+j_{z}(\mathbf{r})\hat{z},$ (27a) has the radial component $\frac{j_{r}(\mathbf{r})}{j_{0}(\mathbf{r})}=\eta\beta kzr\left[1-\tilde{p}_{4}% \left(2+\frac{2+\eta}{2\eta}\tilde{r}^{2}\right)\right]$ (27b) that is conveniently measured in units of an overall scale $j_{0}(\mathbf{r})=\frac{2}{3}\frac{D}{a_{p}}e^{-\beta U(\mathbf{r})}$ (27c) that again reflects diffusion in the potential energy well, as in Eq. (21). The radial component of the current density changes sign as it passes through the equatorial plane, $z=0$. The axial component, $\frac{j_{z}(\mathbf{r})}{j_{0}(\mathbf{r})}=2\tilde{p}_{4}\frac{2+\eta}{\eta^{% 2}}\,\tilde{r}^{4}-4\left(\frac{1}{\eta}+\tilde{p}_{4}\right)\,\tilde{r}^{2}+4% (1-2\tilde{p}_{4}),$ (27d)

varies non-monotonically with distance $r$ from the axis. The net current, $\mathbf{j}(\mathbf{r})=\epsilon\mathbf{j}_{1}(\mathbf{r})$, is proportional to the strength of the non-conservative driving force.

Streamlines of $\mathbf{j}(\mathbf{r})$ plotted in Fig. 2(d) through 2(f) reveal circulatory flows that agree well with the simulated results for the same system under the same conditions. Specifically, the first-order result for the strongly confined system in Fig. 2(d) shows a single roll consistent in position, extent and speed with the simulation result in Fig. 2(a). The double-roll structure in Fig. 2(e) similarly is consistent with that in Fig. 2(b), although quantitative details of the circulation differ, particularly near the optical axis. The single-roll structure in Fig. 2(f) again qualitatively resembles that in Fig. 2(c), although the centers of circulation appear at different radial positions. The first-order perturbation theory developed in Eq. (13) thus captures the essential features of general Brownian vortex circulation in this model system, albeit with quantitative discrepancies. Trends in the analytical results therefore offer useful insights into the nature of the phenomenon.

# III.3 Topological transition and flux reversal

The probability current, $\mathbf{j}(\mathbf{r})$, vanishes at the cores of toroidal vortexes. Solutions of $\mathbf{j}_{1}(\mathbf{r})=0$ take the form of circles in the plane $z=0$ centered on the axis at $r=0$. Two vortex cores exist if $p_{4}$ is sufficiently small, at radii $r_{+}$ and $r_{-}$ that satisfy

 $\displaystyle r_{\pm}^{2}$ $\displaystyle=r_{0}^{2}\,(1\pm\Delta^{2}),\quad\text{where}$ (28a) $\displaystyle r_{0}^{2}$ $\displaystyle=\frac{1}{\beta k}\frac{2}{2+\eta}\left(\frac{1}{\tilde{p}_{4}}+% \eta\right)\quad\text{and}$ (28b) $\displaystyle\Delta^{2}$ $\displaystyle=\frac{\sqrt{(8+4\eta+\eta^{2})\tilde{p}_{4}^{2}-4\tilde{p}_{4}+1% }}{1+\eta\tilde{p}_{4}}.$ (28c)

According to Eq. (28), the vortex at $r_{+}$ is present for all temperatures greater than zero and moves outward as the temperature increases. The other vortex at $r_{-}$ moves inward, and ceases to exist when it reaches $r_{-}=0$, which occurs when $\tilde{p}_{4}=\frac{1}{2}$. This condition corresponds to a threshold temperature

 $k_{B}T_{c}=\frac{4+\eta}{16}\,\frac{k}{p_{4}}$ (29)

below which the probability current consists of two counter-rotating toroidal vortexes and above which only a single roll remains. The threshold temperature is proportional to the light’s intensity through the trap stiffness, $k$, and depends on details of the particle-light interaction through $p_{4}$. Interestingly, it does not depend on the strength, $\epsilon$, of the non-conservative force. In most optical trapping experiments, the temperature remains constant while $k$ is varied. Equation (29) for the threshold temperature then may be recast into an equivalent condition for the trap stiffness.

The system’s behavior for $T corresponds to strong confinement by the harmonic well. The outer toroidal roll therefore may not be perceptible in an experiment of finite duration because the particle spends comparatively little time exploring the outermost reaches of the force landscape. The apparent transition from one vortex in Fig. 2(a) to two in Fig. 2(b) in fact results from increasing the occupation of the outer roll and does not constitute a topological transition. The appearance of a second concentric vortex in previous optical trapping studies ((1); (2)) almost certainly corresponds to the same statistical sampling considerations.

Increasing the temperature beyond $T_{c}$, or equivalently reducing the stiffness of the trapping potential, causes an actual topological transition by eliminating the inner roll at $r_{-}$. In this case, the single remaining roll circulates upstream along the optical axis. The behavior in Fig. 2(f) therefore reflects a flux reversal relative to Fig. 2(d) as well as a topological transition. This more dramatic change also has been observed in experimental studies ((1)).

# III.4 Size-dependent crossovers

The nature of the Brownian vortex behavior in an optically trapped colloidal sphere depends qualitatively on the form of the nonconservative force, $F_{z}(r)$. Neither the concentric roll structure nor the high-temperature topological transition appear if this force is quadratic rather than quartic. What dynamic patterns emerge therefore depends sensitively on the size and composition of the sphere and the structure of the light.

The data in Fig. 3 show how the controlling parameters $\epsilon$ and $p_{4}$ vary with sphere radius $a_{p}$ for a colloidal silica sphere trapped in water by an optical tweezer of vacuum wavelength 532 nm brought to a focus by a lens of numerical aperture 1.4. Figures 3(b), (c) and (d) show plots of $F_{z}(r)$ for the representative radii indicated by vertical dashed lines in Fig. 3(a). As in Fig. 1, these results were obtained numerically using Eq. (43).

In the Rayleigh limit, for $a_{p}\ll\lambda$, the non-conservative force is peaked at the optical axis and falls off very nearly quadratically in $r$. For such particles, $\epsilon>0$. Figure 3(b) is representative of this range of conditions. The quartic contribution being weak, Rayleigh particles enter into single-roll Brownian vortexes circulating down the optical axis in the direction of the light’s propagation. This is the mode of operation that first was identified in Ref. ((3)). The corresponding set conditions is shaded red in Fig. 3(a).

Particles that are much larger than the wavelength of light, $a_{p}\gg\lambda$, have force profiles $F_{z}(r)$ that are peaked far enough from the optical axis that they also appear to be quadratic over the accessible range, but with $\epsilon<0$. Such particles circulate in a single toroidal vortex, but in the sense opposite to that adopted by Rayleigh particles. The retrograde circulation of larger spheres was pointed out in a ray-optics analysis of optically trapped colloidal spheres ((28)) and subsequently was observed experimentally ((1)).

Retrograde circulation in a single roll also can arise for particles that are intermediate in size between the Rayleigh range and the ray-optics regime. The example force field depicted in Fig. 3(c) has this property, and is characterized by $\epsilon<0$. Even though the quartic term is sizable under these conditions, the curvature of $F_{z}(r)$ has the same sign over the entire accessible range of radii, and only a single toroidal roll can be populated. The domain of such behavior is shaded blue in Fig. 3(a).

Topological transition and flux reversals are only possible if $\partial_{r}F_{z}(r)$ changes sign in an accessible part of the force landscape. That occurs for conditions such as those in Fig. 3(d), and corresponds to the green-shaded regions in Fig. 3(a). In such cases, the inner roll circulates in the same sense as the single roll in the Rayleigh regime, and the outer roll rotates in the opposite sense. Interestingly, there appears to be no set of conditions that favor a retrograde double-roll or the corresponding topological transition to a retrograde single roll.

Different patterns of size-dependent crossovers arise for spheres of different materials, or in beams of different wavelengths or focusing properties. Figure 3 makes clear that small variations in properties can have a large influence on Brownian vortex circulation, not simply changing the rate of circulation, but rather reversing the direction of circulation. The present work has focused on the circulating currents’ topology. Its results also could be used to address questions about drift rate and circulation frequency that have been raised in previous studies ((3); (1)).

# IV Conclusion

Brownian vortexes should be generic features of all probability-conserving stochastic systems subject to time-independent driving by non-conservative forces. Equation (13) constitutes the leading-order description of a weakly-driven Brownian vortex. This level of approximation already captures the topological transitions and flux reversals that have been reported for archetypal Brownian vortex circulation in optically trapped colloids. Analytical results for topological transitions in the current density go beyond reproducing experimental results by clarifying their nature and providing a unified explanation for the system’s behavior. The idealized model of an optical tweezer as a harmonic well subject to a steady non-conservative force is likely to serve as a useful model other systems as well.

It is noteworthy that the steady-state probability distribution in a Brownian vortex,

 $\rho(\mathbf{r})=e^{-\beta\phi(\mathbf{r})},$ (30)

depends on characteristics of the force field $\mathbf{F}(\mathbf{r})$, but not on the particle’s mobility $\mu$. The particle’s transport properties therefore do not dictate how the probability redistributes as the system is driven out of equilibrium. In this sense, the general relationship between $\mathbf{F}(\mathbf{r})$ and the effective potential $\phi(\mathbf{r})$ described in Appendix A constitute an analog for the Boltzmann relation for this class of nonequilibrium systems.

The generic nature of the model discussed in Sec. III suggests that this might be a common factor for steady-state circulation in a broad range of nonequilibrium systems. The formalism developed here also should apply in other contexts such as circulatory and oscillatory flows in social networks, financial systems, and chemical networks. It also would be interesting to extend this formalism to more general systems such as systems with multiple fixed points, and systems of multiple interacting particles.

# V Acknowledgments

This work was supported by the National Science Foundation through Grant Number DMR-135875.

# Appendix A Perturbation Theory

The effective potential, $\phi(\mathbf{r})$, that determines the steady-state probability distribution, $\rho(\mathbf{r})$, differs from the imposed potential, $U(\mathbf{r})$, by an amount

 $\psi(\mathbf{r})=\phi(\mathbf{r})-U(\mathbf{r})$ (31)

that vanishes if $\nabla\times\mathbf{A}(\mathbf{r})=0$. We assume therefore that $\psi(\mathbf{r})$ is characterized by the same small parameter, $\epsilon$, that characterizes $\nabla\times\mathbf{A}(\mathbf{r})$. Imposing conservation of probability through Eq. (8) yields a differential equation for $\psi(\mathbf{r})$,

 $\left[\nabla^{2}\psi(\mathbf{r})-\beta\nabla\psi(\mathbf{r})\cdot\nabla U(% \mathbf{r})\right]\\ -\left[\beta\left|\nabla\psi(\mathbf{r})\right|^{2}+\beta\nabla\psi(\mathbf{r}% )\cdot\nabla\times\mathbf{A}(\mathbf{r})\right]\\ =\beta p(\mathbf{r}),$ (32)

that depends on the temperature and characteristics of the force field, but not on the diffusing particle’s mobility. To first order in $\epsilon$, this reduces to

 $\nabla^{2}\psi(\mathbf{r})-\beta\nabla\psi(\mathbf{r})\cdot\nabla U(\mathbf{r}% )=\beta p(\mathbf{r}).$ (33)

The associated field

 $\chi(\mathbf{r})=e^{-\frac{1}{2}\beta U(\mathbf{r})}\psi(\mathbf{r})$ (34)

then may be obtained as an expansion in eigenfunctions $\chi_{n}(\mathbf{r})$ of the Hermitian operator

 $\hat{H}^{\prime}_{U}=\nabla^{2}+\left[\frac{1}{2}\beta\nabla^{2}U(\mathbf{r})-% \frac{1}{4}\beta^{2}\left|\nabla U(\mathbf{r})\right|^{2}\right].$ (35)

Specifically, solutions of

 $\hat{H}^{\prime}_{U}\chi_{n}(\mathbf{r})=\lambda_{n}\chi_{n}(\mathbf{r})$ (36)

form a complete set of basis functions labeled by index $n$ with eigenvalues $\lambda_{n}$. When appropriately normalized they satisfy the orthogonality condition

 $\int\chi_{m}(\mathbf{r})\chi_{n}(\mathbf{r})\,d^{3}r=\delta_{mn}.$ (37)

Equation (13d) follows from Eq. (37).

Imposing conservation of probability at each order of the perturbation expansion ensures that solutions reflect steady state behavior. Higher-order corrections obtained by incorporating the second-order terms in Eq. (32) will redistribute the probability distribution, but are not likely to eliminate qualitative features of the steady-state circulation that arise at first order.

# Appendix B Generalized Lorenz-Mie theory

An optical tweezer can be modeled as a strongly focused Gaussian beam, and its field can be expanded as a series in vector spherical harmonics,

 $\mathbf{E}_{i}(\mathbf{r})=E_{0}\sum_{n=1}^{\infty}\sum_{m=-n}^{n}\left[a_{mn}% \mathbf{M}^{(1)}_{mn}(k\mathbf{r})+b_{mn}\mathbf{N}^{(1)}_{mn}(k\mathbf{r})% \right].$ (38)

The coefficients $a_{mn}$ and $b_{mn}$ have been previously reported for a converging Gaussian beam ((29); (30); (31)). and depend on the numerical aperture of the lens that brings the beam to a focus.

The incident beam illuminates a particle located at $\mathbf{r}_{p}$, which gives rise to a scattered wave $\mathbf{E}_{s}(\mathbf{r}-\mathbf{r}_{p})$ that propagates to position $\mathbf{r}$. The corresponding expansion for the scattered field,

 $\mathbf{E}_{s}(\mathbf{r})=E_{0}\sum_{n=1}^{\infty}\sum_{m=-n}^{n}\left[r_{mn}% \mathbf{M}^{(3)}_{mn}(k\mathbf{r})+s_{mn}\mathbf{N}^{(3)}_{mn}(k\mathbf{r})% \right],$ (39)

has expansion coefficients,

 $\begin{split}r_{mn}&=-b_{n}\ a_{mn}\ \ \text{and}\\ s_{mn}&=-a_{n}\ b_{mn},\end{split}$ (40)

that are related to the incident beam’s coefficients by the particle’s Lorenz-Mie scattering coefficients, $a_{n}$ and $b_{n}$ ((29); (30)). For the particular case of scattering by a sphere, the Lorenz-Mie coefficients depend on the sphere’s radius and refractive index relative to the medium.

The combined incident and scattered fields,

 $\mathbf{E}(\mathbf{r})=\mathbf{E}_{i}(\mathbf{r})+\mathbf{E}_{s}(\mathbf{r}-% \mathbf{r}_{p}),$ (41)

can be used to calculate the Maxwell stress tensor, $\mathbf{T}(\mathbf{r})$, with components

 $T_{ij}(\mathbf{r})=\epsilon_{m}E_{i}(\mathbf{r})E_{j}(\mathbf{r})-\frac{1}{2}% \epsilon_{m}E^{2}(\mathbf{r})\,\delta_{ij},$ (42)

where $\epsilon_{m}$ is the dielectric constant of the medium. The optically induced force is then obtained by integrating the Maxwell stress tensor over a surface $S$ that encloses the sphere,

 $\mathbf{F(\mathbf{r})}=\oint_{S}\hat{n}\cdot\mathbf{T}(\mathbf{r^{\prime}})\,d% \mathbf{r^{\prime}},$ (43)

where $\hat{n}$ is the unit vector normal to $S$. Techniques to perform this integration have been previously reported ((32); (33); (34); (35)).

Figure 1 shows the computed force field acting on a 0.5 $\mathrm{\upmu}$$\mathrm{m}$-radius silica sphere in water in an optical tweezer with vacuum wavelength $\lambda=532~{}\mathrm{n}\mathrm{m}$ projected by a lens with numerical aperture 1.4 (convergence angle 60${}^{\circ}$). We define the particle’s point of mechanical equilibrium to be the origin of the coordinate system in Fig. 1.

# References

• (1) B. Sun, J. Lin, E. Darby, A. Y. Grosberg and D. G. Grier. “Brownian vortexes.” Phys. Rev. E80, 010401(R) (2009).
• (2) B. Sun, D. G. Grier and A. Y. Grosberg. “Minimal model for Brownian vortexes.” Phys. Rev. E82, 021123 (2010).
• (3) Y. Roichman, B. Sun, A. Stolarski and D. G. Grier. “Influence of non-conservative optical forces on the dynamics of optically trapped colloidal spheres: The fountain of probability.” Phys. Rev. Lett.101, 128301 (2008).
• (4) M. de Messieres, N. A. Denesyuk and A. La Porta. “Noise associated with nonconservative forces in optical traps.” Phys. Rev. E84, 031108 (2011).
• (5) P. Y. Wu, R. X. Huang, C. Tischer, A. Jonas and E. L. Florin. “Direct measurement of the nonconservative force field generated by optical tweezers.” Phys. Rev. Lett.103, 108101 (2009).
• (6) R. D. Astumian and M. Bier. “Fluctuation driven ratchets – molecular motors.” Phys. Rev. Lett.72, 1766–1769 (1994).
• (7) J. Prost, J. F. Chauwin, L. Peliti and A. Ajdari. “Asymmetric pumping of particles.” Phys. Rev. Lett.72, 2652–2655 (1994).
• (8) H. Linke. “Ratchets and Brownian motors: Basics, experiments and applications.” Appl. Phys. A75, 167 (2002).
• (9) R. D. Astumian and P. Hänggi. “Brownian motors.” Physics Today55, 33–39 (2002).
• (10) P. Reimann. “Brownian motors: Noisy transport far from equilibrium.” Phys. Rep.361, 57–265 (2002).
• (11) P. Reimann, R. Bartussek, R. Haussler and P. Hänggi. “Brownian motors driven by temperature oscillations.” Phys. Lett. A215, 26–31 (1996).
• (12) Y. Roichman, B. Sun, Y. Roichman, J. Amato-Grill and D. G. Grier. “Optical forces arising from phase gradients.” Phys. Rev. Lett.100, 013602 (2008).
• (13) S. H. Simpson and S. Hanna. “First-order non-conservative motion of optically trapped nonspherical particles.” Phys. Rev. E82, 031141 (2010).
• (14) J. Bammert and W. Zimmermann. “Probability distribution of a trapped Brownian particle in plane shear flows.” Phys. Rev. E82, 052102 (2010).
• (15) M. Khan and A. K. Sood. “Tunable Brownian vortex at the interface.” Phys. Rev. E83, 041408 (2011).
• (16) S. Bradde and G. Biroli. “The generalized Arrhenius law in out of equilibrium systems.” (2012).
• (17) M. Parker, A. Kamenev and B. Meerson. “Noise-induced stabilization in population dynamics.” Phys. Rev. Lett.107, 180603 (2011).
• (18) D. I. Russell and R. A. Blythe. “Macroscopically observable probability currents in finite populations.” Journal of Statistical Mechanics: Theory and Experiment2013, P06008 (2013).
• (19) C. Kwon, P. Ao and D. J. Thouless. “Structure of stochastic dynamics near fixed points.” Proc. Nat. Acad. Sci.102, 13029–13033 (2005).
• (20) C. Kwon and P. Ao. “Nonequilibrium steady state of a stochastic system driven by a nonlinear drift force.” Phys. Rev. E84, 061106 (2011).
• (21) D. B. Ruffner and D. G. Grier. “Optical forces and torques in non-uniform beams of light.” Phys. Rev. Lett.108, 173602 (2012).
• (22) D. Frenkel and B. Smit. Understanding Molecular Simulation (Academic Press, San Diego, 2001).
• (23) B. W. Silverman. Density Estimation for Statistics and Data Analysis (Chapman & Hall, New York, 1992).
• (24) E. L. Florin, A. Pralle, E. H. K. Stelzer and J. K. H. Horber. “Photonic force microscope calibration by thermal noise analysis.” Appl. Phys. A66, S75–S78 (1998).
• (25) M. Polin, K. Ladavac, S.-H. Lee, Y. Roichman and D. G. Grier. “Optimized holographic optical traps.” Opt. Express13, 5831–5845 (2005).
• (26) B. Lukić, S. Jeney, Ž. Sviben, A. J. Kulik, E.-L. Florin and L. Forró. “Motion of a colloidal particle in an optical trap.” Phys. Rev. E76, 011112 (2007).
• (27) G. Pesce, G. Volpe, A. C. de Luca, G. Rusciano and G. Volpe. “Quantitative assessment of non-conservative radiation forces in an optical trap.” EPL86, 38002 (2009).
• (28) R. Huang, P. Wu and E.-L. Florin. “Comment on ”Influence of non-conservative optical forces on the dynamics of optically trapped colloidal spheres: The fountain of probability”, cond-mat.soft 0804.0730v1.”
• (29) C. F. Bohren and D. R. Huffman. Absorption and Scattering of Light by Small Particles (Wiley Interscience, New York, 1983).
• (30) M. I. Mishchenko, L. D. Travis and A. A. Lacis. Scattering, Absorption and Emission of Light by Small Particles (Cambridge University Press, Cambridge, 2001).
• (31) J. M. Taylor and G. D. Love. “Multipole expansion of Bessel and Gaussian beams for Mie scattering calculations.” J. Opt. Soc. Am. A26, 278–282 (2009).
• (32) S. Alexander, L. Hellemans, O. Marti, J. Schneir, V. Ellings, P. K. Hansma, M. Longmire and J. Gurley. “An atomic-resolution atomic-force microscope implemented using an optical lever.” J. Appl. Phys.65, 164–167 (1989).
• (33) T. A. Nieminen, V. L. Y. Loke, A. B. Stilgoe, G. Knoner, A. M. Branczyk, N. R. Heckenberg and H. Rubinsztein-Dunlop. “Optical tweezers computational toolbox.” J. Opt. A9, S196–S203 (2007).
• (34) B. Sun, Y. Roichman and D. G. Grier. “Theory of holographic optical trapping.” Opt. Express16, 15765–15776 (2008).
• (35) D. B. Ruffner and D. G. Grier. “Universal, strong and long-ranged trapping by optical conveyors.” Opt. Express22, 26834–26853 (2014).