Unconventional bound states in the continuum from metamaterial-induced electron acoustic waves

Abstract. Photonic bound states in the continuum (BICs) are spatially localized modes with infinitely long lifetimes, which exist within a radiation continuum at discrete energy levels. These states have been explored in various systems, including photonic and phononic crystal slabs, metasurfaces, waveguides, and integrated circuits. Robustness and availability of the BICs are important aspects for fully taming the BICs toward practical applications. Here, we propose a generic mechanism to realize BICs that exist by first principles free of fine parameter tuning based on non-Maxwellian double-net metamaterials (DNMs). An ideal warm hydrodynamic double plasma (HDP) fluid model provides a homogenized description of DNMs and explains the robustness of the BICs found herein. In the HDP model, these are standing wave formations made of electron acoustic waves (EAWs), which are pure charge oscillations with vanishing electromagnetic fields. EAW BICs have various advantages, such as (i) frequency-comb-like collection of BICs free from normal resonances; (ii) robustness to symmetry-breaking perturbations and formation of quasi-BICs with an ultrahigh Q-factor even if subject to disorder; and (iii) giving rise to subwavelength microcavity resonators hosting quasi-BIC modes with an ultrahigh Q-factor.

To date, photonic BICs involve two mechanisms: (i) symmetry guaranteed mode mismatching and (ii) destructive interference of a topological origin [3].These BICs, while fundamentally interesting, are of limited use in practice as they are readily destroyed by small perturbations [20,21] and coupling to higher diffraction orders [22].Furthermore, a large number of non-BIC modes with a high density of states exists spectrally close to these BICs [4], which considerably hinder the broadband exploitation of BICs.
Here, we propose a generic hydrodynamic double fluid plasma (HDP) effective medium approach to generate new types of photonic BICs in double-net metamaterials (DNMs) consisting of two intertwining, spatially separated metallic networks [23,24].While the discussion in this paper is restricted to Cartesian wire meshes as shown Figure 1, we show that the HDP model equally applies to other DNM geometries in Appendix A (Figure A1 and Figure A2).In contrast to existing designs, the proposed BICs are formed by pure charge fluctuations of the elec-tron acoustic wave (EAW)-like modes in DNMs, which are strikingly robust against a large class of perturbations.They exist within a broad spectral range free of non-BIC resonances and their frequency spacing can be freely tuned by the thickness of the structure, reaching spectrally dense BICs for an optically thick DNM slab.Additionally, new Zak-phase protected topological suface bound states in the continuum (TSBICs) are discovered in DNMs.Both bulk BICs and TSBICs exist at the Γpoint of the in-plane Brillouin zone within the light-cone for DNM slabs, yet they are fully decoupled from radiation.

II. ELECTRONIC ACOUSTIC WAVE IN DOUBLE-NET METAMATERIALS
The new types of BICs originate from a pure charge wave, which does not produce electromagnetic fields.Although impossible in conventional optical materials, this property can be found in non-Maxwellian materials [25], for example in an HDP.The hydrodynamic motion of two charge carriers in combination with Maxwell's equations describes the fluid-like dynamics of an effective double plasma in an electromagnetic field.The full HDP model is derived in Appendix B. In short, we restrict the discussion to a non-dissipative model with only two free parameters: the plasma frequency ω pi := n 0i q 2 i /m i and the thermal pressure parameter κ i := γp i /(n 0i m i c 2 ) for the two charge carriers labelled i = 1, 2. These depend on the equilibrium charge carrier density n 0i , the charges q i , the effective masses m i , the equilibrium pressure p i , the adiabatic exponent γ and the speed of light in vacuum c.
Through a plane-wave ansatz and by linearizing the HDP model [26,27], the generally complex problem is transformed into a 14-dimensional linear eigenproblem for the eigenfrequency ω.As shown in Figure 1(b), from high to low frequencies, a two-fold degenerate electromagnetic transverse band and a longitudinal Langmuir band are obtained, both of which are hyperbolic starting at ω = (ω 2 p1 + ω 2 p2 ) 1/2 for wave number k = 0, and a linear electron acoustic wave (EAW) band emanating from ω = 0.The slope of the EAW band is ± c √ 2 κ 2 1 + κ 2 2 , as analytically shown by k•p theory [28] (Appendix C).Regardless of ω p1,2 , when κ 1 = κ 2 , the EAW consists of two charge waves drifting in opposite direction, with current densitiesj 1 = −j 2 k and vanishing electromagnetic fields E = H = 0 (Figure 1(a)).Due to its zero electromagnetic nature and the resulting decoupling from the electromagnetic vacuum modes, the EAW naturally forms BICs.
To circumvent the difficulties in creating a natural double-plasma [29][30][31], we instead employ a finite DNM lattice, whose fully connected wire morphology effectively provides a freely moving double-plasma with finite equilibrium charge carrier density and Coulomb pressure, leading to a behavior that resembles a plasma with an exact electronic thermal pressure.Note that in DNMs, κ 1,2 are largely determined by the connectivity of the wires, which leads to κ 1,2 ≈ 1 √ 3 (Appendix B 2).The DNM structures discussed in this paper are composed of two single nets of fully connected metal cylinders along the cubic axes with radii r 1,2 , which may be different for the two nets.We use the notation P := (r 1 /a, r 2 /a, R) to represent the network parameters, normalized by the lattice constant a. R is the offset of the two networks that are shifted by S := (1, 1, 1) Ra.Two parameter sets are discussed in this manuscript, The numerically simulated band structure of the P 1 DNM is illustrated in Figure 1(b), featuring a linearly dispersed fundamental band (red dotted line) that emanates from the Γ-point at zero frequency.This result not only qualitatively matches the EAW band in the HDP model, but the dispersion is indeed in good agreement with the analytical prediction from Appendix C ω = 1 √ 3 ck (red dashed line).In the HDP model, the slope of the EAW band is largely independent of the plasma frequencies ω p1,2 (Appendix B), which applies to a variety of DNM realizations (Appendix D), as confirmed by numerical simulations (Figure A6, A7).Analogously to the EAW in the HDP model, the two networks in the DNM carry opposing currents, indicated by arrows in Figure 1(a,b).Simulated surface current densities J are shown in Figure 2(a,b) for DNMs with equivalent radii and different offsets (P 1 and P 2 (1)).They indeed exhibit almost exactly opposite J x in the two networks at k x = π/a and k y,z = 0.For different radii (P 2 (4), Figure 2(c)), a weaker current density is observed on the larger radius network, whose integrated magnitude equals the opposing current density of the lower-radius network.This implies that the total current remains zero, irrespective of the network radii, consistent with the HDP model.
To investigate the nature of the microscopic electromagnetic fields of the EAW-like modes in DNMs, we compute the unit-cell-averaged electric field components of the fundamental mode, normalized by the averaged electric field intensity.The averaged transverse components vanish over the entire Brillouin zone (Figure A5), while the longitudinal component quadratically increases with the wave vector when moving away from Γ (Figure 2(d)).The fundamental modes in DNMs are therefore quasi-longitudinal with E k.We conclude that the HDP model is a valid effective medium description for DNMs near the static limit and that DNMs provide a metamaterial realization of HDPs that require extremely stringent restrictions in natural plasmas, impeding their laboratory observation.We finally note that due to the non-Maxwellian nature of the EAW that lacks electromagnetic fields, the introduction of an effective permittivity seems misleading.In this context, for the P 2 (1) network, a previously reported homogenization approach [32] leading to an effective permittivity ε h for the DNM only accounts for the Langmuir mode residing at high frequencies.In contrast to standard metamaterial homogenization approaches, the HDP model developed here yields a fully consistent effective medium description of the fundamental DNM band and provides an accurate prediction of both the phase velocity and the associated fields.

III. BOUND STATES IN THE CONTINUUM IN DOUBLE NET METAMATERIALS
Here, we show that the EAW character of the bulk modes in DNMs naturally leads to the formation of a new type of BIC, if it is embedded in vacuum in a slab con- figuration, as illustrated in Figure 3(a).In the HDP picture, a standing wave solution exists at certain frequencies where the EAW solutions can be superposed so that the normal component of j i vanishes at the slab boundaries.As the EAW fields satisfy E k, they perfectly decouple from the radiative vacuum modes [33].To illustrate the above prediction, we simulate a 5-layer slab of a P 1 DNM (Figure A18).We apply Bloch-periodic bound-ary conditions in the x−y-plane and scattering boundary conditions in the z-direction, and solve for the complexvalued frequencies.As predicted, multiple BICs with infinite lifetimes are found, located at the Γ-point, marked by red pentagrams in Figure 3(b).The frequencies of the BICs are well approximated by the standing wave solution of the EAW with a predicted phase velocity of 1 √ 3 c and hard wall boundary conditions (j i = 0).The EAW band is well below the first Wood anomaly, so that the higher Bragg orders in the vacuum Rayleigh basis are evanescent and non-radiative [34,35].Because the EAW bulk band is spectrally isolated from the other nonevanescent bulk bands at low frequencies, the BICs are free of non-BIC slab modes in the entire spectral range of the EAW band.In contrast, BICs in photonic crystals are usually mixed with non-BIC modes [3].
At finite wave vectors away from the Γ-point, quasi-BICs are formed since the EAW fields are not strictly orthogonal to the vacuum plane waves, which results in a weak coupling to the vacuum states.These quasi-BICs thus have a very large, but finite, quality factor.The corresponding Q-factor values of the multiple quasi-BIC solutions in the vicinity of the Γ-point are shown in Figure 3(c).The physics is reminiscent of a Fabry-Pérot resonator with mirrors that are perfectly reflecting at normal incidence and otherwise weakly transparent.In scattering experiments, this system yields a simple transmittance spectrum with a Lorentzian line shape at finite wave vectors close to the BIC frequency.The centre frequency forms a band that emanates from the BICs at the Γ-point and hyperbolically blue-shifts with increasing wave vector, while the width of the Lorentzian gradually broadens.The simulated transmittance spectrum along the Γ-X direction is shown in Figure 3(d).The polarization of the quasi-BIC far-field is oriented in the direction of the wave vector, forming a 'star' pattern in 2D momentum space, as shown in Figure 3(e).At the Γ-point where the BICs reside, this generates a singular point where the polarization cannot be defined since the far-field vanishes.Around this point, the polarization pattern evidently exhibits a non-trivial topological winding number [8].In our case, this topological vortex structure of the polarization field is a direct consequence of the quasi-longitudinal nature of the EAW modes and is reminiscent of the electric field lines created by a point charge.
In photonic crystal slabs and metasurfaces, BICs are essentially symmetry protected and can be destroyed by infinitesimal symmetry-breaking perturbations [20,21].Strikingly, the new BICs formed in DNMs are immune to global symmetry breaking and are only prone to perturbations in the individual networks.The DNM slab with P 1 has indeed only one point symmetry, the (110) mirror symmetry, which is insufficient to create symmetryprotected BICs in photonic crystal slabs that require at least two mirror planes.The stability of the DNMbased BICs can be attributed to the surface current densities J on the metallic wires.For the EAW mode at k = π a (0.05, 0, 0) , the currents J x on the two metallic cylinders along the x-direction are exactly opposite (Figure 2(a-c)).J y and J z on each single net resemble electrical quadrupoles, whose form is determined exclusively by the the single net geometry, and which are found irrespective of the global symmetry (Appendix E, Figure A8).Since the interaction between the currents on the two nets is negligibly small in the static limit, the total radiated field from the quadrupole-like current vanishes, irrespective of their relative positions, radii, etc.This ultimately leads to BICs that are robust against global symmetry breaking.This robustness of the BIC modes is demonstrated in Appendix F through simulations of DNM slabs with various parameters (Figure A9).As shown in Appendix G, the BICs can, however, be lifted by breaking the local C 4v symmetry of individual nets (Figure A10), spawning circularly polarized points in the far-field (Figure A11).
Finally, since the BICs originate from the standing wave solutions of the EAWs in the DNM slabs, the number of BIC modes and their frequency spacing are exclusively determined by the thickness of the DNM slab.A full theoretical model that predicts both the quasi-BIC bands including their quality factors (Figure A13) and the transmission spectrum (Figure A12) is presented in Appendix H.In particular, the model yields BIC modes that are equidistantly spaced at frequencies of ωa c = l κπ N , where N a is the thickness of the DNM slab and N > l ∈ N. Extrapolating the linear dispersion relation of the HDP model to the Brillouin zone boundary, we expect N −1 BIC modes over the EAW band, which perfectly matches the full wave simulation results of the P 1 DNM in Figure 3(f).Generally, the theoretical results are in good agreement with the simulations even though the equidistant frequency spacing is of course perturbed in real DNMs.

IV. TOPOLOGICAL SURFACE BOUND STATES IN THE CONTINUUM
Our theory builds on the HDP homogenization model that, while a powerful tool to understand the bulk modes and the formation of BICs, cannot capture the discrete periodicity of DNMs.This discrete lattice structure, however, bears additional features.For example, as predicted in [23], a primitive lattice translation of the P 2 (1) network with BCC symmetry interchanges the two individual networks and results in an EAW band emanating from the corner of the Brillouin zone at the H-point (k = (2π/a, 0, 0)), as shown in Figure A14 [36].A lattice symmetry with an underlying highly symmetric isogonal point group can often be associated with non-trivial topological properties [37][38][39][40].As shown below, the lattice symmetry of the DNMs gives rise to a non-trivial Zak phase and topological surface bound states in the continuum (TSBICs).TSBICs exist at the interface between the DNM slab and vacuum within a radiation continuum, in contrast to typical topological surface states, which reside in bulk bandgaps of two neighbouring domains [41].
The emergence of TSBICs can be understood by examining the family of P 2 (ψ) DNMs.A schematic diagram of a P 2 (4) DNM is shown in Figure 4(a), with the corresponding bandstructure in Figure A15.For ψ = 1 (i.e.equal network radii), the structure has the full body centered cubic (BCC) symmetry.As a result, the EAW and the Langmuir bands meet at X, since they stem from one band that is back-folded half-way between the Γ and the H-point of the BCC Brillouin zone [42].For ψ = 1, a bandgap opens as predicted by the geometrical perturbation theory described in [40].In the presence of inversion and time reversal symmetry, the band topology can be characterized by a Z 2 (binary) topological index, called the Zak phase [43], which is best known for the quantized bulk polarization of electrons in solids.The Zak phase is the well established Berry phase [44] evaluated across a 1D Brillouin zone, which is topologically a circle and therefore closed.It is generally a continuous number on [0, 2π) that depends on the choice of the unit cell.The Zak phase is, however, Z 2 -quantized, with values of 0 and π, if the unit cell has a parity symmetry in its center [45].Therefore, an offset of R = 1  2 between the two nets is necessary to preserve the Cartesian mirror planes in the O h point symmetry of Im3m, allowing a Zak phase classification.
In our case, the Zak phase is obtained by examining parity of the electromagnetic Bloch fields at the high symmetry points of the 1D Brillouin zone [46].We make use of the (001) mirror plane (operator ς z ) in the centre of unit cell in Figure 4(a).The fields must have even (1) or odd (−1) parity with regard to ς z at the Γ and Xpoints, where the wave vector is invariant under ς z .The Zak phase of the band has been shown to be the complex argument of the product of both parities [46].A topological phase transition from trivial to non-trivial Zak states is observed in the projected band structure shown in Figure 4(b) at ψ = 1, where the network has the higher BCC symmetry and the bandgap closes at the X-point of the simple cubic Brillouin zone.The field parities for ψ = 4 are illustrated by the z-components of the periodic part of the electric Bloch field of the EAW mode close to Γ [47] in Figure 4(e) and at X in Figure 4(f).For a DNM slab terminated at its thin wires, we therefore predict, and indeed see, topological surface states within the bandgap, shown as red lines in Figure 4(c).Note that the DNM slab is truncated slightly away from the mirror plane through the thin wires as shown in Figure A18(b), to avoid slicing the thin wire.As discussed in Appendix J, an energy difference between the two surface states is observed, which is caused by hybridization of the surface states on either side of the slab through its finite thickness (Figure A16).The interaction is weakened by increasing the slab thickness, as shown in Figure A15(b).
The topological surface bands in Figure 4(c) are evidently associated with high Q-factor quasi-TSBIC bands revealed by the reflectivity plot in Figure 4(d).TSBICs are, in stark contrast to conventional photonic topological surface states [41], confined through the EAW mechanism and do not require a surrounding topologically trivial band-gap material.When different terminations are applied to the top and bottom surfaces, so that only one surface obtains TSBICs, no resonance is found in the reflection/transmission amplitudes, since no hybridization, which provides a direct transmission channel, occurs.The presence of a single quasi-TSBIC band manifests itself, however, in a resonance in the reflection phase (Figure A17).Despite earlier reports of surface BICs in vi terms of the intrinsic continuum in the material [6], our TSBICs constitute the first example of BICs of topological origin which are decoupled from free space radiation.Controlling the coupling between the topological surface states with free space light could pave the way towards both fundamental scattering properties of topological states and their application, for example in topological antenna arrays [48].

V. CONCLUSIONS
We have demonstrated a new deterministic route towards photonic bound states in the continuum in doublenet metamaterials.A hydrodynamic double fluid plasma serves as an effective medium model for these materials.The hydrodynamic double plasma yields an analytical solution with a new type of quasi-longitudinal EAW band, which gives rise to photonic bound states in the continuum that are robust even against strong geometrical perturbations and symmetry breakings.
Moreover, topological phase transitions are found in double-net metamaterials, accompanied by topological surface bound states in the continuum, protected by a Zak phase topological index.The marriage between the bound-state behavior and the non-trivial topological phase enables the lossless existence of TSBICs even for open boundaries.Our findings endow using doublenet metamaterials as a new platform for the realization of photonic bound states in the continuum by employing a new topological strategy.Advances in 3D fabrication, particularly the 3D printing of metals, enables the manufacture of the proposed double-net metamaterials.While an experimental verification of our predictions is thus readily available at microwave frequencies, we intend to extend our findings to higher frequencies, allowing to, for example, facilitate quantum emission from the coherent spontaneous emission regime, over photonic Bose-Einstein condensation to lasing applications.

ACKNOWLEDGMENTS
We would like to thank Ullrich Steiner for fruitful discussions and general support.We acknowledge funding by the Swiss National Science Foundation through the Spark Grant 190467 and the Project Grant 188647.
Appendix A: Different double network geometries While our theory above has been rigorously derived for the pcu double network only, the two main physical features are inherent to the plasmonic double network topology and symmetry rather than the particular geometrical realization.These features are the slope of the EAW band in the low frequency limit and the longitudinal nature of the homogenized electric field.We here demonstrate that both features are indeed be observed for two other established double net geometries: the balanced double gyroid (srs-c, Ia3d symmetry) and double diamond (dia-c, P n3m symmetry) morphologies (abbreviations as in [49]).
Let us first re-examine the longitudinal nature from a symmetry perspective.In the quasi-static limit, a nontrivial solution to the Poisson equation requires the two nets to be on different potentials [23].Three situations can generally be distinguished: 1. the two nets are only interchanged by one or more space group elements S with point symmetry part that is not the identity [50], 2. the two nets are interchanged by a primitive lattice translation [51], 3. the two nets are no symmetry copy of one another.This case is always possible for unbalanced double nets as for example the qtz-qzd double net [52].
The cubic examples discussed here all belong to the symmetric cases 1 and 2. In these cases, the square S 2 of the symmetry operation that exchanges nets evidently maps the two nets onto themselves.The potential is therefore unchanged by S 2 as the individual nets are short-cut in the quasi-static limit.As a consequence, S yields a multiplication of the potential by −1 and the two nets must hence be on opposite potential.In case 1, the Bloch character is trivial and the EAW band emanates at zero frequency at the Γ-point.In case 2, the Bloch character is −1 with respect to a primitive lattice translation and the EAW hence emanates at zero frequency at the corner of the Brillouin zone.With these preliminary considerations in mind, we first revisit the pcu-c double net.As the pcu-c evidently belongs to case 2 above, with a primitive body-centered cubic lattice translation interchanging nets, the EAW band emanates at zero frequency at the H point [42].Since H lies along the cubic 100 direction, any state at the center of the surface Brillouin zone (with vanishing lateral k) of a (100) inclinated slab is therefore made by a superposition of counter-propagating EAW modes along Γ-H.The behavior in the quasi-static limit now further yields the irreducible representation, or irrep, of the electric field with respect to the group of the wave vector [28], which is the C 4v point group in this case.Since all elements of C 4v do not interchange networks and the individual networks are on constant potential in the quasistatic limit, the EAW irrep is trivial.Therefore, the homogenized field cannot have a component perpendicular to the wave vector direction, as these would transform with a non-trivial 2D E irrep.The fields must be longitudinal from a symmetry perspective over the whole EAW and Langmuir bands, which are glued together and form one band.
The srs-c is of case 1 above, and the EAW band emanates at the Γ-point of the body-centered cubic Brillouin zone at zero frequency.The bandstructure along Γ−H (k = k e x with k ∈ [0, 2π/a]) for wire radius of 0.02a is shown in Figure A1 (a).It agrees perfectly with the HDP prediction for the corrected pressure parameter of κ≈0.675 for the wire radius of 0.02a over the whole band.This can be understood as follows: The group of the wave vector is again C 4v , which is isomorphic to the abstract group G 4  8 in [53].Considering the quasi-static potentials, the EAW band must be even with respect to the C 4 2 screw rotation and the C 2 rotation (which map each single net onto itself), and odd with respect to the glide mirror planes (which interchange networks).It therefore has the 1D representation R 2 in G 4  8 .At the H point, the EAW joins a 6-fold degenerate point, which transforms according to the R 14 representation in G 4  96 , which splits equally into the 4 1D and one 2D representations along Γ-H: R 14 (G 4 96 ) = 8 ).This allows the R 2 representation to form a pair with its time reversal symmetric R 4 representation to form a back-folding at H with finite and opposite slope of the two bands.
As in the pcu-c scenario, the longitudinal EAW behavior is protected by symmetry.Since the group of the wave vector contains a C 42 screw axis and diagonal glide mirrors that both contain a shift along the propagation direction, only the non-symmorphic wallpaper subgroup p2gg (number 8 in [54]), with isogonal C 2v point group, can be employed to make a general prediction regarding the homogenized fields.Considering the quasi-static potentials again, we obtain an even symmetry classification with respect to the C 2 symmetry in p2gg.This classification evidently does not permit homogeneous field components perpendicular to the direction of the wave vector.The odd symmetry behavior with respect to the glide mirror symmetry in p2gg on the other hand prohibits a viii homogenized field component along k.For the srs-c, we therefore expect a truly vanishing homogenized electric field that extends over the whole EAW band.This vanishing field is demonstrated in Figure A1 (b) within the numerical precision of the simulations (note the logarithmic y-axis).
Similar to the pcu-c, the dia-c network with simple cubic P n3m symmetry belongs to case 2 above.The EAW band therefore emanates at the simple cubic Rpoint [42].Consequently, the EAW band lies outside of the light cone along the cubic 100 and 110 directions and cannot lead to BICs.For a slab with (111) inclination, however, the EAW band along the path connecting R with Γ (k = ( a ]) lies at the center of the surface Brillouin zone.The low frequency dispersion agrees again perfectly with the HDP prediction with corrected κ≈0.675 as shown in Figure A2 (a).It therefore yields BICs by the same mechanism as exploited in the pcu-c case with slab inclination of (100), assuming a longitudinal homogenized electric field.Such a field is guaranteed by the trivial symmetry classification of the mode with respect to the C 3v group of the wave vector [55], which only allows a homogenized field parallel to k.This longitudinal nature is once again demonstrated through full wave simulations in Figure A2 (b).We begin with a single plasma fluid model for a warm plasma [56].This model consists of a set of three equations describing the charge carrier (CC) dynamics, Eq. (B1a) is the continuity equation of CC conservation, where n(r, t) is the CC volume density field and u(r, t) is the CC velocity field.Eq. ( B1b) is the Navier-Stokes equation for the electrically excited CC liquid, where we ignore shear forces, the non-linear magnetic field contribution in the Lorentz force density, and gravitational contributions.In this equation, ⊗ is the outer product, m the effective mass of the charge carriers, P (r, t) the CC liquid pressure field, and ρ(r, t) := qn(r, t) the CC charge density.We note that a kinetic pressure term is present in any liquid.Here, P should be understood as the total mesoscopic pressure based on the microscopic statistical ensemble of CCs, including inter-particle Coulomb interactions.Finally, Eq. ( B1c) is a local adiabatic equation with the adiabatic index γ, which can be rigorously derived from an energy balance equation in the limit where the flow velocity is much larger than vicious and thermal diffusion velocities [56].These fluid equations are closed with the Maxwell curl equations, which describe the electro-magnetic field dynamics.In Lorentz-Heaviside units, with E(r, t) and H(r, t) the electric and magnetic fields, respectively, c the speed of light in vacuum, and j(r, t) := ρ(r, t)u(r, t) the current density of the plasma charge carriers.
We now linearize Eq. (B1) assuming small variations in the density and pressure n(r, t) := n 0 + δn(r, t) P (r, t) := P 0 + δP (r, t) and thus small velocities u(r, t) to arrive at Since δP = 0 if δn = 0, Eq. (B3c) is solved by δP = γP 0 δn/n 0 .Finally, since Eq.(B2) and Eq.(B3) are now both spatio-temporally homogeneous and linear in all fields, we can make a plane-wave ansatz n(r, t) := n exp{ı(k•r−ωt)}, etc., such that the field symbols represent their constant coefficients from this point on, to obtain We now introduce the plasma frequency ω 2 p := n 0 q 2 /m, the plasma wave number k p := ω p /c, and the dimensionless pressure parameter κ 2 := γP 0 /(n 0 mc 2 ), to obtain a k-family of linear Hermitian eigenproblems for the eigenvalue ω/ω p .The eigenvector is defined as The dimensionless matrix operator which we from now on refer to as Hamiltonian, has the sub-blocks and Here, we have used the Pauli spin matrix σ y , the vector product matrix K × w := k×w, the Kronecker (tensor) product ⊗, the zero vector 0 and matrix 0 of dimension 3, and the identity matrix 1 d of dimension d; A † denotes the Hermitian adjoint of A.
The problem is of course isotropic, so that we can choose the coordinate frame such that k = k e z without loss of generality.The eigenproblem Eq. (B5) predicts a four-fold degenerate zero-frequency band ω(k) = 0.The associated eigenspace is spanned by two longitudinal spurious modes, which violate the Maxwell divergence equations, and the static magnetic field solutions generated by a transverse current density at constant charge δ n = 0 and vanishing electric field.
Note that the system is particle-hole symmetric [57], that is, for every eigenvector v of H(k) with energy ω, there is one particle-hole symmetric eigenvector v * with energy −ω.Thus we just show the positive frequency bands in Figure A3.They follow the longitudinal mode Langmuir dispersion [58] with fields E = e z , H = 0, δn = ık l q , and j = ıω l e z ; and the two-fold degenerate transverse light mode dispersion with fields E ± = e x ±ıe y , H ± = ckt ωt (e y ∓ıe x ), δn = 0, and j = 2ıω 2 p ωt E ± .

The metal pcu wire-mesh
At low frequencies in the microwave regime, most metals act like prefect electrical conductors (PECs), which do not support propagating modes in the bulk.Employing plasmonic metal wires, we can however reduce the CC density to introduce an effective plasma freqeuncy for excitations in the wire direction, while maintaining the lossless PEC character [59].A 3D metallic pcu [49] wire-mesh with sub-wavelength period therefore effectively generates a lossless plasma at microwave excitations.Based on [60], we here formally show that the electro-dynamic theory of the metallic pcu network in the long wavelength limit, is indeed equivalent to the single plasma model with a constant pressure parameter κ = 1/ √ 3. We use the macroscopic constitutive relations from [60] that have been derived for a square lattice of aligned wires (along e z ) in the homogenization picture.We start with the current I(z) and the electric potential V (z) between neighbouring wires.A comparison of the magnetic flux through a rectangle of infinitesimal height with the corresponding boundary integral (Ampére's law in integral form) yields (assuming a plane-wave form) [61] with the self-inductance per unit length of the wires in the mesh that is a simple consequence of the Biot-Savart law for the individual wires and only depends on the lattice constant (nearest neighbour distance) a and the wire radius r.The potential between neighbouring wires is on the other hand related to the linear charge density λ(z) via the linear relation λ = CV with the effective capacitance per unit length C = 1/(c 2 L) that is obtained using Gauss' law.Combining this relation with charge conservation k z I = ωλ leads  Eq. (B9) and Eq.(B10) can be readily generalized to a 3D connected pcu wire-mesh: Eq. (B9) simply acquires full vector form with k z → k, I → a 2 j , and E z → E .The linear charge density is effectively reduced by a factor of 1/3 as the charges spread over three wires in each direction in the unit cell and the pcu network is fully connected.A more rigorous treatment for wires of finite radius reveals that this factor is indeed better approximated by [62] with the plasma wavenumber In practice, the wire-mesh pressure parameter falls in a range between κ ≈ 2/3 for a wire radius of r = a/100 and κ ≈ 3/4 for r = a/5.Rearranging, we obtain Defining the plasma frequency as the cut-off frequency of the aligned wire mesh [59] k 2 p a 2 := C, which is indeed a good approximation of the definition in Eq. (B11b), and the field vector we formally arrive at the above eigenproblem Eq. (B5) with a κ that is confined in a small interval and only weakly depends on the geometrical parameters of the wire-mesh.We finally note that the scaled potential κ kp V equals the scaled CC density κ kp qδn fom Section B 1, if we associate the linear charge density along the PEC wires per unit cell area λ/a 2 with the charge density ρ of the plasma, so that the field vector v is identical in both models.Figure A4 shows the bandstructure obtained through a full-wave simulation, which agrees well with the single plasma fluid model shown in Figure A3.The two modes in Figure A4(b) correspond to the longitudinal (LMW) and transverse (EMW) modes in Figure A3.

Hydrodynamical double plasma fluid model and plasmonic double pcu wire-mesh
We now consider a plasma consisting of two CC species of respective charge q i , mass m i , and equilibrium CC density n (i) 0 and pressure P (i) 0 (i = 1, 2).The associated plasma frequencies, plasma wave numbers and pressure parameters areω pi , k pi and κ i .We further assume that CC interact only amongst their own species, which is physically of course only possible in a universe where particles of a different type do not interact (via Coulomb forces or otherwise).The theory is, however, approximately valid if the interaction is weak or if the two charge carriers are separated on a mesoscopic scale.
In the prescribed situation, the corresponding fluid equations Eq. (B4c) and Eq.(B4d) separate into two individual equations, while both current densities contribute to the Ampère-Maxwell law Eq.( B4b), As we have shown in Section B 2, a PEC pcu wiremesh acts like a single plasma in the homogenization regime, that is if the Bloch wavevector of the mode satisfies |k| 2π/a.We use the parameter vector P := (r 1 /a, r 2 /a, R/a) to define a pcu double network, where r i are the individual network radii and a the lattice constant, and R is the offset in each Cartesian direction, so that the two networks are separated by a shift vector S := (1, 1, 1) R. Assuming that the networks are well separated, the DNM fields therefore obey Eq. (B13) with the substitutions discussed in Section B 2.
The two networks generally exhibit different plasma frequencies ω pi that depend on r i and lie approximately between k p a/(2π) = 1/5 for r = 1/100 and k p (2π) = 3/5 for r = 1/5.The pressure parameters are generally also different, but contained in a very small intervall between κ = 2/3 and κ = 3/4, so that it is useful to introduce κ 1 := κ−δκ and κ 2 := κ+δκ with δκ κ.The canonical field vector of the double wire-mesh is therefore .
The eigenproblem can be written as with the associated Hamiltonian with the sub-blocks With k = ke z w.l.o.g. because of the isotropy of the problem, this double-net eigenproblem yields the longitudinal band with fields E e z , H = 0 , and j i = e z ; and the twofold degenerate transverse light band with E ± = ωt c (e x ±ıe y ) , H ± = k t (e y ∓ıe x ) , V i± = 0, and j i± = ık pi (e x ±ıe y ).Similar to Section B 1, the zero-frequency band is 6-fold degenerate, with two spuriouos longitudinal solutions and now four transverse current modes with and all other fields vanishing.
xii Importantly, a new fundamental band appears below the plasma frequency.It emanates at the Γ-point at ω = 0 with constant slope as shown in the main manuscript in Figure 1(b).The associated modes are reminiscent of acoustic waves in a conventional fluid, but with two species of charges fluctuating in opposite directions.The dispersion relation is with fields (in leading order in δκ) E = δκke z , H = 0, As we can see, the magnetic field of the EAW is strictly zero, while the electric field is longitudinal and linearly depends on both δκ and k.In the real DNM structures, since the HDP model is only approximately valid, the electric field of the low frequency mode obtains and additional longitudinal term that is zero order in δκ and quadratic in k.We demonstrate this behavior in  b), the x-component increases quadratically when moving away from the Γ-point, while the y and z components vanish for all k e x .In Figure A5 (c), there is an additional small linear dependence as δκ/κ ≈ 0.03.In any case, the EAW modes in DNMs are evidently quasi-longitudinal.

Appendix C: Effective Hamiltonian
The solution of the eigenproblem of the HDP associated with the Hamiltonian in Eq. (B14) can be obtained algebraically using any standard computer algebra tool.A manual solution can, however, be more conveniently and elegantly obtained through the application of k • p perturbation theory at the Γ-point.Note that this procedure is exact when applied to the EAW band, as we shall see below.
Let us briefly summarize the general approach [28].We first expand the eigenstates at some k = k 0 +δk (at the Γpoint we have k 0 = 0) in reciprocal space on the basis of those at k 0 , where δk has infinitesimal length.Consider the solution of the Hamiltonian at k 0 to be known for the eigenspace U belonging to some eigenvalue λ 0 := ω(k 0 )/c.We thus define the orthonormalized eigenvectors u α that span U, such that If the components of the Hamiltonian are analytical at k 0 , we obtain in first order in δk, At the perturbed point k 0 +δk, the eigenvectors v n are an element of U in zero order perturbation theory.This immediately yields the identity (using Einstein notation), As a result, the eigenvalue equation at k 0 +δk becomes Upon testing with u α , we obtain its weak form, which is an algebraic eigenvalue equation that has the dimension of U, In other words, the deviation of the dim(U) eigenvalues λ n from λ 0 is given by the eigenvalues of the effective Hamiltonian H αβ , whose entries are the matrix elements of the first order perturbation δH in U. Since H is linear in k for the HDP Hamiltonian Eq. (B15), the first order expansion is exact, that is H DN = H 0 +δH with Since we are interested in the eigenspace of H 0 at vanishing frequency, we are searching for its kernel, which is 8-dimensional and contains arbitrary H, arbitrary V 1 or V 2 , and arbitrary j 1 = −j 2 , with all other fields vanishing.Labelling the field components in u α with i and disregarding the trivial electrostatic solution with H = 0, we thus obtain the 5-dimensional eigenspace matrix with orthogonormalized eigenvectors as rows The effective Hamiltonian is thus where Three of the eigenvalues of the perturbation Hamiltonian in Eq. (C2) are zero.In addition we reproduce the EAW band in Eq. (B19) and its chirally symmetry counterpart: Evidently, the EAW slope only weakly depends on the difference in plasma frequencies of the two networks in first order δκ.Particularly, if the two networks have note that the dispersion of EAW band in the low frequency limit is robust against changes in the DNM parameters such as the radius of the metal wires and the relative offset between the two nets.This assumption is verified by simulating DNMs with various parameters, shown in Figure A6 and discussed in Section C. (h) Bandstructure for the slab shown in (g).As theoretically predicted, 9 BIC modes are found wires oriented perpendicularly to the propagation direction, which we refer to as x w.l.o.g., at finite frequencies.The pcu-c DNM with body centered cubic symmetry has a global C 4v symmetry through the center of each of the wires.Since the currents in the x wires are constant along the wire surface, they transform trivially under the corresponding C 4v symmetry.As a consequence, the weak EAW currents in the wires in the y-z plane must transform trivially as well, and therefore either all flow towards the nodes or away from them in a 4-fold quadrupole-like fashion.This expected behavior is demonstrated through the surface current densities obtained from full wave simulations, whose y component is shown in Figure A8 (a) and the z component in Figure A8 (d).
Notably, the in-plane geometry of the individual nets dominates the electromagnetic distribution, since the currents on the two nets experience very weak intercoupling at low frequencies.The local C 4v symmetry of each single net thus guarantees the formation of the aforementioned quadrupole-like current distributions, even in the absence of a global C 4v symmetry.We demonstrate this behavior through the current distributions shown in Figure A8(b) and (e) for DNMs with a non-symmetric shift, and in Figure A8(c) and (f) for distinct wire radii.To a very good approximation, these fields retain their trivial character with respect to the individual net's C 4v symmetry.Importantly, this trivial character is evidently incompatible with the free space photon modes, which  carry the non-trivial character of the E ± (right and left circular polarization) irreducible representation [63].
Appendix F: BICs in the DNM immune to global symmetry breaking The in-plane cross-shaped geometry which intersects perpendicularly by identical metal pillar of each individual wire mesh keeps the C 4v symmetry.The cross-shaped metallic pillar guide the surface density current, forming two quadruples in the space.The coupling between the two electrical quadruples is weak which is also the mechanism causing the BICs modes in DNM slab to be immune to the global symmetry.In the main text, we show the DNM slab with the parameters R(0.4,0.4) and offset (a/3, a/3,a/3).To further illustrate robustness of BICs in DNM, we show numerical results of DNM slabs with various parameters in Figure A9, all of which host BIC states.

Appendix G: Lifting BICs by breaking local symmetry
Throughout this work, we have considered the metals constituting the DNMs as perfect electric conductors, so that surface currents flow on the metallic wires and determine the microscopic electromagnetic modes.We have shown in Section F that the surface current density forms two quadrupole-like distributions, which are robust against global symmetry breaking.These quadrupole currents are in the 1D representation of the C 4v point group labelled B 1 in [28], and thus cannot couple to vacuum radiation, which is in the 2D E representation.However, the BICs can be destroyed by breaking the local C 4v symmetry of the individual nets in the DNM.We demonstrate this by introducing an thick metal cylinder segment on the lateral wires along the x and y directions normal to the slab direction z (Figure A10(a)).These additional cylinders are arranged in such a way that the C 4v symmetry is reduced to a single mirror plane on the x−y diagonal.For simplicity, we leave the other net unperturbed.The broken local symmetry, reduced to merely a mirror group, enables couplings of the EAW mode to vacuum modes.When compared to the DNM with two identical unperturbed wire meshes, the slab formed by the perturbed DNM (Figure A10(c)) degrades the BIC modes at the Γ-point to quasi-BICs with finite Q-factors (Figure A10(d)).It is, however, worth noting that the perturbed DNM slab can still be well described by the HDP theory, as illustrated by the bulk bandstructure along the out-of-plane z direction in Figure A10(b), which retains the typical EAW band.
The nature of the quasi-BICs, on the other hand, necessarily imposes changes in polarization structure in the far field.Compared with the 'star' polarization pattern of the BICs in the unperturbed structure (Figure A11(a)), xvii the perturbed DNM slab exhibits a richer polarization distribution, forming two patches of elliptically polarized states with opposite handedness, separated by a diagonal linear polarization line (Figure A11(c)).The position of this line on the diagonal is evident from the fact that the in-plane components of both the dominating longitudinal currents and the quadrupole currents are even with respect to the remaining mirror if the mirror coincides with the plane of incidence, that is, if the wave vector is invariant.Thus, the far-field must obey the same symmetry classification, that is, it must be p polarised, with the electric field confined to the plane of incidence.In all other directions, the EAW field is a combination of even (quadrupole currents) and odd (in-plane longitudinal currents) contributions, and thus the far-field is generally elliptically polarized.
The mirror plane divides the x−y plane into two halfplanes of opposite elliptical polarizations, which is immediately evident from the action of mirror symmetry σ d on the modes on one half-plane.The vortex at the origin splits into two circular polarization points on either side.This mechanism is similar to the one reported in [20]; it is of topological origin.The appearance of the circular polarization points becomes evident from tracing the polarization on various momentum space circles of different raddi centered around the Γ-point (Figure A11(c)).The polarization path projected on the Poincaré sphere goes twice through the linear polarization state on the mirror plane and unfolds from a small figure-of-eight (turquoise, small radius) towards a double-loop around the equator (grey, large radius), and thus necessarily passes through the north/south poles (Figure A11(d)).
Figure A11 further indicates that the circular polarized points lie on the axis perpendicular to the remaining mirror plane.We can understand this behavior from a symmetry perspective.Consider the action of the timereversal operator τ , which transforms the slab mode with wave vector k into one with −k and complex conjugated field, that is τ E k = E * −k .The mirror symmetry on the other hand interchanges the x and y components of the field and also maps k to its negative.The combined operation τ σ d thus leaves k invariant on the line perpendicular to the mirror plane.We thus obtain E x != E * y , so that the normalized in-plane electric far-field depends on one free parameter ϕ only, e ıϕ e −ıϕ .
In other words, the field stays on a line of constant longitude on the Poincaré sphere, connecting zero linear polarization along the mirror plane (at Γ) and perpendicular to it (away from Γ for small perturbations).It must therefore pass the circular polarized poles (ϕ = π/4) on either side of the mirror plane.We here devise an analytical model that reveals the connection between the EAW bulk picture and the slab modes.It generates a good prediction for the quasi-BIC modes and the transmission spectra.
Let us first solve the scattering problem using two counter-propagating EAW slab modes only.We first define the two vacuum domains I and III, semi-infinite in z-direction and infinite in x and y, and the DNM slab domain II.As the electromagnetic fields vanish in the HDP model, we need to extract fields at the domain interfaces from the simulations.The EAW fields in the DNM exhibit longitudinal electric fields and vanishing magnetic fields, as we have already seen in Figure A5.To solve the homogenized Maxwell equations (which is equivalent to cutting the Rayleigh series at the fundamental Bragg order) at the slab surface, however, we need to average the fields over the this surface rather than the whole unit cell [34].Considering an in-plane wave vector q along the x direction, we obtain a homogenized electric field pointing along x, and a magnetic field pointing along y.This behaviour is consistent with the symmetry classification for the mirror symmetric DNM, and means that the EAWs only couple to p polarised waves in the vacuum domains.One further obtains that the surface impedance Z := E x /H y is largely independent of the frequency, and, consistent with optical reciprocity, quadratically depending on q in good approximation.For an EAW wave propagating in positive z direction, we obtain The vacuum impedance in contrast is of course frequency dependent and finite at vanishing q, with for a wave propagating in positive z direction.We assume a slab of thickness d = N a and the wave numbers in propagation direction are in vacuum k 1 := (ω/c) 2 − q 2 , and in the slab k 2 := [ω/(κc)] 2 − q 2 , using Eq.(B19) for equal wire thickness of 0.08a (κ ≈ 0.72).The non-vanishing field components for the scattering problem can thus be expressed as E I = Z 1 e ık1(z+d/2) − re −ık1(z+d/2) e ıqx , (H1b) The homogenized Maxwell equations are satisfied if the fields match at both interfaces at z = ±d/2, yielding a standard scattering problem Here, the transmission is overestimated significantly, as the surface field homogenization breaks down, particularly at large frequencies.The same fundamental idea can be employed to calculate the quasi-BIC bands.For this, we consider only outgoing waves in the vacuum domains, which transport energy away from the slab, that is, with {k 1 }<0 [>0] in region I [III].We can characterize the quasi-normal modes as even or odd regarding the z → −z mirror symmetry of the homogenized slab, and therefore substantially reduce the complexity of the problem.In particular, one only needs to consider the fields in one vacuum domain and the even or odd superposition of EAW modes in the slab.The fields are thus expressed by H II = e ık2z ± e −ık2z e ıqx , (H3c) For convenience, we have here normalized the eigensolutions such that the positive EAW wave in the slab carries a coefficient of 1.The matching conditions at the slab surface thus yield At the Γ-point (q = 0), we immediately obtain A = 0 from Eq. (H4b).This implies a modal solution with vanishing homogenized field outside of the slab and thus infinite quality factor, in other words BICs.Eq. (H4a) further yields p 2 2 = ∓ 1, which reproduces the BIC frequencies obtained from the HDP model under hard wall boundary conditions, that is, For arbitrary q, Eq. (H4a) and Eq.(H4b) yield a finite A and are transcendental in the complex frequency.They therefore need to be solved numerically.We divide Eq. (H4b) by Z 1 =0 and add it to Eq. (H4a) to obtain Since f (ω) is holomorphic, except when crossing the light line, we employ a standard Newton method using ω n as starting guess for q n+1 = q n +δq (n ∈ N) to obtain the quasi-BIC bands and corresponding quality factors Q := {ω}/(2 {ω}) in Figure A13.As expected, the bandstructure and quality factors predict the position and width of the transmission bands in Figure A12, respectively.

Appendix I: Negative refraction in DNMs with BCC symmetry
The DNM modes are naturally Bloch waves due to the discrete periodicity of the geometry.Crystal symmetries  in the DNM play an important role and generate additional physics that cannot be predicted by a homogenization theory such as the HDP model.The crystal symmetry can for example lead to band degeneracies at high symmetry points in the Brillouin zone, where the wave vector is mapped onto itself or a reciprocal lattice equivalent under point group operations.If a non-primitive unit cell is chosen, artificial band folding, or better a band translation into the artificially decreased Brillouin zone, occurs.This can have a significant effect on the physical interpretation of the modes as this back-translation can ostensibly change the sign and magnitude of the phase velocity, while keeping the slope of the band and thus the group velocity of the mode unchanged.One might therefore expect left-handed propagation where the mode is right-handed and vice versa.
Consider a DNM with a body-centered cubic (BCC) symmetry, which is also called the pcu-c wire-mesh with an offset of a/2/ (1, 1, 1) between the two pcu networks.We here set the radius of each individual net to 0.08 a.A primitive cell, a parallelepiped spanned by the BCC basis vectors, is illustrated in Figure A14(a).The canonical Brillouin Zone, a rhombic dodecahedron, is shown in Figure A14(b).We emphasize the two high symmetry points, invariant under all octahedral (O h ) symmetries, Γ = (0, 0, 0) and H = (2π/a, 0, 0), connected along the cu- bic [100] direction.In the HDP model, we derived a positive slope of EAW mode, emanating from Γ in the long wavelength limit, consistent with the picture in the nonprimitive simple cubic Brillouin zone (Section 2 in the main text).However, the band really emanates at the H point, which is mapped to the Γ point in the simple cubic Brillouin zone.The correct BCC symmetry classification, with the corresponding bandstructure shown in Figure A14(c), thus reveals the left-handed nature of the EAW modes and predicts negative refraction for the [100] inclinated slabs in question.We note that the lefthanded nature of the EAW modes persists in perturbed DNMs of lower simple cubic symmetry with different network offsets or wire radii, even though the bandstructure does not reveal it in these cases.The underlying reason is rooted in the fact that the dominating Bragg order of the Bloch wave remains outside of the first Brillouin zone close to (±2π/a, 0, 0) for these perturbed structures.We further note that the weak coupling of the quasilongitudinal EAW modes to vacuum radiation makes it hard for them to be observed experimentally.Topological surface bound states emerge in the bulk bandgap due to the nontrivial topological index and a decoupling from vacuum radiation.In Figure A15 (b), we show the TSBIC frequencies on the Γ-point as a function of slab thickness, for integer numbers of unit cells N , for which the DNMs keep the geometrically equivalent boundary at the top and bottom surfaces.For small slab thicknesses, the two surface states on the top and the bottom of the slab are not spatially separated, leading to a frequency splitting of the two hybridized modes, which are even/odd with respect to the mirror symmetry at the center of the slab.Evidently larger thickness of the the DNM slab cause weaker coupling between the two surface modes, hence a reduced energy splitting that converges to zero.In Figure A16, we show the dominating surface currents J sz of the even/odd hybridized TSBIC states with two different thicknesses of N = 5 and 10.

FIG. 1 . 1 √ 3
FIG. 1.(a) Schematic diagram of an EAW mode in a HDP fluid.The blue and yellow counter-propagating sinusoids arise from two separate charge carriers, resulting in a vanishing field.The DNM consists of two interwoven percolating metallic nets, supporting the two counter-propagating currents.In the electrostatic limit, the DNM is equivalent to an HDP, resulting in a homogenized electromagnetic field that vanishes at the Γ point.(b) Dispersion relation of an HDP with κ = 1 √ 3 on the left, featuring a light-like electromagnetic wave (cyan), a Langmuir wave (black), and an electron-acoustic wave (red).The corresponding simulated band structure of an unperturbed DNM is shown as dots on the right.The three modes resemble those of the HDP.The slope of the fundamental band is in excellent agreement with the EAW band of HDP model (red dashed line).

FIG. 3 .
FIG. 3. (a) Sketch of a DNM slab embedded in vacuum.(Quasi-) BICs on (emerging from) the Γ-point do not (weakly) radiate energy.(b) Quasi-normal mode bandstructure of an N = 5 unit-cell thick P1 DNM slab.The black quasi-BIC bands emanate from the red BIC points at Γ and consist of two EAW bulk modes.The two green bands correspond to spoof-plasmon surface modes.The region outside of the light cone is shaded in grey.(c) The Q-factors of the four quasi-BIC modes are large at finite k and diverge at Γ.(d) Transmittance of an incoming plane wave within the light cone.The transmission is 0, except close to the quasi-BIC bands, on which it is 1.(e) Linear polarization direction of the far-field in k-space.(f) Frequency of the BIC modes for N unit cell thick slabs.The simulated frequencies (black dots) match the analytical prediction from the HDP model (blue) well.

FIG. 4 .
FIG. 4. (a) Geometry of a P2(4) DNM consisting of a thick (blue) and a thin (yellow) wire network.(b) Projected surface bands at Γ as a function of ψ = r1/r2.A topological phase transition occurs for ψ = 1, where the bandgap closes.(c) Quasi-normal slab modes as in Figure 3(b).Additional to the black quasi-BICs and the green surface bands, two TS-BIC bands are found inside the projected bulk gap (red).(d) Transmittance of an incoming plane wave within the light cone as in Figure 3(c).The quasi-TSBICs give rise to high-Q transmission similar to the bulk BIC bands.(e,f) P2(4) EAW electric field component Ex across the (100) plane in the center of the unit cell for kx = ky = 0, and (e) kz = 0.05 π a and (f) kz = π a .The parity change, even in (e) and odd in (f), yields a Zak phase of π for the EAW band.
FIG. A1.EAW mode characterization for the srs-c double net along the Γ-H path (parametrized by k = kx ex).(a) Dispersion relation compared to the HDP prediction with limiting (1/ √ 3) and radius corrected κ = 0.675.(b) Averaged field components normalized by averaged field intensity parallel to k ([100]) and along perpendicular high symmetry directions.

Appendix B: Plasma and pcu network modes 1 .
Single plasma fluid model FIG. A3.Dispersion relation of a single plasma fluid with thermal coefficient κ = 0.4.It has a longitudinal Langmuir (LMW) and a doubly-degenerate electromagnetic wave (EMW) band.
tok z I = ωCV .(B10) FIG. A5.Absolute value of the electric field, averaged over one unit cell, for each component of the EAW bands along the high symmetry path Γ−X along the x direction.The parameter choices are (a) P2(1);(b) P1; and (c) P2(4).
FIG. A9.(a) Geometry of a DNM slab consisting of two identical nets with P = (2/25, 2/25, 1/2).The slab thickness is 2a.(b) Bandstructure of the slab shown in (a).The black band emanating at zero frequency is made by an in-plane EAW mode, while one quasi-BIC band emanates from the BIC mode indicated by a red star.The green bands are spoof plasmon surface modes confined to the slab surface.(c) Slab with tickness 2a and P = (2/25, 2/25, 1/3).(d) Bandstructure for the slab shown in (c).The BIC mode and the associated quasi-BIC band are conserved under the symmetry break.(e) Slab with thickness 2a and P = (2/25, 1/50, 1/3).(f) Bandstructure for the slab shown in (e).(g) Slab with thickness 10a and P = (2/25, 2/25, 1/3).(h) Bandstructure for the slab shown in (g).As theoretically predicted, 9 BIC modes are found xvi FIG. A10.(a) Schematic diagram of a DMN with broken local C4v symmetry; the radius of the added thick metal cylinder is 0.12a, while the radius of the unperturbed net is 0.02a.(b) Corresponding dispersion relation of (a); the black dots are from a full wave simulation, the blue line is the 1/ √ 3 slope.(c) DMN slab with (left) and without (right) breaking the local symmetry.(d) Q-factor of the nets shown in (c).The Q-factor diverges at the Γ-point for the unperturbed network (red line), which is indicative of the existence of BICs.In contrast, the Q-factor stays finite for the network with strong local symmetry breaking (blue), showing quasi-BIC behavior only.
Appendix H: Analytical quasi-BICs and transmission spectra FIG. A12.Analytical transmissivity through an N = 5 DNM slab using the HDP dispersion Eq. (B19) and a homogenized DNM impedance.

H2) with p 2 :
= exp{−ık 2 d/2}.The resulting transmissivity T = |t| 2 is shown in Figure A12.It agrees well with the numerical full wave results in Fig. 3 (d) in the main manuscript, except for a small region close to the light line.
FIG. A13.Bandstructure and quality factor of the quasi-BIC bands for a N = 5 DNM slab using the HDP dispersion Eq. (B19) and a homogenized DNM impedance. xix FIG. A14.(a) A primitive unit cell of the pcu-c DNM with BCC symmetry.The radius of the two nets is r1 = r2 = 0.08a.(b) BCC Brillouin zone.(c) Band diagram of the DNM in (a).The black dots are calculated by full wave simulations, the blue solid line is the 1/ √ 3 HDP predication.The dashed line illustrates the band folding into the simple cubic Brillouin zone.
FIG. A16.Jsz plots of the hybridized modes formed by two TSBICs.Even (a) and odd (b) modes for a slab thickness of N = 5; odd (c) and even (d) modes for a thickness of N = 10.
Appendix J: Topological surface bound states in the continuum In Figure 4(b) of the main text, a topological phase transition is observed that leads to topological surface bound states in the continuum (TSBICs), with the cor-responding TSBIC bands shown in Figure 4(c) and (d) in the main text.In the topologically non-trivial region, the distinct parity of the field on the two high symmetry points X (Figure 4(e)) and Γ (Figure 4(f)) leads to a π Zak phase in the EAW band.The corresponding bulk band structure is shown in Figure A15 (a).
FIG. A17.Reflection amplitude (dashed blue) and phase (red) of a symmetric DNM slab with TSBICs on both top and bottom surfaces (a), and a slab with TSBICs only on the reflection side (b).(a) Critical coupling with vanishing reflection and a π phase jump can be observed at the resonance frequency.(b) The reflection is uniformly 1 across the resonance due to the lack of a transmission channel, while the phase shows a 2π phase jump indicating an overcoupling regime.
FIG. A18.(a) Geometry of DNM slab with P1 parameter; (b) Geometry of DNM slab with different network radii, giving rise to topological surface states.