Skip to content
Quantum.world

Research

QuantumBilliards.jl: An extensible, user friendly, high-performance Julia library for the large scale computation of spectra and eigenstates of quantum billiards and related problems

Overview Research area: Quantum chaos and computational physics — specifically the numerical solution of the two-dimensional interior Helmholtz eigenvalue problem with Dirichlet boundary conditions (q

QuantumBilliards.jl: An extensible, user friendly, high-performance Julia library for the large scale computation of spectra and eigenstates of quantum billiards and related problems
arXiv
2610.32605
Published
2026-10-23
Authors
Matic Orel, Črt Lozej

AI summary

Overview

Research area: Quantum chaos and computational physics — specifically the numerical solution of the two-dimensional interior Helmholtz eigenvalue problem with Dirichlet boundary conditions (quantum billiards), packaged as scientific software.

Technical level: Advanced. The paper assumes familiarity with the Helmholtz/Schrödinger eigenvalue problem, boundary-integral operators, Green's functions, and nonlinear matrix eigenvalue problems. The library itself is presented as user friendly, but the paper's content is aimed at researchers and numerical-methods practitioners.

Scope (one sentence): The paper documents QuantumBilliards.jl, an open-source Julia library that unifies geometry construction, symmetry reduction, several sweep and accelerated eigensolvers, eigenfunction reconstruction, and Poincaré–Husimi phase-space analysis for two-dimensional quantum billiards.

What This Paper Is About

Quantum billiards — a point particle moving freely in a two-dimensional domain and reflecting specularly off the walls — become, when quantized, the Dirichlet Helmholtz eigenproblem. Solving this problem accurately at high excitation (large wavenumber k) is numerically expensive, and different numerical formulations work better in different regimes. The paper presents a single Julia framework in which several independent quantization methods can be applied to the same geometry and their results cross-validated, covering the whole workflow from geometry definition to eigenstate and phase-space analysis.

Key Contributions

  1. A unified workflow framework. QuantumBilliards.jl combines basis-expansion methods (real plane waves, Fourier–Bessel, corner-adapted bases) and boundary-integral formulations (double-layer potential, combined field-integral equation) behind a common interface, so the same billiard geometry can be attacked by several independent methods.

  2. Four accelerated eigensolvers in one package. The library implements the Vergini–Saraceno scaling method, the expanded boundary integral method (EBIM), Beyn's contour-integral method, and a Compact Rational Krylov (CORK) method for the nonlinear boundary-integral eigenvalue problem, alongside conventional wavenumber-sweep solvers (DM, PSM, DLP, CFIE).

  3. Broad geometric coverage. The framework supports smooth, cornered, symmetry-reduced, and multiply connected domains, using Kress product quadrature with graded parametrizations near true geometric corners, and a composite solver that assembles self- and mutual interactions of separate boundary components into one globally coupled operator.

  4. Post-processing built in. The library reconstructs eigenfunctions from the boundary normal derivative (the boundary function) and computes Poincaré–Husimi (PH) phase-space representations, connecting spectra directly to eigenstate and phase-space analysis.

Main Findings

  • Independent solvers agree on the same spectrum. Figure 3 shows a direct spectral sweep for the stadium billiard in the (-1,-1) reflection-symmetry sector over k ∈ [20, 25]: the DM, PSM, DLP, and CFIE tensions have different absolute scales (they are defined by different scalar tension functions) but their minima occur at the same physical eigenwavenumbers.

  • Direct sweeps are asymptotically costly. With N = O(k) boundary degrees of freedom at fixed points per wavelength, a dense factorization costs O(N³) = O(k³), and a direct sweep incurs this repeatedly, including at wavenumbers that are not solutions. This is the stated motivation for the accelerated solvers.

  • Accelerated methods return several eigenvalues per local problem. Vergini–Saraceno converts the basis problem into a generalized eigenvalue problem around a reference wavenumber k₀ using a dilation-weighted boundary norm; EBIM expands the Fredholm matrix to second order, with first neglected term O(ε³) where ε = k − k₀.

  • Accuracy control is heuristic for Vergini–Saraceno. The paper states that no strict error bounds for the eigenvalues are proven for this method; in practice the error is approximated by the tension, written as Δk_n ≤ C τ_n(k).

  • Vergini–Saraceno has a geometric restriction. Formally, with a single expansion center, the method is only applicable to star-shaped billiards.

  • Basis choice is geometry dependent. Real plane waves are fast to evaluate but, due to finite-precision limitations, their convergence is in practice limited to mostly smooth and convex geometries, and convergence is not guaranteed even there. Fourier–Bessel and corner-adapted bases (with ν = π/α for interior corner angle α) incorporate boundary conditions directly; additional evanescent plane waves may be required near parts of the boundary that are difficult to resolve.

  • CFIE trades cost for robustness. The combined-field formulation with coupling parameter η = k (following Kress's near-optimal choice) avoids the spurious resonances of the pure double-layer formulation, at the cost of additional special-function evaluations.

  • Scarred eigenstates illustrated. Figure 1 shows scarred eigenstates of the stadium billiard with half width equal to 1, at wavenumbers k = 74.1838018, 74.22509066, 74.26339429.

  • Benchmarking is qualitative in the presented text. The available content compares solver formulations conceptually and via the sweep in Figure 3; specific runtime, memory, or accuracy benchmark tables are not reported in the supplied text.

Methodology in Plain English

The authors built a software library around two complementary ways of turning the Helmholtz equation into a finite numerical problem.

In the basis approach, the wavefunction is written as a sum of functions that already satisfy the Helmholtz equation in free space (plane waves, or Bessel functions adapted to corners). The only thing left to enforce is that the sum vanishes on the boundary. This is turned into a matrix problem whose smallest "tension" — a measure of how badly the boundary condition is violated — signals an eigenvalue. To avoid scanning wavenumber by wavenumber, the Vergini–Saraceno method Taylor-expands a specially weighted version of this tension around a reference wavenumber, so one matrix problem yields several nearby eigenvalues at once.

In the boundary-integral approach, Green's identities reduce the two-dimensional problem to an equation living only on the one-dimensional boundary. After discretization this is a matrix that depends nonlinearly on k, and eigenvalues occur where the matrix loses invertibility — monitored by its smallest singular value. The boundary integrals are evaluated with Kress product quadrature, which splits off the logarithmic singularity and integrates it with a rule derived from trigonometric interpolation, while the smooth remainder uses the periodic trapezoidal rule. Corners are handled by clustering quadrature nodes near them via a graded parametrization. EBIM applies the same Taylor-expansion idea to this boundary matrix; Beyn's method instead extracts all eigenvalues inside a contour in the complex wavenumber plane; CORK attacks the nonlinear eigenvalue problem with a compact rational Krylov approach.

Discrete symmetries are handled by restricting the calculation to a chosen symmetry sector (for example the (-1,-1) sector of the stadium), which the basis choice can encode directly through sines and cosines. Once eigenvalues are found, eigenfunctions are reconstructed from the boundary normal derivative using the Rellich identity, which converts the interior normalization integral into a boundary integral, and Poincaré–Husimi functions are computed for phase-space analysis.

Why This Matters

Quantum billiards are a cornerstone model in quantum chaos: because the classical dynamics is energy-independent, the semiclassical limit k → ∞ is unusually clean, making billiards an ideal testing ground for theoretical conjectures about spectral statistics, quantum ergodicity, scars, and eigenstate localization. Many contemporary studies need large ensembles of highly excited eigenstates with accurate boundary functions and phase-space representations, not just a handful of low-lying levels. Having several independent solvers in one framework matters because agreement between independent formulations is a practical validation tool, and because different formulations are advantageous in different parameter regimes.

Real-world applications (as discussed in the paper):

  • Microwave cavities: flat metallic and superconducting cavities are direct analogue billiards, since the Helmholtz and Schrödinger equations coincide in two dimensions, allowing spectra and field patterns to be measured directly.
  • Acoustic and elastic resonators: machined metal blocks provide a similar wave-based route to billiard spectra.
  • Electronic billiards: surface-state electrons confined by atomic corrals, or ballistic cavities in semiconductor heterostructures, with Coulomb blockade as a further probe.
  • Optical and cold-atom systems: shaped dielectric microcavities and microlasers act as open billiards, and light-sheet walls have been used to build billiards for cold atoms.

Industry relevance: the numerical machinery — high-frequency Helmholtz solvers, boundary-integral quadrature for cornered geometries, microcavity and microlaser modelling, and microwave cavity design — is relevant to photonics, resonator and antenna engineering, and semiconductor device modelling. The library's role as an extensible open-source ground layer for Julia scientific computing also has relevance for reproducible computational research and for groups building custom solvers on top of it.

Future Directions

  • Spectral-flow formulations. High-frequency spectral-flow methods based on Neumann-to-Dirichlet maps are named as an unexplored direction.
  • Different Green functions and physics. Extensions to billiards on the hyperbolic metric, billiards in a homogeneous magnetic field (discussed in an appendix), and relativistic neutrino billiards.
  • Better corner and singular-quadrature treatment. Recursively compressed inverse preconditioning for corners, charge bases added near corners, and alternative singular and near-singular quadratures such as QBX and QB2X.
  • Broader physical scope. The package is designed as a ground layer for related problems with different boundary conditions — for instance dielectric microcavities — as well as billiard state time evolution, which is described as planned for the future.

Target Audience

Researchers and graduate students in quantum chaos, semicondustor-style mesoscopic physics, wave physics, and computational acoustics/electromagnetics who need large batches of accurate Helmholtz eigenvalues and eigenstates in two-dimensional domains. It is also aimed at numerical analysts interested in nonlinear eigenvalue problems and boundary-integral methods, and at scientific-software developers who want an extensible Julia foundation for related boundary-value problems. Readers seeking a drop-in tutorial will find the paper's framing more methodological than introductory; readers seeking quantitative performance comparisons will need the benchmarks that are not present in the supplied text.

Authors’ abstract

We present QuantumBilliards.jl, an open-source Julia library providing a suite of numerical solvers for the two-dimensional interior Helmholtz eigenproblem with Dirichlet boundary conditions, a paradigmatic model in quantum chaos. The library implements conventional wavenumber-sweep techniques together with several accelerated eigensolvers: the Vergini--Saraceno scaling method, the expanded boundary integral method (EBIM), Beyn's contour-integral method, and the Compact Rational Krylov (CORK) method. It provides a common framework for constructing billiard geometries, exploiting discrete symmetries, computing eigenvalues and eigenstates, and treating smooth, cornered, and multiply connected domains. The library further includes reconstruction of eigenfunctions from boundary normal derivatives and the computation of Poincaré--Husimi (PH) phase-space representations. Available at: https://github.com/Quantum-Chaos-Julia/QuantumBilliards.jl.

Read the original paper