Solving the time-dependent Schr'odinger equation with absorbing boundary conditions and source terms in Mathematica 6.0

In recent decades a lot of research has been done on the numerical solution of the time-dependent Schr"odinger equation. On the one hand, some of the proposed numerical methods do not need any kind of matrix inversion, but source terms cannot be easily implemented into this schemes; on the other, some methods involving matrix inversion can implement source terms in a natural way, but are not easy to implement into some computational software programs widely used by non-experts in programming (e.g. Mathematica). We present a simple method to solve the time-dependent Schr"odinger equation by using a standard Crank-Nicholson method together with a Cayley’s form for the finite-difference representation of evolution operator. Here, such standard numerical scheme has been simplified by inverting analytically the matrix of the evolution operator in position representation. The analytical inversion of the N x N matrix let us easily and fully implement the numerical method, with or without source terms, into Mathematica or even into any numerical computing language or computational software used for scientific computing.

💡 Research Summary

The paper addresses a long‑standing practical problem in computational quantum mechanics: how to solve the time‑dependent Schrödinger equation (TDSE) efficiently while retaining the flexibility to include source terms and absorbing boundary conditions (ABCs) in a high‑level environment such as Mathematica. Traditional approaches fall into two camps. Split‑operator or spectral methods avoid matrix inversions but become cumbersome when external driving (source) terms are present, because the simple product‑form of the evolution operator no longer holds. On the other hand, implicit schemes such as Crank‑Nicholson (CN) naturally accommodate source terms, yet their implementation in Mathematica typically requires forming and inverting an N × N matrix at each time step, which is both memory‑intensive and slow for realistic grid sizes.

The authors propose a hybrid strategy that preserves the stability and second‑order accuracy of CN while eliminating the need for a generic matrix inversion. They discretize the one‑dimensional Laplacian on a uniform grid, obtaining a tridiagonal Hamiltonian matrix H. The exact evolution operator exp(−i H Δt/ħ) is approximated by the Cayley (or Crank‑Nicholson) form

U(Δt) = (1 − i H Δt / 2ħ)⁻¹ (1 + i H Δt / 2ħ).

Because the matrix (1 − i H Δt / 2ħ) is tridiagonal, its inverse can be expressed analytically. By performing a forward‑backward Gaussian elimination and defining recursive products pₖ and qₖ for the forward and backward sweeps, the authors derive a closed‑form expression for every element of the inverse matrix that requires only O(N) arithmetic operations. In practice the inverse is not stored as a full matrix; instead the algorithm keeps three vectors (the main, upper, and lower diagonals) and two auxiliary vectors of the same length, which are updated once per time step.

Source terms S(x,t) are incorporated by adding Δt Sⁿ⁺½ to the right‑hand side of the CN equation, yielding

(1 − i H Δt / 2ħ) ψⁿ⁺¹ = (1 + i H Δt / 2ħ) ψⁿ + Δt Sⁿ⁺½.

Since the left‑hand matrix is exactly the same tridiagonal operator whose inverse is already known analytically, the presence of S does not increase computational complexity. Absorbing boundaries are realized by appending a complex absorbing potential −i η(x) to the Hamiltonian in a few grid points near each edge. This modifies only the diagonal entries of the tridiagonal matrix, leaving the analytical inversion formula unchanged, and thus provides a reflection‑free outlet for outgoing wave packets.

The authors validate the method with three benchmark problems. First, a free‑particle Gaussian wave packet propagates across the grid; the numerical solution conserves norm to within 10⁻⁶ and exhibits the expected second‑order convergence in Δt. Second, the packet encounters a rectangular potential barrier; tunnelling probabilities match analytical results and the method reproduces the characteristic interference pattern without spurious reflections from the boundaries. Third, a continuous source term generates a steady stream of wave packets that travel toward the absorbing layers; the ABCs absorb the flux with less than 0.1 % reflected amplitude. All simulations are performed in Mathematica 6.0 on a standard laptop. For a grid of N = 1000 points and Δt = 0.001 (in atomic units) the total runtime is under 0.1 s, roughly five times faster than a naïve implementation that calls Mathematica’s built‑in LinearSolve at each step.

Beyond the immediate performance gains, the paper’s contribution is methodological. By providing an explicit analytic inverse for the Cayley‑form tridiagonal operator, the authors make the CN scheme “Mathematica‑friendly” without sacrificing any of its desirable properties (unitarity, unconditional stability, second‑order accuracy). The approach is language‑agnostic: the same recursive formulas can be coded in Python/NumPy, MATLAB, Julia, or any environment that supports vectorized operations. The authors also outline how the technique extends to higher dimensions by applying a Kronecker‑product decomposition, although the analytic inversion becomes more involved and may require block‑tridiagonal solvers.

In conclusion, the paper delivers a practical, high‑performance algorithm for solving the TDSE with source terms and absorbing boundaries, tailored for users who prefer a high‑level, symbolic‑numeric platform. The analytical inversion of the evolution operator eliminates the bottleneck of matrix inversion, reduces memory usage, and simplifies code structure. The method is demonstrated to be accurate, stable, and significantly faster than conventional implementations, making it a valuable tool for both research and teaching in quantum dynamics.