As early as 1929, von Neumann and Wigner constructed spatially localized bound states in the continuum (BICs) as solutions of the single particle Schrödinger equation with energies above the associated potential.1^,2 BICs are, however, general wave phenomenona and are therefore not restricted to quantum mechanics and have recently been demonstrated in a variety of classical wave systems.3^–14 Among these different demonstrations, light–matter interactions endowed by BICs in periodic systems have recently gained strong attention and have been exploited for applications in biosensing,15^,16 vortex beam generation,17^,18 high-Q resonators, and lasing.9^,19^–22 To date, several mechanisms for generating BICs in the periodic systems have been reported, including (i) symmetry-guaranteed mode mismatching, (ii) destructive interference with a topological origin,3 and (iii) the Friedrich–Wintgen mechanism.23^–25 Recent evidence, however, suggests that currently available BICs provide unreliable lifetimes. For example, mechanism (i) yields BICs that are readily destroyed by small symmetry-breaking perturbations20^,26^,27 and couple to higher diffraction orders.28 Such features lead to quasi-BIC states with manageable lifetimes. On the other hand, mechanism (ii) gives rise to BICs that are easily shifted in momentum space (excitation angle). Furthermore, non-BIC modes typically exist spectrally close to these BICs,4 which cannot easily be eliminated through existing design principles. BICs created by mechanism (iii) are, for example, realized by the destructive interference between two optical modes satisfying the Friedrich–Wintgen condition and are, therefore, always accompanied by a constructive interference background mode.

Heuristically, an alternative approach to generating photonic BICs is to engineer a material that supports a pure charge wave with a vanishing electromagnetic (EM) field. While such a wave cannot exist in all-dielectric materials, it is naturally obtained in a non-Maxwellian hydrodynamic double plasma (HDP). We show that such an HDP can be approximately realized through double-net metamaterials (DNMs). Metamaterials are artificial structures, typically made of periodically arranged elements, so-called meta-atoms, on a lattice significantly smaller than the wavelength of interest.29 These structures locally mimic the quasi-homogeneous light–matter interaction of ordinary materials but give access to optical properties that cannot be found in nature, including negative refraction,30^,31 near-zero refractive index,32^,33 and strong chiro-optical effects.34^–36 Fully connected metallic network-like metamaterials mimic the EM response of bulk metals. In contrast to bulk metals, whose plasma frequency is typically in the ultraviolet, the effective plasma frequency of metallic network-like metamaterials can be tuned to a specific wavelength range of interest.37^,38 Network metamaterials can be manufactured for visible wavelength operation via self-assembly, where they have been associated with several interesting effects, such as linear and circular dichroism.36^,39^,40 DNMs, also referred to as interlaced wire media, consist of two intertwining, spatially separated metallic network domains and offer the generation of a fundamentally new physics. DNMs have only recently gained theoretical interest in the context of negative index ellipsoids,41 a light tunneling anomaly,42 and longitudinal mode propagation.43 An experimental realization of an optical gyroid DNM dates back to 2011,44 while a microwave DNM has been reported recently.45

Here, we propose a generic effective medium approach to generate a new type of HDP photonic BICs in DNMs. In contrast to existing designs, the proposed BICs are formed by pure charge fluctuations with a vanishing EM field, also known as electron acoustic wave (EAW) in HDPs. These BICs are spectrally isolated from normal resonances, robust against a large class of perturbations. Most surprisingly, they retain an ultrahigh Q-factor within a subwavelength DNM microcavity resonator, even below the diffraction limit. They exist within a broad spectral range, free of non-BIC resonances and their frequency spacing can be freely tuned by the size of the structure, reaching spectrally dense BICs for an optically thick DNM slab in the absence of non-BIC modes.

To better understand the formation of the new type of BIC modes in DNMs, we first derive and discuss the nature of the bulk states, particularly the new EAWs, in the context of metamaterials. We start with the HDP model. The HDP is formally described by (i) particle conservation to describe the individual charge carrier dynamics; (ii) the hydrodynamic Navier–Stokes equations and the Lorentz force to model particle–particle and light–particle interactions; and (iii) the Maxwell equations, which describe the light-field dynamics. A detailed derivation of the HDP model can be found in Note S1 in the Supplementary Material.

For convenience, we restrict the discussion to a nondissipative model and therefore, to low frequencies, where the constituting metals are well approximated by a perfect electric conductor. We further introduce the plasma frequency ωpi≔n0iqi2/mi and the nondimensional thermal pressure parameter κi≔γpi/(n0imic2) for the two charge carriers labeled i=1,2, which remain the only free parameters in the model. These two parameters depend on the equilibrium charge carrier density n0i, the charge qi, the effective mass mi, the equilibrium pressure pi, the adiabatic exponent γ, and the speed of light in vacuum c. Through a plane-wave approach and by linearizing the HDP model,46^,47 the generally complex problem is transformed into a 14-dimensional linear eigenproblem for the eigenfrequency ω (see Note S1 in the Supplementary Material for a detailed derivation). The components of the eigenvector in this low-dimensional problem with eigenfrequency ω are the electric and magnetic field amplitudes E and H and the amplitude of the charge carrier density fluctuations δni and the corresponding current density ji for each plasma (labeled i=1,2).

The band structure of the linearized HDP, obtained by solving the eigenvalue problem, is shown on the left in Fig. 1(b) (note that we only show the physically relevant positive frequencies). From high to low frequencies, we obtain the two-fold degenerate EM transverse band (EMW) (cyan), the longitudinal Langmuir band (LMW) (black), both of which are hyperbolic, emerging from ω=(ωp12+ωp22)1/2 for wavenumber k=0, and the linear EAW band emanating from ω=0. The application of k·p theory48 to the Γ-point to construct an effective Hamiltonian for the EAW band predicts a slope of ωa2=(m2+κ22κ12)1m2+1κ12c2ka2, where m=ωp1ωp2 (Note S1.3 in the Supplementary Material). This reduces to the approximated form ωa=cκ1ka for equal plasma parameters. For equal pressure parameters κ1=κ2, the EAW modes consist of oppositely drifting charge waves with δn1=−δn2 and j1=−j2∥k and vanishing EM fields E=H=0 [Fig. 1(a)]. Due to this zero EM nature and the resulting decoupling from the EM vacuum modes, the EAW naturally forms BICs in the presence of a photonic radiation continuum. In a natural HDP, the two charge carrier species of the plasma, usually two frictionless electrons and ions, interact so that the model described above only applies in a limit, and extremely strict conditions are necessary to construct an observable EAW system.49 Even though the EAW may play an important role in space physics, due to broadband electrostatic noise observed in the Earth’s auroral region, there has been little experimental work on this mode in the laboratory.50^–52

Concurrently, single 3D metallic network metamaterials (Note S1.2 in the Supplementary Material) have been shown to act as an effective realization of a single species plasma with ultralow effective plasma frequency in the gigahertz range.53 An effective hydrodynamic behavior beyond the standard constitutive equations that leads to a nonlocal response can be understood through a charge accumulation and release effect.37^,38 To circumvent the difficulties in creating a natural double-plasma,49^–51 we employ a 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. The double-plasma [Eq. (S14) in the Supplementary Material] and the plasmonic network eigenproblem [Eq. (S20) in the Supplementary Material] are indeed formally identical, as we have shown in Note S1 in the Supplementary Material. Note that the effective thermal pressures are mainly determined by the connectivity of the wires, leading to κ1,2≈13, due to the vertically intersecting morphology in the three axes. For finite radii, κ1,2 are adjusted through the radius of the networks,54 giving more accurate values of the EAW band’s slope (see Note S1.4 in the Supplementary Material for details and Note S3 in the Supplementary Material for an illustration of this effect). 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. A quasi-static circuitry theory of multinet metamaterials has been proposed to predict the low-frequency dispersion,55 but not the EAW nature. The DNM structures discussed in this paper are composed of two single nets of fully connected metal cylinders along the cubic axes with radii r1,2, which may differ between the two nets. We use the notation P≔(r1,r2,d) to represent the network parameters. d is the offset of the two networks that are shifted by S=(1,1,1)⊺d. The double network with d=a2 is the pcu-c net with Im3¯m symmetry. Throughout this paper, the lattice constant is a=5mm.

The numerically simulated band structure of a DNM with P=(0.4,0.4,a3) is illustrated in Fig. 1(b) and features a fundamental band with linear slope (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 approximation from Note S1 in the Supplementary Material in the limit of thin wires ω=13ck (red dashed line). As expected from the HDP model, this slope mainly depends on κ and only weakly depends on the quotient of the effective plasma frequencies. Our theory is valid within a large variety of parameters of DNMs, confirmed by numerical simulations for various values of P (Note S2 and Fig. S3 in the Supplementary Material). For larger wave vectors closer to the boundary of the Brillouin zone, the long wavelength limit breaks down, and the slope deviates and becomes zero at the X-point [k=(π/a,0,0)] due to time-reversal symmetry.48 Analogously to the EAW in the HDP model, the two networks in the DNM carry opposing currents, indicated by arrows in Fig. 1(a). The simulated surface current densities J for DNMs with equivalent radii and different offsets [P=(0.4,0.4,a3) and P=(0.4,0.4,a2)] are shown in Fig. 2(a). They indeed exhibit almost exactly opposite Jx in the two networks at a small k (kx=0.05π/a and ky,z=0). For different radii [P=(0.4,0.1,a3)], 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 microscopic EM fields of the EAW-like modes in DNMs, we compute their 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 (Figs. S4 and S5 in the Supplementary Material), while the longitudinal component quadratically increases with the length of the wave vector [Fig. 2(a)]. Therefore, the fundamental modes in DNMs are quasi-longitudinal with E∥k. We conclude that the HDP model is a valid effective medium description for DNMs in the static limit. DNMs provide a metamaterial realization of HDPs, which require extremely stringent restrictions in natural plasmas for their laboratory observation. We finally note that the notion of an effective permittivity seems misleading due to the non-Maxwellian nature of the EAW mode, meaning the lack of EM fields. In this context, a previously reported homogenization approach56 leading to an effective permittivity εh for the DNM only accurately describes the Langmuir mode residing above the effective plasma frequencies. In contrast with standard metamaterial homogenization approaches, the HDP model developed here yields a fully consistent, effective medium description of the fundamental DNM band and accurately predicts both the phase velocity and the associated fields. Even though the discussion in this paper is restricted to Cartesian wire meshes as shown Fig. 1, we show that the HDP model applies equally to other DNM geometries in Note S3 in the Supplementary Material (Figs. S4 and S5 in the Supplementary Material).

Our DNM theory, presented in Sec. 2, reveals a significant finding concerning a spectrally isolated fundamental band that supports quasi-longitudinal EAW-like bulk modes. Unlike those found in single-network metamaterials or meta-atom-based systems, these modes possess a remarkable characteristic: the vanishing of their EM fields. This distinctive property leads us to observe a new type of BIC that naturally arises from the EAW-like behavior of the bulk modes within a vacuum slab configuration. Figures 2(b)–2(e) visually highlight this relationship between the EAW modes and the emergence of these BICs.

In the HDP picture, a standing wave slab solution exists at frequencies where the EAW solutions can be superposed such that the normal component of ji vanishes at the slab boundaries. As the EAW fields are perpendicular to the vacuum plane waves, they perfectly decouple from radiative modes at the Γ-point of the surface Brillouin zone.57 At the Γ-point, we thus obtain zero radiation and infinite lifetime [or Q-factor, here we define Q=R(f)2I(f)] of the slab modes, which are therefore BICs in the center of the light cone. Note that the full nonapproximate scattering problem, while more complex and involving evanescent Floquet modes,58 yields the same qualitative result, since the higher Bragg orders in the vacuum Rayleigh basis59 are evanescent and nonradiative as the EAW band is below the first Wood anomaly; see Fig. 2(a).

To illustrate the above prediction in practice, we simulate a DNM slab with P=(0.4,0.4,a2) and N=5 unit cell thickness. The geometry is shown in Fig. 2(b). We apply Bloch-periodic boundary conditions in the x-y plane and scattering boundary conditions (numerical Sommerfeld conditions60) in the z direction and solve for the complex-valued frequencies (that is, we compute the quasi-normal slab modes). As predicted, multiple BICs with infinite lifetimes are located at the Γ-point and marked by red pentagrams in Fig. 2(c). The frequencies of the BICs can be well approximated by the standing wave solution of the EAW with the predicted phase velocity of 13c and hard wall boundary conditions (ji=0). Since the EAW bulk band is spectrally isolated from the other nonevanescent bulk bands at low frequencies, the BICs are fundamentally prohibited from interfering with any non-BIC slab modes in the whole spectral range of the EAW band, in sharp contrast to BICs in photonic crystals.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 Q-factor. The corresponding Q-factor values of the multiple quasi-BIC solutions near the Γ-point are shown in Figs. 2(c) and 2(d). The physics is reminiscent of a Fabry–Pérot resonator with mirrors that are perfectly reflective at normal incidence and are otherwise weakly transparent, but without using any mirrors. This physics yields a simple signature in scattering experiments: the transmittance resembles a Lorentzian line shape at a finite wave vector close to the BIC frequency. The center frequency forms a band that emanates from the BICs at the Γ-point and hyperbolically blueshifts with increasing wave vectors while the width of the Lorentzian gradually broadens. The simulated reflectance spectrum along the Γ-X direction is shown in Fig. 2(d) and agrees well with our theoretical modeling (see Note S4 and Figs. S6 and S7 in the Supplementary Material). The reflection spectrum, accounting for metal losses, is provided in Note S7 in the Supplementary Material.

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 Fig. S9 in the Supplementary Material (see Note S5 in the Supplementary Material). This generates a singular point at the Γ-point where the BICs reside, where the polarization cannot be defined because the far-field vanishes. Around this point, the polarization pattern evidently exhibits a nontrivial 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. It is reminiscent of the electric field lines created by a point charge.

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, assuming the HDP in a hard wall cavity, is given in Note S4 in the Supplementary Material. The theoretical model predicts the quasi-BIC bands, including their Q-factors (Fig. S7 in the Supplementary Material) and the transmission spectrum (Fig. S6 in the Supplementary Material). In particular, the model yields BIC modes that are equidistantly spaced at frequencies of ωac=lκπN, where Na is the thickness of the DNM slab and l∈N. The frequencies at which BIC modes occur on the Γ-point, and thus also the number of BIC modes that exist over the spectral range of the EAW band, are a function of the number of unit cells in the slab N (which can be a noninteger number). 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=(0.4,0.4,a2) DNM in Fig. 2(e). Generally, the theoretical results agree with the simulations, even though the equidistant frequency spacing is perturbed in real DNMs.

BICs residing at high symmetry points in the Brillouin zone in photonic crystal slabs and metasurfaces are typically protected by in-plane inversion (C2) symmetry but can be disrupted by infinitesimal symmetry-breaking perturbations, as well as couplings to high diffraction orders.26^,27 Interestingly, despite being susceptible to symmetry breaking, the newly formed BICs in DNMs show exceptional robustness, retaining ultrahigh Q-factors. Our findings indicate that these new BICs persist even in the presence of symmetry breaking, disorder, and scattering into high diffraction orders. The exceptionally high robustness is linked to the HDP-effective medium nature of the DNMs.

To elucidate, we first demonstrate variations to the BICs imposed by local and periodic symmetry breaking. Two representative schemes were studied. (i) Deviations of the nets from their centrosymmetric positions [bottom left corner of Fig. 3(a)]. For example, a DNM with P=(0.4,0.1,R) has only one in-plane mirror symmetry normal to the [11¯0] direction, which is insufficient for creating symmetry-protected BICs.8 (ii) Breaking the C4v symmetry of a single net, while keeping the global mirror symmetry along [11¯0], by adding thicker metallic cylinders between the nodes, as shown in the top-right corner of Fig. 3(a). In both cases, the Q-factor obeys an approximately inverse quadratic relationship Q(α)=Q0α−2 wherein α is the asymmetry parameter,20 illustrated in Fig. 3(a). Q0 is a constant determined by the specific network design.

Both cases possessing Q0 higher than 2×104 were obtained, leading to high Q-factors even with considerable α-values. For example, Q=106 for a relatively large α=10−1. Note that the value of Q0 alone, without consistently normalizing α, is arbitrary. However, if the perturbation is geometrical, α can be normalized by typical feature sizes (object dimension or unit cell size), while α is internally normalized if the perturbation is angular. The Q factor is typically higher than 106, even when structural changes occur at a scale comparable to the unit-cell size. This showcases the remarkable resilience of the novel BIC modes in DNMs. This excellent stability of the DNM-based BICs compared to previous findings20 seems to stem from the macroscopic HDP effective medium nature of the DNMs. Microscopically, the surface current densities J on the metallic wires that are orthogonal to the propagation direction resemble an electrical quadruple that is symmetry-incompatible with free-space radiation. In the static limit, interactions between the currents on the two nets are negligible, so the total EM radiation vanishes irrespective of the networks’ relative position, radii, etc. Finally, the symmetry-breaking schemes mentioned above, (i) and (ii), give rise to a linear radial vortex polarization around the Γ-point converting into linear polarization on the diagonal in wave-vector space, separating two regions with right- and left-turning elliptical polarization, respectively (Note S5 in the Supplementary Material). There is one spot on either side in wave-vector space where the polarization of the quasi-BICs is circularly polarized so that perfect circular dichroism in reflectance can be achieved with oblique angle of incidence at the circular polarization points at the quasi-BIC frequency.

Due to the successful application of the HDP model as an effective medium theory on DNMs and the fact that EAW modes are primarily determined by the morphology rather than geometrical parameters, such as a and P, it is reasonable to speculate that the new BICs are resilient to generic disorders. To illustrate this, we demonstrate high-Q states through numerical simulations, even when scattering into nonzero diffraction orders is possible. We specifically create a DNM slab with a large supercell of eight unit cells in the y direction, with a randomly generated inharmonious disorder in the spacing between the two nets while maintaining a fixed offset in the x and z directions of a2. We fixed the thickness of this slab to N=2. For the unperturbed DNMs, BICs appear at a frequency of 9.55 GHz at the Γ point, as shown in Fig. 2(e). A comparison of the band diagrams between the two structures is shown in Fig. 3(b). The perturbed supercell DNM slab retains the eigenmodes of the perfectly ordered DNM, weakly affecting their frequency position. The relative difference between the positions of the EAW standing waves at the Γ point is 0.05 GHz [see Fig. 3(b)]. Since the bands are back-folded into the Brillouin zone of the supercell geometry, we also fold back the bands of the unperturbed slab, shown as dashed lines. The Q-factors are represented by the color values of the solid lines for the perturbed case. The enlarged lattice plays a crucial role in allowing more diffraction channels, which permits the BIC around the frequency of 9.55 GHz, that is, k0=200m−1 to radiate into the first Bragg order at G≈160m−1 for the supercell lattice constant of 8a=40mm. In contrast, the unperturbed case is below the first Wood anomaly with the first Bragg order, G=πa≈1250m−1. The enlarged supercell DNMs transform BICs into quasi-BICs, maintaining an extremely high Q-factor of 23,427 even above the Wood anomaly, where the first Bragg order in the y direction constitutes an additional radiation channel for which the longitudinal EAW-like modes are not perfectly decoupled from vacuum radiation. The electric field distribution corresponding to the supercell DNM slab with random variations in the y direction is depicted in Fig. 3(c). The fact that the field is tightly confined within the DNM slab confirms its weak coupling with the surrounding environment. More interestingly, since additional branches of the in-plane EAW wave are now back-folded in the continuum [see Fig. 3(c)], they also give rise to a quasi-BIC band with an even higher Q-factor of approximately 4×106 for the kz=0 band. This is because the back-folded modes inherit their main intensity from the unperturbed DNM, which lies in the first Bragg order of the supercell outside of the light cone (G≈160m−1) at the relevant frequencies (k0⪅100m−1).

In summary, the BICs discovered in DNMs exhibit remarkable stability and robustness, even in the presence of extreme symmetry breaking and disorder above the Wood anomaly, which can be attributed to the HDP’s effective-medium nature of DNMs. Consequently, these new BIC modes in DNMs hold immense potential for various optical applications, including optoelectronic devices, integrated photonic circuits, and waveguide-based devices. By leveraging their exceptional stability and robustness, DNMs can pave the way for developing next-generation ultracompact and highly efficient optical systems.

Benefiting from the peculiar non-Maxwellian HDP effective medium nature of the DNMs with approximately vanishing overall EM fields, they might be used to build ultrahigh-Q individual quasi-BIC resonators with subwavelength dimensions. To elaborate on this potential, we study quasi-BIC resonances in a cubic DNM microresonator [Fig. 4(a)]. In the HDP effective medium picture, the eigenmodes in the individual cubic structure obey the 3D hard-boundary condition; that is, the normal current component of both plasmas must vanish at the cubic boundaries. This leads to eigenfrequencies ωmnl=vgπL|n| with n=(m,n,l)∈N3. The (m,n,l) indices label the allowed Cartesian momentum components in the cubic cavity responsible for the resonances. The lowest three resonance frequencies involve momentums marked as colored spheres in Fig. 4(b). The lowest mode is threefold degenerate (|n|2=1, red spheres, T1u irreducible representation in the Oh group48). The second mode is also threefold degenerate (|n|2=2, blue spheres, T2g irrep), while the third mode is nondegenerate (|n|2=3, yellow spheres, A2u irrep). Increasing the box length, L, simply redshifts the modes, as shown in Fig. 4(c) (dashed lines). For the DNMs, similar behavior is expected because of the HDP effective medium approximation. By increasing the DNM cube edge length, L, from 10 to 20 mm (two unit cells to four unit cells along the Cartesian axes), we find approximately the same trend in the frequency shifts, as in the HDP model. On the other hand, the DNM cavities exhibit a finite lifetime. The frequencies of the lowest three modes of the DNM cubic resonator are displayed in Fig. 4(c) as solid-hued lines, whose color represents the Q-factors. Note that we chose the resonator geometry such that the inversion center of the double net is in the center. This implies that individual nets are cut along their long axis at the boundary of the resonator at 15 mm. We avoid such cuttings, explaining the discontinuity of the solid lines in Fig. 4(c). The graph confirms the validity of the effective medium description of the DNMs. It becomes, however, clear that the microscopic boundary configuration plays a crucial role in determining the Q-factor, which ranges from 102 to 107. Similar to the slab configuration, the quasi-BIC states manifest themselves as resonances in the scattering cross sections (SCs). Figure 4(d) shows a slightly asymmetric Fano resonance dip in the SC at the fundamental quasi-BIC mode at 7.54 GHz, whose Q-factor is as high as 2400 for the cube’s edge length of 12.5 mm. The appearance of an SC dip rather than a peak is owed to the destructive interference with the background scatterings, whose asymmetry is more pronounced for lower Q-factor resonances. The fundamental quasi-BICs mode is located in frequency space at approximately the boundary between the Rayleigh scattering and the Mie scattering regimes without high-index dielectric filling within the structure. It is, however, possible to drive the quasi-BICs to operate within the Rayleigh scattering regime by introducing high-index dielectrics. The total field intensity in Fig. 4(e) shows the reduced disturbance to the incident plane wave and the pronounced field build-up in the DNMs on the resonator frequency at 7.541 GHz, when compared to the off-resonance frequency at 7.55 GHz.

The extremely small physical size of the resonator (on the order of a few unit cells) is in stark contrast to previously reported purely dielectric quasi-BIC resonators. The reduced SC and thus perturbation to the incident EM field and the gigantic field build-up in the DNM cube could give rise to stealth and highly sensitive microwave probing antenna opportunities.

In summary, this study presents an approach to achieving photonic BICs in DNMs. The HDP is employed as an effective medium model for DNMs. An exact analytical solution for the refractive index of EAWs is derived in the long wavelength limit. The non-Maxwellian HDP leads to a unique type of quasi-longitudinal EAW band, which results in the emergence of BICs at arbitrary frequencies within the frequency regime of the EAW wave, free of non-BIC resonances and shows remarkable stability and robustness even in the presence of strong geometrical perturbations, symmetry breaking, general disorder, and being above the Wood anomaly, where coupling into higher Bragg orders is possible. Importantly, the BIC modes can maintain high Q-factors even if the structure is contained within a subwavelength scale, as demonstrated for an ultrasmall cubic DNM resonator. The discovery of these resilient BIC modes in DNMs presents exciting opportunities for various optical applications, including optoelectronic devices, integrated photonic circuits, and waveguide-based devices. It also holds the potential for developing innovative ultracompact and highly efficient optical systems. This study provides valuable insights into optical metamaterials and their practical applications. While the theoretical predictions, based on a lossless plasma and a corresponding perfectly conducting double net, can be readily explored in the microwave regime, an extension of the theory in future work would be desirable to quantify the nonradiative contribution to the quasi-BIC lifetimes at infrared frequencies.

Wenhui Wang achieved her PhD from the University of Science and Technology of China, in 2016. Afterward, she became a member of the Meta Materials Group at Birmingham University, United Kingdom, and subsequently joined the Research Laboratory for Quantum Materials at Singapore University of Technology and Design. The research presented here was conducted during her time as a senior researcher within the Soft Matter Physics Group at the Adolphe Merkle Institute in Switzerland. Her research interests span metamaterials, plasma dynamics, and topological photonics.

Matthias Saba received a PhD in physics in 2015 from the University of Erlangen–Nuremberg (Germany). He subsequently joined the Condensed Matter Theory Group at Imperial College London, where he studied 3D topological photonic materials and metamaterials. His current research interests as a group leader in the Soft Matter Physics Group at the Adolphe Merkle Institute (Switzerland) include double-net metamaterials, chiral thermal plasmonic crystals, and the theory of photonic states in finite periodic and disordered structures.

Biographies of the other authors are not available.