The role of non-conservative scattering forces and damping on Brownian particles in optical traps

TL;DR

This paper investigates how non-conservative forces and damping influence the steady state and spectral properties of Brownian particles in optical traps, revealing strong signatures especially in underdamped regimes.

Contribution

It provides a comprehensive theoretical analysis of non-conservative forces on trapped particles, including probability distributions, phase space currents, and spectral densities across all dynamical regimes.

Findings

Non-conservative forces significantly alter steady state distributions.

Spectral density shows strong signatures of non-conservative forces in underdamped regimes.

Perturbative expressions are derived for all dynamical regimes.

Abstract

We consider a model of a particle trapped in a harmonic optical trap but with the addition of a non-conservative radiation induced force. This model is known to correctly describe experimentally observed trapped particle statistics for a wide range of physical parameters such as temperature and pressure. We theoretically analyse the effect of non-conservative force on the underlying steady state distribution as well as the power spectrum for the particle position. We compute perturbatively the probability distribution of the resulting non-equilibrium steady states for all dynamical regimes, underdamped through to overdamped and give expressions for the associated currents in phase space (position and velocity). We also give the spectral density of the trapped particle's position in all dynamical regimes and for any value of the non-conservative force. Signatures of the presence non-conservative forces are shown to be particularly strong for in the underdamped regime at low frequencies.

Figures

Equations

$$m\frac{d^2}{d{t}^2}\mathbf{x}(t)=-\gamma \frac{d}{dt}\mathbf{x}(t)+\mathbf{F}(\mathbf{x})+\sqrt{2T\gamma }\mathbf{\xi }(t),$$

$$⟨{\xi }_i(t){\xi }_j(t^{\prime })⟩={\delta }_{ij}\delta (t-t^{\prime }),$$

$$\gamma =6\pi {\eta }_{atm}{R}_p\frac{0.619}{0.619+Kn}(1+c_K)$$

$$\mathbf{F}=(\mathbf{p}\cdot \nabla )\mathbf{E}+\frac{d\mathbf{p}}{dt}\times \mathbf{B}.$$

$$\mathbf{p}(t)={\int }_{-\infty }^{\infty }d{t}^{\prime }\alpha (t-t^{\prime })\mathbf{E}(\mathbf{r},t).$$

$$\mathbf{E}(\mathbf{r},t)=E_0(\mathbf{r})\cos (\varphi (\mathbf{r})-\omega t){\mathbf{e}}_x,$$

$$\mathbf{B}=\frac{1}{\omega }\nabla \times [E_0(\mathbf{r})\sin (\varphi (\mathbf{r})-\omega t){\mathbf{e}}_x]$$

$$\mathbf{p}(t)$$

$$\mathrm{Re}\mathbf{E}(\mathbf{r},t)\alpha (\omega )$$

$$\mathbf{F}=\frac{1}{4}{\alpha }^{\prime }\nabla I(\mathbf{r})+\frac{1}{2}{\alpha }^{\prime \prime }I(\mathbf{r})\nabla \varphi (\mathbf{r}),$$

$$V(\mathbf{r})=\frac{{\alpha }^{\prime }}{4}I(\mathbf{r}).$$

$$\mathbf{E}(\mathbf{r})=E_0{\left[1+\frac{z^2}{z_0^2}\right]}^{-\frac{1}{2}}\exp \left[-\frac{x^2}{w_x^2(z)}-\frac{y^2}{w_y^2(z)}+i\varphi (\mathbf{r})\right]{\mathbf{e}}_{\mathbf{x}},$$

$$w_{\alpha }(z)=w_{\alpha }{\left[1+\frac{z^2}{z_0^2}\right]}^{\frac{1}{2}},$$

$$\varphi (\mathbf{r})=kz-\arctan \left(\frac{z}{z_0}\right)+\frac{k(x^2+y^2)}{2z(1+\frac{z_0^2}{z^2})},$$

$$\alpha =\frac{{\alpha }_0}{1-i\frac{k^3{\alpha }_0}{6\pi {\epsilon }_0}},$$

$${\alpha }_0=4\pi {ϵ}_0{R}_p^3\frac{ϵ-1}{ϵ+2},$$

$$V(\mathbf{r})={\kappa }_x\frac{x^2}{2}+{\kappa }_y\frac{y^2}{2}+{\kappa }_z\frac{z^2}{2}$$

$$\varphi (\mathbf{r})=z\left(k-\frac{1}{z_0}\right)+\frac{z}{z_0^2}\left[\frac{k}{2}(x^2+y^2)+\frac{z^2}{3z_0}\right].$$

$$E_x(\mathbf{r})\approx E_0\left\{1-\frac{x^2}{w_x^2}-\frac{y^2}{w_y^2}-\frac{z^2}{2}\left(\frac{2}{z_0^2}+k^2-\frac{2k}{z_0}\right)+i\left[z\left(k-\frac{1}{z_0}\right)+\frac{kz(x^2+y^2)}{2z_0^2}+\frac{z^3}{3z_0^3}-\frac{z^3}{6}{\left(k-\frac{1}{z_0}\right)}^3\right]\right\}$$

$$\frac{1}{z_0^2}-\frac{k}{z_0}+\frac{1}{w_x^2}+\frac{1}{w_y^2}=0,$$

$$z_0=\frac{2}{k\left[1-{\left(1-\frac{4}{k^2{w}_x^2}-\frac{4}{k^2{w}_y^2}\right)}^{\frac{1}{2}}\right]},$$

$$z_0=\frac{\lambda }{\pi }{\left[1-{\left(1-\frac{{\lambda }^2}{{\pi }^2{w}_x^2}-\frac{{\lambda }^2}{{\pi }^2{w}_y^2}\right)}^{\frac{1}{2}}\right]}^{-1}.$$

$${\mathbf{F}}_{\mathrm{scat}}=\frac{{\alpha }^{\prime \prime }}{2}{E}_0^2\left\{{\mathbf{e}}_x\frac{kxz}{z_0^2}+{\mathbf{e}}_y\frac{kyz}{z_0^2}+{\mathbf{e}}_z\left[(k-\frac{1}{z_0})\left(1-\frac{2x^2}{w_x^2}-\frac{2y^2}{w_y^2}-\frac{z^2}{z_0^2}\right)+k\frac{x^2+y^2}{2z_0^2}+\frac{z^2}{z_0^3}\right]\right\},$$

$${\mathbf{F}}_{\mathrm{scat}}=\frac{{\alpha }^{\prime \prime }}{{\alpha }^{\prime }}{\kappa }_z\left\{{\mathbf{e}}_{x}kxz+{\mathbf{e}}_{y}kyz+{\mathbf{e}}_z\left[z_0(z_{0}k-1)\left(1-\frac{2x^2}{w_x^2}-\frac{2y^2}{w_y^2}-\frac{z^2}{z_0^2}\right)+k\frac{x^2+y^2}{2}+\frac{z^2}{z_0}\right]\right\},$$

$${\mathbf{F}}_{\mathrm{scat}}=\frac{{\alpha }^{\prime \prime }}{{\alpha }^{\prime }}{\kappa }_z\left\{{\mathbf{e}}_{x}kxz+{\mathbf{e}}_{y}kyz+{\mathbf{e}}_z\left[z_0(z_{0}k-1)-x^2\left(\frac{2z_0(z_{0}k-1)}{w_x^2}-\frac{k}{2}\right)-y^2\left(\frac{2z_0(z_{0}k-1)}{w_y^2}-\frac{k}{2}\right)-z^2\left(k-\frac{2}{z_0}\right)\right]\right\}.$$

$${\nabla }^2{E}_x-[\nabla \varphi {]}^2{E}_x+k^2{E}_x+i\left(2\nabla E_x\cdot \nabla \varphi +E_x{\nabla }^2\varphi \right)=0,$$

$${\mathbf{F}}_{\mathrm{scat}}=\frac{{\alpha }^{\prime \prime }}{2}I\nabla \varphi =\frac{{\alpha }^{\prime \prime }}{2}{E}_x^2\nabla \varphi$$

$$\nabla \cdot {\mathbf{F}}_{\mathrm{scat}}=\frac{{\alpha }^{\prime \prime }}{2}{E}_x\left(2\nabla E_x\cdot \nabla \varphi +E{\nabla }^2\varphi \right)=0,$$

$${\mathbf{F}}_{\mathrm{scat}}={\mathbf{F}}_{\mathrm{scat},1}+{\mathbf{F}}_{\mathrm{scat},2},$$

$${\mathbf{F}}_{\mathrm{scat},1}=\frac{{\alpha }^{\prime \prime }}{{\alpha }^{\prime }}{\kappa }_z{z}_0(z_{0}k-1)\left[1-\frac{2x^2}{w_x^2}-\frac{2y^2}{w_y^2}\right]{\mathbf{e}}_z.$$

$${\mathbf{F}}_{\mathrm{scat},2}=\frac{{\alpha }^{\prime \prime }}{{\alpha }^{\prime }}k{\kappa }_z\left[{\mathbf{e}}_{x}xz+{\mathbf{e}}_{y}yz-{\mathbf{e}}_z{z}^2\right].$$

Full text

I Introduction

The seminal work of Ashkin in the 1970s instigated the development of optical traps ash . Optical traps or tweezers are an invaluable tool for studying statistical mechanics, both in an out of equilibrium, at the level of single particles. Moving the center of the optical trap allows one to measure very small forces in a wide range of systems, including colloids colloid , biological systems (for instance viruses, proteins and biopolymers) bio , as well as dielectric and metallic nano-particles nano . The same kind of perturbation can be used to analyze non-equilibrium processes, for instance optical tweezers have been used to demonstrate various non-equilibrium fluctuation theorems fluct . However in such studies the optical trap does not simply furnish a simple potential, optical scattering forces actually generate a non-conservative component to the overall force exerted by the trap. Because of the presence of these scattering forces, even static optical traps represent non-equilibrium systems as the particles trapped within them are driven by the non-conservative forces and thus are not described by a Gibbs-Boltzmann probability distribution. For instance, it has been demonstrated grier2008 ; grier2009 ; grier2010 ; grier2015 both theoretically and experimentally that scattering forces induce non-equilibrium currents, so called Brownian vortices, for trapped overdamped Brownian particles. At atmospheric pressures and near room temperatures, the underlying dynamics of a trapped particle can be taken to be overdamped Brownian. However optical traps operated at low pressure can probe dynamics where inertial effects become important and cross over to the underdamped regime. Clearly the current state of the art of the theory has to be modified to take into account the effects of inertia and this is the goal of this paper. Using the model of a Gaussian optical trap in the presence of a simplified non-conservative force we examine the impact of inertial effects across the whole regime between overdamped and underdamped Brownian motion.

First a general theory for the probability distribution of the perturbed steady state due to a small non-conservative force is developed for general damped Brownian dynamics. The modification of the steady state probability distribution function for the full phase space in the optical model of grier2015 is found and the resulting marginal probability distributions for the particle position and velocity is computed. We demonstrate that the Brownian vortices found in the current of the particle position persists and remains geometrically the same as that found in grier2015 . However the amplitude of the vortices depends on the damping parameters, we will show that this damping parameter can be tuned in order to optimize the non-equilibrium current. In addition, we show that there is a non-equilibrium current in the space of velocities, the vortex here turns out to be remarkably similar to the vortex in position space.

For the optical trap model of grier2015 we also compute the spectral density of the particle position for general damping. Our result recovers the computation of laporta2011 for overdamped Brownian motion and generalizes it to all damping regimes. Here we see that the non-conservative force generates very strong changes in the spectral density at low frequencies and that the amplitude of the low frequency spectral density can be tuned by varying the friction coefficient $\gamma$ and that there is a critical value of $\gamma$ where it is minimized. Finally, we also derive a function correlation which is a signal of the breaking of time reversal symmetry due to the non-conservative optical scattering force.

Some of the results derived in this paper are used in a accompanying letter, where the role of inertia in the presence of non-conservative forces is demonstrated in the experimental context letter . In this paper we concentrate on the presence of a nonconservative force in addition to a purely harmonic trap. This simplified model allows the derivation of a large number of analytical results. However, to fully explain the results of letter , where small anharmonicities can give rise to large effects in the underdamped regime due to near resonance phenomena, quartic terms in the potential must be considered to fully reproduce the experimental results. However the results derived here can be adapted to the experimental context by using effective harmonic parameters describing the nonlinear system Renz1 . This therefore yields an analytic theory sufficient to completely describe the experiments and fully non-linear numerical simulations letter .

II Modeling particles in optical traps

II.1 Langevin Dynamics

In this section we first discuss the underlying Langevin dynamics used to model the trapped particle. We then discuss the underlying model for trapping, aimed at experts in non-equilibrium statistical mechanics who are perhaps not familiar with the physics of optical trapping.

In this paper we will consider the dynamics of Brownian particles in an optical trap. If ${\bf x}(t)$ denotes the particle position, the underlying Langevin equation is

[TABLE]

where $m$ is the particle mass, $\gamma$ the friction coefficient and ${\bf F}({\bf x})$ is a time independent force generated by the optical trap. The noise term is assumed to have Gaussian white noise correlations in time, this noise has zero mean and correlation function

[TABLE]

in the amplitude term $T$ is the temperature measured in units where $k_{B}=1$. When the trap is applied in a gas phase the friction coefficient can be significantly modified and one can go from the underdamped to overdamped regimes with the same experimental set up if one can control the pressure in the trapping cell.

A semi-phenomenological relationship the friction coefficient in air as a function of the pressure is given by beresnev1990 ; Li2013

[TABLE]

where in the first, Stokes drag-like, term $R_{p}$ is the radius of the particle and $\eta_{\rm atm}$ the viscosity coefficient of air at atmospheric pressure. In addition, $Kn$ is the Knudsen number, accounting from the deviation from continuum behavior, $Kn=\lambda_{\rm mfp}/R_{p}$ where $\lambda_{\rm mfp}$ is the mean free path of air molecules and is given by $\lambda_{\rm mfp}=(68{\rm nm})P_{\rm atm}/P$, where $P$ is the applied pressure in the trap and $P_{\rm atm}$ is the atmospheric pressure. Finally $c_{K}=(0.31Kn)/(0.785+1.152Kn+Kn^{2})$. We thus see that the value of the friction coefficient can be tuned by modifying the pressure of the trap.

II.2 Model of the optical trap

In the experiment described in the accompanying paper letter , the wave length of the laser is $\lambda=1064{\rm nm}$, while the radius of the trapped (fused silica) particle is $R_{p}=68{\rm nm}$. As $R_{p}\ll\lambda$, we are in the Rayleigh regime, where the particle can be treated as a point dipole.

Here for completeness we review the model of the optical trap used in quid2013 as well as reviewing the underlying physics giving rise to both conservative and non-conservative forces in the dipole approximation hecht2006 . The Lorentz force acting on an electric dipole ${\bf p}$ is given by

[TABLE]

The response of the dipole to the applied electric field is taken to be linear

[TABLE]

For a mono-chromatic electromagnetic field of angular frequency $\omega$ polarized in the $x$ direction we can write

[TABLE]

where ${\bf e}_{x}$ denotes the unit vector in the direction ${\bf x}$ and $\phi({\bf r})$ is the phase. Using the Maxwell equation $\partial{\bf B}/\partial t=-\nabla\times{\bf E}$ we can write

[TABLE]

Using linear response we can write

[TABLE]

so $\alpha(\omega)$ is the polarizability at frequency $\omega$ and $\alpha^{\prime}$ and $\alpha^{\prime\prime}$ denote the real (dispersive) and imaginary (dissipative) parts respectively. Now evaluating the force and replacing the terms $\cos^{2}(\phi({\bf r})-\omega t)$ and $\sin^{2}(\phi({\bf r})-\omega t)$ by their temporal average $1/2$ and cross terms $\cos(\phi({\bf r})-\omega t)$ and $\sin(\phi({\bf r})-\omega t)$ by their temporal average [math] we obtain

[TABLE]

where $I({\bf r})=|{\bf E}({\bf r})|^{2}$ is the field intensity. The first term above is clearly a conservative force generated by the potential

[TABLE]

The second term contains in general a non-conservative component. Its underlying origin comes from the interaction of the magnetic field with the rate of change of the dipole moment.

The electric field produced by the laser is approximated by the Gaussian form

[TABLE]

where the field is polarized in the direction ${\bf e}_{x}$ and for $\alpha=x,\ y$ we have

[TABLE]

is the $z$ dependent beam radius. The phase is given by

[TABLE]

where $k=2\pi/\lambda$ is the wave vector. The parameter $z_{0}$ is the Rayleigh range and is normally used as a fitting parameter, however it is predicted to be given by $z_{0}\sim\pi w_{\alpha}^{2}/\lambda$. Finally, the particle is assumed to be spherical and the polarizability $\alpha$ is taken to be hecht2006

[TABLE]

where $\alpha_{0}$ is the zero-frequency polarizability given by

[TABLE]

and $\epsilon_{0}$ is the dielectric constant of the vaccum while $\epsilon$ is the relative dielectric constant of the particle. The form of the electric field given by Eq. (11) arises from the paraxial approximation to Maxwell’s equation for the electric field in the vacuum, and in this case the corresponding Helmholtz equation. The paraxial approximation assumes that for the whole electric field function denoted by $E$ (polarized in the direction $x$) one has that $|\partial^{2}E/\partial z^{2}|\ll|k\partial E/\partial z|$. Inspection of Eqs. (11) and (13) show that the paraxial approximation is valid when $z/z^{2}_{0}\ll k$, and one can expand the force perturbatively for small $z$ when $z/z_{0}\ll 1$. The first inequality implies that we must have $z\ll 2\pi^{3}\omega^{4}_{\alpha}/\lambda^{3}$.

In the accompanying letter letter the laser waist parameters are estimated in the high pressure regime where both the effects of trap nonlinearities and the scattering forces are weak. It is found that $w_{x}=0.915\mu m$, $w_{y}=1.034\mu m$ and $z_{0}=w_{z}=2.966\mu m$ (in the notation of letter ) and so we must have $z\ll 0.29\ {\rm\mu m}$. In addition, one may also Taylor expand the forces when $x/w_{x}$ and $y/w_{z}$ are small, experimental results show that these two values are both of the order of $0.1$, so while non-linear effects are weak they are not completely negligible.

Under the above assumptions, to lowest order, the potential part of the force is generated by a quadratic potential

[TABLE]

with spring constants given by $\kappa_{x}=E_{0}^{2}\alpha^{\prime}/w_{x}^{2}$, $\kappa_{y}=E_{0}^{2}\alpha^{\prime}/w_{y}^{2}$ and $\kappa_{z}=E_{0}^{2}\alpha^{\prime}/2z^{2}_{0}$.

The phase factor $\phi$ is given, expanding to first and next order of the particle displacement, by

[TABLE]

As mentioned above, the term $z_{0}$ is often treated as a fitting parameter as the solution Eq. (11) is an approximate solution. However an alternative physical approach is to insist that the approximated solution obeys Maxwell’s equations exactly near the center of the trap. Expanding the solution for small $x$, $y$ and $z$ up to cubic order we find that

[TABLE]

Now using the Helmholtz equation $\nabla^{2}E_{x}+k^{2}E_{x}=0$ near the origin we find the ${\cal O}(0)$ equation, coming from the real term,

[TABLE]

giving an exact relation between $z_{0}$ and the two waists $w_{x}$ and $w_{y}$ in terms of the wave vector $k$. In the case where $w_{x}=w_{y}=w$ and when $k\gg 1/z_{0}$ (so we assume that the paraxial approximation is valid) we obtain the standard formula $z_{0}=w^{2}k/2=\pi w^{2}/\lambda$. However this formula is only valid if $2\pi^{2}w^{2}/\lambda^{2}\gg 1$. Taking the root of Eq. (19) compatible with the standard formula we find

[TABLE]

or in terms of the wavelength $\lambda$ we have

[TABLE]

Using this to determine $z_{0}$, with the values of $w_{x}$ and $w_{y}$ determined above yields $z_{0}=2.781\mu m$, which is different from the fully fitted value of $z_{0}$ by about $6\%$. This estimate thus provides a useful starting point to determine the trap parameters in the fitting procedure.

The scattering part of the force is given to leading, quadratic, order by

[TABLE]

which can be written in terms of the spring constant $\kappa_{z}$ as

[TABLE]

which gives

[TABLE]

The lowest order approximation that takes into account both conservative/potential and scattering forces is thus the expressions for ${\bf F}_{p}=-\nabla V$ (${\cal O}(1)$) and ${\bf F}_{\rm scat}$ (${\cal O}(0)+{\cal O}(2)$). The next order corrections are ${\cal O}(3)$ terms to ${\bf F}_{p}$ coming from anharmonic contributions.

Using the form Eq. (11) in the Helmholtz equation without making any approximations yields

[TABLE]

The real and imaginary parts of the above give two equations. The scattering force is given by

[TABLE]

and taking the divergence of this yields

[TABLE]

from Eq. (25). However the result Eq. (24) gives $\nabla\cdot{\bf F}_{\rm scat}=4z\alpha^{\prime\prime}\kappa_{z}/\alpha^{\prime}z_{0}$. In the limit $z_{0}\to\infty$, where the paraxial approximation becomes exact we find that this term is zero, but in general it is not. This actually mean that the very last term in Eq. (24) should not be there. This can be regarded as being due to an error in the phase term.

Taking this into account we can naturally divide the scattering force into two components

[TABLE]

where

[TABLE]

and

[TABLE]

For the values of $z_{0}$ obtained here the dominant scattering force is ${\bf F}_{\rm scat,1}$, however as $z$ scales as $z_{0}$ when just the harmonic confinement is taken into account it is not evident that the last term in the direction ${\bf e}_{z}$ can be neglected. However we find that numerically this turns out to be the case.

III Perturbative calculation of the non-equilibrium steady state probability density function and current

In the standard paradigm of statistical physics, Brownian particles are subjected to forces generated by potentials and this means that the underlying steady state distribution is given by the Gibbs-Boltzmann distribution. The crucial point here is that the steady state in the presence of a nonconservative force has an associated current. We note that Brownian particles subject to linear potentials in unbounded systems also have currents, however because the system is unbounded the probability distribution spreads with time and can be associated with effective transport coefficients dea2007 . However in periodic systems the particle position modulo the period of the system does have a steady state which carries a current, and this system can be analysed for overdamped Brownian particles reim2001 . In general very few results exist on the steady states of systems with nonconservative forces or in the presence of driving, even for trapped mono-particle systems where only one body forces are present. Here we will resort to the use of perturbation theory, which can be justified when the nonconservative component of the force is weak in comparison to the conservative component. Indeed, perturbation theory was the approach employed in grier2015 to analyse a trapped overdamped Brownian particle.

We will consider systems which are driven by underlying white noise, and thus include both underdamped and overdamped Brownian motion. The probability distribution of the particle’s position in phase space is thus described by the Fokker-Planck equation

[TABLE]

where ${\bf y}$ denote the phase-space coordinates and $H$ is the forward Fokker-Planck operator. Before perturbation by a non-conservative force we assume that the system obeys detailed balance and has an underlying Gibbs-Boltzmann distribution denoted by $P_{0}$, obeying

[TABLE]

Now we consider what happens when the system is perturbed by an additional force which is in terms of the Fokker-Planck operator denoted by the perturbation $\Delta H$. The question we address is what is the new steady state distribution $P_{s}$ defined by

[TABLE]

To first order in perturbation theory we write $P_{s}=P_{0}+\Delta P_{0}$ where $\Delta P_{0}$ is given by

[TABLE]

This now has the formal solution

[TABLE]

Here $H^{-1}$ is the pseudo-Green’s function for the operator $H$ as the solution $\Delta HP_{0}$ must obey the normalization

[TABLE]

The transition probability $P({\bf y}|{\bf y}^{\prime};t)$ can formally be written as

[TABLE]

and the pseudo-Green’s function thus as

[TABLE]

This yields

[TABLE]

Now using the time reversal symmetry of the equilibrium state

[TABLE]

where ∗ indicates time reversal of the coordinates, this means that coordinates with an odd number of temporal derivatives change sign. Taking $t=0$ we see that time reversal symmetry also implies that $P_{0}({\bf y}^{*})=P_{0}({\bf y})$. For overdamped Brownian particles this has no effect as the position vector obeys ${\bf X}={\bf X}^{*}$. However if ${\bf Y}$ represents the position and velocity ${\bf Y}=({\bf X},{\bf V})$ we have ${\bf Y}^{*}=({\bf X},-{\bf V})$. We thus find

[TABLE]

If the stochastic process represented by the Fokker-Planck equation is denoted by ${\bf Y}_{t}$ and we use the notation

[TABLE]

i.e. the average value of $f({\bf Y}_{t})$ at time $t$ given that ${\bf Y}_{0}={\bf y}$, and

[TABLE]

the equilibrium value of $f({\bf Y})$. We thus change the integration variable to write

[TABLE]

where we have used ${\bf y}^{**}={\bf y}$. This then yields

[TABLE]

The change in the equilibrium distribution at the point ${\bf y}$ in phase space is thus given by a Kubo-type formula kubo . This result is general and applies to any system described by a Fokker-Planck equation.

III.1 Overdamped Brownian motion

Here we consider a process characterized solely by the position ${\bf X}$, that is to say the overdamped regimes. Here we assume that the unperturbed Fokker Planck operator, defined via its action on a test function $f$, is given by

[TABLE]

i.e. it describes a particle with diffusivity $D$ in a potential $V$ at temperature $T$ (where for convenience we use units where $k_{B}=1$, the results for $T$ measured in Kelvin are thus obtained by replacing $T$ by $k_{B}T$ in all the formulas which follow) and $\beta=1/T$.

Here the steady state is given by the Gibbs-Boltzmann equilibrium distribution

[TABLE]

The perturbation of the Fokker-Planck operator due to the presence of an extra force ${\bf F}$ is given by

[TABLE]

From this we find

[TABLE]

This now gives

[TABLE]

The case where the force is derived from a potential, ${\bf F}=-\nabla U$, should of course be simple as we know the new equilibrium measure should be

[TABLE]

and for small $U$ should find

[TABLE]

It is illuminating, as well as a useful check, to obtain this from the general formalism derived previously. To do this, we use the following result of stochastic calculus osk for an arbitrary function $U$:

[TABLE]

Using this and the relation ${\bf F}=-\nabla U$ in Eq. (50) now gives

[TABLE]

Clearly we have

[TABLE]

and also we have by definition that

[TABLE]

The time integral in Eq. (54) can now be simply carried out to yield

[TABLE]

in agreement with Eq. (52).

In general the perturbing force can be decomposed into a conservative and nonconservative part via the Helmholtz decomposition where one writes

[TABLE]

where $U$ is determined from the equation

[TABLE]

and the remaining term can be similarly evaluated to write

[TABLE]

where the nonconservative part of the force ${\bf F}_{n}$ is given by

[TABLE]

however, the solution for $U$ must be chosen so that $\nabla\cdot{\bf F}_{n}=0$. Using this decomposition then yields

[TABLE]

The modification of the steady state probability density function by the non-conservative force is thus seen to be related to the work done by the non-conservative force in the direction of the original potential $V$. This formula is not particularly illuminating from a physical point of view for overdamped Brownian motion, when we study the case of general damping we will find a much more explicit interpretation.

III.1.1 Application to a model for an optical trap

Here we study the model proposed by Grier et al grier2015 , where

[TABLE]

so this represents an anisotropic harmonic trap unless $\eta=1$. We have seen that the potential part of the scattering force in Eq. (16) in general also has $\kappa_{x}\neq\kappa_{y}$ as the difference between $\kappa_{x}$ and $\kappa_{y}$ is numerically small we carry out the analysis for this case to simplify the notation and the following analysis. In grier2015 the nonconservative component of the force generated by the optical trap is taken to be

[TABLE]

$a$ is proportional to the waist $w_{x}\approx w_{y}$ (we use the same notation as grier2015 for ease of comparison). The part of the scattering force used is thus ${\bf F}_{\rm scat,1}$, expressed in Eq. (29) which is that which is numerically the most important. However the nature of our perturbative computation means that the corrections to $\Delta P_{0}$ from each component of the non-conservative force is additive, and so we can simply compute the perturbations due to each component and then sum them. The perturbing parameter $\varepsilon$ is dimensionless and we will compute the modified steady state probability density function to first order in $\varepsilon$.

The decomposition of the perturbation in terms of a conservative and nonconservative contribution now yields the potential

[TABLE]

and the nonconservative component of the force

[TABLE]

Using this we find

[TABLE]

The evaluations of the expectation values required to compute the renormalization of the steady state probability distribution can be carried out by explicitly solving the Langevin equation for the unperturbed process, the corresponding equations are

[TABLE]

where $\xi_{x}(t)$, $\xi_{y}(t)$ and $\xi_{z}(t)$ are independent Gaussian white noises. The above equations can be solved in terms of the particle positions at $t=0$, denoted by $x$, $y$ and $z$ to give

[TABLE]

Using these solutions we find

[TABLE]

and we also see that

[TABLE]

Putting all this together (and noting that $\langle U({\bf X})\rangle_{0}=0$) we obtain

[TABLE]

where $\Delta\phi$ is interpreted as the renormalization of the trapping potential. This agrees with the result found in grier2015 for the same trapping potential and nonconservative perturbation considered here.

III.1.2 Alternative formulation for overdamped particles

The probabilistically based solution for the perturbed steady state distribution gives us an interpretation of the perturbation of the steady state in terms of work done by the perturbing forces. However in the case where the unperturbed solution is Gaussian we see that the solution is particularly simple, the solution having a polynomial form for a polynomial perturbing force. Given the structure of the solution found we make the same ansatz as in grier2015

[TABLE]

this yields

[TABLE]

For the case studied above as defined by the potential in Eq. (63) and the force in Eq. (64) we find

[TABLE]

Inspection of the probabilistic representation of the solution, and from the symmetry of the problem in the $(x,y)$ plane, now tells us that the solution must have the form

[TABLE]

this yields

[TABLE]

and comparing the coefficients of the polynomials now gives

[TABLE]

in agreement with the earlier calculation leading to Eq. (76) and thus the result given in grier2015 .

Here we also consider the perturbation to a force of type ${\bf F}_{\rm scat,2}$ derived in Eq. (30) which we write as

[TABLE]

as it has no potential part. Here we find

[TABLE]

where the ′ on $K$ denotes that it is the form of $K$ for the perturbation ${\bf F}^{\prime}$. Here we see that the solution must be of the form

[TABLE]

This gives

[TABLE]

and thus

[TABLE]

This therefore gives a contribution to the steady state distribution

[TABLE]

III.2 Underdamped Brownian particles

Here the phase space corresponds to ${\bf Y}=({\bf X},{\bf V})$ where ${\bf X}$ and ${\bf V}$ denote respectively the particle position and its velocity. Here the Fokker-Planck operator is given by

[TABLE]

where $\gamma$ is the friction coefficient, $m$ is the particle mass, and ${\bf F}_{t}$ is the total force on the particle which is assumed to depend only on the position ${\bf x}$. For the unperturbed system we take ${\bf F}_{t}({\bf x})=-\nabla V({\bf x})$. The equilibrium distribution is then given by

[TABLE]

In the perturbed system we denote the additional force simply by ${\bf F}({\bf x})$. This yields

[TABLE]

and thus

[TABLE]

where we have used ${\bf V}^{*}=-{\bf V}$. The result in Eq. (95) immediately gives us a physical interpretation of the modification of the steady state probability density function. We see that

[TABLE]

where $\overline{W}({\bf x},-{\bf v})$ is the total work done, after a very long time $t^{*}$, by the non-conservative force on a particle started at the posiiton ${\bf x}$ but and with initial velocity $-{\bf v}$ , whereas $\overline{W}_{0}$ is the total work done by the non-conservative force, at the same late time $t^{*}$, for a particle whose initial velocity and position are taken from the unperturbed equilibrium probability density function. Both of these total works diverge as $t^{*}\to\infty$ but their difference is constant. When the perturbing force is purely conservative, i.e. ${\bf F}=-\nabla U$ we immediately obtain, much more transparently that in the Brownian case,

[TABLE]

and as ${\bf V}=d{\bf X}/dt$ we find, using the same arguments as for the overdamped case

[TABLE]

thus recovering the corresponding modified Gibbs-Boltzmann distribution.

III.3 Application to under generally damped particles in optical trap model

In principle, for a harmonic trap all the average values in Eq. (95) can be computed. However, as for the overdamped case, the form of the solution suggests writing

[TABLE]

This then yields

[TABLE]

There is a distinct advantage, at least for complex perturbations by non-conservative forces, of using this algebraic formulation as opposed to the probabilistic representation. However the probabilistic representation allows one to predict the polynomial form of $K$ and only the coefficients of this polynomial need be determined after substituting the correct form of $K$ into Eq. (100). This is because in the expression $\langle{\bf V}_{t}\cdot{\bf F}({\bf X}_{t})\rangle_{{\bf x},-{\bf v}}$ we can see explicitly what terms in the initial position ${\bf x}$ and velocity ${\bf v}$ must arise by inspecting the form of ${\bf F}$ given in Eq. (64).

From the symmetry in the $(x,y)$ plane, the solution of Eq. (100) must have the form

[TABLE]

After simplification, the coefficients are given by

[TABLE]

where $S=a[\eta\kappa m+2\gamma^{2}][(\eta-4)^{2}\kappa m+6(2+\eta)\gamma^{2}]$. Using these expressions we have thus computed the perturbed probability density function in the full $({\bf x},{\bf v})$ phase space. However the final expressions are rather unwieldy, we will therefore give explicit expressions for the marginal probability density functions of ${\bf X}$ and then ${\bf V}$.

Denoting the marginal steady state probability density function for the position by $\Delta P_{0}({\bf x})=\int d{\bf v}K({\bf x},{\bf v})P_{0}({\bf x},{\bf v})$, we find from Eq. (93)-(101)

[TABLE]

where $\rho=\sqrt{x^{2}+y^{2}}$. In the underdamped limit $\gamma\rightarrow 0$, we find thus

[TABLE]

and in the overdamped limit $m\rightarrow 0$, we recover (for a third time) the result of grier2015 Eq. (76).

Using the modified probability density function in Eq. (110) can study various moments of the particle position. The average value of the height of the particle is given by

[TABLE]

while the average velocity in the direction $z$ is given by $\langle V_{z}\rangle=0$, after some mathematical cancellations which appear remarkable from a mathematical point of view but obvious from the simple physical insight that the particle is trapped! The above result for $\langle z\rangle$ is thus independent of the details of the dynamics (it is independent of the damping parameters). This surprising result is a consequence of the fact that the average velocity is zero (as obviously must be the case if the particle remains trapped): as the average velocity is zero the mean position along the $z$ axis is given by a force balance equation where the friction term is zero and, as the average thermal force and acceleration are zero, the average force balance equation depend only on the external forces acting on the particle and is thus independent of the damping parameters.

Here we introduce the underlying frequency of the harmonic part of the trap $\omega_{0}^{2}=\kappa/m$ and the damping rate $\Gamma=\gamma/m$. The quality factor is defined as $Q=\omega_{0}/\Gamma=\sqrt{\kappa m/\gamma^{2}}$, it quantifies the degree of damping in the system, large $Q$ corresponding to the underdamped regime and small $Q$ corresponding to the overdamped regime. In terms of these variables, which are often preferred in experimental studies, we find

[TABLE]

Denoting the marignal steady state probability density function for the velocity by $\Delta P_{\rm MB}({\bf v})=\int d{\bf x}K({\bf x},{\bf v})P_{0}({\bf x},{\bf v})$, we find from Eq. (93)-(101)

[TABLE]

where $v_{\rho}=\sqrt{v_{x}^{2}+v_{y}^{2}}$ and

[TABLE]

is the Maxwell-Boltzmann velocity distribution. The marginal probability distributions for $v_{x}$, $v_{y}$ and $v_{z}$ however do not change in first order perturbation theory and are thus given by the Maxwell distribution.

In grier2015 the steady state current was computed for an overdamped Brownian particle. Here we have currents in both position and velocity spaces. From Eq. (92), we see that currents ${\bf J_{x}}$ and ${\bf J_{v}}$ can be read off from

[TABLE]

where

[TABLE]

The effective currents in position space $\overline{\bf J_{x}}({\bf x})$ and velocity space $\overline{\bf J_{v}}({\bf v})$, can be found by integrating Eq. (116) over ${\bf v}$ and ${\bf x}$ respectively to obtain the Fokker-Plank equation for the marginal probability distribution in position, ${\bf x}$ and velocity ${\bf v}$ respectively. From this, and using the divergence theorem, the associated currents can be read off as

[TABLE]

In the non-equilibrium steady state, we find that the currents are given, to first order in the perturbation by the nonconservative force as,

[TABLE]

These expressions can then be rewritten in terms of the function $K({\bf x},{\bf v})$ to give

[TABLE]

Using the form of $K({\bf x},{\bf v})$ derived for the optical trap model we obtain

[TABLE]

where ${\bf e_{\rho}}$ is the radial unit vector in polar coordinates (recall here $\rho=\sqrt{x^{2}+y^{2}}$). In the limit $m\rightarrow 0$, we recover the result of grier2015 for the overdamped case

[TABLE]

where $D$ is the effective diffusion constant, as defined by the Einstein relation, $D=T/\gamma$. In the underdamped limit $\gamma\rightarrow 0$, we obtain

[TABLE]

Remarkably the geometry of the current field $\overline{\bf J_{x}}({\bf x})$ in space is always the same, the damping parameters only appear via an overall space independent amplitude term $A_{x}(\eta,Q)$ such that

[TABLE]

where $\rho_{c}=\sqrt{2/\beta\kappa}$ and

[TABLE]

We see above that when $Q\to 0$ and $Q\to\infty$, $A_{x}\to 0$ and thus the amplitude of the non-equilibrium current can be maximized at a given value of $Q$ when all other physical parameters are fixed, as can be clearly seen in Fig. 1a. We also see that the amplitude increases as $\eta$ decreases, this is to be expected as the trapping in the $z$ direction becomes weaker on reducing $\eta$.

As pointed out in grier2015 , at the point $\rho_{c}$ the direction of the current in the $z$ direction reverses, the value of $\rho_{c}$ is the value where the particles potential energy due to the harmonic trap attains the equilibrium value as predicted by the equipartition of energy. We also see that the current in the $\rho$ direction is positive for $z>0$ and negative for $z<0$. The current in the radial (polar coordinate) direction $\rho$ is in the direction ${\bf e_{\rho}}$ and its sign is proportional to $z$, thus the current flows away from the origin in the $(x,y)$ plane when $z>0$ and towards the origin when $z<0$.

The current in velocity space is given by $\bf v$

[TABLE]

where $v_{\rho}=\sqrt{v_{x}^{2}+v_{y}^{2}}$. In the overdamped limit, $m\rightarrow 0$, one finds

[TABLE]

and the limit $\gamma\rightarrow 0$, we find

[TABLE]

Again we see that the geometric structure in phase space is independent of the damping parameters and that the damping parameters only appear in a velocity independent constant amplitude. Remarkably the geometric structure is basically identical to that of the position current if we reverse the direction and make the correspondence ${\bf v}\equiv{\bf x}$. Here the $Q$-dependent amplitude of the current $A_{v}(\eta,Q)$ is defined such that

[TABLE]

where $v_{c}=\sqrt{2/\beta m}$, and

[TABLE]

Here we find that as $Q\to 0$ $A_{v}(\eta,Q)\to 0$, whereas $A_{v}(\eta,Q)=1/(4-\eta)$ as $Q\to\infty$. This function is monotonically increasing with $Q$, as can be seen in Fig. 2a, it reaches a plateau at large $Q$. In contrast to the amplitude of the spatial current $A_{x}(\eta,Q)$ we see that the amplitude is largest for large values of $\eta$, i.e. for strong trapping in the $z$ direction.

Note that the apparent divergence for $\eta=4$ is removed due to the presence of terms of next order in $\gamma$, however the point $\eta=4$ should present some form of resonance in the current. The value $\eta=4$ is however far from its typical experimental value which is around $0.3$ letter .

III.4 Application to general particles in optical trap model with the ${\bf F}_{scat,2}$ force

We want to find $\Delta P_{s}({\bf x},{\bf v})=P_{0}({\bf x},{\bf v})K^{\prime}({\bf x},{\bf v})$ satisfying the equation (100) with ${\bf F^{\prime}}$ defined in Eq. (84) instead of ${\bf F}$. The solution of this equation must have the form

[TABLE]

After simplification, the coefficients are given by

[TABLE]

where $S=a[\eta\kappa m+2\gamma^{2}][(\eta-4)^{2}\kappa m+6(2+\eta)\gamma^{2}]$. Denoting the marginal steady state probability density function for the position by $\Delta P_{s}({\bf x})=\int d{\bf v}K^{\prime}({\bf x},{\bf v})P_{0}({\bf x},{\bf v})$, we find from Eq. (93) and (135)

[TABLE]

We recover the Eq. (91) taking the limit $m\rightarrow 0$. Denoting the marginal steady state probability density function for the velocity by $\Delta P_{\rm MB}^{\prime}({\bf v})=\int d{\bf x}K^{\prime}({\bf x},{\bf v})P_{0}({\bf x},{\bf v})$, we find from Eq. (93) and (135)

[TABLE]

The effective current in position space is defined as $\overline{\bf J_{x}^{\prime}}({\bf x})=\int d{\bf v}\ {\bf v}P_{0}({\bf x},{\bf v})K^{\prime}({\bf x},{\bf v})$ from Eq. (123). Using the expression of $K^{\prime}$ from Eq. (135), we get

[TABLE]

The effective current in velocity space is defined from Eq. (124) as

[TABLE]

After some simplifications, we find

[TABLE]

We see that the above results currents induced by the force ${\bf F}^{\prime}$ are remarkably similar to those induced by the non- conservative part of the force ${\bf F}_{n}$ of the force ${\bf F}^{\prime}$. This can in fact be seen much more easily by noting that the nonconservative parts of ${\bf F}$ and ${\bf F}^{\prime}$ are proportional. Indeed, we can write

[TABLE]

where $U=-\kappa(\frac{x^{2}+y^{2}}{2}-\frac{z^{3}}{3})$ is a potential difference between the two nonconservative forces. This potential term does not contribute to the current at the order if perturbation theory used here, thus explaining the simple relationship between the two currents found here.

III.5 Comparison with numerical simulation

We will now verify the analytical results for the nonequilibrium steady state probability distribution function and the associated currents, presented above, by comparing them with the results of numerical integration of the corresponding Langevin equations which are given by

[TABLE]

Here, $\xi_{x}$, $\xi_{y}$ and $\xi_{z}$ are independent Gaussian white noises. We consider the system with the following adimensionalized parameters, where $a=1$ and $\kappa_{x}/m=\omega_{0}^{2}=1$ such that

[TABLE]

where $Q=\omega_{0}/\Gamma$ is the quality factor, $\sigma^{2}=(\beta\kappa_{x}a^{2})^{-1}$ and $\eta_{y}=\kappa_{y}/\kappa_{x}$ and $\eta_{z}=\kappa_{z}/\kappa_{x}$. The numerical integration of these equations is carried out using the algorithm of Sivak et al. sivak2014 . The probability density function is estimated using a spatial binning procedure, for boxes of volume $v=\Delta x\Delta y\Delta z$, as

[TABLE]

where $I({\bf X}_{t_{i}},{\bf x})=1$ if ${\bf X}_{t_{i}}$, the particle position at time $t_{i}=i\Delta t$, is in the bin of ${\bf x}$ and ${\cal N}$ is the total number of points in the trajectory. The spatial current is computed from Eq. (119) using the estimator

[TABLE]

where ${\bf V}_{t_{i}}$ is the measured particle velocity at time $t_{i}$.

The numerically obtained value of the steady state probability distribution for the parameters $\varepsilon=0.1$, $\sigma=0.5$, $\eta_{y}=1$, $\eta_{z}=0.2$ and $Q=1$ while the bins are taken with respect to $\Delta x=\Delta y=6\sigma/100$ and $\Delta z=6\sigma_{z}/100$ (where $\sigma_{z}=\sigma/\sqrt{\eta}$). The simulation is carried out using for the integration time step $\Delta t=10^{-2}$ and over a total measurement time $t_{m}=10$ for a total of ${\cal N}=5.10^{9}$ different trajectories. We show in Fig 3(a) the estimated steady state distribution as a grey-scale (color online) plot, along with the vector field plot of the current (shown as arrows). In Fig. 3b we see the corresponding quantities predicted by the perturbation theory employed in this paper. The agreement is excellent, as it should be for the small value of $\varepsilon$ used here.

IV Power spectra for underdamped particles in perturbed Harmonic potentials

The computation of the spectral density function for the particle position in the trap model considered here has been carried out in laporta2011 . Here we generalize this calculation to the case of underdamped Brownian particles and show that significant differences arise. We then verify the calculation by comparing it with numerical simulations.

IV.1 Analytical calculation

Here we compute the power function of an underlying harmonic process subjected to a perturbing non-conservative force. The equations of motion for the components of the particle position are again given by Eqs. (152, 153, 154), and in this particular approximation only the steady state distribution of $z$ is changed (the positions $x$ and $y$ are not coupled to the non-conservative force and the process $z$ is a slave to the positions $(x,y)$ slave ). The addition of the perturbation gives $z$ a non-zero average value given by

[TABLE]

Due to the linearity of Eq. (154) we can decompose $z$ as

[TABLE]

where $z_{ef}$ is the equilibrium fluctuation due to the thermal noise obeying

[TABLE]

and $\zeta_{n}$ is the fluctuation in $z$ due to the fluctuations of the non-conservative force, $\zeta_{n}$ then obeys

[TABLE]

where

[TABLE]

is the fluctuating part of the non-conservative noise in the $z$ direction. Fourier transforming in time then gives

[TABLE]

and

[TABLE]

The power spectrum of the fluctuations $\zeta=\zeta_{ef}+\zeta_{n}$ is then given by

[TABLE]

where

[TABLE]

is just given by the ordinary equilibrium fluctuation spectrum

[TABLE]

The contribution from the non-conservative force to the power spectrum is

[TABLE]

Defining

[TABLE]

then gives

[TABLE]

To proceed with the computation of $C(\omega)$ we use

[TABLE]

and so the associated connected correlation function is given by

[TABLE]

In terms of the equilibrium spectral function for $x$

[TABLE]

we thus find

[TABLE]

Using this we obtain

[TABLE]

the last line holding when the system is isotropic in the $(x,y)$ plane.

In order to compute $C(\omega)$ we note that the poles of $S_{xx-eq}(\nu)$ are given by

[TABLE]

In the underdamped case where $4\kappa m>\gamma^{2}$ the poles can be written as

[TABLE]

where

[TABLE]

In the evaluation of the integral in Eq. (177) the poles that contribute from $S_{xx-eq}(\nu)$ come from the upper complex plane and thus only $\nu=\nu_{\pm}$ contribute. The poles of $S_{xx-eq}(\omega-\nu)$ can be written to be at $\nu-\omega=\nu_{p}$ and thus $\nu=\nu_{p}+\omega$. Therefore, only the poles $\nu=\nu_{\pm}+\omega$ contribute from $S(\omega-\nu)$. In terms of these poles we can write

[TABLE]

Evaluating the residues in the upper complex plane we thus find

[TABLE]

We now proceed by writing

[TABLE]

where $\Gamma=\gamma/m$ and $\Omega=\sqrt{4\kappa m-\gamma^{2}}/m$. Using this yields

[TABLE]

This finally gives

[TABLE]

Substituting this into Eq. (172) now gives

[TABLE]

which can be written as

[TABLE]

In terms of, nearly, all the physical parameters we have

[TABLE]

where

[TABLE]

This can also be written as

[TABLE]

In the overdamped limit, where we set $m=0$, one finds

[TABLE]

which agrees with the calculation for an overdamped Brownian motion in laporta2011 .

In the case where the trap is also anisotropic in the $(x,y)$ plane we write the stiffnesses in the directions $x$, $y$ and $z$ as $\kappa_{x}$, $\kappa_{y}$ and $\kappa_{z}$, and use the corresponding values of $w_{x}$ and $w_{y}$ for the fluctuating nonconservative force which we write from Eq. (29) as

[TABLE]

where

[TABLE]

with $a=\sqrt{a_{x}a_{y}}$ where $a_{x}=w_{x}/\sqrt{2}$ and $a_{y}=w_{y}/\sqrt{2}$. Here $\kappa_{g}=\sqrt{\kappa_{x}\kappa_{y}}$ is the geometric mean of the stiffness in the directions $x$ and $y$. We can thus write

[TABLE]

where $\varepsilon_{x}=\varepsilon a/a_{x}$ and $\varepsilon_{y}=\varepsilon a/a_{y}$. We can now simply use the results obtained above (as the contributions from ${\bf F}_{x}$ and ${\bf F}_{y}$ are additive) to obtain

[TABLE]

If we denote the geometric mean of the two waists in the $x$ and $y$ directions by $w_{g}=\sqrt{w_{x}w_{y}}$, we find that in terms of the waste variables we have

[TABLE]

If we consider the low frequency limit $\omega\to 0$ and consider the simplest case where the system is isotropic in the $(x,y)$ plane, i.e. $w_{x}=w_{y}$, we find

[TABLE]

which has the form

[TABLE]

where $A$ and $B$ are damping independent constants given by

[TABLE]

This predicts that $S_{z}(0)$ attains a minimal value at $\gamma_{c}=\sqrt{B/A}$.

IV.2 Signatures of time reversal symmetry breaking

The crucial difference between equilibrium and non-equilibrium systems stems from the breaking of time reversal symmetry in the latter systems. It is this symmetry breaking that leads to the appearance of currents in non-equilibrium steady states kubo . Here we show time reversal symmetry breaking can be inferred from temporal measurements of the particle positions.

Consider two observables $A$ and $B$ in a steady state, their correlation function is defined by

[TABLE]

where the last step above uses the fact that any steady state by definition must be invariant by time translations kubo .

The time reversal symmetry is broken when

[TABLE]

so that in this case the Onsager reciprocal relation kubo ,

[TABLE]

does not hold.

For an equilibrium system, as $C_{AB}(t)=C_{BA}(t)$, we find that $\tilde{C}_{AB}(\omega)$ is thus real. In the problem studied here we need to find two operators $A$ and $B$, with non-zero correlation to demonstrate the violation of time reversal symmetry. The obvious choice $A=z$ and $B=x$ (or $B=y$) is not useful as the corresponding correlators are zero. However choosing $B=x^{2}$ (or $B=y^{2}$) does yield a non-zero correlation function. However, to simplify the calculations which follow, we will choose $B=f_{n}(t)$, the fluctuating part of the component of the nonconservative force in the $z$ direction, which, see Eq. (164), is a function of $x^{2}$ and $y^{2}$ with zero mean, and from which the aforementioned correlators can trivially extracted. We define the correlation function

[TABLE]

which simplifies to

[TABLE]

because $\zeta_{ef}$ is independent of $f_{n}$ and $f_{n}$ is chosen to have zero mean. This now gives

[TABLE]

where $C(\omega)$ is given by Eq.(185). As $C(\omega)$ is real, we see that $\tilde{C}_{zf_{n}}(\omega)$ is not real and has an imaginary component

[TABLE]

Therefore if one can experimentally measure the imaginary component of this correlation function, one has another method of probing the non-equilibrium steady state of trapped particles in optical traps.

IV.3 Numerical simulations

Following the equations (155)-(157), where we adimensionalized the space and time variables by writing them in units of $a$ and $\omega_{0}=\sqrt{\kappa_{x}/m}$ respectively. We compute the correlation function $\langle X(\omega)Y(\omega^{\prime})\rangle=2\pi\delta(\omega+\omega^{\prime})S_{XY}(\omega)$ from a Fast Fourier Transform of the time variable functions $X(t)$ and $Y(t)$ over $2^{27}\sim 10^{8}$ points and $10^{3}$ stochastic averages where after discretization $S_{XY}(\omega_{i})=\langle X(\omega_{i})Y^{*}(\omega_{i})\rangle$ (where $.^{*}$ represents the complex conjugate). We verify here the exact equation for $S_{zz}$ as predicted in Eq. (197) and the time reversal breaking predicted in via the non-zero value of $C^{\prime\prime}_{zf}$ function as is given in Eq. (207).

V Conclusion

Optical trapping can be carried out in both liquid and gas phases. In the gas phase, controlling the pressure can be used to modify the friction coefficient of the trapped particle as the gas viscosity changes. The fact that the system is subject to a non-conservative force means that the resulting steady state depends on the precise form of the dynamics, in this case the damping coefficient $\gamma$. We have shown that this dependence on $\gamma$ can be demonstrated by measuring currents in the steady state, notably the currents associated with the marginal probability densities in both position and velocity space. Futhermore, the signature of optical scattering can be seen in the behavior of the spectral density function of the position and is particularly strong at low frequencies. We have also derived a correlation function between the particle position and the non-conservative force (equivalent to examining correlations between the position $Z(t)$ and the variable $X^{2}(t)$) which has a component which does not respect time reversal symmetry, another measurement of non-equilibrium behavior.

Many of the non-equilibrium results derived here have a non-monotonic dependence on the particle’s friction coefficient and are thus susceptible to being optimized to maximize their potential experimental signals, pointing to directions for future experimental studies.

A promising future direction of this study is to study the full nonlinear model for the optical trap to see what modifications are engendered by non-linear effects. However our numerical studies of the non-linear model (see accompanying paper letter ) suggest that the effects predicted by the harmonic model here persist, and are thus robust. Finally, even within the harmonic model here, the non-linear form of the conservative force means that we have only computed the first two moments of the temporal correlation functions and our results for the modified steady state distribution are perturbative. However it would of be interest to compute the full distribution function even for this simple harmonic model.

Acknowledgements.

This work was partially funded by the Bordeaux IdEX program - LAPHIA (ANR-10-IDEX-03-02) and the ANR project RaMaTRaF.