Quickstart: Your First SBE Simulation

What you will do

  1. Set up a 1D spatial grid along a quantum wire

  2. Define an 800 nm Gaussian laser pulse (10 fs pulsewidth)

  3. Initialize the full SBE solver (Coulomb, phonons, dephasing, transport, emission)

  4. Time-evolve for 3000 steps (30 fs) and watch the polarization response

  5. Plot the electric field and induced polarization

Prerequisites

This example requires the parameter files params/qw.params (material properties) and params/mb.params (many-body physics flags) to exist in the working directory. These ship with the examples and configure a GaAs quantum wire in an AlAs host with all major physics enabled.

See 02_architecture for a detailed description of every parameter.

1. Setup

import os
import numpy as np
import matplotlib.pyplot as plt
from scipy.constants import c as c0_SI

from pulsesuite.PSTD3D.SBEs import InitializeSBE, QWCalculator
from pulsesuite.PSTD3D.typespace import GetKArray, GetSpaceArray

# Ensure output directories exist (InitializeSBE writes diagnostic files here)
for d in ['dataQW/Wire/C', 'dataQW/Wire/D', 'dataQW/Wire/Ee',
          'dataQW/Wire/Eh', 'dataQW/Wire/P', 'dataQW/Wire/Win',
          'dataQW/Wire/Wout', 'dataQW/Wire/Xqw', 'dataQW/Wire/info',
          'dataQW/Wire/ne', 'dataQW/Wire/nh', 'output']:
    os.makedirs(d, exist_ok=True)

print("Setup complete.")

Setup complete.

2. Spatial Grid

The quantum wire lives on a 1D real-space grid. We use 50 pixels at 10 nm spacing (500 nm total length). GetSpaceArray and GetKArray create the matched real-space and momentum-space grids that the SBE solver needs.

Nr = 50             # Number of grid points along the wire
drr = 10e-9         # Pixel size: 10 nm

rr = GetSpaceArray(Nr, (Nr - 1) * drr)   # Real-space grid (m)
qrr = GetKArray(Nr, Nr * drr)            # Momentum-space grid (rad/m)

print(f"Grid: {Nr} points, {drr*1e9:.0f} nm spacing")
print(f"Wire length: {(Nr-1)*drr*1e9:.0f} nm")
print(f"Momentum range: [{qrr.min():.2e}, {qrr.max():.2e}] rad/m")

Grid: 50 points, 10 nm spacing
Wire length: 490 nm
Momentum range: [-3.14e+08, 3.02e+08] rad/m

3. Laser Pulse Parameters

We drive the quantum wire with a linearly-polarized Gaussian pulse at 800 nm (near the GaAs band gap at ~1.5 eV). The pulse has a 10 fs full-width and peaks at 15 fs into the simulation.

The pulse shape includes a super-Gaussian envelope (\(e^{-x^{20}}\)) that gives a sharp turn-on/turn-off, preventing long tails from wasting computation.

# Physical constants
c0 = c0_SI
twopi = 2.0 * np.pi

# Pulse parameters
E0x = 1e7           # Peak E-field: 10 MV/m (V/m)
lam = 800e-9        # Wavelength: 800 nm
tw = 10e-15         # Pulsewidth: 10 fs (FWHM of Gaussian envelope)
tp = 15e-15         # Pulse peak time: 15 fs

# Derived quantities
w0 = twopi * c0 / lam   # Angular frequency (rad/s)
Tc = lam / c0            # Optical cycle period (s)

# Time stepping
Nt = 3000           # Number of time steps
dt = 10e-18         # Time step: 10 attoseconds (s)
T_total = Nt * dt   # Total simulation time

print(f"Wavelength: {lam*1e9:.0f} nm")
print(f"Photon energy: {twopi * c0 / lam * 1.055e-34 / 1.6e-19:.3f} eV")
print(f"Optical cycle: {Tc*1e15:.2f} fs")
print(f"Pulse peak: {tp*1e15:.0f} fs, width: {tw*1e15:.0f} fs")
print(f"Simulation: {Nt} steps x {dt*1e18:.0f} as = {T_total*1e15:.0f} fs")

Wavelength: 800 nm
Photon energy: 1.553 eV
Optical cycle: 2.67 fs
Pulse peak: 15 fs, width: 10 fs
Simulation: 3000 steps x 10 as = 30 fs

4. Initialize the SBE Solver

InitializeSBE is the central entry point. It:

  • Reads material parameters from params/qw.params (band gap, masses, dephasing rates)

  • Reads physics flags from params/mb.params (Coulomb, phonons, transport, etc.)

  • Computes the momentum grid and energy bands

  • Initializes all subsystem modules (Coulomb matrices, phonon interactions, etc.)

  • Sets up initial Fermi-Dirac carrier distributions at 77 K

We configure 2 quantum wire subbands (Nw=2).

Emax0 = E0x  # Peak field magnitude for the wireoff guard

# Initialize: this reads params files and sets up everything
InitializeSBE(qrr, rr, 0.0, Emax0, lam, 2, True)

print("SBE solver initialized with all many-body physics.")

dcv / e0 = 3.3538811268940077e-10
alphae = 303090660.34386384, alphah = 768475083.4742767
1/qc = 2.3003106022655407e-09, sqrt(2) / sqrt(alphae² + alphah²) = 1.7119449377756965e-09
ehint = 0.8262109841666643
Wire Radius = 2.3003106022655407e-09
Wire sqrt(area) = 3.8950888292173374e-09
Wire Thickness = 5e-09

Calculating Coulomb Arrays
  Progress: 0/230 (0.0%)

Vint: Using JIT-compiled version
  Progress: 10/230 (4.3%)
  Progress: 20/230 (8.7%)
  Progress: 30/230 (13.0%)
  Progress: 40/230 (17.4%)
  Progress: 50/230 (21.7%)
  Progress: 60/230 (26.1%)
  Progress: 70/230 (30.4%)
  Progress: 80/230 (34.8%)
  Progress: 90/230 (39.1%)
  Progress: 100/230 (43.5%)
  Progress: 110/230 (47.8%)
  Progress: 120/230 (52.2%)
  Progress: 130/230 (56.5%)
  Progress: 140/230 (60.9%)
  Progress: 150/230 (65.2%)
  Progress: 160/230 (69.6%)
  Progress: 170/230 (73.9%)
  Progress: 180/230 (78.3%)
  Progress: 190/230 (82.6%)
  Progress: 200/230 (87.0%)
  Progress: 210/230 (91.3%)
  Progress: 220/230 (95.7%)
Finished Calculating Unscreened Coulomb Arrays

Quantum Wire Linear Chi = (-0.13219666733079763+0.14243061062572296j)

InitializeSBE?
SBE solver initialized with all many-body physics.

5. Time Evolution

Now we run the simulation. At each time step:

  1. Compute the laser E-field at the current time

  2. Call QWCalculator for each subband – this solves the SBEs for one time step

  3. Average the polarization from both subbands

  4. Record the field and polarization at the wire midpoint

# Allocate field and polarization arrays
Exx = np.zeros(Nr, dtype=np.complex128)   # E-field along wire (V/m)
Eyy = np.zeros(Nr, dtype=np.complex128)
Ezz = np.zeros(Nr, dtype=np.complex128)
Vrr_pot = np.zeros(Nr, dtype=np.complex128)  # Potential (V)

# Subband 1 polarizations
Pxx1 = np.zeros(Nr, dtype=np.complex128)
Pyy1 = np.zeros(Nr, dtype=np.complex128)
Pzz1 = np.zeros(Nr, dtype=np.complex128)

# Subband 2 polarizations
Pxx2 = np.zeros(Nr, dtype=np.complex128)
Pyy2 = np.zeros(Nr, dtype=np.complex128)
Pzz2 = np.zeros(Nr, dtype=np.complex128)

Rho = np.zeros(Nr, dtype=np.complex128)   # Charge density

# Boolean flags for QWCalculator (lists for in-place modification)
boolT = [True]
boolF = [False]

# Storage for time-series at wire midpoint
mid = Nr // 2
t_arr = np.zeros(Nt)
Ex_arr = np.zeros(Nt)
Px_arr = np.zeros(Nt)

# Time evolution loop
t = 0.0
for n in range(Nt):
    # Gaussian pulse with super-Gaussian cutoff
    phase = w0 * (t - tp)
    Exx[:] = (E0x
              * np.exp(-(phase**2) / (w0 * tw)**2)
              * np.cos(phase)
              * np.exp(-(phase**20) / (2 * tw * w0)**20))

    # Evolve SBEs for subband 1
    QWCalculator(Exx, Eyy, Ezz, Vrr_pot, rr, qrr, dt, 1,
                 Pxx1, Pyy1, Pzz1, Rho, boolT, boolF)
    # Evolve SBEs for subband 2
    QWCalculator(Exx, Eyy, Ezz, Vrr_pot, rr, qrr, dt, 2,
                 Pxx2, Pyy2, Pzz2, Rho, boolT, boolF)

    # Record midpoint values
    t_arr[n] = t
    Ex_arr[n] = np.real(Exx[mid])
    Px_arr[n] = np.real((Pxx1[mid] + Pxx2[mid]) / 2)

    t += dt

    # Progress indicator
    if (n + 1) % 1000 == 0:
        print(f"  Step {n+1}/{Nt} ({(n+1)*dt*1e15:.0f} fs)")

print(f"Done. Simulated {T_total*1e15:.0f} fs of quantum wire dynamics.")

/home/docs/checkouts/readthedocs.org/user_builds/pulsesuite/envs/latest/lib/python3.10/site-packages/numba/core/typed_passes.py:338: NumbaPerformanceWarning: 
The keyword argument 'parallel=True' was specified but no transformation for parallel execution was possible.

To find out why, try turning on parallel diagnostics, see https://numba.readthedocs.io/en/stable/user/parallel.html#diagnostics for help.

File "../../../../../envs/latest/lib/python3.10/site-packages/pulsesuite/PSTD3D/coulomb.py", line 957:
@jit(nopython=True, parallel=True, fastmath=True, cache=True)
def _MBCE_core(ne, nh, Veh2, Vee2, Win, Wout, k3, Ceh, Cee):
^

  warnings.warn(errors.NumbaPerformanceWarning(msg,

/home/docs/checkouts/readthedocs.org/user_builds/pulsesuite/envs/latest/lib/python3.10/site-packages/numba/core/typed_passes.py:338: NumbaPerformanceWarning: 
The keyword argument 'parallel=True' was specified but no transformation for parallel execution was possible.

To find out why, try turning on parallel diagnostics, see https://numba.readthedocs.io/en/stable/user/parallel.html#diagnostics for help.

File "../../../../../envs/latest/lib/python3.10/site-packages/pulsesuite/PSTD3D/coulomb.py", line 1016:
@jit(nopython=True, parallel=True, fastmath=True, cache=True)
def _MBCH_core(ne, nh, Veh2, Vhh2, Win, Wout, k3, Ceh, Chh):
^

  warnings.warn(errors.NumbaPerformanceWarning(msg,

/home/docs/checkouts/readthedocs.org/user_builds/pulsesuite/envs/latest/lib/python3.10/site-packages/pulsesuite/PSTD3D/SBEs.py:5769: ComplexWarning: Casting complex values to real discards the imaginary part
  self.I0[w - 1] = CalcI0(

     500 Elec       1.304783e-04       2.638160e-06       6.582184e-13      -5.192877e+02      -3.311274e-29
     500 Hole       1.304783e-04       2.622998e-06       2.084388e-08      -2.594662e+02      -1.063610e-28

    1000 Elec       1.458617e-03       3.463166e-05       2.831190e-10      -4.859135e+02      -3.098461e-29
    1000 Hole       1.458617e-03       3.421813e-05       1.950913e-07      -2.859514e+02      -1.172178e-28

  Step 1000/3000 (10 fs)

    1500 Elec       5.411518e-03       1.448956e-04       1.706842e-09      -3.421537e+03      -2.181766e-28
    1500 Hole       5.411518e-03       1.416956e-04       5.535595e-07      -3.508286e+01      -1.438125e-29

    2000 Elec       1.059943e-02       3.236310e-04       6.214206e-10      -7.229744e+03      -4.610096e-28
    2000 Hole       1.059943e-02       3.149357e-04       4.130049e-07       2.214505e+02       9.077750e-29

  Step 2000/3000 (20 fs)

---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
Cell In[5], line 40
     34 Exx[:] = (E0x
     35           * np.exp(-(phase**2) / (w0 * tw)**2)
     36           * np.cos(phase)
     37           * np.exp(-(phase**20) / (2 * tw * w0)**20))
     39 # Evolve SBEs for subband 1
---> 40 QWCalculator(Exx, Eyy, Ezz, Vrr_pot, rr, qrr, dt, 1,
     41              Pxx1, Pyy1, Pzz1, Rho, boolT, boolF)
     42 # Evolve SBEs for subband 2
     43 QWCalculator(Exx, Eyy, Ezz, Vrr_pot, rr, qrr, dt, 2,
     44              Pxx2, Pyy2, Pzz2, Rho, boolT, boolF)

File ~/checkouts/readthedocs.org/user_builds/pulsesuite/envs/latest/lib/python3.10/site-packages/pulsesuite/PSTD3D/SBEs.py:6021, in QWCalculator(Exx, Eyy, Ezz, Vrr, rr, q, dt, w, Pxx, Pyy, Pzz, Rho, DoQWP, DoQWDl)
   6017 def QWCalculator(
   6018     Exx, Eyy, Ezz, Vrr, rr, q, dt, w, Pxx, Pyy, Pzz, Rho, DoQWP, DoQWDl  # noqa: F811
   6019 ):
   6020     """Backward-compatible shim — delegates to _default_solver."""
-> 6021     _default_solver.QWCalculator(
   6022         Exx, Eyy, Ezz, Vrr, rr, q, dt, w, Pxx, Pyy, Pzz, Rho, DoQWP, DoQWDl
   6023     )

File ~/checkouts/readthedocs.org/user_builds/pulsesuite/envs/latest/lib/python3.10/site-packages/pulsesuite/PSTD3D/SBEs.py:5888, in SBESolver.QWCalculator(self, Exx, Eyy, Ezz, Vrr, rr, q, dt, w, Pxx, Pyy, Pzz, Rho, DoQWP, DoQWDl)
   5886     rh[:] = 0.0
   5887 else:
-> 5888     self._sbe_calculator(Ex, Ey, Ez, Vr, dt, Px, Py, Pz, re, rh, WriteFields, w)
   5890 RhoE = np.zeros(len(rr), dtype=complex)
   5891 RhoH = np.zeros(len(rr), dtype=complex)

File ~/checkouts/readthedocs.org/user_builds/pulsesuite/envs/latest/lib/python3.10/site-packages/pulsesuite/PSTD3D/SBEs.py:5681, in SBESolver._sbe_calculator(self, Ex, Ey, Ez, Vr, dt, Px, Py, Pz, Re, Rh, WriteFields, w)
   5679     ne3[k] = C3[k, k]
   5680     nh3[k] = D3[k, k]
-> 5681 self._relaxation(ne3, nh3, VC, E1D, Rsp, dt, w, WriteFields)
   5682 for k in range(self.Nk):
   5683     C3[k, k] = ne3[k]

File ~/checkouts/readthedocs.org/user_builds/pulsesuite/envs/latest/lib/python3.10/site-packages/pulsesuite/PSTD3D/SBEs.py:5279, in SBESolver._relaxation(self, ne, nh, VC, E1D, Rsp, dt, w, WriteFields)
   5277 """Apply phonon and carrier-carrier relaxation."""
   5278 if self.EHs or self.Phonon:
-> 5279     WinE, WoutE = self._relaxation_E(ne, nh, VC, E1D)
   5280     WinH, WoutH = self._relaxation_H(ne, nh, VC, E1D)
   5282     if WriteFields:

File ~/checkouts/readthedocs.org/user_builds/pulsesuite/envs/latest/lib/python3.10/site-packages/pulsesuite/PSTD3D/SBEs.py:5238, in SBESolver._relaxation_E(self, ne, nh, VC, E1D)
   5236 WoutE = np.zeros(Nk, dtype=float)
   5237 if self.EHs:
-> 5238     MBCE(
   5239         np.real(ne),
   5240         np.real(nh),
   5241         self.kr,
   5242         self.Ee,
   5243         self.Eh,
   5244         VC,
   5245         self.gam_eh,
   5246         self.gam_e,
   5247         WinE,
   5248         WoutE,
   5249     )
   5250 if self.Phonon:
   5251     MBPE(np.real(ne), VC, E1D, WinE, WoutE)

File ~/checkouts/readthedocs.org/user_builds/pulsesuite/envs/latest/lib/python3.10/site-packages/pulsesuite/PSTD3D/coulomb.py:1419, in MBCE(ne0, nh0, ky, Ee, Eh, VC, geh, ge, Win, Wout, k3, Ceh, Cee)
   1418 def MBCE(ne0, nh0, ky, Ee, Eh, VC, geh, ge, Win, Wout, k3=None, Ceh=None, Cee=None):
-> 1419     _require_instance().MBCE(ne0, nh0, ky, Ee, Eh, VC, geh, ge, Win, Wout)

File ~/checkouts/readthedocs.org/user_builds/pulsesuite/envs/latest/lib/python3.10/site-packages/pulsesuite/PSTD3D/coulomb.py:1311, in CoulombModule.MBCE(self, ne0, nh0, ky, Ee, Eh, VC, geh, ge, Win, Wout)
   1309 ne[1 : Nk + 1] = np.abs(ne0)
   1310 nh[1 : Nk + 1] = np.abs(nh0)
-> 1311 _MBCE_core(ne, nh, Veh2, Vee2, Win, Wout, self.k3, self.Ceh, self.Cee)

KeyboardInterrupt:

6. Results

The top panel shows the driving laser E-field – a few-cycle pulse at 800 nm. The bottom panel shows the polarization induced in the quantum wire. The polarization responds to the field but with a phase lag and amplitude set by the quantum mechanical response of the electron-hole system.

t_fs = t_arr * 1e15  # Convert to femtoseconds

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 6), sharex=True)

# E-field
ax1.plot(t_fs, Ex_arr * 1e-6, 'b-', linewidth=0.8)
ax1.set_ylabel('E-field (MV/m)')
ax1.set_title('Laser Pulse Driving a GaAs Quantum Wire')
ax1.axhline(0, color='gray', linewidth=0.5)

# Polarization
ax2.plot(t_fs, Px_arr, 'r-', linewidth=0.8)
ax2.set_ylabel('Polarization (C/m$^2$)')
ax2.set_xlabel('Time (fs)')
ax2.axhline(0, color='gray', linewidth=0.5)

plt.tight_layout()
plt.show()

print(f"Peak E-field: {np.max(np.abs(Ex_arr))*1e-6:.1f} MV/m")
print(f"Peak polarization: {np.max(np.abs(Px_arr)):.2e} C/m^2")

What just happened

Behind the scenes, InitializeSBE set up and QWCalculator solved:

  • Coulomb interactions: Screened electron-hole, electron-electron, and hole-hole scattering matrices (excitonic correlations, band-gap renormalization)

  • Phonon scattering: LO phonon emission and absorption at 36 meV (GaAs)

  • Dephasing: Diagonal and off-diagonal dephasing of the density matrix

  • DC transport: Drift of carriers under any applied DC field

  • Spontaneous emission: Radiative recombination of electron-hole pairs

  • Optical coupling: Dipole matrix elements connecting valence and conduction bands

All of these ran self-consistently at every time step, controlled by the flags in params/mb.params. This is the same physics that produces the published results in Gulley & Huang, Opt. Express 27, 17154 (2019).

Next steps

  • 02_architecture – understand the module structure and parameter files

  • 03_coulomb – inspect the Coulomb matrices computed during initialization

  • 04_phonons – explore phonon scattering rates and temperature dependence

  • 05_optics – see how fields are projected between propagation and QW spaces

References

  1. J. R. Gulley and D. Huang, “Self-consistent quantum-kinetic theory for interplay between pulsed-laser excitation and nonlinear carrier transport in a quantum-wire array,” Opt. Express 27, 17154-17185 (2019).