Unitarity

Unitarity guarantees us that the sum of probabilities of a measurement is constant, and that constant is 1. This provides a second verification of our FDTD code. Moreover, unitarity holds even when the wave function interacts with a potential, so it is not as rigorous of a check, but it is more broadly applicable.

Unitarity gives us $$\forall t\in ℝ,{\int }_{-\infty }^{\infty }{\psi }^{*}(x,t)\psi (x,t)dx=1$$ However, we don't have a continuous wave function defined over all space. We only have a numeric approximation to $\psi (x,t)$ on a finite grid. By design, the grid includes the area where $\psi$ has a significant value, so we make a fast and simple approximation of our integral with the trapezoidal rule.

$${\int }_{-\infty }^{\infty }{\psi }^{*}(x,t)\psi (x,t)dx\approx \sum _{i=0}^{N-2}({{\Psi }^{*}}_i{\Psi }_i+{{\Psi }^{*}}_{i+1}{\Psi }_{i+1})\mathrm{\Delta x}/2.0$$

where $\mathrm{\Delta x}$ is the spacing between grid points.

We expect the simulation to closely follow unitarity, so the integral will be close to 1. Rather than display the value of the integral, we focus on its variance from one by displaying $$1-{\int }_{-\infty }^{\infty }{\psi }^{*}\psi dx$$

We collect data for two cases, for a free particle with a constant $V=0$, and for a particle interacting with a potential barrier.

The difference between our time evolving wave function and the ideal unitary result.

Frame	Steps	$1-{\int }_{-\infty }^{\infty }{\psi }^{*}\psi dx$
		$V=0$	Barrier
0	0	2.966e-7	2.966e-7
200	100200	1.045e-7	1.855e-7
400	200400	3.937e-7	8.744e-7
600	300600	1.838e-6	2.210e-6
800	400800	2.471e-6	2.199e-6
1000	501000	3.859e-6	3.515e-6
1200	601200	5.211e-6	5.387e-6

The difference from unitarity slowly increases, but it does this is in line with the nature of our approximations. Most importantly, in both scenarios, the free particle, and a particle interacting with a potential barrier, there are no rapid changes in our integral.

Implementation

As you might expect, we borrow a lot from the previous verification step. One again, every 200^th frame we leverage the DumpSchrodinger class to fetch the wave function.

    const computed = await dumpSchrodinger.dumpWavefunction(schrodinger.getWaveFunctionBuffer0());
          

We also added method to the DumpSchrodinger class to compute the integral for us.

    /**
     * Simple trapezoidal integration of the wave function. The expectation is that the value will be roughly 1,
     * and roughly constant over time.
     *
     * @param {Float32Array} wavefunction An array containing the values of the wave function.
     * @param {Number}       length       The physical length of the simulation.
     * @returns {Number} The integral of the wave function over the length of the simulation.
     */
    integrateWaveFunction(wavefunction, length)
    {
      // Remember that each Ψ value contains both a real and imaginary part.
      const nIntervals = wavefunction.length/2.0 - 1.0;
      const dx = length / nIntervals;
      let sum = 0;
      for (let i=0; i<nIntervals; i++) {
        sum += ((wavefunction[2*i]*wavefunction[2*i] + wavefunction[2*i+1]*wavefunction[2*i+1]  // Ψᵢ
                  + wavefunction[2*i+2]*wavefunction[2*i+2] + wavefunction[2*i+3]*wavefunction[2*i+3]) // Ψᵢ₊₁
                /2.0) * dx;
      }
      return sum;
    }
          

Invoke this method with a wave function, such as the one we just retrieved, along with the physical length of the grid.

    const totalProbability = dumpSchrodinger.integrateWaveFunction(computed, length);
          

And then we log the results to the console where you can see the results for your run. We polish the output a little bit by using exponential notation for a number we expect to be small.

    console.log("<tr>");
    console.log("<td>" + nframes + "</td>");
    console.log("<td>" + nframes*nsteps+ "</td>")
    console.log("<td>" + (1.0-totalProbability).toExponential(3) + "</td>");
    console.log("</tr>");
          

Interestingly, if you take a closer look at the DumpSchrodinger class you might notice the simpson method. This is a more complex and more accurate integration method. In this case it produced results within a few percent of the simpler trapezoid rule, so I stuck with the simpler approach.

• Back