Implementing the Dirac Equation in Processing: Visualizing Zitterbewegung and Linear Dispersion

In the previous post, as an introduction, I made a declaration:

I will try to implement the Dirac equation.

This time is the first step.

I will numerically implement the 1D Dirac equation and visualize the time evolution of a wave packet.

Don't just read it.
Run it.

---

1. 1D Dirac Equation (Minimal Form)

The 1+1 dimensional Dirac equation can be written in the following form:

$$i \partial_t \psi = (-i \alpha \partial_x + \beta m)\psi$$

Where:

• (\psi) is a 2-component spinor
• (m) is the mass
• (\alpha, \beta) are Pauli matrices

Specifically, we take:

$$\alpha = \sigma_x,\quad \beta = \sigma_z$$

Then the Hamiltonian becomes:

$$H = -i\sigma_x \partial_x + m\sigma_z$$

The important points are:

It has 2 components
The spatial derivative is linear

This is the decisive difference from the Schrödinger equation.

---

2. Discretization

Divide the space into a grid:

• Number of grid points: N
• Spatial step: dx
• Time step: dt

The spatial derivative is approximated by the central difference:

$$\partial_x \psi(x_i) \approx \frac{\psi_{i+1} - \psi_{i-1}}{2dx}$$

For time evolution, we use the Euler method (simply, for now):

$$\psi(t+dt) = \psi(t) - i H \psi(t) dt$$

*Prioritizing "seeing the movement" over strict rigor.

---

3. Processing Implementation Code

Below is a minimal configuration code example.

int N = 400;
float dx = 0.05;
float dt = 0.01;
float m = 1.0;

float[][] psiRe = new float[N][2];
float[][] psiIm = new float[N][2];

void setup() {
  size(800, 400);
  initializeWavePacket();
}

void draw() {
  background(0);
  updateDirac();
  drawWave();
}

void initializeWavePacket() {
  for (int i = 0; i < N; i++) {
    float x = (i - N/2) * dx;
    float gauss = exp(-x*x * 2.0);
    psiRe[i][0] = gauss;
    psiRe[i][1] = 0;
    psiIm[i][0] = 0;
    psiIm[i][1] = 0;
  }
}

void updateDirac() {
  float[][] newRe = new float[N][2];
  float[][] newIm = new float[N][2];

  for (int i = 1; i < N-1; i++) {

    for (int s = 0; s < 2; s++) {
      newRe[i][s] = psiRe[i][s];
      newIm[i][s] = psiIm[i][s];
    }

    // Spatial derivative term
    float dRe0 = (psiRe[i+1][1] - psiRe[i-1][1]) / (2*dx);
    float dIm0 = (psiIm[i+1][1] - psiIm[i-1][1]) / (2*dx);

    float dRe1 = (psiRe[i+1][0] - psiRe[i-1][0]) / (2*dx);
    float dIm1 = (psiIm[i+1][0] - psiIm[i-1][0]) / (2*dx);

    // σx ∂x term
    newRe[i][0] += dt * dIm0;
    newIm[i][0] -= dt * dRe0;

    newRe[i][1] += dt * dIm1;
    newIm[i][1] -= dt * dRe1;

    // m σz term
    newRe[i][0] -= dt * m * psiIm[i][0];
    newIm[i][0] += dt * m * psiRe[i][0];

    newRe[i][1] += dt * m * psiIm[i][1];
    newIm[i][1] -= dt * m * psiRe[i][1];
  }

  psiRe = newRe;
  psiIm = newIm;
}

void drawWave() {
  stroke(0, 255, 255);
  for (int i = 0; i < N-1; i++) {
    float prob = psiRe[i][0]*psiRe[i][0] + psiIm[i][0]*psiIm[i][0]
               + psiRe[i][1]*psiRe[i][1] + psiIm[i][1]*psiIm[i][1];
    line(i*2, height/2 - prob*200, (i+1)*2, height/2 - prob*200);
  }
}

---

4. What do you see?

① Linear Dispersion

In the Schrödinger equation,

$$E \propto k^2$$

But in Dirac,

$$E \propto k$$

The way the wave packet spreads is different.

This is part of the reason why,

graphene is said to be "Dirac-like."

---

② Zitterbewegung (Trembling motion)

When it has mass,

the center of the wave packet vibrates slightly.

This is a phenomenon caused by the interference between:

• Positive energy components
• Negative energy components

Concepts that were once abstract

now tremble before your eyes.

---

5. What was seen here

The Dirac equation is

not some obscure formula for a 4-component spinor.

• 2 components
• Linear derivative
• Mass term

With just these,

• Relativistic dispersion
• Positive-negative energy symmetry
• Spinor structure

all emerge.

This is the structure.

---

🔍 What is Happening

Looking at the screenshot:

• Initial Gaussian wave packet in the center
• It collapses into "grains" over time
• The amplitude seems to scatter away

This is typically due to:
① Numerical instability (limitations of the Euler method)

The current code is simple forward Euler.
The essence of the Dirac equation is unitary time evolution.
However, the Euler method:

• Does not preserve the norm
• Accumulates phase errors
• Causes high-frequency components to go wild

That is why it "granulates" and disappears.
This isn't wrong, per se; it's just numerically crude.

🧠 This part is actually important

The fact that it becomes grainy means:

• Positive and negative energy components are mixing
• High wave-number modes are being excited
• Lattice effects (Fermion doubling problem) might be starting to appear

In other words, the Dirac-like structure is emerging.

🔬 Is it correct?

YES. However:

• It's not perfect physical evolution
• It's numerically coarse
• The norm is not preserved

This is the current state. But as a first stage, it is correct.

🔧 What happens if you improve it?

If you switch to:

• Crank–Nicolson method
• Split-operator method
• Unitary updates using identity matrix approximation

Then you will see:

• Norm preservation
• Clean Zitterbewegung
• Clear linear dispersion

The current stage is "it moved." Next, we can proceed to "physically correct time evolution."

6 Since it is currently using the forward Euler method, norm preservation is not guaranteed.

🌟 The Most Interesting Point

This "granulation" is also a visualization of:

• High wave-number modes
• Lattice artifacts
• Numerical instability

In other words, it is evidence that the Dirac equation has a numerically delicate structure.
This serves as a perfect foreshadowing for the FPGA chapter.

🔁 Supplementary Experiment: "Seeing" Mass and Norm Preservation

In the previous section, we implemented the 1D Dirac equation in Processing and observed the time evolution of a wave packet.

Here, we take it one step further.

We will observe:

• Preservation of the norm (probability)
• Time step dependency
• Differences in Zitterbewegung due to the mass term

---

🧠 Observation Points

The Dirac equation inherently possesses unitary time evolution.

In other words,

$$\int |\psi|^2 dx$$

remains constant regardless of time.

However, in a simple forward Euler method, it is not perfectly preserved.

Let's "see" this difference.

---

🧪 Improved Processing Code

Below is the complete code for the supplementary experiment.

You can run it by replacing the previous code entirely.

// ===============================
// 1D Dirac Equation (1+1D)
// Norm visualization & mass test
// ===============================

int N = 400;
float dx = 0.05;

// ▼▼▼ Test Parameters ▼▼▼
// Reducing dt improves numerical stability
// First run at 0.01 → then change to 0.002 for comparison
float dt = 0.01;

// m = 1.0 → Zitterbewegung (trembling motion) occurs
// m = 0.0 → Zitterbewegung disappears (confirm linear dispersion)
float m = 1.0;
// ▲▲▲ Test Parameters ▲▲▲

float[][] psiRe = new float[N][2];
float[][] psiIm = new float[N][2];

void setup() {
  size(900, 500);
  initializeWavePacket();
}

void draw() {
  background(0);
  updateDirac();
  drawWave();
  displayNorm();
}

void initializeWavePacket() {
  for (int i = 0; i < N; i++) {
    float x = (i - N/2) * dx;
    float gauss = exp(-x*x * 2.0);
    psiRe[i][0] = gauss;
    psiRe[i][1] = 0;
    psiIm[i][0] = 0;
    psiIm[i][1] = 0;
  }
}

void updateDirac() {

  float[][] newRe = new float[N][2];
  float[][] newIm = new float[N][2];

  for (int i = 1; i < N-1; i++) {

    // Copy current values
    for (int s = 0; s < 2; s++) {
      newRe[i][s] = psiRe[i][s];
      newIm[i][s] = psiIm[i][s];
    }

    // Central difference (spatial derivative)
    float dRe0 = (psiRe[i+1][1] - psiRe[i-1][1]) / (2*dx);
    float dIm0 = (psiIm[i+1][1] - psiIm[i-1][1]) / (2*dx);

    float dRe1 = (psiRe[i+1][0] - psiRe[i-1][0]) / (2*dx);
    float dIm1 = (psiIm[i+1][0] - psiIm[i-1][0]) / (2*dx);

    // σx ∂x term
    newRe[i][0] += dt * dIm0;
    newIm[i][0] -= dt * dRe0;

    newRe[i][1] += dt * dIm1;
    newIm[i][1] -= dt * dRe1;

    // m σz term (mass term)
    newRe[i][0] -= dt * m * psiIm[i][0];
    newIm[i][0] += dt * m * psiRe[i][0];

    newRe[i][1] += dt * m * psiIm[i][1];
    newIm[i][1] -= dt * m * psiRe[i][1];
  }

  psiRe = newRe;
  psiIm = newIm;
}

void drawWave() {

  stroke(0, 255, 255);
  noFill();

  beginShape();
  for (int i = 0; i < N; i++) {

    float prob =
      psiRe[i][0]*psiRe[i][0] + psiIm[i][0]*psiIm[i][0] +
      psiRe[i][1]*psiRe[i][1] + psiIm[i][1]*psiIm[i][1];

    vertex(i * width / N, height/2 - prob * 250);
  }
  endShape();
}

float totalNorm() {
  float sum = 0;
  for (int i = 0; i < N; i++) {
    sum += psiRe[i][0]*psiRe[i][0] + psiIm[i][0]*psiIm[i][0]
         + psiRe[i][1]*psiRe[i][1] + psiIm[i][1]*psiIm[i][1];
  }
  return sum * dx;
}

void displayNorm() {
  fill(255);
  textSize(16);
  text("Norm = " + nf(totalNorm(), 1, 4), 20, 30);
  text("m = " + m + "   dt = " + dt, 20, 50);
}

---

🔬 Experiment Guide

① Compare by changing dt

• dt = 0.01
• dt = 0.002

→ You will see how the norm collapses differently.
→ It becomes clear that the Dirac equation is numerically delicate.

---

② Compare by changing m

• m = 1.0
• m = 0.0

→ Trembling motion (Zitterbewegung) appears when m = 1.0.
→ When m = 0.0, the trembling disappears and the wave packet spreads linearly.

Here, you can visually confirm that:

0 Zitterbewegung arises from the interference between the mass term and the positive/negative energy components.

---

🧩 Emerging Structure

• Linear dispersion (E ∝ k)
• 2-component spinor structure
• Interference effects due to mass
• Numerical instability

All of these are written in the equations from the beginning.

However,

only by running it does it become a "realization."

---

• Trembling state with m=1.0
  
• Straight-moving state with m=0.0, dt=0.001
  

Next Time

We will translate this structure into the physical constraints of an FPGA.

Can an abstraction
be turned into a logic circuit?

Next is "Physicalization."

---