Classical Potential Barriers and their Resultant Forces

[ZombieKitty]

So, I'm writing a code to solve Hamilton's equations for a free classical particle incident on different potential barriers in 1D. I'm eventually going to be using this for some sort of semi-classical analysis, but for now have to get this working as expected.

I first tested it using a Gaussian barrier, which works.. I'm now trying to incorporate a rectangular barrier, but I know I'm getting something wrong with the physics.

The square barrier is just a piecewise defined as follows:
U = Uo for -a<=x<=a
= 0 elsewhere

The Hamiltonian is as simple as it gets:
H = p^2/(2m) + U

So, dH/dp = p/m = dx/dt and dH/dx = dU/dx = -dp/dt
These are the two equations that I propagate simultaneously using an RK4 algorithm.

My problem is how to handle the square potential numerically. Since dU/dx is infinite at the edges, there's basically an infinite force acting over an infinitesimally small time dt. Any hints as to how to code this into my program? My thoughts are that it has to do with the timestep... because any arbitrarily large slope that I code in is acting over one timestep, dt, which is too large.

I'm trying to think about the change in momentum analytically. I basically have the following:

Δp = ∫F dt = ∫-dU/dx dt

So, at the edges, dU/dx -> ∞ and dt -> 0

I know that analytically, if the incident particle does not have E > V, then it will reflect with equal and opposite momentum.

So Δp = -2p = ∫-dU/dx dt

I'm not quite sure how to resolve -2p in the limiting case analytically.. I think if I can figure this out, it will point me in the right direction as far as the coding goes. Any help would be greatly appreciated.

 

[ZombieKitty]

To clarify I'm looking for how to handle dt, in the limiting case:

lim dU/dx → ∞
Δp = -2p = ∫-dU/dx dt

 

[Khashishi]

You probably need to write a special case for the edge and handle that separately. This is a collision detection problem. There are numerous resources in computer game programming on this. You need to find the normal vector of the barrier and reflect the momentum over the plane perpendicular to the normal of the barrier.

You can't apply a force over a small time step because the force is infinite. You could treat it as an impulse, which is just an instantaneous change in momentum.

 

[ZombieKitty]

^^
Thanks for the reply.

I can see how I could do the former as you said, treating the edges as a special case. It's rather simple in 1D as well.. since the reflection of the momentum over the plane perpendicular to the normal = -pi (incident momentum)

I was, however, hoping to find a more elegant mathematical way of doing this. I was thinking of it in terms of an impulse, but I'm not sure how to handle it. I was thinking that this expression would lead me in the right direction:

I = Δp = ∫-dU/dx dt

If the incident momentum is p, then I know intuitively that I = -2p, since it will reflect with the same momentum in the opposite direction. I'm not sure how to resolve this analytically though.. Any ideas on this "proof":

lim dU/dx → ∞
I = ∫-dU/dx dt
************
Therefore, I=-2p

 

[Khashishi]

Consider if the barrier is not infinitely steep, but still very steep. Then the particle will penetrate the "edge" of the barrier (where it starts to increase), slow to a stop where T-V=0, and turn around and come right back out.
I = ∫-dU/dx dt {integral over t:-inf to inf}
Let t=0 be when T-V=0
I = ∫-dU/dx dt {integral over t:-inf to 0} + ∫-dU/dx dt {integral over t:0 to inf}
By symmetry these need to be equal
I = 2∫-dU/dx dt {integral over t:-inf to 0}
But p = 0 at t=0. So if p = p0 before it reaches the barrier, then I = -2*p0.

 

[ZombieKitty]

^^
Cute, thanks!