Unified Strong-Field Dynamics Simulations from Atoms to Heterostructures

AI summary

Overview

Research area: Computational quantum physics — strong-field light–matter interaction, attosecond physics, and semiconductor heterostructure dynamics.

Technical level: Advanced. The paper assumes familiarity with the time-dependent Schrödinger equation (TDSE), B-spline and finite-element/isogeometric discretisations, sparse parallel linear algebra (PETSc/SLEPc), and strong-field concepts such as high-harmonic generation (HHG).

Scope (one sentence): The paper presents and validates TDSE-Z, an open-source parallel framework that solves the single-active-electron TDSE with the BenDaniel–Duke position-dependent-mass kinetic operator on adaptive B-spline meshes, spanning atomic and semiconductor heterostructure systems in one to three dimensions.

What This Paper Is About

Most strong-field simulation codes assume a constant electron mass, which is exact for dilute atomic gases but wrong for semiconductor heterostructures, quantum wells, and quantum dots where the effective mass varies in space and can jump sharply at material interfaces. The standard kinetic operator then fails to conserve probability current at those interfaces, producing spurious reflections that can masquerade as physical effects. The paper's goal is a single framework that rigorously handles the BenDaniel–Duke (BDD) operator for position-dependent mass while remaining efficient enough for demanding time-dependent strong-field calculations, and to demonstrate that it reproduces known atomic benchmarks and heterostructure physics.

Key Contributions

1. A weak-form Galerkin discretisation of the BenDaniel–Duke operator on geometry-adapted B-spline meshes. The framework supports arbitrary potentials and customisable laser configurations in one to three dimensions, and enforces self-adjointness and continuity of the wavefunction and probability flux at mass interfaces without extra interface boundary conditions.

2. An integrated, open-source simulation platform (TDSE-Z) built on PETSc, SLEPc, and PetIGA. It combines parallel element assembly, a shift-and-invert Krylov–Schur eigensolver with optional GPU offloading for initial-state preparation, Crank–Nicolson time propagation with a CPU-optimised MPI-parallel solver, and structured HDF5 output of spectra, populations, dipoles, energies, currents, and wavefunction snapshots.

3. A two-tier validation programme. A static position-dependent-mass eigensolver checked against the analytical Quesne model and a GaAs/Al₀.₃Ga₀.₇As double quantum well, plus a time-propagation engine checked against hydrogenic Rydberg series, the HHG cutoff law, harmonic-oscillator analytics, length–velocity equivalence, and angular-momentum commutator diagnostics.

4. Six non-uniform knot-grading strategies and demonstrated parallel scalability. The "knot engineering" schemes (uniform, symexp, symtan, symtann, symtanu, hydrogenic) target different physical regimes and are reported to reduce eigensolver and propagator cost by over an order of magnitude relative to uniform meshes of comparable resolution.

Main Findings

• Analytical PDM benchmark passed at machine precision: the weak-form BDD eigensolver reproduces the analytical Quesne position-dependent-mass model to machine precision.
• Heterostructure tunnel splittings follow WKB: for the GaAs/Al₀.₃Ga₀.₇As double quantum well, the computed tunnel splittings replicate the Wentzel–Kramers–Brillouin exponential decay over more than two orders of magnitude, with agreement at the sub-percent level.
• Atomic benchmarks reproduced: the time-propagation engine accurately reproduces HHG spectra in atomic benchmarks and confirms that dimensionality matters for fully capturing the strong light–matter interaction.
• Long-time norm and energy conservation: Crank–Nicolson propagation conserves the discrete M-norm to solver tolerance, and benchmarks against exact quantum harmonic oscillator solutions keep the energy error at machine precision.
• Strong scaling above 80%: the framework exceeds 80% strong-scaling efficiency up to 128 cores in the reported tests, and the abstract describes maintained performance across hundreds of CPU cores. The scaling study used a 3D hydrogen potential with 97 B-splines per axis, degree p = 3, a hydrogenic knot sequence over [−25, 25]³ a.u., and an AMD Zen 4 cluster with 32 MPI ranks per node; the ground-state energy on this grid was E₀ = −0.499995 a.u. without regularisation.
• Adaptive meshes reduce cost: non-uniform knot sequences place points near the Coulomb cusp and bound-state region while remaining sparse in the continuum, cutting the cost of the eigensolver and time propagator by over an order of magnitude compared with uniform meshes of comparable resolution.
• Gauge consistency verified: the code checks the length–velocity relation v⁽ᵅ⁾ₘₙ = i(εₙ − εₘ)d⁽ᵅ⁾ₘₙ for every non-degenerate state pair, and reports the Frobenius norm of the [H, L̂_z] commutator, which should vanish to machine precision in a spherically symmetric potential.
• Stable under extreme fields: for field strengths exceeding E₀ ≈ 0.5 a.u., time steps below Δt < 0.01 a.u. are adopted to resolve rapidly oscillating dipole dynamics.

Methodology in Plain English

The wavefunction is written as a sum of local piecewise-polynomial basis functions (B-splines) laid out on a mesh whose points are deliberately clustered where physics happens — near the atomic core, around interfaces, and in regions where the wavefunction oscillates quickly. Instead of evaluating the kinetic term directly, the equation is recast in "weak" integral form, which is the natural language of finite-element methods. This reformulation means a jump in the effective mass at a material boundary is handled automatically by the integrals, rather than requiring special interface conditions that, if omitted, cause artificial reflections.

With that discretisation, the time-dependent Schrödinger equation becomes a sparse matrix equation. The initial stationary state is found by solving a generalised eigenvalue problem with a Krylov–Schur method and shift-and-invert spectral transformation, using preconditioned iterative linear solvers and optionally a direct LU solver. The state is then marched forward in time with the Crank–Nicolson scheme, which is second-order accurate, unconditionally stable, and norm-conserving; each step solves a sparse linear system with GMRES preconditioned by additive Schwarz with overlap 1, ILU(0) on each subdomain, and tolerances of 10⁻¹².

During propagation the code computes the dipole moment directly from the coefficients and the dipole acceleration from Ehrenfest's theorem, including an extra term that accounts for the spatially varying mass and vanishes in the constant-mass limit. A smooth absorbing potential in the boundary layer prevents outgoing waves from reflecting back, and a Kaiser window suppresses edge artefacts before Fourier transforming to obtain HHG spectra. The dipole velocity is decomposed into intraband, interband, and bound–continuum contributions, with the interband part tied directly to quantum coherences through an optical Bloch equation.

Why This Matters

Impact on research: The work fills a specific software gap. Existing public codes — QPROP for B-spline atomic calculations, PCTDSE and 3D-GTDSE for Cartesian split-operator propagation, SCID-TDSE for spherical coordinates — all assume a constant mass in the kinetic term. TDSE-Z combines B-spline flexibility with the BDD formalism, giving the condensed-matter and attosecond communities a tool that treats mass discontinuities correctly rather than approximately. The validation against an analytical PDM model and against WKB tunnel splitting provides confidence that interface physics is not being contaminated by numerical artefacts.

Real-world applications:

• Semiconductor heterostructure and quantum-well HHG: interpreting and predicting harmonic emission from engineered semiconductor structures, where the effective-mass profile and sharp interfaces shape the ultrafast electron response.
• Attosecond optoelectronics: designing next-generation devices whose operational limits are set by ultrafast electron dynamics in nanostructures.
• Quantum dots and solid-state quantum information: modelling the strong-field response of platforms used for quantum information processing and quantum photonic devices.
• Moiré heterostructures and band-resolved dynamics: a stated application area, including electron–hole dynamics and intra-/inter-band current decomposition.
• Molecular attosecond spectroscopy: high-fidelity simulations to help decode complex electronic and nuclear dynamics observed in experiments.

Industry relevance: Semiconductor device and photonics designers benefit from a solver that respects interface boundary conditions; the quantum computing and quantum photonics sectors benefit from heterostructure-level dynamical modelling; and the HPC/scientific-software community gets a PETSc/SLEPc/PetIGA-based code with GPU offloading for the eigensolver and demonstrated scaling on hundreds of cores.

Future Directions

1. Extending GPU acceleration to time propagation. The static eigensolver already supports optional GPU offloading, but the propagation engine is described as CPU-optimised; the modular architecture is explicitly intended for future expansion toward exascale quantum dynamics.

2. Applying the framework to the full heterostructure programme. The paper lists semiconductor heterostructure harmonic response, electron–hole dynamics in Moiré heterostructures, and band-resolved intra-/inter-band dynamics as target applications that the BDD operator now makes accessible.

3. Benchmarking against solid-state HHG experiments. Recent experiments on semiconductor heterostructures and doped quantum wells are cited as motivation; quantitative comparison against such measurements is the natural next validation tier beyond the analytic and atomic benchmarks presented here.

4. Scaling studies beyond the reported configuration. The scalability demonstration used a 3D hydrogen potential on a specific grid and cluster topology; extending strong-scaling characterisation to production heterostructure runs, larger basis sets, and exascale-class machines remains open work.

Target Audience

Computational physicists working in atomic, molecular, and optical physics and in condensed-matter ultrafast dynamics; method developers interested in finite-element/isogeometric discretisations of the Schrödinger equation with position-dependent mass; researchers modelling semiconductor heterostructures, quantum wells, quantum dots, and Moiré systems under intense laser fields; and high-performance computing specialists who use or extend PETSc/SLEPc-based scientific software. Readers without a background in TDSE numerics or sparse parallel solvers will find the formalism and implementation sections demanding.

Authors’ abstract

We present \textsc{TDSE-Z}, a high-performance open-source framework for strong-field quantum dynamics in atomic, molecular, and semiconductor effective-mass systems. The core engine implements a weak-form Galerkin discretisation of the Hermitian BenDaniel-Duke operator, $\hat{T}_{\mathrm{BDD}} = -\frac{1}{2}\nabla\cdot(m^{-1}(\mathbf{r})\nabla)$, on geometry-adapted B-spline meshes, supporting arbitrary potentials and customisable laser configurations in one to three dimensions. We validate the static position-dependent-mass (PDM) eigensolver through two stringent benchmarks: a comparison to the analytical Quesne PDM model and a $\text{GaAs/Al}_{0.3}\text{Ga}_{0.7}\text{As}$ double quantum well, where the exponential decay of computed tunnel splittings follows Wentzel-Kramers-Brillouin (WKB) theory at the sub-percent level. We further demonstrate the time-propagation engine on constant-mass systems, accurately reproducing high-harmonic generation (HHG) spectra in atomic benchmarks and confirming the importance of dimensionality in fully capturing the strong light-matter interaction. Our implementation demonstrates robust strong-scaling efficiency, maintaining performance across hundreds of CPU cores. While the static eigensolver currently supports optional GPU offloading, the time-propagation engine is CPU-optimised, providing a modular architecture for future expansion toward exascale quantum dynamics.