# Schrödinger equation with a potential step

I want to simulate a Gaussian wave packet colliding with a potential step. I tried this (it indicates what do I want to solve, however I know it is naive):

a = 200;
b = 10.;
sol = NDSolve[{I D[u[t, x], t] == (-1/2) D[u[t, x], {x, 2}] +
HeavisideTheta[x] u[t, x],
u[0., x] == Exp[-((x + 0.7*a)^2./(2*b^2))] Exp[I x/4.],
u[t, a] == 0, u[t, -a] == 0},
u, {t, 0, 4000}, {x, -a, a},
AccuracyGoal -> 4, PrecisionGoal -> 8];

Animate[
Plot[Evaluate[Abs[u[t, x] /. First[sol]]^2], {x, -a, a}, PlotRange -> {0, 1}],
{t, 0, 413, 0.01}]


This code without the Heaviside term gives a freely moving wave-packet to the right.

I don't know how to incorporate the boundary conditions in this problem, that must be continuity of the solution and of its first derivative at the potential step, and finiteness of the solution as $x \rightarrow \infty$.

The solution of this problem is known exactly from many quantum mechanics textbooks (see for example Cohen-Tannoudji, Diu, Laloe, Quantum Mechanics, vol I, complement JI, pg. 79).

Post edition: apparently -after studying the answers given in this post- the problem is with the Heaviside Theta function, which is similar but not totally equivalent to the UnitStep. In fact, the Heaviside is defined as $=0$ if $x<0$; or $=1$ if $x>0$, while the UnitStep is defined as $=0$ if $x<0$; or $=1$ if $x \ge 0$. Could this simple difference be the source of the error I get with the original code?. I think that was the source of the error messages I got with this code since the HeavisdeTheta remains unevaluated at $x=0$.

You can use LogisticSigmoid[c x] with a large enough $c \gg 1$ to approximate the square potential with a smooth step. For example

a = 200;
b = 10.;
sol = NDSolve[{I D[u[t, x], t] == (-1/2) D[u[t, x], {x, 2}] +
LogisticSigmoid[20 x] u[t, x],
u[0., x] == Exp[-((x + 0.7*a)^2./(2*b^2))] Exp[I x/4.],
u[t, 2 a] == 0, u[t, -2 a] == 0}, u, {t, 0, 4000}, {x, -2 a, 2 a},
AccuracyGoal -> 4, PrecisionGoal -> 8];

Animate[Plot[Evaluate[Abs[u[t, x] /. First[sol]]^2], {x, -2 a, 2 a},
PlotRange -> {0, 1}], {t, 0, 4000, 0.01}]


gives a decent result.

• thank you, I tried with (Tanh[x]+1) which is close to your solution, however the solution classically rebounds at the step, and there are not a transmitted wave to the right. Apr 24, 2018 at 2:19
• the solution must be two gaussians separating themselves from the discontinuity. Apr 24, 2018 at 2:24
• I think to get the behavior in the book your initial wave packet needs to be of the form Eq.(4) on page 80; in other words a superposition with a Gaussian envelope in momentum space, rather than an actual Gaussian in x-space.
– ulvi
Apr 24, 2018 at 7:10
• No, that formula is for total reflection $E<V$. The correct initial wave packet is given in eq 16-a, page 64, which is I used in the post. The answer given above is not the solution for the problem. I think appropriate boundary conditions at the jump must be given explicitly for NDSolve to be capable of handling the problem, or we must use another method. Apr 24, 2018 at 15:36
• Your wave packet doesn't have enough energy to penetrate the barrier. Try lowering the barrier height. a = 200; b = 10.; sol = NDSolve[{I D[u[t, x], t] == (-1/2) D[u[t, x], {x, 2}] + 0.025 LogisticSigmoid[20 x] u[t, x], u[0., x] == Exp[-((x + 0.7*a)^2./(2*b^2))] Exp[I x/4.0], u[t, 3 a] == 0, u[t, - 3 a] == 0}, u, {t, 0, 4000}, {x, - 3 a, 3 a}, AccuracyGoal -> 4, PrecisionGoal -> 8]; Animate[Plot[Evaluate[Abs[u[t, x] /. First[sol]]^2], {x, -3 a, 3 a}, PlotRange -> {0, 1}], {t, 0, 4000, 0.01}] gives an outcome close to half transmission / half reflection.
– ulvi
Apr 25, 2018 at 0:15

Just play with parameters and replace HeavisideTheta with UnitStep

a = 200;
b = 10.;
sol = NDSolve[{I D[u[t, x], t] == (-1/2) D[u[t, x], {x, 2}] +
0.15 UnitStep[x] u[t, x],
u[0., x] == Exp[-((x + 0.7*a)^2./(2*b^2))] Exp[I x/2.],
u[t, a] == 0, u[t, -a] == 0}, u, {t, 0, 4000}, {x, -a, a},
AccuracyGoal -> 4, PrecisionGoal -> 8]

Animate[Plot[Evaluate[Abs[u[t, x]]^2 /. First[sol]], {x, -a, a},
PlotRange -> {0, 1}], {t, 0, 750, 0.005}]


At t=400 we have both incident and reflected wave packets:

In your code I increased the initial momentum by a factor of 2 and reduced the height of the potential barrier by a factor of 0.15. The penetration probability strongly depends on these parameters. This is all!

Notice that a canonical way of solving this kind of problems is via the split-operator technique employing FFT to switch between the real (for potential energy) and the momentum spaces (for kinetic part of the Hamiltonian).

• Thank you, we already get the same solution with the help of ulvi. I didn't experience changing the methods yet. I suppose the splitting method must be the most accurate, I'll check it one of these days. Apr 25, 2018 at 14:53