Introduction

Originally discovered in the context of electronic systems, the concept of topological insulators has been generalized to photonic structures thanks to the platform-independent framework of topological band theory used to classify such systems. Subsequently, over the past decade, there has been substantial interest in photonic topological insulators due to their potential to yield next-generation optical devices based on their topologically protected edge states1. For example, non-reciprocal waveguiding modes can be achieved in photonic crystals exhibiting non-trivial topology from broken time-reversal symmetry, which can be realized using gyro-electric or gyro-magnetic materials2,3 as well as in driven nonlinear systems4. Similar waveguiding modes can also be found in a variety of metamaterials, such as those based on shifted ring-resonator arrays5,6, helical waveguide arrays7, or that use synthetic dimensions8,9,10,11. Moreover, using solely the crystalline symmetries of the photonic crystal, different classes of non-trivial topology can also be attained, leading to robust waveguiding states along bends that preserve the crystalline symmetry12,13,14,15,16 or robust cavity-like states for enhanced light-matter interactions17,18,19.

Material topology in electronic systems is identified through the system’s band structure and Bloch eigenstates using invariants defined on the system’s Brillouin zone; similarly, topological robustness is defined in terms of the system’s bulk band gap20,21. Traditionally, these same invariants have been used to classify topology in photonic systems as well. However, the analogy between topological insulators in electronic systems and in photonic systems is not always exact. For realizations of photonic topological insulators operating at longer wavelengths2,22,23, the system can be bounded by materials that provide excellent approximations of perfect electric conductors24, yielding similar open boundary conditions to those that appear in electronic systems. In contrast, photonic topological insulators operating at technologically relevant wavelengths and length-scales are generally based on photonic crystal slab motifs, and border-free space on at least one surface, which is gapless above the light cone [Fig. 1]. This gapless radiative environment has two relevant consequences for the classification of photonic topology: (1) For those wavevectors above the light-cone, the photonic crystal slab does not possess a true band structure consisting entirely of bound state solutions with real frequencies and instead only exhibits a resonance band structure characterized by complex frequencies. As such, it is not known whether standard topological invariants, defined by integrating over a system’s occupied bands, can be meaningfully applied to the resonances of such photonic slab systems. (2) Similarly, the full system of photonic crystal slab plus surrounding environment is gapless above the light-cone, meaning that according to band theory the system’s topological protection is not well-defined. Thus, as already noted by Raghu and Haldane25, the radiative environment has traditionally limited the notion of topological protection—as defined by topological band theory—in photonic topological insulators. Therefore, although topological band theory can be applied to photonic systems to yield some insight into their topological properties, this cannot be viewed as a complete picture as the setting goes beyond the scope of band theoretic approaches; any complete picture of topological photonic crystal slabs must account for their inherent out-of-plane radiative losses, and provide a definition of topological protection despite these losses.

Fig. 1: Gapless environment for photonic crystal slabs.
figure 1

a Schematic of a free-standing photonic crystal slab in a three-dimensional (3D) geometry. Photonic structures are inherently 3D and are usually surrounded by an homogeneous material that features a light cone. As such, out-of-plane radiative losses, depicted by the red arrows, are inherent to such structures. The photonic crystal is a triangular lattice with lattice constant a composed of dielectric rods, \({\bar{\epsilon }}_{jj}=14\) for j = x, y, z, of radius r = 0.37a and height t = 0.5a embedded in a gyro-electric material slab, \({\bar{\epsilon }}_{jj}=1\) and \({\bar{\epsilon }}_{xy}=-0.4i\), of thickness t = 0.5a. b Band structure of the photonic crystal slab in (a) over the first Brillouin zone, for a transverse magnetic-like polarization. The photonic band gap at around ω = 0.42[2πc/a] is the “topological band gap” known from extrapolation from the two-dimensional photonic crystal approximation. The shaded region depicts those frequencies and wavevectors that are at, or above, the light line of the surrounding air. The red line depicts the light line, ω = ck. Above the light line, photonic crystal slabs generally exhibit resonances, not bound states.

Here, we present a general framework for classifying topology in realistic three-dimensional (3D) photonic systems that directly accounts for both in-plane and out-of-plane radiative losses. To do so, we solve two interconnected challenges, demonstrating how to perform dimensional reduction to calculate invariants of 2D systems, such as Chern numbers, for photonic crystal slab systems that are inherently 3D, and showing how this generalizes to non-Hermitian systems so that radiative losses can be properly accounted for. The framework we develop is rooted in the spectral localizer, which combines the position matrices and the Hamiltonian matrix of the system to draw a local picture of a system’s topology26. Using the operator-based approach of the spectral localizer, we construct a spectral FEM-localizer built on an effective Hamiltonian directly derived from the full-wave Maxwell equations via finite-element method (FEM). We demonstrate the utility of the FEM-localizer framework through two fundamental examples in topological photonics: a 2D photonic Chern insulator and a 2D photonic Chern quasicrystal, proving that the spectral FEM-localizer approach is applicable to aperiodic structures that lie beyond the scope of topological band theory. Finally, we apply the spectral FEM-localizer to study the topology of a photonic Chern slab, showing how to generate a local strong 2D invariant for the 3D system while directly accounting for the out-of-plane radiative losses. As a significant portion of photonic topological insulators are metasurfaces, the proposed classification method will benefit to the characterization of topological metasurfaces where topological boundary modes are used to control radiation, scattering, and emission27,28,29,30,31. Given the wide variety of systems that can be faithfully approximated by FEMs, we anticipate that our framework has broad applicability both within photonics and beyond, and will be useful for the study of the topology in complex physical systems where a band theoretic picture is not available.

Results

Overview of the spectral localizer

Despite being a real-space approach to topology, the spectral localizer shares some conceptual similarities with topological band theory. A modern understanding of traditional band theoretic approaches to topology can be built from the concept of atomic limits—the limit in which the couplings between adjacent atoms, molecules, or structural decorations are turned off, and the system’s band structure becomes completely flat. Atomic limits are topologically trivial as they always possess a complete symmetry-preserving localized Wannier basis32,33. Moreover, different systems are topologically equivalent if one system can be smoothly deformed (i.e., path continued) into the other system without closing the relevant bulk band gap (i.e., the band gap at the frequency of interest) or breaking any necessary symmetry. Thus, from this perspective, the question of material topology becomes one of whether a system can be path continued to an atomic limit; if so, it is trivial. Band theoretic approaches offer a few different methods for making this determination, either using standard topological invariants21, or by comparing a system’s band structure against the possible elementary band representations34.

In contrast, the spectral localizer framework has emerged as a method to diagnose a system’s topology from the real-space perspective of the atomic limit. The key idea is that just as the wavevector-space description of atomic limits is as a material with completely flat bands35, a real-space description of atomic limits can be understood in terms of the system’s position matrices Xj and Hamiltonian H via the commutation relations

$$\left[{X}_{j}^{({{{\rm{AL}}}})},{H}^{({{{\rm{AL}}}})}\right]=0,\quad j=1,\ldots ,d,$$
(1)

with d the dimension of the system. In other words, in an atomic limit, H(AL) is block diagonal, with each block corresponding to each decoupled atom, molecule, or structural element. Likewise, the position of each such interior degree of freedom is condensed into a single location. Hence, \({X}_{j}^{({{{\rm{AL}}}})}\) and H(AL) commute for atomic limits.

From this real-space perspective, the question of topology then becomes one of understanding whether a system’s non-commuting H and Xj matrices are nevertheless path continuable to commuting matrices, while preserving any necessary symmetry and maintaining the relevant bulk spectral gap. To perform this assessment directly in real-space for a d-dimensional Hermitian system in any of the ten Altland-Zirnbauer symmetry classes35,36,37,38, one first forms the system’s spectral localizer by combining Xj and H using a Clifford representation26,

$$\begin{array}{ll}{L}_{({x}_{1},\ldots ,{x}_{d},E)}({X}_{1},\ldots ,{X}_{d},H) =\\ \mathop{\sum }\limits_{j=1}^{d}\kappa ({X}_{j}-{x}_{j}I)\otimes {\Gamma }_{j}+(H-EI)\otimes {\Gamma }_{d+1}.\end{array}$$
(2)

Here, Γ1, …, Γd+1 are (d + 1)-dimensional Clifford representation satisfying \({\Gamma }_{j}^{{\dagger} }={\Gamma }_{j}\), \({\Gamma }_{j}^{2}=I\), and ΓjΓl = − ΓlΓj for j ≠ l, while I is the identity matrix. The spectral localizer is an inherently local approach to material topology; in Eq. (2), (x1, …, xd, E) ≡ λ is the spatial coordinate (x1, …, xd) and energy E where the system’s topology is being probed. Finally, κ is a hyperparameter chosen to make the units consistent between the position and Hamiltonian matrices, and to additionally balance the emphasis on the system’s position information relative to its Hamiltonian. Typically, for gapped systems, it suffices to choose κ ~ 2g/L39,40,41 where g is the width of the bulk band gap and L is the length of the finite system considered. Extensive studies have demonstrated the spectral localizer’s versatile applicability across a broad range of topological systems, even beyond the traditional topological band theory26,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53.

Using results from the study of C*-algebras26,39,43, the spectrum of the spectral localizer has been proven to be connected to local topological markers for every discrete symmetry class (i.e., Altland-Zirnbauer class35,36,37,38) for every physical dimension. For a 2D Hermitian system in class A, namely a system lacking any discrete symmetries, the appropriate local topological invariant is the local Chern number, defined as

$${C}_{(x,y,E)}^{{{{\rm{L}}}}}(X,Y,H)=\frac{1}{2}{{{\rm{sig}}}}\left[{L}_{(x,y,E)}(X,Y,H)\right],$$
(3)

where sig denotes the matrix’s signature, the difference between its total number of positive and negative eigenvalues. A non-zero local Chern number at (x, y, E) indicates that the re-centered position (X − xI), (Y − yI), and Hamiltonian (H − EI) matrices cannot be path continued to be commuting, meaning that the system cannot be deformed to the atomic limit. In other words, a non-zero local Chern number, \({C}_{(x,y,E)}^{{{{\rm{L}}}}}\,\ne\, 0\) tells us that the system is topologically distinct from the atomic limit, and is locally topologically non-trivial at the spatial and energy coordinate (x, y, E).

Separate from its ability to identify material topology, the spectral localizer can be seen as a tool to identify localized states by looking at an approximate state ψ of both the position and Hamiltonian matrices of the system

$${X}_{j}{{{\boldsymbol{\psi }}}}\,\approx\, {x}_{j}{{{\boldsymbol{\psi }}}}\quad \,{{\mbox{and}}}\,\quad H{{{\boldsymbol{\psi }}}}\,\approx\, E{{{\boldsymbol{\psi }}}},$$
(4)

In particular, the localizer gap, defined as the smallest singular value of the spectral localizer

$$\begin{array}{ll}{\mu }_{({x}_{1},\ldots ,{x}_{d},E)}^{{{{\rm{C}}}}}({X}_{1},\ldots ,{X}_{d},H)=\\ \min \left\{\left|\sigma \left({L}_{({x}_{1},\ldots ,{x}_{d},E)}({X}_{1},\ldots ,{X}_{d},H)\right)\right|\right\},\end{array}$$
(5)

with σ(L) being the set of eigenvalues of L, gives a measure about large of a modification of the system is needed to realize such state ψ. As such, a localizer gap closing \({\mu }_{({x}_{1},\ldots ,{x}_{d},E)}^{{{{\rm{C}}}}}=0\) indicates that there exists such an approximate state ψ at energy E that is localized at spatial position (x1, …, xd).

Altogether, the local topological markers and the localizer gap give a consistent picture of the topology locally in space and energy. The local topological markers are used to probe the topology locally at space-energy coordinate (x1, …, xd, E), and cannot change as long as the localizer gap does not close (\({\mu }_{({x}_{1},\ldots ,{x}_{d},E)}^{{{{\rm{C}}}}}\,\ne\, 0\)). When the local markers change across some spatial or energy path from one topological phase to another topological phase, the localizer gap must close, resulting in a localized state [Eq. (4)] at the interface between the two topological distinct phases: this is precisely bulk-edge correspondence.

The spectral localizer is not the only real-space theory of material topology that provides local topological markers54,55,56,57,58,59,60,61,62,63,64. However, it is the only currently known theory of topology that preserves system sparsity, i.e., if H is sparse, \({L}_{({x}_{1},\ldots ,{x}_{d},E)}\) is sparse. Thus, unlike other local markers that typically require projecting into an occupied subspace, yielding still-relatively-large dense matrices, the calculation of local markers using the spectral localizer can leverage advances in sparse matrix algorithms to realize substantial numerical speedups. For example, finding a sparse matrix’s signature does not require finding any of its eigenvalues, and can instead be found using Sylvester’s law of inertia65,66.

Building the spectral FEM-localizer

The spectral localizer is used to study a system’s topology directly from its equations of motion, i.e., its master equation. For example, it is straightforward to apply the spectral localizer framework to tight-binding models26,42,44,46,47,48,49,50,51,52,53 where the position operators are diagonal matrices with entries being the spatial position (x1, …, xd) of the model’s sites, and H the tight-binding Hamiltonian. However, the spectral localizer is not limited to such approximate descriptions of physical system. Instead, so long as a system admits a discretization of some form (or some other method for generating a bounded matrix description of the system), the matrix of its discretized master equation can be inserted into the spectral localizer [Eq. (2)], where practical parameters and geometry can be used. With such a choice of discretization, the entries of the position matrices Xj can then be chosen as the grid mesh positions from discretization of the master equation and the Hamiltonian matrix H can be chosen as being any matrix whose eigenvector is a solution of the master equation, namely a matrix whose eigenproblem is physically meaningful in order to have a relevant joint-spectrum problem. As such, for the study of the topology in physical system within the spectral localizer framework, no further approximation is required beyond the discretization of the master equation, yielding a potentially more accurate description of a system’s topology.

The photonic master equation

The master equations for photonic systems are given by Maxwell’s equations. Assuming that the materials used are linear and time-independent, the electromagnetic fields can be written in a time-harmonic form eiωt such that a photonic system can be described by the source-free Maxwell’s equations

$$\nabla \times {{{\boldsymbol{E}}}}({{{\boldsymbol{x}}}})=i\omega \bar{\mu }({{{\boldsymbol{x}}}}){{{\boldsymbol{H}}}}({{{\boldsymbol{x}}}}),$$
(6)
$$\nabla \times {{{\boldsymbol{H}}}}({{{\boldsymbol{x}}}})=-i\omega \bar{\epsilon }({{{\boldsymbol{x}}}}){{{\boldsymbol{E}}}}({{{\boldsymbol{x}}}}),$$
(7)
$$\nabla \cdot [\bar{\epsilon }({{{\boldsymbol{x}}}}){{{\boldsymbol{E}}}}({{{\boldsymbol{x}}}})]=0,$$
(8)
$$\nabla \cdot [\bar{\mu }({{{\boldsymbol{x}}}}){{{\boldsymbol{H}}}}({{{\boldsymbol{x}}}})]=0,$$
(9)

where ω is the angular frequency, E(x) and H(x) are the electric and magnetic fields, and \(\bar{\epsilon }({{{\boldsymbol{x}}}})\) and \(\bar{\mu }({{{\boldsymbol{x}}}})\) are the permittivity and permeability tensors. Previously, the spectral localizer framework has been applied to describe the topology of 2D photonic crystals40,41,45. These prior studies took advantage of finite-difference discretizations of Maxwell equations [Eqs. (6) and (7)] to automatically satisfy the divergence-free condition [Eqs. (8) and (9)]. However, this approach is not easily scalable, as it requires a uniform mesh and therefore requires very large matrices for the simulation of realistic designs of 3D systems.

Instead, here, we use the finite-element method (FEM) for a more general approach to any physical system described by a master equation. The key advantage of FEM-localizer framework is that it allows for a much coarser meshing and therefore reduces the size of the matrices involved. Moreover, the method easily incorporates physical processes with different characteristic length scales, where additional equations of motion can be included to described the relevant processes and their couplings. For example, dispersion can be described by introducing auxiliary equations for the material’s internal degrees of freedom responsible for this effect67 that are then coupled to Maxwell’s equations [Eqs. (6)–(9)]. To best make use of the FEM approach for photonic systems, here we use the Helmholtz equation for the electric field E(x), derived from Eqs. (6)–(9)

$$\nabla \times \left({\bar{\mu }}^{-1}({{{\boldsymbol{x}}}})\nabla \times {{{\boldsymbol{E}}}}({{{\boldsymbol{x}}}})\right)-{\omega }^{2}\bar{\epsilon }({{{\boldsymbol{x}}}}){{{\boldsymbol{E}}}}({{{\boldsymbol{x}}}})=0,$$
(10)

as the master equation to probe the topology in photonic systems.

Overview of the finite-element method

In this section, we briefly outline FEM discretizations, with a focus on the relevant information for then incorporating the resultant equations of motion into the spectral localizer framework. A more detailed description of FEM can be found in FEM textbooks24.

The FEM discretization starts by reformulating the master equation into its weak form, where the differential equation is no longer satisfied exactly at every point of the mesh. The weak form of the master equation consists in turning the differential equation into an integral equation, via the method of integration by parts, to improve numerical stability. Indeed, the differentiation of the solutions of the master equation may be limited at the boundaries of the simulation domain or at some material interfaces, where a jump in their values can be observed. The weak form is therefore obtained by multiplying the master equation with a weight function α(x) and by evaluating the overlap integral over the simulation region Ω. The weak form for the master equation Eq. (10) is then

$$\begin{array}{ll}{\int_{\Omega }}\left({\bar{\mu }}^{-1}({{{\boldsymbol{x}}}})\nabla \times {{{\boldsymbol{E}}}}({{{\boldsymbol{x}}}})\right)\cdot \left(\nabla \times {{{\boldsymbol{\alpha }}}}({{{\boldsymbol{x}}}})\right)d{{{\boldsymbol{x}}}}\\ -{\omega }^{2}{\int_{\Omega }}\bar{\epsilon }({{{\boldsymbol{x}}}}){{{\boldsymbol{E}}}}({{{\boldsymbol{x}}}})\cdot {{{\boldsymbol{\alpha }}}}({{{\boldsymbol{x}}}})d{{{\boldsymbol{x}}}}=0.\end{array}$$
(11)

The solution vector, E(x) for Eq. (10), is decomposed in terms of the shape functions wn

$${{{\boldsymbol{E}}}}({{{\boldsymbol{x}}}})=\mathop{\sum}\limits_{n}{E}_{n}{{{{\boldsymbol{w}}}}}_{n}({{{\boldsymbol{x}}}}),$$
(12)

with unknown weights En located along the extended mesh nodes [see for example the green crosses in Fig. 2]. Importantly, these shape functions are chosen to satisfy the desired properties of the system’s equation. For instance, when considering Maxwell’s equations, it is imperative to fulfill both the divergence-free conditions and interface conditions. In this context, the so-called curl elements68 can be chosen as shape functions. These curl elements inherently satisfy the divergence condition [Eqs. (9) and (8)] as they are divergence-free functions. Furthermore, the curl elements enforce the criterion that the tangential component of the electric field must be continuous while the normal component can be discontinuous across interfaces, in accordance with the interface conditions set by the Maxwell’s equations.

Fig. 2: Position coordinate of the degrees of freedom used for constructing the position matrices within the finite-element method discretization.
figure 2

a Schematic of the discretization into finite elements of the two-dimensional (2D) simulation domain composed of a single dielectric rod (red shaded region) at the center of a triangular unit cell. b Zoom-in of (a) where the mesh nodes and extended mesh nodes are depicted using blue and green markers, respectively. In this context, the 2D structure is solved for transverse magnetic polarization (Hz ≠ 0) and the shape function used are the curl elements. The position of the extended mesh nodes corresponding to the unknown weighting coefficient En are located at (xn, yn).

Altogether, the solution to the master equation is obtained by solving the FEM-discretized master equation after performing numerical integration, which can be written as a system of linear equations in the matrix-form

$${H}_{{{{\rm{eff}}}}}(\omega ){{{\boldsymbol{\Psi }}}}={{{\boldsymbol{0}}}},$$
(13)

where \({{{\boldsymbol{\Psi }}}}={(\ldots ,{E}_{n},\ldots )}^{{{{\rm{T}}}}}\) is the solution vector, the electric field E(x), at frequency ω. The magnetic field H(x) can then be derived using Eq. (6). Notably, Heff(ω) is a good candidate to be an effective Hamiltonian for insertion into the spectral localizer as the eigenvector with zero eigenvalue is a solution to the master equation.

Incorporating boundary conditions

Although there are many available implementations of FEM, here we describe the methodology based on the FEM discretization from the commercial software COMSOL MULTIPHYSICS69 as this is a widely used software across many different physical platforms. Using details from this specific FEM implementation, we then show how to develop a Hamiltonian that incorporates the system’s boundary conditions.

Using the Eigenvalue Solver algorithm in COMSOL, the FEM discretization leads to solving the following set of equations for the solution vector in the extended mesh Ψ,

$${H}_{{{{\rm{eff}}}}}(\omega ){{{\boldsymbol{\Psi }}}}+{N}_{F}{{{\boldsymbol{\Lambda }}}}={{{\boldsymbol{0}}}},$$
(14)
$$N{{{\boldsymbol{\Psi }}}}={{{\boldsymbol{0}}}},$$
(15)

where NF and N are respectively the constraint force Jacobian matrix and the constraint Jacobian matrix, and Λ is a vector made of the Lagrange multipliers for the boundary conditions and fictitious degrees of freedom. In Eqs. (14) and (15), Heff is a matrix-valued function

$${H}_{{{{\rm{eff}}}}}:\omega \to {(-i\omega )}^{2}M-(-i\omega )C+K,$$
(16)

where \(\omega \in {\mathbb{C}}\), M is the mass matrix, C is the damping matrix, and K is the stiffness matrix. The real part of ω corresponds to the solution’s angular frequency while its imaginary part is its decay rate. In sum, Eqs. (14) and (15) contain the discretized weak form of the master equation [see the term Heff(ω)Ψ in Eq. (14)], as well as additional terms that account for the boundary conditions.

In order to remove the additional terms that incorporate the boundary conditions and solve an equation that resembles Eq. (13), the constraint equation is first solved. Namely, Eq. (15) is solved as

$$N{{{{\boldsymbol{\Psi }}}}}_{d}={{{\boldsymbol{0}}}}$$
(17)

with Ψd = 0. Then

$${{{\boldsymbol{\Psi }}}}=\left({{{{\boldsymbol{\Psi }}}}}_{d}+\,{{\mbox{Null}}}\,{{{{\boldsymbol{\Psi }}}}}_{c}\right)$$
(18)

is a solution of Eqs. (14) and (15). Equations (14) and (15) are therefore reduced to find Ψc, the solution of the eliminated matrix equation

$${H}_{{{{\rm{eff}}}},c}(\omega ){{{{\boldsymbol{\Psi }}}}}_{c}={{{\boldsymbol{0}}}}$$
(19)

where Heff,c is the matrix-valued function defined as

$${H}_{{{{\rm{eff}}}},c}(\omega )={{{\mbox{Nullf}}}}^{T}{H}_{{{{\rm{eff}}}}}(\omega )\,{{\mbox{Null}}}\,$$
(20)

with Null and Nullf being composed of basis vectors spanning the null space of N and \({N}_{F}^{{{{\rm{T}}}}}\), respectively,

$$N\,{{\mbox{Null}}}\,={{{\boldsymbol{0}}}},$$
(21)
$${{{\mbox{Nullf}}}}^{T}{N}_{F}={{{\boldsymbol{0}}}}.$$
(22)

Physically, the eliminated matrix equation [Eq. (19)] considers the eliminated effective Hamiltonian Heff,c(ω) where all the degrees of freedom involved in the boundary conditions have been accounted for and removed. Consequently, Heff,c(ω) is the effective Hamiltonian compatible with the spectral localizer framework, and not Heff(ω).

The spectral localizer via finite-element method

As the eigenvector of Heff,c(ω) with eigenvalue zero is physically meaningful because it corresponds to a solution of eliminated Maxwell’s equations, Heff,c(ω) can be chosen to be the effective Hamiltonian matrix in the spectral localizer framework. Accordingly, the position matrices Xj,c need to be constructed in the same vector space as Heff,c(ω). In the non-eliminated vector space, the position matrices Xj are diagonal matrices with entries corresponding to the j-th coordinate position of the n-th FEM degrees of freedom xj,n, i.e., the positions of the extended mesh nodes [see for example the green crosses in Fig. 2],

$${X}_{j}=\left(\begin{array}{ccc}\ddots &&\\ &{x}_{j,n}&\\ &&\ddots \end{array}\right).$$
(23)

The position matrices in the eliminated space Xj,c are therefore obtained by projecting Xj onto the eliminated space via Null and Nullf as

$${X}_{j,c}={{{\mbox{Nullf}}}}^{T}{X}_{j}{{\mbox{Null}}}\,.$$
(24)

Given the unique mathematical and computational complexities associated with using a FEM, there are several modifications that must be made to the spectral localizer framework [Eq. (2)] to probe the systems local topology. First, to probe the topology at spatial position coordinate (…, xj, …), xj needs to be expressed in the eliminated space in accordance to the position matrix Xj,c. This can be implemented by directly projecting (Xj − xjI), appearing in Eq. (2), onto the eliminated space, \({({X}_{j}-{x}_{j}I)}_{c}={{{\mbox{Nullf}}}}^{T}({X}_{j}-{x}_{j}I){{\mbox{Null}}}\,\), and by probing at spatial position zero namely λ = (0, …, 0, E). Second, as only the zero eigenvalue of Heff,c corresponds to a physically meaningful solution of the Maxwell’s equations, the localizer has to be probed at λ = (…, 0), namely the effective Hamiltonian Heff(ω) itself carries the information about the frequency ω at which the system’s topology is being classified. Thus, for shorthand notation, the topology at location (x1, …, xd, ω) is probed using the spectral FEM-localizer defined as

$$\begin{array}{ll}{\hat{L}}_{({x}_{1},\ldots ,{x}_{d},\omega )}({X}_{1,c},\ldots ,{X}_{d,c},{H}_{{{{\rm{eff}}}},c})=\\ {L}_{(0,\ldots ,0,0)}({X}_{1,c}-{x}_{1}{I}_{c},\ldots ,{X}_{d,c}-{x}_{d}{I}_{c},{H}_{{{{\rm{eff}}}},c}(\omega )),\end{array}$$
(25)

with

$${I}_{c}={{{\mbox{Nullf}}}}^{T}{{\mbox{Null}}}\,.$$
(26)

For example, for a 2D Hermitian system, the spectral FEM-localizer can be written using the Pauli spin matrices as the choice of Clifford representation, yielding

$$\begin{array}{l}{\hat{L}}_{(x,y,\omega )}({X}_{c},{Y}_{c},{H}_{{{{\rm{eff}}}},c})=\\ \left(\begin{array}{cc}{H}_{{{{\rm{eff}}}},c}(\omega )&\kappa ({X}_{c}-x{I}_{c})-i\kappa ({Y}_{c}-y{I}_{c})\\ \kappa ({X}_{c}-x{I}_{c})+i\kappa ({Y}_{c}-y{I}_{c})&-{H}_{{{{\rm{eff}}}},c}(\omega )\end{array}\right),\end{array}$$
(27)

and the local topological marker, for class A systems, at position (x, y) and angular frequency ω is obtained from

$$\begin{array}{ll}{C}_{(x,y,\omega )}^{{{{\rm{L}}}}}({X}_{c},{Y}_{c},{H}_{{{{\rm{eff}}}},c})=\\ \displaystyle\frac{1}{2}{{{\rm{sig}}}}\left[{\hat{L}}_{(x,y,\omega )}({X}_{c},{Y}_{c},{H}_{{{{\rm{eff}}}},c})\right].\end{array}$$
(28)

Altogether, the effective Hamiltonian Heff,c(ω) solves the system’s master equation rigorously with the only approximation being the discretization, taking into account all the possible processes in the physical system. The retained information in the effective Hamiltonian therefore gives us a more rigorous description of the topology in the physical system. For instance, in photonic systems, this approach can directly incorporate the radiative processes overlooked in the literature by using the non-Hermitian line-gap extension of the spectral localizer49. Instead of using the Hermitian localizer [Eq. (27)], the non-Hermitian spectral localizer for classifying 2D non-Hermitian (lossy) systems is now written as

$$\begin{array}{ll}{\hat{L}}_{(x,y,\omega )}^{({{{\rm{NH}}}})}({X}_{c},{Y}_{c},{H}_{{{{\rm{eff}}}},c})=\\ \left(\begin{array}{cc}{H}_{{{{\rm{eff}}}},c}(\omega )&\kappa ({X}_{c}-x{I}_{c})-i\kappa ({Y}_{c}-y{I}_{c})\\ \kappa ({X}_{c}-x{I}_{c})+i\kappa ({Y}_{c}-y{I}_{c})&-{H}_{{{{\rm{eff}}}},c}{(\omega )}^{{\dagger} }\end{array}\right),\end{array}$$
(29)

and the local topological marker at position (x, y) and angular frequency ω is obtained from

$$\begin{array}{ll}{C}_{(x,y,\omega )}^{{{{\rm{L}}}},(\,{{\mbox{NH}}})}({X}_{c},{Y}_{c},{H}_{{{{\rm{eff}}}},c})=\\ \displaystyle\frac{1}{2}{{{{\rm{sig}}}}}_{{\mathbb{R}}}\left[{\hat{L}}_{(x,y,\omega )}^{({{{\rm{NH}}}})}({X}_{c},{Y}_{c},{H}_{{{{\rm{eff}}}},c})\right].\end{array}$$
(30)

where \({{{\mbox{sig}}}}_{{\mathbb{R}}}\) denotes the matrix’s difference between its number of positive and negative eigenvalues with respect to their real part. Additionally, the localizer gap becomes

$$\begin{array}{ll} \mu^{{{\rm{C}}}, ({{\text{NH}}})}_{(x,y,\omega)}(X_c, Y_c, H_{{{\rm{eff}}},c}) = \\ \min \left\{ \left| {{\text{Re}}} \left[\sigma \left({{\hat{L}}}^{({{\rm{NH}}})}_{(x, y, \omega)}(X_c, Y_c, H_{{{\rm{eff}}},c}) \right) \right] \right| \right\}. \end{array}$$
(31)

Notably, ω can be complex in the lossy system and can include the damping term for calculating the local markers. However, by using the line-gap extension of the spectral localizer for 2D class A systems, the imaginary part of ω should not matter in the calculation of the topology49.

Finally, it is emphasized here that the matrices for constructing the spectral localizer can be readily obtained from COMSOL once the Eigenvalue Solver study has been run, regardless of the module used: the matrices M, C, K, Null, and Nullf for determining Heff,c and the spatial position (…, xj, …) of the extended mesh nodes can be directly accessed from the COMSOL functions in MATLAB70. The spectral FEM-localizer framework can therefore be immediately applied to the wealth of examples and designs considered by the COMSOL community.

Example of 2D photonic Chern structures

As an initial demonstration of the versatility of the spectral FEM-localizer framework, we consider two fundamental examples in topological photonics: a 2D Haldane photonic crystal heterostructure25 that is the canonical photonic Chern insulator [Fig. 3], and a 2D photonic quasicrystal71 that is an aperiodic system where topological band theory cannot be applied [Fig. 4]. In both cases, the photonic system is first described by the full-wave Helmholtz equation (using the Wave Optics module) and then discretized using FEM via COMSOL MULTIPHYSICS. Finally, the spectral FEM-localizer is constructed from the Eigenvalue Solver study and used to classify the systems’ topology.

Fig. 3: Probing of the local topology in a two-dimensional photonic Chern crystal system.
figure 3

a Design of the 2D photonic Chern heterostructure. The inner parallelogram is a non-trivial topological lattice with lattice constant a made of dielectric rods, \({\bar{\epsilon }}_{jj}=14\) for j = x, y, z, with radius r = 0.37a embedded in a gyro-electric background, \({\bar{\epsilon }}_{jj}=1\), \({\bar{\epsilon }}_{xy}=-0.4i\), that breaks time-reversal symmetry. The outer lattice is a topologically trivial lattice with the same lattice constant a composed of air rods, \({\bar{\epsilon }}_{jj}=1\), with radius r = 0.35a in a dielectric background with \({\bar{\epsilon }}_{jj}=5.5\). b Ribbon band structure along the Γ−K direction for the transverse magnetic modes (Hz ≠ 0), with a topological edge mode present around ω0 = 0.37[2πc/a]. c Spectrum of the FEM-localizer \(\sigma ({\hat{L}}_{{{{\boldsymbol{\lambda }}}} = (x,{y}_{0},{\omega }_{0})})\) normalized by 10−4Heff,c(ω0) and the local Chern number \({C}_{{{{\boldsymbol{\lambda }}}} = (x,{y}_{0},{\omega }_{0})}^{{{{\rm{L}}}}}\) along the green line in (a) at y0 = 0 and frequency ω0 = 0.37[2πc/a], with \(\kappa =1.5 \left[1{0}^{-4}\parallel \! {H}_{{{{\rm{eff}}}},c}({\omega }_{0}) \! \parallel / \parallel \! {X}_{c} \! \parallel \right]\). d Local density of states (LDOS) for the Hz component of the field at ω0 = 0.37[2πc/a].

Fig. 4: Probing of the local topology in a two-dimensional photonic Chern quasicrystal system.
figure 4

a Design of the photonic Chern quasicrystal based on a Penrose tiling. The rhombuses composing the Penrose tiling have sides of length a, and the dielectric rods are located at the vertices of the Penrose tiling in a gyro-electric background, \({\bar{\epsilon }}_{jj}=1\), \({\bar{\epsilon }}_{xy}=-0.4i\), with radius r = 0.18a. b Spectra of the photonic system along the real axis, with continuous periodic (PBC) and with perfect electric conductor (PEC) boundary conditions to simulate the system without and with edges, respectively. c Truncated design of the photonic Chern quasicrystal in (a) with perfect electric conductors (PEC) boundary conditions. d Local density of states (LDOS) for the Hz component of the field at ω0 = 0.37[2πc/a] with the geometry in panel (c). e Spectrum of the FEM-localizer \(\sigma \left({\hat{L}}_{{{{\boldsymbol{\lambda }}}} = (x,{y}_{0},{\omega }_{0})} \right)\) normalized by 10−4Heff,c(ω0) and the local Chern number \({C}_{{{{\boldsymbol{\lambda }}}} = (x,{y}_{0},{\omega }_{0})}^{{{{\rm{L}}}}}\) along the green line in (a) at y0 = 0 and frequency ω0 = 0.37[2πc/a], with \(\kappa =1 \left[1{0}^{-4}\parallel \! {H}_{{{{\rm{eff}}}},c}({\omega }_{0}) \! \parallel /\parallel \! {X}_{c} \! \parallel \right]\). f Local density of states (LDOS) for the Hz component of the field at ω0 = 0.37[2πc/a] with the geometry in panel (a).

2D Photonic Chern crystal

For a first example, we focus on the 2D Haldane photonic heterostructure25 shown in Fig. 3a. Here, the Haldane heterostructure is made of two triangular lattices of lattice constant a with perfect electric conductor (PEC) boundary conditions, and we are considering the topology of the transverse magnetic modes (Hz ≠ 0). The inner lattice is a topologically non-trivial insulator composed of dielectric rods, \({\bar{\epsilon }}_{jj}=14\) for j = x, y, z, with radius r = 0.37a embedded in a gyro-electric background, \({\bar{\epsilon }}_{jj}=1,{\bar{\epsilon }}_{xy}=-0.4i\), to break time-reversal symmetry. The outer lattice is a topological trivial insulator composed of air rods, \({\bar{\epsilon }}_{jj}=1\), of radius r = 0.35a in a dielectric background with \({\bar{\epsilon }}_{jj}=5.5\). The interface between the two topological distinct lattices yields a topological edge state that can be seen in the ribbon band structure [Fig. 3b] and from the local density of states at ω0 = 0.37[2πc/a] [Fig. 3d].

The spectral FEM-localizer and local Chern marker in Eq. (27) and Eq. (28), respectively, can be used to diagnose the topology of this lossless 2D system. However, one should note that the eliminated effective Hamiltonian Heff,c is non-Hermitian due to the projection onto the eliminated space. Nevertheless, the non-Hermitian part is found to be negligible and only the Hermitian part is kept in Heff,c as Hermitian materials and lossless boundary conditions have been used. The spectrum of the spectral FEM-localizer \(\sigma ({\hat{L}}_{(x,{y}_{0},{\omega }_{0})})\) and the local Chern number \({C}_{(x,{y}_{0},{\omega }_{0})}^{{{{\rm{L}}}}}\) along the path depicted by the green line in Fig. 3a at y0 = 0 and frequency ω0 = 0.37[2πc/a] are shown in Fig. 3c, demonstrating the local topological picture of the heterostructure. As expected from topological band theory, inside the topological band gap at around ω0 = 0.37[2πc/a], the inner lattice is topological with CL = 1 while the outer lattice is trivial with CL = 0. Therefore, the spectral FEM-localizer correctly captures the change of topology as demonstrated by the eigenvalue crossing with respect to zero of the spectrum of L near the heterostructure’s interface.

2D Photonic Chern quasicrystal

As a second example, we investigate the topology of a magnetooptic 2D photonic quasicrystal surrounded by a homogeneous material, again for transverse magnetic modes (Hz ≠ 0). Notably, this example cannot be classified using topological band theory as the system is not periodic and therefore does not possess a band structure. Moreover, the quasicrystal is surrounded by a homogeneous material that is gapless rather than gapped (or insulating), yielding an ill-defined notion of bulk topological invariant and topological robustness. The topological quasicrystal is constructed from a Penrose tiling72, with the rhombuses having sides of length a and where the dielectric rods are positioned on the vertices of the tiling71, as shown in Fig. 4a. The photonic quasicrystal is composed of dielectric rods with permittivity \({\bar{\epsilon }}_{jj}=14\) and radius r = 0.18a, embedded in a gyro-electric background, \({\bar{\epsilon }}_{jj}=1,{\bar{\epsilon }}_{xy}=-0.4i\) to break time-reversal symmetry, and PEC boundary conditions are used. In this quasicrystal heterostructure [Fig. 4a], the identification of any edge state resonances via the local density of states [Fig. 4f] is obscured by the presence of the bulk states present in the homogeneous material surrounding the quasicrystal. Instead, the edge states can be identified from the local density of states only when the surrounding homogeneous material is truncated, as shown Fig. 4c and d. The edge and bulk spectra of the photonic quasicrystal [Fig. 4b] can also be obtained using, respectively, PEC and (continuous) periodic boundary conditions in a polygonal geometry to simulate the photonic system both with and without boundaries (see Supplemental Material).

Similar to the crystalline example, the spectral FEM-localizer and local marker in Eq. (27) and Eq. (28) are used to probe the topology of this photonic system once the non-Hermitian part is removed from the eliminated effective Hamiltonian Heff,c. Figure 4e shows the spectrum of the FEM-localizer \(\sigma ({\hat{L}}_{(x,{y}_{0},{\omega }_{0})})\) and the local Chern number \({C}_{(x,{y}_{0},{\omega }_{0})}^{{{{\rm{L}}}}}\) as the probe location is varied, revealing the topology despite the system’s aperiodicity and the lack of a surrounding insulator. In particular, there is a crossing of the spectrum of the localizer near the boundary of the structural interface, and thus a change of the local Chern number as the probe coordinate is moved from the trivial homogeneous material to the center of the quasicrystal along the green line depicted in Fig. 4a, and at ω0 = 0.37[2πc/a]. The location where the local index changes, associated to a vanishing localizer gap, is also an indication of the location of the topological edge state [Eq. 4], which cannot be resolved from the system’s local density of states [Fig. 4f]. A similar plot for the spectrum of the spectral localizer and the local Chern number can be realized along the frequency axis, as shown in the Supplemental Material, demonstrating some range of frequencies ω for which the quasicrystal is topologically non-trivial. As such, this example highlights how the spectral FEM-localizer is capable to identify the topology of the system without the need of a band structure or a bulk band gap, and how the location of the topological edge states can be explicitly determined from the vanishing of the localizer gap.

Photonic Chern slab

Photonic structures are inherently 3D, and as such the 2D photonic designs extensively studied in the literature and described in the previous section are only approximations to actual photonic slabs realized experimentally. While there are standard methods for developing 3D photonic crystal slabs whose resonance band structures are quantitatively similar to the band structures of 2D photonic crystals73, these methods only approximate the Hermitian portion of the band structure, not the radiative portion. Moreover, approximating radiative losses due to a system’s environment as material absorption within the system is uncontrolled, as this is fundamentally relocating degrees of freedom that were outside of a structure to be inside of it. Indeed, this is a particularly problematic approximation for topological systems, whose primary features are boundary-localized states—relocating degrees of freedom changes what it means for a state or resonance to be localized. Therefore, out-of-plane radiative losses, an inherent aspect of photonic slabs, have been neglected in previous theoretical treatments of topological photonics, as the radiative loss cannot be accounted for using a band theoretic approach. Nevertheless, as we have already seen in the case of the quasicrystal, the spectral FEM-localizer allows to diagnose the topology beyond the scope of topological band theory. Here, we show that the spectral FEM-localizer can be used to classify the topology of photonic crystal slabs and directly incorporate out-of-plane radiative losses using radiative boundary conditions.

As an example illustrating the probing of the topology in photonic slab while accounting for the out-of-plane radiative loss, we consider a free-standing photonic crystal slab embedded in air, as shown in Fig. 5a. This is the 3D slab version of the Haldane photonic heterostructure studied in Fig. 3a, with PEC on the x- and y-boundaries and a radiative boundary condition implemented through perfectly matched layers (PML) on the z-boundaries. The parameters used are the same as in the Chern heterostructure example [Fig. 3a] except that now the background is a slab of thickness t = 0.5a embedded in air, and the rods have the same finite height t = 0.5a. The topological edge state can be seen in the ribbon band structure [Fig. 5b] and from the local density of states at z = 0 and at ω0 = 0.42[2πc/a] [Fig. 5d].

Fig. 5: Probing of the local topology in a photonic Chern slab system.
figure 5

a Design of the photonic Chern slab. The inner parallelogram is a non-trivial topological lattice with lattice constant a made of dielectric rods, \({\bar{\epsilon }}_{jj}=14\) for j = x, y, z, with radius r = 0.37a and height t = 0.5a embedded in a gyro-electric slab, \({\bar{\epsilon }}_{jj}=1\), \({\bar{\epsilon }}_{xy}=-0.4i\), with thickness t = 0.5a. The outer lattice is a topologically trivial lattice with the same lattice constant a, composed of air rods, \({\bar{\epsilon }}_{jj}=1\), with radius r = 0.35a and height t = 0.5a in a dielectric slab, \({\bar{\epsilon }}_{jj}=5.5\), with thickness t. b Ribbon band structure along the ΓK direction for the transverse magnetic modes (Hz ≠ 0), with the topological mode present around ω0 = 0.42[2πc/a]. The shaded region depicts those frequencies and wavevectors that are at, or above, the light line of the surrounding air. The red line depicts the light line, ω = ck. c Spectrum of the FEM-localizer \(\sigma \left({\hat{L}}_{{{{\boldsymbol{\lambda }}}} = (x,{y}_{0},{\omega }_{0})}^{(NH)}\right)\) normalized by 10−4Heff,c(ω0) along the green line in (a) at y0 = 0 and frequency ω0 = 0.42[2πc/a], with \(\kappa =1.5 \left[1{0}^{-4}\parallel \! {H}_{{{{\rm{eff}}}},c}({\omega }_{0}) \! \parallel /\parallel \! {X}_{c} \! \parallel \right]\). d Local density of states (LDOS) for the Hz component of the field at z = 0 and at ω0 = 0.42[2πc/a].

As the system is 3D and in class A in the Altland-Zirnbauer classification35,36,37,38, the topology for the 2D topological edge state in the slab can be classified using an integer invariant such as the Chern number. Within the spectral localizer framework, the topology is diagnosed by disregarding the z-direction (i.e., all mesh vertices are retained, but their coordinates is reduced (x, y, z) → (x, y)), performing a version of dimensional reduction to enable the calculation of a strong 2D invariant of a 3D system. Physically, this is equivalent as looking at the change of topology as we move in the (x, y)-plane, irrespective of the z-coordinate. Despite the removal of the z-direction, all the information from the 3D geometry, including the out-of-plane radiative loss, is retained for the assessment of the topology as the effective Hamiltonian Heff,c is derived from the full 3D geometry. In other words, even though the system is 3D, Eqs. (29) and (30) can still be used— the Xc, Yc, and Heff,c matrices all contain information about the full 3D system, and the Chern number does not depend on the position matrix Z or the z-coordinates. This is because the 2D Chern number is still a strong topological invariant in a 3D system; and calculating it using Xc and Yc means we are looking at edge states around the boundary in the (x, y)-plane. Thus, the 2D non-Hermitian FEM-localizer in Eq. (29) can be used to identify the topology in the slab with radiative boundary conditions, giving an accurate probing of the topology in the photonic slab where all the possible processes are accounted for.

The topology is studied by looking at the spectrum of the spectral FEM-localizer \(\sigma ({\hat{L}}_{(x,{y}_{0},{\omega }_{0})}^{(NH)})\) along the green path in Fig. 5a at y0 = 0, and at the (incomplete) band gap around ω0 = 0.42[2πc/a], as shown in Fig. 5b. The plot demonstrates a net crossing in the spectrum with respect to zero [see red arrow in Fig. 5b], indicating a change of the signature of the spectral FEM-localizer [Eq. (30)], hence a change of topology, near the boundary between the outer and inner lattices, similar to what is observed in Fig. 3c. However, the topological protection given by the localizer gap now takes into account the radiative losses of the topological edge slab state, and as such the protection is weaker than would be predicted from the 2D band structure.

Discussion

Using the operator-based approach of the spectral localizer, we have developed a general framework for studying the topology in realistic photonic structures directly from the discretized master equations of the system using finite-element methods (FEM). In particular, we studied the topology in photonic systems derived directly from the full-wave Maxwell equations. Using the photonic Chern insulators and the photonic Chern quasicrystal, we have demonstrated the ability of the proposed spectral FEM-localizer framework to correctly capture the local topology in photonic topological materials. Moreover, the framework have been applied to a photonic Chern slab predicting genuine topological protection of the topological edge slab state when taking into account possible radiative loss of the slab state. As the radiative feature of optical systems play a significant role for a broad range of photonics applications, we expect that the spectral FEM-localizer’s ability to classify the topology of photonic systems will be useful for developing next-generation devices, including topological metasurfaces which uses topological boundary modes to control radiation, scattering and emission27,28,29,30,31. Looking forward, we anticipate the generality of the framework to be of practical use for systems in any of the ten Altland-Zirnbauer symmetry classes35,36,37,38 as well as in topological crystalline insulators41,45, and for tackling topological problems in other complex physical platforms such as in acoustic systems74, plasmonic systems67,75,76, and in polaritonic systems77,78,79.