# how to simulate a steep potential barrier in langevin equation

When simulating a Langevin equation, how is a vertical potential barrier handled?

I have the time overdamped evolution of the position $x$, described by

$\gamma\frac{dx}{dt}=-V'(x)+\eta(t)$

where $V(x)$ is the potential, $\gamma$ the drag, and $\eta(t)$ the thermal noise. When I want to simulate a particular trajectory, I discretize the time in steps of a small finite $\Delta t$, and I use simply iterations with

$x_i = x_{i-1} - \frac{1}{\gamma}V'(x_{i-1})\Delta t +\sqrt{2D\Delta t} \, \eta_i$

where $D=K_BT/\gamma$ is the diffusion coefficient, and $\eta_i$ a random number.

Now, imagine $V(x)$ has an infinitely high barrier, whose slope can be tuned (and made vertical). I was wondering: if the particle by chance hits the barrier, there the value of the derivative $V'$ can be indefinitely large. In the simulation, this means that $x_i$ can be seriously "kicked" away by the term $-V'(x_{i-1})$. It seems unphysical.

I guess reducing the $\Delta t$ used helps, but is there a standard way to handle this? Or too much vertical barriers cannot be simulated?

• Are you using something like $V \propto tanh(x)$? If not, that might be the way to go because you can "tune" your gradient but still keep it from diverging. You could also use the wave solutions for wave bores to approximate a steep $V(x)$, but those have downstream oscillations that you may not want... – honeste_vivere Nov 18 '15 at 1:12
• @honeste_vivere : I'm not using $tanh(x)$ in particular, but yes, I have a "tunable" $V'(x)$ which I can prevent from diverging. Still, I wonder how to deal with the problem of a too steep slope. – scrx2 Nov 18 '15 at 11:59
• If you know the analytical/functional form for V(x) then you should be able to determine whether the gradient scale length is approaching the grid cell size of your numerical solving scheme, correct? If you do this in Mathematica, that software has several options in NDSolve that can deal with this issue and prevent diverging results (so long as you are a Mathematica-whisperer)... – honeste_vivere Nov 18 '15 at 14:11
• @fpdx "Now, imagine V(x) has an infinitely high barrier, whose slope can be tuned (and made vertical)." is not correct what you meant? I think you meant "infinitely steep barrier", otherwise Tom-Tom's answer should solve your problem. – chau Nov 19 '15 at 11:03

A infinite barrier of potential reflects the fact the particle cannot enter a certain region of space. Solving the Langevin equation with such a barrier means that you have to find a way to state that the particle cannot enter the domain, but you also have to describe what happens at the boundary, because several scenarios are possible :

1. the particle stops when it reaches the boundary (absorption)

2. the particle is reflected at the boudary, with opposite velocity

3. the particle stops during a waiting time $\tau$ (determined or random) before moving back

4. the particle is reflected with a speed change $\delta v$

5. the particle follows any scenario 1.-4. with certain probabilities

The behaviour at the boundary is extremely important.

If you run a numerical simulation

1. is easy to encode with a position test

2. can be easily mimicked in one dimension by extending the domain and folding it. For instance if you have a barrier at $x=0$, extend the domain to $\mathbb R$ and then add the probabilities of presence at $x$ and $-x$. The velocity distribution is obtained by substraction of course.

3. is obtained by adding the waiting time $\tau$ to the simulation time using the Scenario 2. trick.

4. also uses Scenario 2. trick, the velocity is modified when crossing the boundary as in Scenario 3.

5. is also quite easy to code in a numerical simulation using the above techniques.

I Hope is helps !

Here may not be a complete answer, but it may give you (and me) some hints and perhaps some alternative solution.

1. Euler's scheme should not work, even for the deterministic equation where noise is set to zero. This is because in Euler's scheme, one always requires $\Delta t$ small so that the $\Delta x$ is small. When $\Delta x$ is large, one runs into the numerical instability of the solution. This problem is well-known for ODEs (see, e.g., Numerical Recipes, section on "Stiff systems".)

2. I think one way to avoid this is to use discretised-space simulation. Here you discritise the space into discrete steps and simulate the discrete Markov process with transition rates between neighbouring sites proportional to the difference in potential between the sites.

3. Another option is to use under-damped simulation. Note that the problem arises because your over-damped equation is first order. [It is not clear to me that over-dampped approximation is valid for such a infinitely strong force; although the over-damped equation seems to be mathematically well-defined.] For the under-damped equation, you can simulate $\Delta x$ smoothly, while $\Delta v$ jumps whenever the particle cross the step of the potential.

• Thanks. About your point 1: in this answer, Euler's scheme is said to be the best. – scrx2 Nov 19 '15 at 8:35
• About your point 3: the simulation is about the movement of a micron scale particle in water (with a tunable external potential). Overdamping comes from the Reynold's number of ~1e-5. – scrx2 Nov 19 '15 at 8:39
• About your point 2: interesting. Do you have any ref I can check? – scrx2 Nov 19 '15 at 8:40
• @fpdx 1. Euler's is better than Runge-Kutta's and the like (not the best). It works well when $\Delta t$ is small. Suppose you solve $x'=f(x)$. You then update $x_{n+1}=x_{n}+f(x_n) \Delta t$ assuming that $f(x_n) \Delta t$ is small. However in your case $f(x) \Delta t$ is never small, it is in fact infinite. More worse when $f(x)$ is concentrated in a small region (delta function): in your simulation, you will never have $x_n$ exactly at the delta function. Therefore effectively your simulation does not see any step potential. – chau Nov 19 '15 at 11:10
• 2. No, I don't. It is rather based on physical intuition; that is why I say it was not a complete answer. For sure it gives the correct dynamics when the particle is away from the step (because it is simply the Brownian motion) as well as the stationary distribution (because it satisfies the same detailed balance). For the dynamics near the barrier, I don't know if it is 100% correct or one has to choose the transition rates more wisely than Glauber's. – chau Nov 19 '15 at 13:07