A new method to find full complex roots of a complex dispersion equation for light propagation

A new numerical method is presented to find full complex roots of a complex dispersion equation. For the application of the solution, the complex dispersion equation of a cylindrical metallic nanowire is investigated. By using this method, locus of Brewster angle, complex dispersion curves of Surface Plasmon Polaritons (SPPs) and complex bulk modes can be obtained in once calculation. Approximate analytical solution to the complex dispersion equation has also been derived to verify our method.

💡 Research Summary

**
The paper introduces a novel numerical algorithm designed to locate all complex roots of a complex dispersion equation, a task that is central to the analysis of light propagation in plasmonic and other nanophotonic structures. Traditional approaches—exhaustive grid searches and Newton‑Raphson iteration—suffer from prohibitive computational cost or from missing roots when multiple solutions or strong losses are present. The authors overcome these limitations by combining a two‑stage local approximation (linear and quadratic) with rigorous criteria that determine whether a given mesh cell contains a root.

In the first stage, the complex equation f(x, y) = 0 (where x and y represent the real and imaginary parts of either frequency ω or wave‑vector k) is sampled at the four corners of a rectangular mesh cell A. By applying a first‑order Taylor expansion, the partial derivatives ∂f/∂x and ∂f/∂y are expressed as finite differences of the corner values. The algorithm then checks whether the derived increments Δa₁, Δa₂ (for x) and Δb₁, Δb₂ (for y) are positive and nearly equal within a user‑defined tolerance δ. If these conditions hold, cell A is flagged as containing a root. This linear test is computationally cheap and can be applied independently to every cell, guaranteeing that all possible roots are at least identified as candidate regions.

When the linear test fails because one of the partial derivatives is close to zero—i.e., the function exhibits a flat or saddle‑like behavior—the method proceeds to the second stage. A larger mesh cell B, comprising eight sampling points, is used to construct a second‑order Taylor expansion. From the eight values the second‑order partial derivatives fₓₓ, f_yy, and f_xy are obtained via finite‑difference formulas. Solving the resulting linear system yields refined estimates of the increments Δc and Δd (the corrections in x and y). Additional consistency checks (e.g., fₓₓ ≈ f₂ₓₓ, f_yy ≈ f₂_yy) are performed with the same tolerance δ. If these quadratic criteria are satisfied, cell B is accepted as a root‑containing region. This quadratic stage captures roots that would be missed by the linear approximation, such as those lying on inflection points or near branch points of the complex function.

The overall algorithm proceeds as follows: (1) define the search domain and discretize it with mesh sizes Sₓ and S_y; (2) apply the linear test to every cell, marking candidates; (3) for cells where the linear test is inconclusive, apply the quadratic test; (4) optionally refine the identified cells by sub‑gridding or by selecting the point with minimal |f| as the final root estimate. Because each cell is processed independently, the method is naturally parallelizable and scales well to high‑dimensional parameter sweeps.

To demonstrate the method, the authors apply it to the complex‑ω dispersion relation of surface plasmon polaritons (SPPs) on a cylindrical metallic nanowire. The dispersion equation (Eq. 9) is derived from Maxwell’s equations with appropriate boundary conditions, involving Bessel and Hankel functions of integer order n. The metal’s dielectric response follows a Drude model with a finite relaxation time τ, introducing loss and making the equation fully complex. The authors normalize frequencies by the bulk plasma frequency ω_p and lengths by c/ω_p, and they explore modes with azimuthal numbers n = 0, 1, 2, 3. Using mesh sizes Sₓ = S_y = 0.01 and a tolerance δ = 0.1, they obtain the full set of complex‑ω solutions in a single computation.

The results reveal two distinct branches for n = 0 and n = 1: (i) an asymptotic SPP branch approaching the surface‑plasmon frequency ω_sp ≈ 0.8 (in normalized units), and (ii) a higher‑frequency branch (ω > 1) identified as the Brewster‑angle locus, where the reflection coefficient vanishes for a specific incidence angle. Higher‑order modes (n = 2, 3) display only the SPP branch, consistent with stronger confinement of higher‑order azimuthal fields. Importantly, the complex‑ω curves show no back‑bending, confirming that they represent temporal decay (Im ω < 0) rather than spatial attenuation, which would appear in a complex‑k analysis. The authors also derive an approximate analytical solution to the dispersion equation and verify that it matches the numerical results within a few percent, providing an additional validation of the algorithm’s accuracy.

In summary, the paper delivers a robust, efficient, and easily parallelizable technique for locating all complex roots of transcendental dispersion equations. By avoiding the pitfalls of exhaustive grids and Newton‑Raphson iteration, it enables rapid, comprehensive mapping of plasmonic dispersion relations, including subtle features such as Brewster‑angle loci and multi‑mode interactions. The method is readily extensible to other systems where complex eigenvalues arise, such as metamaterials, photonic crystals with loss, and waveguides with gain or non‑linear response, making it a valuable tool for the broader nanophotonics and computational electromagnetics community.