Matrix-based computational photonics with applications in beam propagation and Mach–Zehnder interferometric design

Computational photonics demands efficient numerical techniques to handle high-dimensional field evolution with strong coupling. This work presents a hierarchical matrix-pairing framework for N-dimensional numerical integration, combined with a branching-matrix coarse-to-fine optimization strategy, applied to beam propagation method (BPM) simulations of tapered waveguides and Mach–Zehnder interferometer (MZI) design. In BPM simulations, replacing element-wise Python loops with vectorized NumPy array operations for the transverse-plane power and mode-overlap integrals yields a per-simulation analysis-stage speedup of approximately 91 × on a 128 × 128 grid with 50 propagation steps (Python 3.11/NumPy 1.26/Intel Core i7-12700H CPU). This speedup is specific to the pure-Python execution environment; JIT-compiled or compiled-language baselines would reduce the relative gain. Both implementations apply the same trapezoidal quadrature rule to identical discrete field arrays and agree to within 2.5 × 10⁻ 13 in the power and overlap integrals; these values represent mutual agreement between the two implementations and are not errors relative to an independent reference or ground truth. In MZI optimization, the branching-matrix coarse-to-fine search recovers the target differential arm length ΔL in under 60 iterations for the single test case demonstrated; a systematic comparison with gradient-based and derivative-free optimizers shows that its advantage lies in global-search reliability on a multimodal cost surface rather than in speed. These results demonstrate that vectorized integration provides a computationally efficient alternative to loop-based integration in Python-based BPM workflows, with performance gains that are environment- and grid-size-specific.