1
$\begingroup$

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$.

Improve this question
$\endgroup$

2 Answers 2

Reset to default
4
$\begingroup$

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.

Improve this answer
$\endgroup$
6
  • $\begingroup$ 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. $\endgroup$
    – Gluoncito
    Apr 24, 2018 at 2:19
  • $\begingroup$ the solution must be two gaussians separating themselves from the discontinuity. $\endgroup$
    – Gluoncito
    Apr 24, 2018 at 2:24
  • $\begingroup$ 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. $\endgroup$
    – ulvi
    Apr 24, 2018 at 7:10
  • $\begingroup$ 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. $\endgroup$
    – Gluoncito
    Apr 24, 2018 at 15:36
  • $\begingroup$ 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. $\endgroup$
    – ulvi
    Apr 25, 2018 at 0:15
1
$\begingroup$

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:

enter image description here

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).

Improve this answer
$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$
    – Gluoncito
    Apr 25, 2018 at 14:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.