previous up next ','..','$myPermit') ?> SYLLABUS  Previous: 2.5.1 Advection-diffusion equation  Up: 2.5 Implicit Crank-Nicholson  Next: 2.6 Computer quiz


2.5.2 Schrödinger equation

In quantum mechanics, the evolution of a particle is generally modelled with a complex wave function $ \vert\psi>=\psi(x,t)$ . If the particle does't decay with time, the total probability of finding this particle somewhere in space of course remains exactly unity at all times $ <\psi\vert\psi>=\int \vert\psi\vert^2 dx = 1 \;\;\forall t$ . Solving the Schrödinger equation normalized to have Planck's constant $ \hbar=1$ , the particle mass $ m=1/2$ and the static potential $ V(x)$

$\displaystyle i\frac{\partial\psi}{\partial t}=H(x)\psi$ $\displaystyle \hspace{5mm}\textrm{with}\hspace{5mm}$ $\displaystyle H(x)=-\frac{\partial^2}{\partial x^2}+V(x)$ (2.5.2#eq.13)

it is clearly an advantage to formulate a numerical scheme that conserves the probability by keeping the evolution operator unitary even as it is discretized. This can be achieved by writing a formal solution $ \psi(x,t)=\exp(-iHt)\psi(x,0)$ in the so-called Cayley form with
$\displaystyle \exp(-iH\Delta t)$ $\displaystyle \simeq$ $\displaystyle \frac{1-\frac{i}{2}H\Delta t}{1+\frac{i}{2}H\Delta t}$ (2.5.2#eq.14)

This approximation is unitary, accurate to second-order in time $ \mathcal{O}(\Delta t^2)$ . It suggests an evolution of the form
$\displaystyle (1+\frac{i}{2}H\Delta t)\psi^{t+\Delta t}$ $\displaystyle =$ $\displaystyle (1-\frac{i}{2}H\Delta t)\psi^t$ (2.5.2#eq.15)

Replacing the second order derivative in the Hamiltonian operator $ H(x)$ with finite differences centered in space $ \mathcal{O}(\Delta x^2)$ , one obtains a scheme that is stable, unitary and is in fact the same Crank-Nicholson method applied here for the Schrödinger equation:
$\displaystyle \psi_j^{t+\Delta t}$ $\displaystyle +$ $\displaystyle \frac{i\Delta t}{2}\left( -\frac{ \psi_{j-1}^{t+\Delta t} -2\psi_... ...Delta t} +\psi_{j+1}^{t+\Delta t}}{\Delta x^2} +V_j\psi_j^{t+\Delta t}\right) =$  
$\displaystyle = \psi_j^t$ $\displaystyle +$ $\displaystyle \frac{i\Delta t}{2}\left( -\frac{ \psi_{j-1}^t -2\psi_{j }^t +\psi_{j+1}^t}{\Delta x^2} +V_j\psi_j^t\right)$ (2.5.2#eq.16)

which is finally cast into the linear system

$\displaystyle \left(\begin{array}{c} -\frac{i\Delta t}{2\Delta x^2} \\ 1 +\frac... ...ft(\begin{array}{c} \psi_{j-1}^t \\ \psi_j^t \\ \psi_{j+1}^t \end{array}\right)$ (2.5.2#eq.17)

Exploiting the tri-diagonal structure of the matrix, the scheme has been implemented in JBONE using complex arithmetic 
      BandMatrixC a = new BandMatrixC(3, h.length);  //Complex objects
      BandMatrixC b = new BandMatrixC(3, h.length);
      Complex[]   c = new Complex[h.length];
      Complex     z = new Complex();
      double[] V = physData.getPotential();          //Heavyside(x-L/2) times
      double  scale  = 10.*velocity;                 //  an arb. scaling factor
      double  dtodx2 = timeStep/(dx[0]*dx[0]);
      Complex pih  = new Complex(0., 0.5*dtodx2);
      Complex mih  = new Complex(0.,-0.5*dtodx2);
      Complex pip1 = new Complex(1.,     dtodx2);
      Complex mip1 = new Complex(1.,    -dtodx2);

      for (int j=0; j<=n; j++) {
        z = new Complex(0.,0.5*scale*timeStep*V[j]);
        a.setL(j,mih);                               //Matrix elements
        a.setD(j,pip1.add(z));
        a.setR(j,mih);
        b.setL(j,pih);                               //Right hand side
        b.setD(j,mip1.sub(z));
        b.setR(j,pih);
      }
      c=b.dot(h);                                    //Right hand side with
      c[0]=c[0].add(b.getL(0).mul(h[n]));            // with periodicity
      c[n]=c[n].add(b.getR(n).mul(h[0]));            
      hp=a.solve3(c);                                //Solve linear problem

      for (int j=0; j<=n; j++){                      //Plot norm & real part
        fp[j]=hp[j].norm();
        gp[j]=hp[j].re();
      }
The code again relies on the complex BandMatrixC.solve3() method to solve the linear system efficiently with $ \mathcal{O}(N)$ complex operations using the standard LU-factorization that will be studied in 3.

The applet below illustrates the scheme with a calculation of the scattering of a wave packet on a square potential barrier.
JBONE applet:  press Start/Stop to compute the scattering of a wave packet on a square potential barrier with a height specified with PhysDataValue and a width extending throughout the right side of the periodic domain. The real part of the wave function (blue line) starts to oscillates and the evolution shows how the probability density (black line) splits into a reflected and a transmitted part.
"; if ($user_nbr<100) echo " "; echo "

"; ?> The perfect conservation of the lowest order moment shows that the total probability is indeed conserved for all times.



Numerical experiments: Crank-Nicholson for Schrödinger
  1. Modify the energy in the wavepacket ICWavelength and verify that the applet indeed reproduces the reflection / transmission coefficients that you expect.
  2. Increase the value of the TimeStep and explain what you observe.
  3. Extend the RunTime and observe how quasi-modes are formed, with reflections in the periodic domain that remind what happens to an electron in a crystal.

SYLLABUS  Previous: 2.5.1 Advection-diffusion equation  Up: 2.5 Implicit Crank-Nicholson  Next: 2.6 Computer quiz