🔬 mesoscale physics

Non-unitary Time Evolution via the Chebyshev Expansion Method

This paper demonstrates that the Chebyshev expansion method can be effectively extended to non-Hermitian systems with unbounded spectra by identifying numerical rounding errors as the primary limitation and providing an analytic error bound to guide the selection of spectral radius and time steps for accurate simulations, as illustrated by the Hatano-Nelson model.

Original authors: Áron Holló, Dániel Varjas, Cosma Fulga, László Oroszlány, Viktor Könye

Published 2026-09-10

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Áron Holló, Dániel Varjas, Cosma Fulga, László Oroszlány, Viktor Könye

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

In the quantum world, the way a system changes over time is usually governed by a strict rule: the total amount of "stuff" in the system, known as probability, must always add up to one. This is called unitary evolution, and it applies to isolated systems where nothing enters or leaves. To predict how a quantum particle moves or how a wave spreads, scientists use mathematical tools to calculate this change step by step. For decades, one of the most reliable tools for this job has been a method based on a specific set of curves called Chebyshev polynomials. These curves are excellent at approximating functions, but they have a well-known limitation: they were thought to work perfectly only when the system's energy values stayed within a narrow, real-number range, like a ruler marked from negative one to positive one.

However, the real world is often messier than that. Many important physical situations, such as open systems where energy leaks out or materials that behave differently depending on the direction of flow, are described by non-Hermitian mathematics. In these cases, the energy values can be complex numbers, existing in a two-dimensional plane rather than just on a line. For a long time, physicists believed the powerful Chebyshev method would fail in these complex scenarios, forcing them to use slower or less accurate techniques. A team of researchers has now challenged this assumption, showing that the method works just as well in the complex realm as it does on the real line, provided the calculations are handled with a specific kind of care.

The researchers, working across institutions in Hungary and Germany, set out to test the limits of this mathematical tool. They focused on a specific model known as the Hatano-Nelson model, which describes a chain of sites where a particle can hop from one to the next, but with a twist: the particle finds it easier to hop in one direction than the other. This creates a non-reciprocal system, a common feature in non-Hermitian physics. By simulating the movement of a wave packet—a localized group of particles—through this chain, the team demonstrated that the Chebyshev expansion could accurately track the system's evolution over time, even when the underlying math involved complex numbers far outside the traditional safe zone.

The key discovery was not that the method works everywhere without cost, but rather that the errors previously blamed on the method's fundamental incompatibility with complex systems were actually just the result of standard computer rounding. When computers perform calculations, they store numbers with finite precision, much like a ruler that can only measure down to the nearest millimeter. When the researchers pushed the Chebyshev method into the complex plane, they found that the accumulated rounding errors grew larger, eventually making the results look wrong. They realized this was not a flaw in the mathematics itself, but a limitation of the digital tools used to solve it.

To fix this, the team developed a clear, analytical guide for how to use the method safely. They derived a rule that links the size of the time step in a simulation to the spread of the system's energy values. If the energy values are spread out over a large area in the complex plane, the simulation must take smaller, more frequent steps to keep the rounding errors from piling up. Conversely, if the energy values are tightly clustered, larger steps can be taken. This relationship allows scientists to predict exactly how precise their simulation will be before they even run it. By following this guideline, they showed that the Chebyshev method can produce results that are virtually indistinguishable from exact solutions, even for systems that were previously thought to be too difficult for this approach.

The study also highlighted a subtle advantage of this method over older techniques. In some complex systems, the mathematical description of the particle's state can become extremely large or extremely small, causing standard calculation methods to lose precision or crash. The Chebyshev expansion, when tuned correctly, avoids these pitfalls because it builds the solution from a series of stable, bounded steps. The researchers verified this by comparing their results against known analytical solutions for the Hatano-Nelson model under different boundary conditions. In every case, the method held up, confirming that the barrier between Hermitian and non-Hermitian physics is not a wall, but a threshold that can be crossed with the right numerical strategy.

Ultimately, this work expands the toolkit available to physicists studying open and dissipative systems. It suggests that the same efficient algorithms used for simple, closed quantum systems can be applied to the more complex, real-world scenarios where energy and information are constantly being exchanged. The findings do not claim to have solved every problem in non-Hermitian physics, but they remove a significant obstacle, proving that a trusted method can be extended far beyond its original boundaries. By identifying the true source of error and providing a practical way to control it, the researchers have opened the door to more accurate and efficient simulations of the quantum dynamics that govern everything from light in optical fibers to the behavior of electrons in exotic materials.

Technical Summary: Non-unitary Time Evolution via the Chebyshev Expansion Method

Problem Statement
The numerical computation of time evolution in quantum systems is a fundamental challenge, particularly for large systems where direct exponentiation of the Hamiltonian is computationally prohibitive. While the Chebyshev expansion method is a well-established, efficient technique for unitary time evolution in Hermitian systems with bounded spectra (typically mapped to the real interval $[-1, 1]$ ), its applicability to non-Hermitian systems has been considered limited. Physics literature often asserts that Chebyshev expansions fail outside the real $[-1, 1]$ interval due to the behavior of Chebyshev polynomials, leading researchers to adopt alternative methods such as Taylor expansions, Runge-Kutta algorithms, Faber polynomials, or Hamiltonian Hermitization for non-Hermitian dynamics. Conversely, mathematical literature suggests the expansion should be valid for arbitrary complex numbers. This paper addresses the apparent inconsistency between these perspectives and investigates the numerical limitations of applying Chebyshev expansions to non-Hermitian matrices.

Methodology
The authors generalize the Chebyshev expansion method to arbitrary non-Hermitian matrices by analyzing the expansion of the complex exponential function $e^{-itz}$ and extending it to matrix exponentials $e^{-itH}$ .

Theoretical Framework: The method relies on the expansion:
$e^{-itz} = J_0(t) + 2 \sum_{m=1}^{\infty} (-i)^m J_m(t) T_m(z)$
where $J_m$ are Bessel functions of the first kind and $T_m$ are Chebyshev polynomials of the first kind. The authors demonstrate that while the series converges for the entire complex plane, numerical stability is compromised by floating-point rounding errors when $z$ lies outside the standard interval.
Error Analysis: The paper derives an analytic upper bound for the numerical rounding error. It identifies that for $z$ outside $[-1, 1]$ , Chebyshev polynomials $T_m(z)$ grow exponentially (bounded by Bernstein ellipses with radius $\varrho$ ), while Bessel functions $J_m(t)$ decay. The interplay between these terms, combined with machine precision limits ( $\epsilon \approx 10^{-16}$ ), dictates the maximum achievable accuracy.
Algorithmic Implementation: The authors implement an adaptive time-stepping scheme. Instead of a fixed number of expansion terms, the series is truncated when terms fall below a threshold ( $10^{-14}$ ) for five consecutive iterations. The method is applied to the Hatano-Nelson (HN) model, a 1D tight-binding chain with non-reciprocal hopping, under both Periodic Boundary Conditions (PBC) and Open Boundary Conditions (OBC).
Validation: The numerical results obtained via Chebyshev expansion are compared against exact analytic solutions derived from the eigenvalues and eigenvectors of the HN Hamiltonian.

Key Contributions and Results

Extension to Non-Hermitian Systems: The paper establishes that the Chebyshev expansion method is valid for arbitrary non-Hermitian matrices, provided the time step and spectral radius are chosen carefully. The method is not fundamentally restricted to the $[-1, 1]$ interval.
Identification of Error Source: The authors demonstrate that the breakdown of accuracy outside the conventional spectral bounds is not due to a fundamental mathematical limitation of the expansion, but rather a consequence of numerical rounding errors in floating-point arithmetic. This is evidenced by comparisons using integer arithmetic, which maintain accuracy at much larger parameter values where floating-point arithmetic fails.
Analytic Error Bounds: The paper derives a rigorous analytic upper bound for the total rounding error (Eq. 17) and an inverse formula (Eq. 18) to determine the maximum allowable spectral radius $\varrho_{\text{max}}$ for a given time step $t$ and desired tolerance $\Delta_{\text{max}}$ .
$\varrho_{\text{max}} = \frac{2}{t} W\left(\frac{\Delta_{\text{max}}}{4\epsilon}\right)$
where $W$ is the Lambert W function. This provides a practical guideline for selecting simulation parameters.
Performance on the Hatano-Nelson Model:
- PBC: The spectrum lies on an ellipse in the complex plane. The method performs accurately as long as the time step is small enough to keep the effective spectral radius within the error bounds.
- OBC: The spectrum is entirely real. The method performs exceptionally well when the spectrum is within $[-1, 1]$ . For larger hopping parameters where the spectrum exceeds this interval, accuracy is maintained by reducing the time step.
- Comparison: The Chebyshev method yields results comparable in accuracy to Hermitian cases ( $p=0$ ) and outperforms exact diagonalization approaches in regimes where analytic solutions suffer from numerical instability due to exponentially localized eigenstates (specifically in the OBC case with large non-reciprocity).

Significance and Claims
The paper claims that the Chebyshev expansion method is a robust and general tool for computing non-unitary time evolution in non-Hermitian systems, correcting the misconception that it is restricted to real spectra or the $[-1, 1]$ interval. The primary significance lies in providing a quantitative, analytic framework to manage numerical errors. By identifying rounding errors as the limiting factor, the authors offer a clear strategy: by selecting appropriate time steps based on the derived error bounds, one can achieve high-precision simulations of non-Hermitian dynamics without resorting to more complex or less efficient alternative algorithms. The work bridges the gap between mathematical theory (validity on the complex plane) and physical numerical practice (handling floating-point limitations).

More like this