Initial Conditions

Our simulation evolves the Schrödinger equation in time, but without initial conditions it has nothing to evolve, and is a very boring simulation. Let's spice things up a bit and fire in a free particle from the left. Later, we will put the variations in the potential, $V\left(x\right)$ , to the right.

The free particle Schrödinger equation $$i\hslash \frac{\partial }{\partial t}\Psi \left(x,t\right)=\frac{-{\hslash }^2}{2m}\frac{{\partial }^2}{\partial x^2}\Psi \left(x,t\right)$$ has solutions representing a localized free particle as a Gaussian wave packet. A Gaussian wave packet centered around $x=x_0$ with momentum $p_0$ is¹ $$\Psi \left(x,0\right)=\frac{1}{{\left(2\pi {\mathrm{\sigma }}^2\right)}^{\frac{1}{4}}}{e}^{-\frac{{(x-x_0)}^2}{4{\mathrm{\sigma }}^2}}{e}^{i{p}_0\mathrm{x}}$$

The most interesting part is the momentum exponential, where $p=\hslash k$ , but in our code $\hslash =1$ , so we make the expansion $$e^{i{p}_0\mathrm{x}}=\cos kx+i\sin kx$$ This is the only term with an imaginary part, so it determines the division between real and imaginary parts of the wave function, and hence the division between the r and g components of the texture.

There are no derivatives in this expression, so we need only evaluate if for each position on our grid. Luckily that is exactly what we have been doing. Once again we use the standard geometry and texture coordinates which will evaluate the function on exactly the grid points we are using in our calculations.

  // The physical length of the grid in nm.
  uniform float length;
  // center of the wave packet
  uniform float x0;
  // width of the wave packet
  uniform float w;
  // wave number
  uniform float k;
  
  varying vec2 vTextureCoord;

  vec2 computePsi(float s)
  {
    float x = length*s;
    // Generalization of Gaussian width
    float alpha = w*w;
    float deltaX = x-x0;
    // Normalization constant http://quantummechanics.ucsd.edu/ph130a/130_notes/node80.html
    float a = pow(2.0/(PI*alpha), 0.25);
    float gaussian = a*exp(-deltaX*deltaX/alpha);
    vec2  phase = vec2(cos(k*x), sin(k*x));
    return gaussian*phase;
  }

  void main()
  {
    gl_FragColor.rg = computePsi(vTextureCoord.s);
  }
      

To setup this Gaussian wave packet as the initial conditions for our simulation we need to write it into the texture used for input for the first step of our simulation. Luckily, we built a method into the Schrödinger solver that returns a framebuffer for rendering into exactly this texture.

There are some interesting things we can observe in even this initial simulation. The wave packet is not a sharp peak, but is spread out in space. Further, the wave packet spreads out as it propagates so it is wider and shorter at the end of the simulation. That is uncertainty of a particle naturally increases over time even with no interactions.

• Back
• Next