∎

e1e-mail<EMAIL_ADDRESS>
11institutetext: Institut für Theoretische Physik II: Weiche Materie, Heinrich
Heine-Universität Düsseldorf, Universitätsstraße 1, 40225 Düsseldorf, Germany
22institutetext: University of Lille, UGSF CNRS UMR8576, 59000 Lille, France

We discuss the dynamics of a Brownian particle under the influence of a
spatially periodic noise strength in one dimension using analytical theory and
computer simulations. In the absence of a deterministic force, the Langevin
equation can be integrated formally exactly. We determine the short- and long-
time behaviour of the mean displacement (MD) and mean-squared displacement
(MSD). In particular we find a very slow dynamics for the mean displacement,
scaling as $t^{-1/2}$ with time $t$. Placed under an additional external
periodic force near the critical tilt value we compute the stationary current
obtained from the corresponding Fokker-Planck equation and identify an
essential singularity if the minimum of the noise strength is zero. Finally,
in order to further elucidate the effect of the random periodic driving on the
diffusion process, we introduce a phase factor in spatial noise with respect
to the external periodic force and identify the value of the phase shift for
which the random force exerts its strongest effect on the long-time drift
velocity and diffusion coefficient.

# Brownian particles driven by spatially periodic noise

Davide Breonie1,addr1  Ralf Blosseyaddr2  Hartmut Löwenaddr1

(Received: date / Accepted: date)

††journal: Eur. Phys. J. E

##  1 Introduction

Dating back to the important paper by Einstein in the annus mirabilis 1905
einstein_uber_1905 , the dynamics of Brownian particles has been in the focus
of statistical physics for more than 100 years now frey_brownian_2005 . The
constant interest in Brownian particles is basically inspired by the two
facts: First, their stochastic description requires fundamental principles
such as the Langevin or Smoluchowski picture such that they serve as
paradigmatic models which can be made systematically more complex. Second,
there is a variety of excellent realizations of Brownian particles including
mesoscopic colloidal particles in suspension hansen_liquids_1991 , random
walkers in the macroscopic world (such as xiong_pedestrian_2012 ) and in the
microscopic biological context codling_random_2008 , and even elements of the
stock exchange market tsekov_brownian_2013 . This facilitates a direct
comparison of the stochastic averages between the stochastic modelling and
real experimental data.

In its simplest one-dimensional form, the most basic model Langevin equation
for a particle trajectory $x(t)$ as a function of time $t$ is
$\dot{x}(t)=\sqrt{D}\eta(t)$ in which $\eta(t)$ is white noise with zero mean
and variance $\langle\eta(t)\eta(t^{\prime})\rangle=\delta(t-t^{\prime})$ and
$D>0$ is the diffusion constant. Here, $\langle...\rangle$ denotes a noise
average. With the initial position $x(t=0)=x_{0}$, the mean displacement
vanishes due to symmetry, $\langle x(t)-x_{0}\rangle=0$, and the mean-squared
displacement is purely diffusive, $\langle(x(t)-x_{0})^{2}\rangle=2Dt$.
Clearly, this basic equation can be extended towards more complicated
situations including an additional static external force, time-dependent
external forcing, higher spatial dimensions, and many interacting particles,
see hess_generalized_1983 ; nagele_dynamics_1996 ; dhont_introduction_1996 ;
hanggi_artificial_2009  for some reviews.

In particular, the basic equation can be generalized to a situation of
spatially dependent noise strength modelled by a positive function $D(x)$,
i.e. a space-dependent diffusion coefficient, such that the most basic model
for such processes is given by the Langevin equation

$\dot{x}(t)=\sqrt{D(x(t))}\eta(t).$ (1)

The special case of multiplicative noise where $\dot{x}(t)=-\kappa
x(t)\eta(t)$ with positive $\kappa$ volpe_effective_2016 , which is somehow
related to this model, documents already that the spatial dependence of the
noise gives rise to fundamentally new mathematical concepts also known as the
Itô-Stratonovich problem mannella_ito_2012 . The mathematical difficulties
associated with the formal treatment of Eq.(1) are subject to intense
discussion, see, e.g., the recent work by Leibovich and Barkai for the
specific choice of $D(x)$ as a power-law leibovich_infinite_2019  and numerous
other studies lancon_drift_2001 ; pesek_mathematical_2016 ;
farago_connection_2016 ; bray_random_2000 ; pieprzyk_spatially_2016 ;
berezhkovskii_communication_2017 ; kaulakys_modeling_2009 ;
aghion_infinite_2020 ; dos_santos_critical_2020 ; xu_heterogeneous_2020 ;
ray_space-dependent_2020 ; li_particle_2020 ; breoni_active_2021 ;
malgaretti_confined_2013 .

In this paper we consider a variant of this model in the context of the
discussion of particle motion in tilted potentials. There is a large
literature on this topic, see reimann_giant_2001 ; reimann_diffusion_2002 ;
sasaki_diffusion_2005 ; reimann_weak_2008 ; evstigneev_diffusion_2008 ;
evstigneev_interaction-controlled_2009 ; cheng_long_2015 ;
juniper_microscopic_2015 ; guerin_universal_2017 ; bai_diffusion_2018 ;
bialas_colossal_2021 . Our model considers overdamped Brownian particles
subject to an oscillating tilted potential and a space-dependent periodic
noise amplitude with the same wave vector $k$ as the force; furthermore, we
will ultimately also allow a shifted phase $\phi$ in the random force. In its
general form, the model is given by the Langevin equation

$\gamma\dot{x}(t)=-\nabla V(x)+\sqrt{2\gamma k_{B}T(x)}\eta(t),$ (2)

where $V(x):=-F_{0}\left(x+\epsilon\sin(kx)/k\right)$ is the potential,
$T(x):=T_{0}\left(1+\nu\cos(kx+\phi)\right)^{2}$ is the space dependant noise
strength, $\gamma$ is the friction coefficient, $F_{0}$ is the tilting force,
$T_{0}$ is a reference temperature, $\eta(t)$ is a white noise, as introduced
before, and $\epsilon$ and $\nu$ are dimensionless parameters. The critical
tilt in this model arises when $\epsilon=1$. In order to avoid an unphysical
negative noise strength we consider $0\leq\nu\leq 1$. The period of both the
force and the noise will be $L=2\pi/k$. We remark that the case $\nu=1$ plays
a special role insofar as there are special positions at which the noise is
zero. In absence of forces, the particle will therefore never cross these
positions but stay confined within a periodicity length $L$.

Our goal in this paper is to describe the particle dynamics as functions of
$\epsilon$, $\nu$ and $\phi$, either in the vicinity of the critical tilt, or
in the absence of the deterministic force, $F_{0}=0$, i.e. in the purely
spatial random noise case. Among our main results are the very slow dynamics
in the relaxation of the mean displacement (MD) and mean-squared displacement
(MSD) for long times in the $F_{0}=0$ case, an essential singularity in the
stationary current for $F_{0}\neq 0$ and $\epsilon\simeq\nu\simeq 1$ and the
phase value $\phi$ for which we have the largest increase in long-time drift
$v_{L}$ and diffusion constant $D_{L}$ for $\epsilon\neq 1$ and $\nu\neq 0$.
Our results have been obtained both with numerical and analytical methods.

The paper is organized as follows: in the beginning we focus on the free case,
for which we study the short- and long-time behaviour of MD and MSD, then we
proceed with the full model, including the tilted potential, for which we
study the stationary distribution and the dependence of long time diffusion
and drift on $\phi$ and $\nu$.

##  2 Free particle case

In the case of a vanishing external force ($F_{0}=0$), the Langevin equation
(2) now reads as

$\gamma\dot{x}(t)=\sqrt{2\gamma k_{B}T(x)}\eta(t),$ (3)

where we set $\phi=0$ without loss of generality. For a given representation
of the noise, this equation can be solved by direct integration in the
particular case of periodic boundary conditions (PBC) in which we identify
$x(t)\pm L$ with $x(t)$. The PBC correspond to a ring-like geometry of the
one-dimensional system.

$\displaystyle x(t)$ $\displaystyle=$
$\displaystyle\frac{2}{k}\arctan\left[\frac{\sqrt{1-\nu^{2}}}{1-\nu}\tan\left(k\sqrt{\frac{k_{B}T_{0}(1-\nu^{2})}{2\gamma}}\right.\right.$
$\displaystyle\left.\left.\times\int_{0}^{t}\eta(t^{\prime})dt^{\prime}+\arctan\left(\frac{\sqrt{1-\nu^{2}}}{1-\nu}\tan\left(\frac{kx_{0}}{2}\right)\right)\right)\right]$

and the limit of this solution for $\nu=1$ is

$\displaystyle
x(t)=\frac{2}{k}\arctan\left[k\sqrt{\frac{2k_{B}T_{0}}{\gamma}}\int_{0}^{t}\eta(t^{\prime})dt^{\prime}+\tan\left(\frac{kx_{0}}{2}\right)\right]\,.$

We remark here that for open boundary conditions the analysis is harder and we
were not able to find an analytical expression except for the special case
$\nu=1$. In this limit PBC and open boundaries are identical as the particle
can never trespass the points where the noise is zero.

Equations (2) and (2) can be used to express noise-averages of any power of
displacement. For an arbitrary moment
$M_{n}(t):=\langle(x(t)-x_{0})^{n}\rangle$ we obtain

$\displaystyle M_{n}(t)$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}\left\\{\frac{2}{k}\arctan\left[\frac{\sqrt{1-\nu^{2}}}{1-\nu}\right.\right.$
$\displaystyle\times\tan\left(k\sqrt{\frac{k_{B}T_{0}(1-\nu^{2})}{2\gamma}}W\right.$
$\displaystyle\left.\left.+\arctan\left(\frac{\sqrt{1-\nu^{2}}}{1-\nu}\tan\left(\frac{kx_{0}}{2}\right)\right)\right)\right]$
$\displaystyle\left.-x_{0}\right\\}^{n}\frac{\textrm{e}^{-\frac{W^{2}}{2t}}}{\sqrt{2\pi
t}}dW$

for $\nu\neq 1$ and

$\displaystyle M_{n}(t)$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}\left\\{\frac{2}{k}\arctan\left[k\sqrt{\frac{2k_{B}T_{0}}{\gamma}}W+\tan\left(\frac{kx_{0}}{2}\right)\right]\right.$
(7)
$\displaystyle\left.-x_{0}\right\\}^{n}\frac{\textrm{e}^{-\frac{W^{2}}{2t}}}{\sqrt{2\pi
t}}dW$

for $\nu=1$. Since we are going to focus on the mean displacement $\langle
x(t)-x_{0}\rangle$ and the mean-squared displacement
$\langle(x(t)-x_{0})^{2}\rangle$, we write the expressions for these two
moments $(n=1,2)$ explicitly:

$\displaystyle\langle x(t)-x_{0}\rangle$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}\frac{2}{k}\arctan\left[\frac{\sqrt{1-\nu^{2}}}{1-\nu}\right.$
(8)
$\displaystyle\times\tan\left(k\sqrt{\frac{k_{B}T_{0}(1-\nu^{2})}{2\gamma}}W\right.\left.\left.+\arctan\left(\frac{\sqrt{1-\nu^{2}}}{1-\nu}\right.\right.\right.$
$\displaystyle\left.\left.\left.\times\tan\left(\frac{kx_{0}}{2}\right)\right)\right)\right]\frac{\textrm{e}^{-\frac{W^{2}}{2t}}}{\sqrt{2\pi
t}}dW-x_{0}$

and

$\displaystyle\langle(x(t)-x_{0})^{2}\rangle$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}\left\\{\frac{2}{k}\arctan\left[\frac{\sqrt{1-\nu^{2}}}{1-\nu}\right.\right.$
(9)
$\displaystyle\times\tan\left(k\sqrt{\frac{k_{B}T_{0}(1-\nu^{2})}{2\gamma}}W+\arctan\left(\frac{\sqrt{1-\nu^{2}}}{1-\nu}\right.\right.$
$\displaystyle\left.\left.\left.\left.\times\tan\left(\frac{kx_{0}}{2}\right)\right)\right)\right]-x_{0}\right\\}^{2}\frac{\textrm{e}^{-\frac{W^{2}}{2t}}}{\sqrt{2\pi
t}}dW$

for $\nu\neq 1$ and

$\displaystyle\langle x(t)-x_{0}\rangle$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}\frac{2}{k}\arctan\left[k\sqrt{\frac{2k_{B}T_{0}}{\gamma}}W\right.$
$\displaystyle\left.+\tan\left(\frac{kx_{0}}{2}\right)\right]\frac{\textrm{e}^{-\frac{W^{2}}{2t}}}{\sqrt{2\pi
t}}dW-x_{0}$

and

$\displaystyle\langle(x(t)-x_{0})^{2}\rangle$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}\left\\{\frac{2}{k}\arctan\left[k\sqrt{\frac{2k_{B}T_{0}}{\gamma}}W\right.\right.$
$\displaystyle\left.\left.+\tan\left(\frac{kx_{0}}{2}\right)\right]-x_{0}\right\\}^{2}\frac{\textrm{e}^{-\frac{W^{2}}{2t}}}{\sqrt{2\pi
t}}dW$

for $\nu=1$.

###  2.1 Short-time behavior

We can use equations (8-2) to extract the short-time behavior of the MD and
MSD. Expanding the integrand in powers of $t$ using a Taylor series and
integrating the terms separately we obtain for the MD:

$\displaystyle\langle x(t)-x_{0}\rangle$ $\displaystyle=$
$\displaystyle-\frac{kk_{B}T_{0}t}{\gamma}\sin(kx_{0})\nu\left[(1+\nu)\cos^{2}\left(\frac{k}{2}x_{0}\right)\right.$
(12)
$\displaystyle\left.+(1-\nu)\sin^{2}\left(\frac{k}{2}x_{0}\right)\right]+\mathcal{O}\left(t^{2}\right)$

and for the MSD

$\displaystyle\langle(x(t)-x_{0})^{2}\rangle=\frac{2k_{B}T_{0}t}{\gamma}\left[(1+\nu)\cos^{2}\left(\frac{k}{2}x_{0}\right)\right.$
$\displaystyle\left.+(1-\nu)\sin^{2}\left(\frac{k}{2}x_{0}\right)\right]^{2}+\mathcal{O}\left(t^{2}\right).$
(13)

In the special limit $\nu=1$ we also add the second order correction as:

$\displaystyle\langle x(t)-x_{0}\rangle$ $\displaystyle=$
$\displaystyle-\frac{2kk_{B}T_{0}t}{\gamma}\cos^{2}\left(\frac{k}{2}x_{0}\right)\sin(kx_{0})$
(14)
$\displaystyle\times\left[1-\frac{6k^{2}k_{B}T_{0}t}{\gamma}\cos^{2}\left(\frac{k}{2}x_{0}\right)\cos(kx_{0})\right]$
$\displaystyle+\mathcal{O}\left(t^{3}\right)$

and

$\displaystyle\langle(x(t)-x_{0})^{2}\rangle$ $\displaystyle=$
$\displaystyle\frac{8k_{B}T_{0}t}{\gamma}\cos^{4}\left(\frac{k}{2}x_{0}\right)$
$\displaystyle\times\left[1+\frac{k^{2}k_{B}T_{0}t}{\gamma}\cos^{2}\left(\frac{k}{2}x_{0}\right)\right.$
$\displaystyle\left.\times\left(7-11\cos\left(kx_{0}\right)\right)\right]+\mathcal{O}\left(t^{3}\right).$

Clearly, the first-order correction of (14) and (2.1) coincides with equations
(12) and (2.1) in the limit $\nu\to 1$. Moreover for $\nu=0$ we recover the
white noise case solved by Einstein einstein_uber_1905 .

We now define an effective potential of the mean displacement as follows.
First we fix the effective force of the mean displacement such that this force
produces the same drift for a particle with space-independent noise as that
for a force-free particle under space-dependent noise, i.e. this effective
force can be viewed as a substitute source for the drift when only white noise
is considered. Hence we define this force $F_{M}(x)$ up to a friction
coefficient prefactor $\gamma$ as the first coefficient of the short-time
expansion of the MD

$\langle x(t)-x_{0}\rangle=a_{1}(x_{0},\nu)t+\mathcal{O}\left(t^{2}\right)$
(16)

as follows

$\displaystyle F_{M}(x)$ $\displaystyle:=$ $\displaystyle a_{1}(x,\nu)\gamma$
(17) $\displaystyle=$ $\displaystyle-
kk_{B}T_{0}\sin(kx)\nu\left[(1+\nu)\cos^{2}\left(\frac{k}{2}x\right)\right.$
$\displaystyle\left.+(1-\nu)\sin^{2}\left(\frac{k}{2}x\right)\right].$

The effective potential of the mean displacement is then defined by
$V_{M}(x)=-\int_{0}^{x}F_{M}(x^{\prime})dx^{\prime}$ yielding

$\displaystyle V_{M}(x)$ $\displaystyle=$ $\displaystyle
k_{B}T_{0}\nu\left[(1-\nu)\sin^{4}\left(\frac{k}{2}x\right)\right.$ (18)
$\displaystyle\left.-(1+\nu)\cos^{4}\left(\frac{k}{2}x\right)+1+\nu\right].$

This potential is shown in Fig. 1. Even though this potential is defined just
by the short-time expansion of the MD, it is still significant for any finite
time, as the particle is overdamped and feels at every time a short-time drift
depending only on its position.

Figure 1: Effective potential of the mean displacement $V_{M}$, obtained from
the short-time drift of the mean displacement, and space-dependent noise
$T(x)$ for $\nu=0.5$ and $\phi=0$ as functions of space $x$. While the
averaged MD tends to the minima of $V_{M}$, where the noise strength $T(x)$ is
largest, individual trajectories spend most of their time around the maxima of
$V_{M}$.

While the average mean displacement behaves according to $V_{M}$, moving over
time towards the regions where $V_{M}$ is smaller and the noise strength is
larger, we want to stress that individual trajectories will not accumulate in
the minima of $V_{M}$ but will instead freely move over all the domain,
spending most of their time in the maxima of $V_{M}$ instead. This is because
when particles reach such low-noise regions they take a longer time escaping,
as their fluctuations there are severely reduced.

###  2.2 Dynamics for finite and long times

Now we explore the behavior of the MD and MSD for finite and long times. First
we present an asymptotic analysis for the special case $\nu=1$. Then we use a
numerical solution of the integrals in (8) and (9) as well as computer
simulations of the original Langevin equation to obtain data for finite times
and arbitrary $\nu$.

####  2.2.1 Asymptotic analysis for $\nu=1$ for long times

Here we present an asymptotic analysis for the MD and MSD by starting from
Eq.(2) and using the asymptotic approximation

$\arctan(\theta)\simeq\frac{\pi}{2}\textrm{sign}(\theta)-\arctan\left(\frac{1}{\theta}\right),$
(19)

for large $\theta$. We now expand $\arctan\left(\frac{1}{\theta}\right)$ using
Euler’s formula chien-lih_8967_2005

$\arctan\left(\frac{1}{\theta}\right)=\frac{\theta}{\theta^{2}+1}+\mathcal{O}\left(\frac{1}{\theta^{3}}\right)$
(20)

and insert this expansion in Eq.(2) to obtain

$\displaystyle\langle x(t)-x_{0}\rangle$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}\frac{2}{k}\left[\frac{\pi}{2}\textrm{sign}\left(k\sqrt{\frac{2k_{B}T_{0}t}{\gamma}}W\right.\right.$
(21) $\displaystyle+\left.\tan\left(\frac{kx_{0}}{2}\right)\right)$
$\displaystyle\left.-\frac{k\sqrt{\frac{2k_{B}T_{0}t}{\gamma}}W+\tan\left(\frac{kx_{0}}{2}\right)}{1+\left(k\sqrt{\frac{2k_{B}T_{0}t}{\gamma}}W+\tan\left(\frac{kx_{0}}{2}\right)\right)^{2}}\right]$
$\displaystyle\times\frac{\textrm{e}^{-\frac{W^{2}}{2t}}}{\sqrt{2\pi t}}dW-
x_{0}+\mathcal{O}\left(t^{-3/2}\right),$

which yields

$\displaystyle\langle x(t)-x_{0}\rangle$ $\displaystyle=$
$\displaystyle\sqrt{\frac{\pi\gamma}{k_{B}T_{0}t}}\frac{1}{k^{2}}\tan\left(\frac{kx_{0}}{2}\right)$
$\displaystyle\times\left(1-\frac{1}{2k}\sqrt{\frac{\gamma}{k_{B}T_{0}t}}\right)-x_{0}+\mathcal{O}\left(t^{-3/2}\right).$

As a result, the leading asymptotic behavior of $\langle x(t)\rangle$ is
determined by the first term involving a scaling behavior of the MD in
$1/\sqrt{t}$. This is remarkably slow compared to typical behavior of a
Brownian particle in a harmonic potential or of active Brownian motion where
the MD reaches its asymptotic value exponentially in time howse_self-
motile_2007 ; ten_hagen_brownian_2011 ; sprenger_time-dependent_2021  thus
constituting an example of a very slow relaxation as induced by space-
dependent noise.

Likewise an asymptotic analysis for $\nu=1$ yields for the long-time limit of
the MSD

$\lim_{t\rightarrow\infty}\left(\langle(x(t)-x_{0})^{2}\rangle\right)=x_{0}^{2}+\frac{\pi^{2}}{k^{2}}.$
(23)

which represents the degree of smearing of the particle distribution for long
times.

####  2.2.2 Computer simulations

We performed direct Brownian dynamics computer simulations of the original
Langevin equations with a finite time step $\Delta t$ to obtain numerically
results for the MD and MSD at any times. In order to properly simulate a
system with space-dependent noise, we used the order $\mathcal{O}\left(\Delta
t\right)$ Milstein scheme milshtejn_approximate_1975  with a time step of
$\Delta t=10^{-3}\tau$, where $\tau:=\frac{\gamma L^{2}}{k_{B}T_{0}}$ is a
typical Brownian time scale of the system. For each simulation set we fixed
the initial position $x_{0}$ within the first period $[-L/2,L/2]$ and averaged
typically over 200 trajectories of length $\simeq 500\tau$.

####  2.2.3 MD and MSD for finite times

Figure 2: Absolute values of the mean displacement $\langle x(t)-x_{0}\rangle$
(a) and the mean position with a minus sign to ensure positivity $-\langle
x(t)\rangle$ (b) for $\nu=1$ and $x_{0}=-0.2L$ as a function of time $t$. The
numerical evaluation of the integral in Eq. (2) (theory) and its asymptotic
short- and long-time expansions (14) and (21) are shown together with
simulation data. The MD increases linear in time $t$ for short times, while
the decay to its limit scales in a very slow way with
$\mathcal{O}\left(t^{-1/2}\right)$. Figure 3: (a) Probability density function
$p(x,t)$ for the particle position at different times $t$, with $\nu=1$ and
$x_{0}=-0.2L$. Here we averaged over 10000 different trajectories of length
$10\tau$. (b) Hopping rate $\Gamma_{h}$ between the two peaks in the particle
distribution as a function of time $t$ for a space-dependent noise with points
of vanishing noise ($\nu=1$), and for a double-well potential with white
noise. Here we have chosen $x_{0}=-0.2L$. We also show the stationary state
theoretical value of $\Gamma_{h}$ for the double well potential, defined as
the inverse of the mean first passage time $t_{e}$, derived in Eq.(26).

Data for the mean displacement and the mean position as a function of time are
obtained by a numerical evaluation of the integral in Eq. (2) and by computer
simulation. For $\nu=1$ results are presented in Fig. 2 together with the
corresponding short-time and long-time asymptotics (14) and (2.2.1). The
displacement starts linear in time $t$ and saturates for long times. The mean
position approaches zero slowly as a power law in time proportional to
$t^{-1/2}$. For large times the statistical error in the simulation data is
significant but nevertheless these data are compatible with the scaling
prediction of the theory.

In order to understand the very slow behavior of the MD we note that while the
MD tends to zero, i.e. to the point with largest noise, this is just an effect
of averaging over particles spending most of their time at the points with the
smallest noise on both sides of the $x$-axis: $x\simeq-L/2$ and $x\simeq L/2$.
This particular mechanism explains why the MD approaches its final value so
slowly, as the particles have to hop from one side to the other to symmetrize
their distribution. In Fig. 3a this is clearly documented in the time
evolution of the particle distribution function $p(x,t)$, which gives the
probability to find a particle after a time $t$ at position $x$ provided it
started at time $t=0$ at position $x_{0}$. The system evolves from a single-
peaked distribution around $x_{0}$ to a double-peaked distribution in $\pm
L/2$. Near the two points $x=\pm L/2$ of zero noise the peaks are getting
sharper as $t\to\infty$ approaching to $\delta$-peaks such that
$\lim_{t\to\infty}p(x,t)=(\delta(x-L/2)+\delta(x+L/2))/2$. The intuitive
reason for this is that once a particle adsorbs at the points $x=\pm L/2$ of
zero noise it will never return to the region where the noise is finite.
This peculiar behavior is clearly delineated from the relaxation in a
symmetric double-well potential with white noise of strength $T_{0}$. In order
to reveal this, we have performed simulations for a Brownian particle in the
double-well potential with two equal minima

$U(x):=A(x^{4}-Bx^{2}).$ (24)

We set $A:=48k_{B}T_{0}/L^{4}$ and $B:=L^{2}/2$ in order to have the two wells
in $\pm L/2$ such that the energy barrier between the two minima is
$3k_{B}T_{0}$. Our simulation for this white-noise reference case show that
both the MD and the MSD decay exponentially in time $t$ rather than with
$1/\sqrt{t}$, and hence much faster than for our case of space-dependent
noise. We also defined a particle hopping rate $\Gamma_{h}$ between the two
peaks of the distribution as

$\Gamma_{h}(t):=\frac{N_{h}(t)}{t/2},$ (25)

where $N_{h}(t)$ is the number of times a particle hops from one peak to the
other in the time interval $[t/2,t]$. Note that the relevant time window in
which hopping is considered is chosen to be proportional in time in order to
improve the statistics. We have a hop whenever the particle trespasses the
$x=L/4$ or $x=-L/4$ thresholds and previously was respectively in the left or
right peak.
In fact, as we show in Fig. 3b, for the double-well potential, the hopping
rate $\Gamma_{h}(t)$ converges to a constant for long times. This rate is
maintaining the equilibrium state with a symmetrized occupation around the two
minima. The rate saturates for $t\to\infty$ to a value very close to the
inverse of the mean first passage time (see for example malgaretti_active_2021
) in the double-well potential $t_{e}$ caprini_correlated_2021 , which in our
case is given by:

$\displaystyle t_{e}$ $\displaystyle\simeq$
$\displaystyle\frac{2\pi}{\sqrt{\nabla^{2}U(L/2)|\nabla^{2}U(0)|}}\exp{\left(\frac{U(0)-U(L/2)}{k_{B}T_{0}}\right)}$
(26) $\displaystyle\simeq$ $\displaystyle 1.859\tau.$

Conversely, for our case of space-dependent noise, the hopping rate keeps
decreasing as a function of time again with an inverse power law $t^{-1/2}$.
This reflects the fact that the peaks of the space-dependent noise
distribution keep growing indefinitely as the particles get in average closer
to the points of zero noise.

Figure 4: Mean displacement $\langle x(t)-x_{0}\rangle$ (a,b) as a function of
time $t$ for three values of $\nu=0.25,0.5,0.75$ and $x_{0}=-0.2L$ for open
boundary conditions (a) and for periodic boundary conditions (b).

Now in Fig. 4 we explore the MD for the case $\nu\not=1$ where the particle
crosses the position of minimal noise. Here the boundary conditions do matter
and we distinguish here between open boundary conditions (Fig. 4a) with
infinitely many oscillations and periodic boundary conditions of a ring-like
geometry (Fig. 4b). While the short-time behavior is linear in time for both
kind of boundary conditions, the MD saturates for long times to a finite value
depending on $\nu$ and $x_{0}$ for open boundary conditions. This finite value
is $-x_{0}$ for periodic boundary conditions since in this case the mean
position will always end at zero due to symmetry. The asymptotic approach to
zero is exponential in time as in the case of the double-well potential with
noise as the particle stays mobile even when approaching the position where
the noise is minimal. This is in marked contrast to the limit of $\nu=1$ where
the particle gets immobilized at the boundaries.

Now we turn to the MSD, first for the special case $\nu=1$ shown in Fig. 5a
where boundary conditions do not matter. The MSD starts linear in time and
then saturates to its long-time limit $C:=x_{0}^{2}+(\pi/k)^{2}$. Its
asymptotic approach to this saturation value is revealed by plotting the MSD
shifted by $C$ which decays to zero for large times, see Fig. 5b. Similar to
the MD for $\nu=1$, we find that the asymptotic behavior is compatible with a
$1/\sqrt{t}$ scaling.

In Fig. 6 we show the MSD for $\nu\not=1$ for both types of boundary
conditions. Open boundary conditions (see Fig. 6a) lead to a long-time
behavior linear in time $\approx 2D_{L}t$ involving a long-time diffusion
coefficient $D_{L}$. Clearly the latter depends on $\nu$ but not on the
initial position $x_{0}$. This dependence is depicted in the inset of Fig. 6a.
We found the empirical expression $D_{L}(\nu)=D_{0}(1-\nu^{2})$ with
$D_{0}=k_{B}T_{0}/\gamma$ to be a very good fit to the data. This can be
regarded as a parabolic fit which fulfills the inflection symmetry in $\nu$
and the constraint $D_{L}(\nu=1)=0$. The same behavior was recently found in a
similar system caprini_2021 .

Figure 5:  Absolute values of the mean-squared displacement (MSD)
$\langle(x(t)-x_{0})^{2}\rangle$ (a) and shifted MSD$-C$, where $C$ is the
long time limit of the MSD (b) as a function of time $t$ for $\nu=1$ and
$x_{0}=-0.2L$. The numerical evaluation of the integral in Eq. (2) and its
asymptotic short- and long-time expansions are shown together with simulation
data.  Figure 6: (a) Mean-squared displacement (MSD)
$\langle(x(t)-x_{0})^{2}\rangle$ as a function of time $t$ for $\nu\neq 1$ and
$x_{0}=-0.2L$ both for open boundary conditions (a) and for periodic boundary
conditions (b). The inset shows the long-time diffusion constant $D_{L}$ as a
function of $\nu$ for the open boundary conditions. Figure 7: Dynamical
exponents of the MD (a,c) and MSD (b,d) close to 0 (a,b) and their final limit
(c,d) for $\nu=1$ and $x_{0}=-0.2L$ in theory and simulation as functions of
time $t$. As we have already seen in Fig. 2 for the MD and Fig. 4 for the MSD,
both quantities grow initially linearly in time and decay to their final limit
with $1/\sqrt{t}$ for the MD and slower than $1/\sqrt{t}$ for the MSD.

Finally, to better clarify the behaviors of the MD and MSD for $\nu=1$, we
plot the dynamical exponents (Fig.7) that define the scaling regimes for the
MD ($\beta$, $\beta^{\prime}$) and MSD ($\alpha$, $\alpha^{\prime}$) close to
their short-time and long-time limits, respectively:

$\displaystyle\beta$ $\displaystyle:=$ $\displaystyle\frac{d(\log_{10}|\langle
x(t)-x_{0}\rangle|)}{d(\log_{10}(t))}\,,\,\,\,$ $\displaystyle\beta^{\prime}$
$\displaystyle:=$ $\displaystyle\frac{d(\log_{10}|\langle
x(t)\rangle|)}{d(\log_{10}(t))},$ $\displaystyle\alpha$ $\displaystyle:=$
$\displaystyle\frac{d(\log_{10}\langle(x(t)-x_{0})^{2}\rangle)}{d(\log_{10}(t))}\,,\,\,\,$
$\displaystyle\alpha^{\prime}$ $\displaystyle:=$
$\displaystyle\frac{d(\log_{10}|\langle(x(t)-x_{0})^{2}\rangle-C|)}{d(\log_{10}(t))}.$

Both the MD and MSD for short times are linear, while for long times the
scaling of the MD converges clearly to -0.5, that corresponds to $1/\sqrt{t}$.
Within the time window explored the MSD has not yet saturated to an ultimate
dynamical exponent for long times. The asymptotics shown is compatible with a
final scaling exponent of $-1/2$ although the approach to this final exponent
is much slower for the MSD than for the MD where the saturation is clearly
visible.

We remark that an algebraic asymptotic approach in the MSD was also found for
equilibrium Brownian dynamics of repulsive interacting particles. Here the
time-derivative of the time-dependent diffusion coefficient MSD$/t$ scales as
$t^{-d/2}$ in $d$ spatial dimensions cichocki_time-dependent_1991 ;
ackerson_correlations_1982 ; lowen_structure_1992 ; kollmann_single-file_2003
but the physical origins of the algebraic scaling laws are different.

##  3 Tilted potential

In this section we leave the situation in which the Brownian particle is a
free particle only driven by spatially-dependent noise. We now consider the
full model, including the deterministic tilted potential. We first look at the
situation near the critical value of the amplitude $\epsilon=1$, where the
tilted potential develops a plateau. The situation addressed in shown in Fig.
8.

###  3.1 The stationary current

Being weakly confined to a region of the deterministic potential in which the
dynamics can be considered ‘slow’, a quasi-stationary distribution can be
defined guerin_universal_2017 .

Figure 8: Potential (a) and corresponding force (b) near the ‘flat’ regime
$\epsilon=1$ as a function of the spatial coordinate $x$.

The Fokker-Planck equation corresponding to the Langevin equation, Eq. (1)
reads as

$\displaystyle\partial_{t}p(x,t)$ $\displaystyle=$
$\displaystyle-\partial_{x}a(x)p(x,t)+\frac{1}{2}\partial_{x}[b(x)\partial_{y}[b(x)p(x,t)]]$

with $a(x)$ the force and $b(x)$ the noise amplitude,

$a(x):=F_{0}(1+\epsilon\cos(kx))\,,\,\,\,b(x):=\sqrt{2\gamma k_{B}T(x)}\,.$
(29)

Following the discussion in guerin_universal_2017 , the dynamics near the
critical tilt value for $\epsilon\geq 1$ is characterized by a stationary
current given by the one-time integrated FP-equation

$-J_{s}=-a(x)p_{s}(x)+\frac{1}{2}b(x)\partial_{x}[b(x)p_{s}(x)]\,.$ (30)

Defining $(b(x)/2)p_{s}(x)=\widehat{p}_{s}(x)$ we can rewrite the last
expression as

$-\frac{J_{s}}{b(x)}=-R(x)\widehat{p}_{s}(x)+\partial_{x}\widehat{p}_{s}(x)$
(31)

with

$R(x)=\frac{2a(x)}{b^{2}(x)}\,.$ (32)

The equation can be solved with the Ansatz $\widehat{p}_{s}(x)=u(x)\cdot v(x)$
which reduces the problem to two readily integrable first-order ordinary
differential equations for $u(x)$ and $v(x)$. One obtains the final expression

$p_{s}(x)=\frac{2J_{s}}{b(x)}\int_{x}^{\infty}dy\frac{1}{b(y)}\exp{\left(-\int_{x}^{y}dzR(z)\right)}\,,$
(33)

in which the current $J_{s}$ can be obtained from the normalization integral
$\int_{-\infty}^{\infty}dxp_{s}(x)=1$. In the following we take for simplicity
(setting all other constants to one)

$R(z)=\frac{1+\epsilon\cos(z)}{(1+\nu\cos(z))^{2}}\,.$ (34)

Setting $b(z):=\exp\left(-F(z)\right)=\exp(-\ln(1+\nu\cos(z))$, and expanding
both $b(z)$ and $R(z)$ in Taylor series around the center of the flat region
near $z=L/2$, the stationary current $J_{s}$ is given by

$\displaystyle J_{s}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left[\int_{-\infty}^{\infty}dx\exp{\left(-\widehat{F}(x)\right)}\int_{0}^{\infty}dy\exp{\left(-\widehat{F}(y)\right)}\right.$
(35)
$\displaystyle\left.\times\exp{\left(-\int_{x}^{y}dz\widehat{R}(z)\right)}\right]^{-1}\,$

in which the symbol $\widehat{...}$ indicates the Taylor-expanded functions,

$\widehat{F}(x)=\ln(1-\nu)+\frac{1}{2}\left(\frac{1}{1-\nu}\right)(x-L/2)^{2}$
(36)

and

$\widehat{R}(z)=2\frac{1-\epsilon}{(1-\nu)^{2}}+\frac{\epsilon(\nu+1)-2\nu}{(1-\nu)^{3}}(z-L/2)^{2}\,.$
(37)

The integration of $\widehat{R}(z)$ yields a cubic polynomial, but due to
cancellations the resulting expression in the exponential is Gaussian in $x$
and cubic in $y$. The Gaussian integral in $y$ can be calculated exactly,
while the remaining expression in $y$ needs to be evaluated numerically for
each value of $\epsilon$ and $\nu$.

The most interesting behaviour of the stationary current is found in the limit
$\nu\rightarrow 1$, $\epsilon\approx 1$. The fact that the coefficients in
Eqs.(36),(37) are singular in $1/(1-\nu)$ leads to a singular behaviour of
$J_{s}$ in the form

$J_{s}\propto(1-\nu)^{m}\exp\left[-\frac{I(\epsilon,\nu)}{(1-\nu)^{n}}\right]\,,$
(38)

with $n=3$, since the dominant singularity in $\widehat{R}(z)$ is
$\propto(1-\nu)^{-3}$, see Eq.(37). The amplitude is
$I(\epsilon,\nu)=(\epsilon(\nu+1)-2\nu)/4$ and the rational factors combine to
$m=0$. The stationary current thus goes to zero with an essential singularity
in $(1-\nu)$.

###  3.2 Phase difference between noise and potential

For a tilted potential, we now explore the effect of a non-zero phase
$\phi\neq 0$ on the long-time behavior of the particle for different values of
$\nu$ by using computer simulation.
As shown in reimann_giant_2001 , the long-time drift velocity and diffusion
coefficients ($v_{L}$ and $D_{L}$ respectively) can be analytically calculated
for the case $\nu=0$, where we set $V(x)$ as potential.
Here the question is how the mismatch of the periodic noise and external
forcing affects the long-time behavior of the particle. Intuitively one would
expect that overcoming an energetic barrier is best if the maximum of the
noise occurs where the external force is opposing most. Then the noise would
help to bring the particle over the energetic barrier. The position where the
force is opposing most is clearly given for $x=L/2+nL$, where $n$ is an
integer. Then it is expected that mobility gets a maximum if the phase shift
is $\phi=\pi$. This is indeed what we confirm by simulation. We chose
$k_{B}T_{0}=0.01F_{0}L$ and $\epsilon=1.3$. The potential barrier $\Delta E$
is given by

$\displaystyle\Delta E(\epsilon)$ $\displaystyle=$
$\displaystyle\frac{F_{0}L}{\pi}\left(\sqrt{\epsilon^{2}-1}-\textrm{arcsec}(\epsilon)\right),$
(39)

yielding $\Delta E\simeq 0.04F_{0}L>0.01F_{0}L$ for $\epsilon=1.3$.

Given these parameters, we simulated the system for different values of $\phi$
and $\nu$ and results are summarized in Fig. 9. We have compared the
simulation data with a mapping on the analytical results for $v_{L}^{(0)}$ and
$D_{L}^{(0)}$ reimann_giant_2001  which were obtained for a spatially constant
temperature. Since the crucial position to hop over the barrier is at
$x=L/2+nL$ where the opposing force is maximal, this represents the kinetic
bottleneck for the dynamical process. Therefore it is tempting to compare our
simulation results with the analytical ones where this local noise strength
$T(x=L/2)$ is inserted as a homogeneous temperature. This mapping theory
should work best if the particle spends most of its time close to the point
$x=L/2$. In fact, Fig. 9 reveals that this simple mapping theory describes the
simulation data well even for large $\nu$. As a function of the phase mismatch
$\phi$, both $D_{L}$ and $v_{L}$ are enhanced when $\phi$ is between about
$\frac{3}{5}\pi$ and $\frac{7}{5}\pi$. Clearly around the value $\phi=\pi$ we
find the maximal enhancement of both $v_{L}$ and $D_{L}$. In the complementary
case, the noise strength $T(x)$ has its minimum closer to the crucial region
where the opposing force is maximal, and as a result the drift velocity and
diffusion are severely reduced. For $\nu=1$ they are even brought exactly to
zero when $|\phi|<\textrm{arcsec}(\epsilon)$, since the particle is stuck and
there is no systematic external force to drift over the positions of vanishing
noise.

Figure 9: Long time drift (a) and diffusion (b) for $\epsilon=1.3$ as
functions of $\nu$ and $\phi$, simulations results and theory. For $\phi$
between about $\frac{3}{5}\pi$ and $\frac{7}{5}\pi$ both $D_{L}$ and $v_{L}$
are enhanced, having a maximum in $\phi=\pi$, while otherwise they are
reduced.

##  4 Conclusions and outlook

In conclusion we have presented a detailed study of a model for a Brownian
particle moving in a one-dimensional environment with a space-periodic noise
and under an external potential with a tilt near its critical value. In the
free case we calculated the exact solution of the associated Langevin
equation, and further explicitly obtain short- and long-time approximations of
the MD and MSD. These results allow us to characterize the slow decay of these
quantities at long times. Interesting relaxation dynamics occurs around points
of vanishing noise which establish centers of growing peaks in the particle
distribution, as particles slow down significantly in the neighborhoods of
these points. Introducing the tilted periodic potential we first determined
the stationary current for the quasi-stationary state, which for $\epsilon\geq
1$ displays an essential singularity for the maximal strength of the noise
oscillations, $\nu$. Finally, we determined numerically the effects of a
space-periodic noise on the long-time diffusion and drift as functions of the
phase difference between noise and potential $\phi$ and the strength of the
noise oscillations $\nu$, finding the largest enhancements to take place for a
phase of $\phi=\pi$ and the maximal possible noise oscillations for $\nu=1$.

Our one-dimensional model with both periodic and open boundary conditions can
be realized by a colloidal particle confined in a ring or a linear channel by
e.g. optical forces lutz_single-file_2004 ; lutz_diffusion_2004 ;
juniper_microscopic_2015 ; juniper_dynamic_2017 . The space-dependent noise
can be artificially added by programming random magnetic fields to kick the
particles in the one-dimensional direction fernandez-rodriguez_feedback-
controlled_2020 . In this set-up a full function $D(x)$ can be programmed and
realized if the magnetic field strength depends on the spatial coordinate $x$.

##  5 Acknowledgements

We thank L Caprini for confirming the long-time behavior of the MSD described
in Fig. 6a in an analogous system and T Voigtmann for the interesting
discussions. DB is supported by the EU MSCA-ITN ActiveMatter, (proposal No.
812780). RB is grateful to HL for the invitation to a stay at the Heinrich-
Heine University in Düsseldorf where this work was performed. HL was supported
by the DFG project LO 418/25-1 of the SPP 2265.

## Author contribution statement

HL and RB directed the project. DB performed analytic calculations and
numerical simulations. RB contributed analytical results. All authors
discussed the results and wrote the manuscript.

## References

  * (1) Einstein A 1905 Annalen der Physik 322 549–560
  * (2) Frey E and Kroy K 2005 Annalen der Physik 14 20–50
  * (3) Hansen J P, Zinn-Justin J and Levesque D 1991 Liquids, Freezing and Glass Transition: Les Houches Session 51., 3-28 Juillet 1989 (North Holland) in chapter ”Colloidal Suspensions” by Pusey P N
  * (4) Xiong H, Yao L, Tan H and Wang W 2012 Discrete Dynamics in Nature and Society 2012 e405907
  * (5) Codling E A, Plank M J and Benhamou S 2008 Journal of The Royal Society Interface 5 813–834
  * (6) Tsekov R 2013 Chinese Physics Letters 30 088901
  * (7) Hess W and Klein R 1983 Advances in Physics 32 173–283
  * (8) Nägele G 1996 Physics Reports 272 215–372
  * (9) Dhont J K G 1996 An Introduction to Dynamics of Colloids (Amsterdam: Elsevier Science)
  * (10) Hänggi P and Marchesoni F 2009 Reviews of Modern Physics 81 387–442
  * (11) Volpe G and Wehr J 2016 Reports on Progress in Physics 79 053901
  * (12) Mannella R and McClintock P V E 2012 Fluctuation and Noise Letters 11 1240010
  * (13) Leibovich N and Barkai E 2019 Physical Review E 99 042138
  * (14) Lançon P, Batrouni G, Lobry L and Ostrowsky N 2001 EPL (Europhysics Letters) 54 28
  * (15) Pešek J, Baerts P, Smeets B, Maes C and Ramon H 2016 Soft Matter 12 3360–3387
  * (16) Farago O and Grønbech-Jensen N 2016 The Journal of Chemical Physics 144 084102
  * (17) Bray A J 2000 Physical Review E 62 103–112
  * (18) Pieprzyk S, Heyes D M and Brańka A C 2016 Biomicrofluidics 10 054118
  * (19) Berezhkovskii A M and Makarov D E 2017 The Journal of Chemical Physics 147 201102
  * (20) Kaulakys B and Alaburda M 2009 Journal of Statistical Mechanics: Theory and Experiment 2009 P02051
  * (21) Aghion E, Kessler D A and Barkai E 2020 Chaos, Solitons & Fractals 138 109890
  * (22) dos Santos M A F, Dornelas V, Colombo E H and Anteneodo C 2020 Physical Review E 102 042139
  * (23) Xu Y, Liu X, Li Y and Metzler R 2020 Physical Review E 102 062106
  * (24) Ray S 2020 The Journal of Chemical Physics 153 234904
  * (25) Li Y, Mei R, Xu Y, Kurths J, Duan J and Metzler R 2020 New Journal of Physics 22 053016
  * (26) Breoni D, Löwen H and Blossey R 2021 Physical Review E 103 052602
  * (27) Malgaretti P, Pagonabarraga I and Rubi J M 2013 The Journal of Chemical Physics 138 194906
  * (28) Reimann P, Van den Broeck C, Linke H, Hänggi P, Rubi J M and Pérez-Madrid A 2001 Physical Review Letters 87 010602
  * (29) Reimann P, van den Broeck C, Linke H, Hänggi P, Rubi J M and Pérez-Madrid A 2002 Physical Review E 65 031104
  * (30) Sasaki K and Amari S 2005 Journal of the Physical Society of Japan 74 2226–2232
  * (31) Reimann P and Eichhorn R 2008 Physical Review Letters 101 180601
  * (32) Evstigneev M, Zvyagolskaya O, Bleil S, Eichhorn R, Bechinger C and Reimann P 2008 Physical Review E 77 041107
  * (33) Evstigneev M, von Gehlen S and Reimann P 2009 Physical Review E 79 011116
  * (34) Cheng L and Yip N K 2015 Physica D: Nonlinear Phenomena 297 1–32
  * (35) Juniper M P N, Straube A V, Besseling R, Aarts D G A L and Dullens R P A 2015 Nature Communications 6 7187
  * (36) Guérin T and Dean D S 2017 Physical Review E 95 012109
  * (37) Bai Z W and Zhang W 2018 Chemical Physics 500 62–66
  * (38) Białas K and Spiechowicz J 2021 arXiv:2111.07661 [cond-mat]
  * (39) Chien-Lih H 2005 The Mathematical Gazette 89 469–470
  * (40) Howse J R, Jones R A L, Ryan A J, Gough T, Vafabakhsh R and Golestanian R 2007 Physical Review Letters 99 048102
  * (41) ten Hagen B, van Teeffelen S and Löwen H 2011 Journal of Physics: Condensed Matter 23 194119
  * (42) Sprenger A R, Jahanshahi S, Ivlev A V and Löwen H 2021 Physical Review E 103 042601
  * (43) Mil’shtejn G N 1975 Theory of Probability & Its Applications 19 557–562
  * (44) Malgaretti P, Puertas A M and Pagonabarraga I 2021 arXiv:2111.07645 [cond-mat, physics:physics]
  * (45) Caprini L, Cecconi F and Marconi U M B 2021 arXiv:2110.03042 [cond-mat]
  * (46) Caprini L 2021 to be submitted
  * (47) Cichocki B and Felderhof B U 1991 Physical Review A 44 6551–6558
  * (48) Ackerson B J and Fleishman L 1982 The Journal of Chemical Physics 76 2675–2679
  * (49) Löwen H 1992 Journal of Physics: Condensed Matter 4 10105–10116
  * (50) Kollmann M 2003 Physical Review Letters 90 180602
  * (51) Lutz C, Kollmann M and Bechinger C 2004 Physical Review Letters 93 026001
  * (52) Lutz C, Kollmann M, Leiderer P and Bechinger C 2004 Journal of Physics: Condensed Matter 16 S4075–S4083
  * (53) Juniper M P N, Zimmermann U, Straube A V, Besseling R, Aarts D G A L, Löwen H and Dullens R P A 2017 New Journal of Physics 19 013010
  * (54) Fernandez-Rodriguez M A, Grillo F, Alvarez L, Rathlef M, Buttinoni I, Volpe G and Isa L 2020 Nature Communications 11 4223


# Human Motion Tracking by Registering an Articulated Surface to 3-D Points
and Normals

Radu Horaud, Matti Niskanen, Guillaume Dewaele, and Edmond Boyer R. Horaud, G.
Dewaele, and E. Boyer are with INRIA Grenoble Rhône-Alpes, 655 avenue de
l’Europe, 38330 Montbonnot Saint-Martin, France.M. Niskanen is with the
Machine Vision Group, Infotech Oulu, University of Oulu, PO Box 4500, Finland.
M. Niskanen acknowledges financial support from the Seppo Synjkangas
foundation.

###### Abstract

We address the problem of human motion tracking by registering a surface to
3-D data. We propose a method that iteratively computes two things: Maximum
likelihood estimates for both the kinematic and free-motion parameters of a
kinematic human-body representation, as well as probabilities that the data
are assigned either to a body part, or to an outlier cluster. We introduce a
new metric between observed points and normals on one side, and a
parameterized surface on the other side, the latter being defined as a
blending over a set of ellipsoids. We claim that this metric is well suited
when one deals with either visual-hull or visual-shape observations. We
illustrate the method by tracking human motions using sparse visual-shape data
(3-D surface points and normals) gathered from imperfect silhouettes.

###### Index Terms:

model-based tracking, human motion capture, articulated implicit surface,
shape from silhouettes, robust surface registration, expectation-maximization.

##  I Introduction

We address the problem of recovering articulated human-motion parameters using
3-D data gathered from multiple image sequences. We advocate that this type of
data has several advantages over 2-D data: it is less ambiguous and it is less
sensitive to self-occlusions. 3-D features may be obtained by stereo [12, 5,
6]. Alternatively one can capitalize on 3-D shape from silhouettes. In
general, 2-D silhouettes are explicitly associated with a 3-D smooth surface
[12, 9, 10]. Another way to use silhouettes is to infer volumetric
representations and to fit articulated models to the voxels thus obtained [4,
13], or to extract skeletal representations from these voxels [3]. It is also
possible to infer 3-D surfaces from silhouettes, namely the visual hull [2] or
the visual shape [8]. The advantage of surface-from-silhouettes is that it
allows the recovery of both 3-D surface points and surface normals. Moreover,
there is no matching process associated with the reconstruction algorithm.
Visual hull algorithms have been proved to be extremely useful for recovering
3-D meshes which, in turn, are very useful for surface rendering. The drawback
is that they need perfect silhouettes. Alternatively, visual shape methods
(such as the one described in [8]) produce sparse surface descriptions (points
and normals) and can operate on imperfect silhouettes.

In this short paper we present a new method for tracking human motion based on
fitting an articulated implicit surface to 3-D points and normals. There are
two important contributions. First, we introduce a new distance between an
observation (a point and a normal) and an ellipsoid. We show that this can be
used to define an implicit surface as a blending over a set of ellipsoids
which are linked together to from a kinematic chain. Second, we exploit the
analogy between the distance from a set of observations to the implicit
surface and the negative log-likelihood of a mixture of Gaussian
distributions. This allows us to cast the problem of implicit surface fitting
into the problem of maximum likelihood (ML) estimation with hidden variables.
We argue that outliers are best described by a uniform component that is added
to the mixture, and we formally derive the associated EM algorithm.

Casting the data-to-model association problem into ML with hidden variables
has already been addressed in the past within the framework of point
registration [14], [11], [6]. In [5] observations are deterministically and
iteratively assigned to each individual body part. We appear to be the first
to apply a probabilistic data-to-model association framework to the problem of
fitting a blending of ellipsoids to a set of 3-D observations and to
explicitly model outliers within this context.

The remainder of the paper is organized as follows. Section II describes how
to compute a distance between a 3-D observation (point and normal) and an
ellipsoid, and how to build an implicit articulated surface based on this
distance. Finally, it introduces the concept of a probabilistic implicit
surface. Section III describes the formal derivation of the EM algorithm in
the case of implicit surface fitting. Section IV describes experiments with
simulated data and with multiple-camera video data.

##  II Modeling articulated objects

In order to model articulated objects such as human bodies we must define a
number of open kinematic chains that link the various body parts. We will use
ellipsoids for modeling these parts. Since we measure 3-D data (point and
orientation vectors) we must properly define a metric that measures the
discrepancy between the data and the model. This metric will be used to define
a distance function as well as a probabilistic implicit surface.

###  II-A The distance from a 3-D datum to an ellipsoid

One convenient way to describe 3-D ellipsoids is to use an implicit equation
and to embed the 3-D Euclidean space into the 3-D projective space. This
yields a 4$\times$4 full-rank symmetric matrix Q:

$\mbox{\bf Q}=\left[\begin{array}[]{cc}\overline{\mbox{\bf
Q}}&\mbox{\boldmath$q$}\\\
\mbox{\boldmath$q$}^{\top}&q_{44}\end{array}\right]=\left[\begin{array}[]{cc}\mbox{\bf
R}\mbox{\bf D}\mbox{\bf R}^{\top}&-\mbox{\bf R}\mbox{\bf D}\mbox{\bf
R}^{\top}\mbox{\boldmath$t$}\\\ -\mbox{\boldmath$t$}^{\top}\mbox{\bf
R}\mbox{\bf D}\mbox{\bf R}^{\top}&\mbox{\boldmath$t$}^{\top}\mbox{\bf
R}\mbox{\bf D}\mbox{\bf R}^{\top}\mbox{\boldmath$t$}-1\end{array}\right]$ (1)

where $\mbox{\bf D}=\mbox{Diag}\>[a^{-2},b^{-2},c^{-2}]$ is a $3\times$3
diagonal matrix, R is a $3\times 3$ rotation matrix and $t$ is a 3-D
translation vector. In practice $b=c$ and we choose $a\geq b$. We denote by
$X$ the homogeneous coordinates of a point $x$ lying onto the surface of the
ellipsoid, $\mbox{\boldmath$X$}^{\top}\mbox{\bf Q}\mbox{\boldmath$X$}=0$. The
adjoint matrix $\mbox{\bf Q}^{\star}=\mbox{\bf Q}^{-\top}$ defines the dual
ellipsoid. The family $P$ of planes which are tangent to the ellipsoid Q
satisfy the constraint $\mbox{\boldmath$P$}^{\top}\mbox{\bf
Q}^{-1}\mbox{\boldmath$P$}=0$ since $\mbox{\bf Q}^{\top}=\mbox{\bf Q}$. We
denote by $p$ the 3-D vector which is orthogonal to the plane $P$ (therefore
$p$ is normal to the ellipsoid at point $x$):
$\mbox{\boldmath$p$}=\overline{\mbox{\bf
Q}}\mbox{\boldmath$x$}+\mbox{\boldmath$q$}$, where the notations of eq. (1)
are used.

The algebraic distance from a 3-D point $Y$ to the surface of an ellipsoid was
used in [12], defined by
$q(\mbox{\boldmath$Y$})=\mbox{\boldmath$Y$}^{\top}\mbox{\bf
Q}\mbox{\boldmath$Y$}$. The value of $q$ varies from $-1$ at the center of the
ellipsoid, to $0$ on its surface, and then to $+\infty$ outside the ellipsoid
as the point is farther away from the surface. The Euclidean distance from a
point to an ellipsoid requires to solve a six-degree polynomial. In [6] an
approximation of the Euclidean distance is used, i.e, a pseudo-Euclidean
distance, as shown on Figure 1.

An observation will be refered to as a 3-D datum and consists of both a 3-D
point and a 3-D vector. We define a metric between such a 3-D datum and an
ellipsoid as follows. Let
$\mbox{\boldmath$Y$}^{\top}=(\mbox{\boldmath$y$}^{\top}\;1)$ be the
homogeneous coordinates of an observed point, and let $n$ be a 3-D observed
vector. An observation or a 3-D datum is denoted by
${\cal{Y}}=(\mbox{\boldmath$y$},\mbox{\boldmath$n$})$. We seek an ellipsoid
point $\mbox{\boldmath$X$}=(\mbox{\boldmath$x$}^{\top}\;1)$ under the
constraint that the vector $p$ (normal to the ellipsoid at $x$) is aligned
with $n$, e.g. Figure 1. In other words, we seek an association between
${\cal{X}}=(\mbox{\boldmath$x$},\mbox{\boldmath$p$})$ and
${\cal{Y}}=(\mbox{\boldmath$y$},\mbox{\boldmath$n$})$. Figure 2 compares the
distance used in this paper with the Euclidean distance from a point to an
ellipsoid.

Let $d_{{\cal{E}}}({\cal{Y}},{\cal{X}})$ be the Euclidean distance from the
datum-point $y$ to the ellipsoid-point $x$ under the constraint that the
datum-vector $n$ and the ellipsoid-vector $p$ are parallel:

$d_{{\cal{E}}}({\cal{Y}},{\cal{X}})=\|\mbox{\boldmath$x$}-\mbox{\boldmath$y$}\|_{2}\mbox{
with }\mbox{\boldmath$n$}\times\mbox{\boldmath$p$}=0$ (2)

where $\|\mbox{\boldmath$a$}\|_{2}$ denotes the Euclidean norm. We seek a
solution for $x$ under the constraints that $p$ and $n$ are parallel and yield
the same orientation. Using eq. (1) we obtain the following set of
constraints:

$\displaystyle\mbox{\boldmath$x$}^{\top}\overline{\mbox{\bf
Q}}\mbox{\boldmath$x$}+2\mbox{\boldmath$x$}^{\top}\mbox{\boldmath$q$}+q_{44}$
$\displaystyle=$ $\displaystyle 0$ (3) $\displaystyle\overline{\mbox{\bf
Q}}\mbox{\boldmath$x$}+\mbox{\boldmath$q$}$ $\displaystyle=$
$\displaystyle\lambda\mbox{\boldmath$n$}$ (4)

Figure 1: The distance from the datum
${\cal{Y}}=(\mbox{\boldmath$y$},\mbox{\boldmath$n$})$ to the ellipsoid Q is
estimated by seeking the point $\mbox{\boldmath$x$}\in\mbox{\bf Q}$ such that
the normal $p$ at $x$ is aligned with vector $n$.  Figure 2: The classical
Euclidean distance from a point to an ellipsoid (left) does not assign a point
to an ellipsoid in an unambiguous way. The 3-D datum distance (right) assigns
without ambiguity $a$ to A and $b$ to B.

From eq. (4) we obtain $\mbox{\boldmath$x$}=\overline{\mbox{\bf
Q}}\;^{-1}\left(\lambda\mbox{\boldmath$n$}-\mbox{\boldmath$q$}\right)$. By
substitution in eq. (3) we obtain two solutions for $\lambda$. From
$\mbox{\boldmath$p$}^{\top}\mbox{\boldmath$n$}>0$ we have $\lambda>0$ and
$\lambda=(\mbox{\boldmath$n$}^{\top}\mbox{\bf R}\mbox{\bf D}^{-1}\mbox{\bf
R}^{\top}\mbox{\boldmath$n$})^{-1/2}$. Therefore, the point onto the ellipsoid
where its normal $p$ is aligned with $n$ is given by:

$\mbox{\boldmath$x$}=\lambda\;\mbox{\bf R}\mbox{\bf D}^{-1}\mbox{\bf
R}^{\top}\mbox{\boldmath$n$}+\mbox{\boldmath$t$}$ (5)

It will be convenient to use the Mahalanobis distance:

$d^{2}_{\cal
M}(\mbox{\boldmath$y$},\mbox{\boldmath$x$})=(\mbox{\boldmath$y$}-\mbox{\boldmath$x$}(\mbox{\bf
R},\mbox{\boldmath$t$},\mbox{\boldmath$n$}))^{\top}\mbox{\bf$\Sigma$}^{-1}(\mbox{\boldmath$y$}-\mbox{\boldmath$x$}(\mbox{\bf
R},\mbox{\boldmath$t$},\mbox{\boldmath$n$}))$ (6)

###  II-B Kinematic chains and human-body modeling

Articulated motion has a long history in mechanics, biomechanics, robotics,
and computer vision. A human body can be described by a number of open
kinematic chains that share a common root. Such an open chain is composed of a
number of rigid objects and two consecutive rigid objects in the chain are
mechanically linked to form a joint. Rotational (or spherical) joints are the
most convenient representations and they are well suited for human body
modeling. Each such joint may have one, two, or three rotational degrees of
freedom. Therefore, within such a chain, a body part Q is linked to a root
body part $\mbox{\bf Q}_{r}$ through a constrained motion, i.e., a kinematic
chain with a number of rotational degrees of freedom. Since each joint may
have several degrees of freedom, the total number of rotational parameters of
a chain is larger than the number of rigid parts composing the chain.
Moreover, the root body part undergoes a free motion itself, i.e., a rigid
displacement with up to six degrees of freedom: three rotations and three
translations.

Therefore the motion of a body-part (or ellipsoid) Q is composed of the root’s
free motion followed by the chain’s constrained motion. We will denote the
motion of Q by the 4$\times$4 homogeneous matrix T which in turn is
parameterized by the joint and free-motion parameter vector $\Lambda$:

$\mbox{\bf T}(\mbox{{\boldmath$\Lambda$}})=\left[\begin{array}[]{cc}\mbox{\bf
R}(\mbox{{\boldmath$\Lambda$}})&\mbox{\boldmath$t$}(\mbox{{\boldmath$\Lambda$}})\\\
\mbox{\boldmath$0$}&1\end{array}\right]$ (7)

A complete human-body model may be described with five kinematic chains that
share a common root body-part. In this paper we use the following simplified
human-body model. There are 14 body parts and 11 joints (2 ankles, 2 knees, 2
hips, 2 elbows, 2 shoulders, and a neck) with $22$ rotational degrees of
freedom (there are two degrees of freedom for each joint). We also consider 3
rotations and 3 translations for the free motion. Hence, there is a total of
26 degrees of freedom.

As detailed above, body parts are described by one or several ellipsoids: The
feet and the thighs are described by two ellipsoids, the torso is described by
three ellipsoids, and all the other body parts are described by a single
ellipsoid, hence there are 21 ellipsoids, e.g., Figure 3. The body parts are
denoted by $\mbox{\bf Q}_{p},1\leq p\leq P$, where, for convenientce,
$\mbox{\bf Q}_{1}$ corresponds to the common root body part.

|  |  |
---|---|---|---
Figure 3: From left to right: The set of 21 ellipsoids used to model 14 body
parts with 11 joints and 2 rotations per joint. The implicit surface defined
as a blending of these ellipsoids. A set of 3-D “surface” observations (points
and normals) and the articulated implicit surface that has been fitted to
these observations.

###  II-C Articulated implicit surfaces

In addition to using a collection of kinematically linked ellipsoids, we will
fuse them in order to define a smooth surface $\mathcal{S}$ over the entire
body. This surface will be described by the implicit equation
$f(\mbox{\boldmath$y$})=C$ with $C=1$. The contribution of an ellipsoid
$\mbox{\bf Q}_{p}$ is defined by:

$f_{p}(\mbox{\boldmath$y$})=\exp\left(-\frac{d^{2}_{\mathcal{M}}(\mbox{\boldmath$y$},\mbox{\boldmath$x$}_{p})}{\nu_{p}^{2}}\right)$
(8)

where $d_{\mathcal{M}}$ is the Mahalanobis distance from $y$ to the ellipsoid
defined by eq. (6), the point $\mbox{\boldmath$x$}_{p}$ lies onto the
ellipsoid, and $\nu_{p}^{2}$ is a parameter that tunes the spatial influence
of the ellipsoid.

An implicit surface is defined as a level set of the following implicit
function which is the fusion (or blending) of $P$ ellipsoids verifying:

$f(\mbox{\boldmath$y$})=\sum_{p=1}^{P}f_{p}(\mbox{\boldmath$y$})$ (9)

The class of implicit surfaces defined as above, i.e.,
$\mbox{\boldmath$y$}\in\mathcal{S}\Leftrightarrow f(\mbox{\boldmath$y$})=C$,
has successfully been used in computer graphics and in computer vision in
conjunction with the algebraic distance [12] and with the pseudo-Euclidean
distance [6]. Within this paper we extend this concept to the distance defined
above. As it will be detailed below, this is well suited to cast the problem
of implicit surface fitting into the framework of maximum likelihood in the
presence of outliers.

In order to track articulated objects, the task at hand consists of fitting
the articulated implicit surface just described to a set of observations. For
this purpose, we first define a distance from a set of observations to the
implicit surface. We have to solve the equation $f(\mbox{\boldmath$y$})=1$,
where $y$ is, as before, an observed 3-D point. One may notice that the first
order Taylor expansion of $\ln a$ at $a=1$ is: $\ln a=a-1+O(a^{2})$. We retain
the following approximation of the distance from a set of $I$ observations to
the articulated implicit surface formed by $P$ ellipsoids and parameterized by
the kinematic variables $\Lambda$:

$F(\mbox{{\boldmath$\Lambda$}})=-\nu^{2}\sum_{i=1}^{I}\ln\sum_{p=1}^{P}\exp\left(-\frac{d_{ip}^{2}}{\nu^{2}}\right)$
(10)

where:

$d_{ip}^{2}=(\mbox{\boldmath$y$}_{i}-\mbox{\boldmath$x$}_{ip}(\mbox{{\boldmath$\Lambda$}},\mbox{\boldmath$n$}_{i}))^{\top}\mbox{\bf$\Sigma_{p}$}^{-1}(\mbox{\boldmath$y$}_{i}-\mbox{\boldmath$x$}_{ip}(\mbox{{\boldmath$\Lambda$}},\mbox{\boldmath$n$}_{i}))$

For convenience, we set $\nu=\nu_{1}=\ldots=\nu_{P}$. The notation
$\mbox{\boldmath$x$}_{ip}$ means that the 3-D point $x$ lies on ellipsoid $p$
and is associated with observation $i$. It is worthwhile to notice that,
whenever a set of observations is closed to one of the ellipsoids, the
distance function is strictly equal to the sum of Mahalanobis distances from
each such observation to the ellipsoid.

###  II-D Probabilistic implicit surfaces

In this section we introduce a probabilistic interpretation of eq. (10). For
this purpose we denote by $z_{i}$ a random variable that assigns an
observation $i$ to an ellipsoid $p$, namely the notation $z_{i}=p$ means that
the $i^{th}$ observation is assigned to the $p^{th}$ ellipsoid. There are as
many hidden variables $z_{i}$ as observations: $i\in\\{1,\ldots,I\\}$. The set
of all the hidden variables is denoted by
$\mathcal{Z}=\\{z_{1},\ldots,z_{I}\\}$.

The likelihood of an observed 3-D point, given its assignment to an ellipsoid
and given an observed 3-D normal, is drawn from a Gaussian distribution:

$P(\mbox{\boldmath$y$}_{i}|z_{i}=p,\mbox{\boldmath$n$}_{i})=\mathcal{N}(\mbox{\boldmath$y$}_{i}|\mbox{\boldmath$x$}_{ip}(\mbox{{\boldmath$\Lambda$}},\mbox{\boldmath$n$}_{i}),\mbox{\bf$\Sigma$}_{p})$
(11)

In practice the data are corrupted by noise and by errors and therefore there
are observations which should not be assigned to an ellipsoid. For this reason
we introduce an outlier class denoted by $P+1$, and we assume that the
likelihood of an observation given that is classified as an outlier is a
uniform distribution over the volume $V$ of the working space:

$P(\mbox{\boldmath$y$}_{i}|z_{i}=P+1,\mbox{\boldmath$n$}_{i})=\mathcal{U}(\mbox{\boldmath$y$}_{i}|V,0)=\frac{1}{V}$
(12)

Therefore, one can write the likelihood of an observation as a mixture of $P$
Gaussian components and one uniform component:

$P(\mbox{\boldmath$y$}_{i}|\mbox{\boldmath$n$}_{i})=\sum_{p=1}^{P+1}\pi_{p}P(\mbox{\boldmath$y$}_{i}|z_{i}=p,\mbox{\boldmath$n$}_{i})$
(13)

The notation:

$\pi_{p}=P(z_{i}=p|\mbox{\boldmath$n$}_{i})$ (14)

denotes the priors, the proportions, or the mixing parameters, and they obey
the obvious constraint $\sum_{p=1}^{P+1}\pi_{p}=1$. Notice that this prior
probability depends on the observed vector $\mbox{\boldmath$n$}_{i}$. In this
paper we do not treat these observed vectors as random variables. By assuming
independent and identically distributed observations one can write the joint
likelihood of all the observations as:

$P(\mathcal{Y}_{1},\ldots,\mathcal{Y}_{I})=P(\mbox{\boldmath$y$}_{1},\mbox{\boldmath$n$}_{1},\ldots,\mbox{\boldmath$y$}_{I},\mbox{\boldmath$n$}_{I})=\prod_{i=1}^{I}P(\mbox{\boldmath$y$}_{i}|\mbox{\boldmath$n$}_{i})P(\mbox{\boldmath$n$}_{i})$

Using Bayes’ formula and the equations above, the negative log-likelihood
writes:

$\displaystyle-$ $\displaystyle\ln
P_{\mbox{{\boldmath$\Lambda$}}}(\mathcal{Y}_{1},\ldots,\mathcal{Y}_{I})=$ (15)
$\displaystyle-$
$\displaystyle\sum_{i=1}^{I}\ln\left(\sum_{p=1}^{P}\pi_{p}\mathcal{N}(\mbox{\boldmath$y$}_{i}|\mbox{\boldmath$x$}_{ip}(\mbox{{\boldmath$\Lambda$}},\mbox{\boldmath$n$}_{i}),\mbox{\bf$\Sigma$}_{p})\right.$
$\displaystyle+$
$\displaystyle\left.\pi_{P+1}\mathcal{U}(\mbox{\boldmath$y$}_{i}|V,0)\phantom{\sum_{p=1}^{P}}\right)$

Notice that there is a strong analogy between eqs. (10) and (15): the former
is a distance between a set of $I$ observations and an articulated implicit
surface while the latter is the joint likelihood of the same observation set,
where the likelihood is a mixture of $P$ normal distributions plus a uniform
distribution that captures the bad observations. This analogy will be
exploited in the next section in order to cast the estimation of the kinematic
parameters in the framework of maximum likelihood with hidden variables via
the EM algorithm.

##  III Robust tracking with the EM algorithm

Because of the presence of the hidden variables,
$\mathcal{Z}=\\{z_{1},\ldots,z_{I}\\}$, the maximum-likelihood estimation
problem, i.e., eq. (15) does not have a simple solution. The most convenient
way to maximize the likelihood of a mixture of distributions is to use the EM
algorithm. The latter has been thoroughly studied in the context of data
clustering [7]. In this paper we formally derived an expectation-maximization
scheme in the particular case of robustly fitting an implicit surface to a set
of 3-D observations. It is worthwhile to notice that the formulae below are
valid independently of the distance function being used, i.e., Figure 1.

First we derive the posterior class probabilities conditioned by the
observations, namely:

$P(z_{i}=p|\mbox{\boldmath$y$}_{i},\mbox{\boldmath$n$}_{i})=\frac{P(z_{i}=p,\mbox{\boldmath$y$}_{i},\mbox{\boldmath$n$}_{i})}{P(\mbox{\boldmath$y$}_{i},\mbox{\boldmath$n$}_{i})}$

We denote these posteriors by $t_{ip}$ and with the notations introduced in
the prevous section we obviously obtain:

$t_{ip}=\frac{\pi_{p}P(\mbox{\boldmath$y$}_{i}|z_{i}=p,\mbox{\boldmath$n$}_{i})}{P(\mbox{\boldmath$y$}_{i}|\mbox{\boldmath$n$}_{i})}$
(16)

Second we consider the joint probability of the set of observations
$\mathcal{Y}=\\{\mathcal{Y}_{1},\ldots,\mathcal{Y}_{I}\\}$ and of their
assignments $\mathcal{Z}$ which yield the following expression:

$P(\mathcal{Y},\mathcal{Z})=\prod_{i=1}^{I}\prod_{p=1}^{P+1}\big{(}P(\mbox{\boldmath$y$}_{i}|z_{i}=p,\mbox{\boldmath$n$}_{i})P(z_{i}=p|\mbox{\boldmath$n$}_{i})\big{)}^{\delta_{p}(z_{i})}P(\mbox{\boldmath$n$}_{i})$

with the following definition for the function $\delta_{p}(z_{i})$:

$\delta_{p}(z_{i})=\left\\{\begin{array}[]{ccc}1&\text{if}&z_{i}=p\\\
0&&\text{otherwise}\end{array}\right.$

Third we derive the expression of the conditional expectation of the log-
likelihood taken over $\mathcal{Z}$, which in this case yields:

$\displaystyle E[\ln P(\mathcal{Y},\mathcal{Z})|\mathcal{Y}]=$
$\displaystyle\sum_{i=1}^{I}\sum_{p=1}^{P+1}E[\delta_{p}(z_{i})|\mathcal{Y}](\ln
P(\mbox{\boldmath$y$}_{i}|z_{i}=p,\mbox{\boldmath$n$}_{i})$
$\displaystyle+\ln\pi_{p})+(P+1)\sum_{i=1}^{I}\ln P(\mbox{\boldmath$n$}_{i})$

One may notice that:

$E[\delta_{p}(z_{i})|\mathcal{Y}]=\sum_{p=1}^{P+1}\delta_{p}(z_{i}=p)P(z_{i}=p|\mbox{\boldmath$y$}_{i},\mbox{\boldmath$n$}_{i})=t_{ip}$

By using the expressions of the normal and uniform distributions, and by
grouping constant terms we obtain:

$\displaystyle E[\ln P(\mathcal{Y},\mathcal{Z})|\mathcal{Y}]=$ (17)
$\displaystyle-$
$\displaystyle\frac{1}{2}\sum_{i=1}^{I}\left(\sum_{p=1}^{P}t_{ip}\big{(}(\mbox{\boldmath$y$}_{i}-\mbox{\boldmath$x$}_{ip})^{\top}\mbox{\bf$\Sigma_{p}$}^{-1}(\mbox{\boldmath$y$}_{i}-\mbox{\boldmath$x$}_{ip})\right.$
$\displaystyle+$
$\displaystyle\left.\ln\det\mbox{\bf$\Sigma_{p}$}-\ln\pi_{p}\big{)}+t_{iP+1}\ln\pi_{P+1}\vphantom{\sum_{p=1}^{P}}\right)+\text{const}$

The maximization of eq. (17) (or equivalently the minimization of its
negative) will be carried out via the EM algortihm (expectation-maximization).
There are, however three notable differences between the standard EM for
Gaussian mixtures [1] and our formulation:

  * •

We added a uniform-noise component to the mixture. The role of this component
is to “capture” outliers and hence to avoid that they influence the estimation
of the model parameters;

  * •

The means of the Gaussian components, $\mbox{\boldmath$x$}_{ip}$ are
parameterized by the kinematic parameters that control the articulated motion
of each ellipsoid; This has an important consequence because the M-step of the
algorithm will incorporate a non-linear minimization procedure over the
kinematic joints.

  * •

At the start of the algorithm each observation is associated with all the
ellipsoids. As the algorithm proceeds, each observation is eventually
associated with one of the ellipsoids. Due to occlusions, missing data, etc.,
there may be ellipsoids with no associated observation. Therefore, there is a
risk that the corresponding covariance becomes infinitely small. To overcome
this problem we use a unique covariance matrix common to all the densities in
the mixture.

Since we formally derived eq. (17), the EM algorithm outlined below guarantees
likelihood maximization. To summarize, the advantages of this formulation are
(i) fast convergence properties of EM and (ii) the fact that _it minimizes the
negative log-likelihood_ given by eqs. (15) and (17). In practice, the
following EM procedure can be used for robust tracking of an articulated
implicit surface:

  1. 1.

Initialization. Compute the locations of the ellipsoid points
$\mbox{\boldmath$x$}_{ip}^{(q)}$ from the current kinematic parameters
$\mbox{{\boldmath$\Lambda$}}^{(q)}$ using eq. (5). Similarly, initialize the
covariance matrix $\mbox{\bf$\Sigma$}^{(q)}$ common to all the ellipsoids.
Initialize the priors, or the mixing parameters
$\pi_{1}^{(q)}=\ldots=\pi_{P+1}^{(q)}=1/(P+1)$;

  2. 2.

E step. Evaluate the posterior probabilities $t_{ip}^{(q)}$ using the current
parameter values, through eq. (16);

  3. 3.

M step. Estimate new values for the kinematic parameters
$\mbox{{\boldmath$\Lambda$}}^{(q+1)}$:

$\displaystyle\arg\min_{\mbox{{\boldmath$\Lambda$}}}\frac{1}{2}\sum_{i=1}^{I}\sum_{p=1}^{P}t_{ip}^{(q)}(\mbox{\boldmath$y$}_{i}-\mbox{\boldmath$x$}_{ip}(\mbox{{\boldmath$\Lambda$}}))^{\top}{\mbox{\bf$\Sigma$}^{(q)}}^{-1}(\mbox{\boldmath$y$}_{i}-\mbox{\boldmath$x$}_{ip}(\mbox{{\boldmath$\Lambda$}}))$

Update the covariance matrix and the priors:

$\displaystyle\mbox{\bf$\Sigma$}^{(q+1)}$ $\displaystyle=$
$\displaystyle\frac{1}{\sum_{p=1}^{P}T_{p}}\sum_{i=1}^{I}\sum_{p=1}^{P}t_{ip}^{(q)}(\mbox{\boldmath$y$}_{i}-\mbox{\boldmath$x$}_{ip}(\mbox{{\boldmath$\Lambda$}}^{(q+1)}))$
$\displaystyle(\mbox{\boldmath$y$}_{i}-\mbox{\boldmath$x$}_{ip}(\mbox{{\boldmath$\Lambda$}}^{(q+1)}))^{\top}$
$\displaystyle\pi_{p}^{(q+1)}$
$\displaystyle=\frac{1}{I}\sum_{i=1}^{I}t_{ip}^{(q)}$

  4. 4.

Maximum likelihood. Evaluate the log-likelihood, i.e., eq. (15) and check for
convergence.

##  IV Experimental results

The tracking algorithm described in the previous section is not tight to any
particular method for extracting the 3-D data. In practice we used 3-D points
and normals that lie on the visual shape [8]. Notice that the visual-shape
algorithm does not require perfect silhouettes and provides as output a sparse
set of 3-D points and normals, not a 3-D mesh. The visual-shape method uses
multiple-view epipolar geometry in conjunction with the assumption that the
object’s surface is locally continuous and twice differentiable (see [8] for
details).

###  IV-A Experiments with simulated data

We used an animation package to build a human body, to simulate various human
motions, and to render image silhouettes. The simulator uses its own shape
representation, that is different than ours, but it allows the user to define
her/his own kinematic model. Therefore we used the same kinematic model with
the same number of degrees of freedom as the one described in section II-B.
Nevertheless, we have not attempted to finely tune the shape parameters of our
model to the simulated data. We simulated a setup composed of seven calibrated
cameras. Sequences of image silhouettes were generated from the 3-D model and
rendered with the cameras’ parameters. We computed 3-D points and normals from
these silhouettes and then we applied our method to these data sets. The
articulated-motion parameters were recovered using our tracker. In order to
assess the merits of the data-to-model fitting process, we added 20% of
outliers uniformly distributed in the volume of the working space. These
simulations allowed us to (i) assess the quality of the tracker with respect
to ground-truth joint trajectories, (ii) analyse the behavior of the method in
the presence of various perturbations that alter the quality of the data,
(iii) quantify the merit of using both 3-D points and normals instead of using
only points, and (iv) determine the optimal number of observations needed for
reliable estimation of articulated motion.

|
---|---
|
Figure 4: The error between simulated angle values and estimated ones (measured in radians), from left ot right: Ground-truth and measured trajectories over 100 frames. Average angle error as a function of silhouette noise for points and for points and normals. Average angle error as a function of frame rate. Average angle error as a function of the number of observations being used. |  |  |  |  |  |
---|---|---|---|---|---|---
|  |  |  |  |  |
|  |  |  |  |  |
|  |  |  |  |  |
Figure 5: The taekwendo sequence. The images and the associated silhouettes
are from the first and fourth cameras. Imperfect silhouettes generate outlier
data which are properly handled by our method.

We simulated a running sequence that involves both kinematic and free-motion
parameters. The graphs of Figure 4 illustrate the average error between the
simulated angle values and their estimated values (measured in radians) for
the kinematic parameters but not for the free-motion parameters. The first
graph compares the ground-truth (simulated) joint trajectories of the left
knee and of the right elbow (dashed curve) with the trajectories estimated
with our method (solid curve) over 100 frames. The second graph illustrates
the behavior of the method in the presence of silhouette noise. The results of
using both 3-D points and normals (dashed curve) are plotted against the
results obtained using 3-D points and the algebraic distance (solid line). The
relative large error corresponds to the fact that the shape model used by the
animation package is not the same as our shape model. Hence, there is a
systematic offset between the ground-truth kinematic parameters and the
estimated parameters.

The third graph shows the average angle error as a function of the number of
frames per second. The last graph shows the influence of the number of
observations, where the latter varies from 50 to 550. The average angle error
drops as the number of observations increases and our method (dashed) performs
better than using 3-D points alone (solid). From all these experiments one may
conclude that tracking is improved when both points and normals are used
instead of just points.

|  |  |  |  |  |
---|---|---|---|---|---|---
|  |  |  |  |  |
|  |  |  |  |  |
Figure 6: The 3-D visual hull (top), the implicit surface fitted to the taekwendo sequence using points and normals (middle), and using points alone (bottom). |  |  |  |  |  |
---|---|---|---|---|---|---
|  |  |  |  |  |
|  |  |  |  |  |
Figure 7: The leaning sequence: Images and silhouettes associated with the
first camera, and the fitted implicit surface.

###  IV-B Experiments with multiple-video data

The experimental data that we used for validating the human motion tracker was
gathered with six calibrated and finely synchronized cameras. Each camera
delivers 780$\times$580 color images at 28 frames per second with a
synchronization accuracy within 1 mili-second. The figures below show these
image sequences sampled at 14 frames per second.

We applied articulated human motion tracking to two multiple-image sequences,
the taekwendo sequence shown on Figure 5 and the leaning sequence shown of
Figure 7. The first data set is composed of 6$\times$700 frames while the
second one is composed of 6$\times$200 frames. We used the same body-part
dimensions for the two characters. One may notice that the silhouettes have
holes and missing pieces, which results in the presence of 3-D outliers. The
top row in Figure 6 shows the 3-D points and normals that were reconstructed
from the imperfect silhouettes; The middle row shows the articulated implicit
surface resulting from application of our method, while the bottom row shows
the same surface resulting from application of our method in conjunction with
the algebraic distance. Obviously, in this last case, there is a discrepancy
between the data and the fitted model: the recovered motion of the right feet
and the right thigh are incorrect. Similarly, figure 7 shows a sample of the
leaning sequence, the corresponding silhouettes, and the fitted model using
the proposed method.

##  V Conclusions

In this paper we described a method for tracking articulated motion with
several cameras. We introduced a new metric that measures the discrepancy
between observations (composed of both 3-D points and 3-D normals) and an
articulated implicit surface. This metric is more powerful than previously
used distance functions because it allows for less ambiguous associations
between the data and the model. Moreover, it is well suited when one deals
with either visual-hull or visual-shape representations of the data.

We cast the data-to-model fitting process into a robust probabilistic
framework. We showed that there is a strong similarity between the
mathematical representation of an implicit surface and a mixture of Gaussian
distributions. We explored this similarity and we showed that the articulated
motion tracking problem can be formulated as maximum likelihood with hidden
variables. We added a uniform component to the mixture to account for
outliers. We formally derived an algorithm that computes ML estimates for the
motion parameters within the framework of expectation-maximization (EM).
Therefore, the tracker may well be interpreted in the framework of robust data
clustering, where the observations are assigned to one of the ellipsoids, or
to an outlier component.

There are many questions that remain open and that we plan to investigate in
the near future: The algorithm may be trapped in local minima if it is not
properly initialized; There are some similarities between our robust tracker
and the use of M-estimators and these similarities diserve further
investigation. There are other interesting issues such as: A thorough and
quantitative evaluation of the results and their comparison with marker-based
motion capture systems and the possibility to capture several articulated
motions at once.

## References

  * [1] C.M. Bishop.  Pattern Recognition and Machine Learning.  Springer, 2006.
  * [2] E. Boyer and J-S. Franco.  A hybrid approach for computing visual hulls of complex objects.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 695–701, Madison, Wisconsin, USA, June 2003. IEEE Computer Society Press.
  * [3] G.J. Brostow, I.A. Essa, D. Steedly, and V. Kwatra.  Novel skeletal representation for articulated creatures.  In European Conference on Computer Vision, pages Vol III: 66–78, 2004.
  * [4] K.-M. Cheung, S. Baker, and T. Kanade.  Shape-from-silhouette across time part ii:applications to human modeling and markerless motion.  International Journal of Computer Vision, 63(3):225–245, 2005.
  * [5] D. Demirdjian.  Combining geometric- and view-based approaches for articulated pose estimation.  In European Conference on Computer Vision, pages Vol III: 183–194, 2004.
  * [6] G. Dewaele, F. Devernay, R. Horaud, and F. Forbes.  The alignment between 3-d data and articulated shapes with bending surfaces.  In Proceedings of the 9th European Conference on Computer Vision, Graz, Austria, volume III of LNCS, pages 578–591. Springer, May 2006.
  * [7] C. Fraley and A. E. Raftery.  Model-based clustering, discriminant analysis, and density estimation.  Journal of the American Statistical Association, 97:611–631, 2002\.
  * [8] J-S. Franco, M. Lapierre, and E. Boyer.  Visual shapes of silhouette sets.  In Proceedings of the 3rd International Symposium on 3D Data Processing, Visualization and Transmission, Chapel Hill (USA), 2006.
  * [9] S. Ilic, M. Salzmann, and P. Fua.  Implicit meshes for effective silhouette handling.  International Journal of Computer Vision, 72(2):159 – 178, 2007\.
  * [10] D. Knossow, R. Ronfard, and R. Horaud.  Human motion tracking with a kinematic parameterization of extremal contours.  International Journal of Computer Vision, 2008.  To appear.
  * [11] B. Luo and E.R. Hancock.  Structural graph matching using the EM algorithm and singular value decomposition.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 23(2001):1120–1136, October 2001.
  * [12] R. Plankers and P. Fua.  Articulated Soft Objects for Multi-View Shape and Motion Capture.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 25(9):1182–1187, 2003.
  * [13] J. Starck and A. Hilton.  Model-based multiple view reconstruction of people.  In IEEE International Conference on Computer Vision, October 2003\.
  * [14] W.M. Wells III.  Statistical approaches to feature-based object recognition.  International Journal of Computer Vision, 28(1/2):63–98, 1997.