how to simulate a steep potential barrier in langevin equation

Question

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?

Answer 1

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 !

Answer 2

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.