The isofrequency surfaces (IFSs) of a medium play an important role in determining the electromagnetic (EM) wave propagation[1] and light–matter interactions [2, 3] inside the medium. For a natural material, the shape of IFSs depends on its microscopic atomic arrangement and is usually an isotropic sphere or an anisotropic ellipsoid in k space. Metamaterials (MMs) [4–10], artificially synthesized materials with subwavelength microstructures, can exhibit various types of closed/open IFSs with exotic topologies (concave sphere [11, 12], hyperboloid [2, 3, 13, 14], cuboid [15, 16], etc.), which lead to anomalous wave phenomena (negative refraction [1, 11, 12], subwavelength imaging [17–20], self-collimation [15, 16, 21, 22], etc.) alongside numerous engineering applications. However, since their unnatural responses rely on the EM resonance of metallic structures [4, 5, 10], these exotic IFSs and associated phenomena only occur near the resonant frequency, which limits their working bandwidth. Also, when the wave frequency approaches quasistatic limit, the IFSs of these artificial media usually deteriorate to an ordinary ellipsoid centered at a zero k-point.

It would be highly desirable if the shapes of low-frequency IFSs can be customized at will, then we could control the EM wave propagation behavior in a relatively broader bandwidth. It is well known that metallic wire medium [4, 23–40], such as disconnected wire arrays, can exhibit the strong spatial dispersion (nonlocal effect) even in long-wavelength limit [25, 26]. By incorporating two or more sets of metallic meshes, one can engineer the density of states at quasistatic limit [28, 32, 33]. The EM properties of these multiple interlocking meshes are governed by electric circuit equations [28, 29] rather than Maxwell equations [41]; they are, thus, called non-Maxwellian medium. Interestingly, researchers [33] found that the global connectivity of these metallic meshes can introduce a momentum shift of the low-frequency index ellipsoid away from the zero k-point. Recently, wire medium again attracted the attention of researchers, since it could serve as a new design degree of freedom to tailor the low-frequency dispersion [34–37] and associated wave properties [38–40] of MMs.

In this study, based on connection-type MMs [28, 32–34, 36] that support multiple quasistatic modes at a nonzero k-point, we propose a scheme to engineer the shapes of low-frequency IFSs. As frequency increases from zero, isotropic or anisotropic index ellipsoids emerge from these nonzero k-points. By properly shaping the index ellipsoids, one can, in principle, realize arbitrary-shaped isofrequency contours (IFCs) even for the lowest frequency band. The benefit of our method is that it may provide a novel way for manipulating long-wavelength EM wave. To illustrate our idea clearly, we focus on a 2D slab of connection-type MMs, instead of complex 3D MMs. We develop an equivalent circuit model for the connection-type MMs to describe their low-frequency EM behaviors. The model implies that the shape of index ellipsoids can be engineered via the equivalent inductances and capacitances of the metallic meshes, which are determined by the structural details of the meshes. By appropriately adjusting the metallic structure, we can achieve isotropic or anisotropic index ellipsoids with different group velocities, which agree well with our equivalent circuit model.

Connection-type MMs are a type of MMs composed of N sets of interpenetrating metallic meshes, which support (N−1) quasistatic modes lying at a nonzero k-point. Their positions in k space are determined by the global connectivity of the meshes, rather than the structural detail of the meshes. For instance, 2D double-mesh MMs with square lattice [33] can support one quasistatic mode at the Brillouin zone corner (M point), as shown in Figure 1. In addition, this quasistatic mode is not an isolated solution at zero frequency. When frequency increases, the index ellipsoids turn into a circle or ellipse from an isolated point at M, which corresponds to a linear cone dispersion emerging from the M point. While the emerging point can be controlled by designing the mesh connectivity, the shapes of index ellipsoids or the group velocities of linear bands are determined by the structural details. Once we can engineer both the centers and the shapes of the low-frequency index ellipsoids, arbitrary-shaped IFCs can be achieved in principle. For the case in Figure 1, if the index ellipsoids emerging from the M point are a circle, then cruciform-shaped IFCs (dashed line in Figure 1A) would be expected for the midband frequency of band 1. Also, if the index ellipsoids emerging from the M point are an ellipse, bowtie-shaped IFCs (dashed line in Figure 1B) would be expected. Generally speaking, multiple zero-frequency modes would appear at multiple k-points when considering more sets of interlaced wire meshes. By engineering the shape of each ellipse via the equivalent inductances and capacitances between multiple meshes, arbitrary-shaped IFCs can, in principle, be realized in connection-type MMs.

To analyze the low-frequency index ellipsoids, we establish an equivalent circuit model which could serve as a guideline in designing connection-type MMs. Here, we consider a 2D double mesh arranged in a square lattice. Centered in its unit cell is a metallic block (height d and width w), whose top (bottom) end is connected to the adjacent unit cell via the metallic bar in the u ̂ 1 ( u ̂ 2 ) direction (Figure 2A). When the unit cell is repeated in x and y directions (3 × 3 supercell in Figure 2B), it is actually composed of two disconnected meshes (colored in blue and red) arranged in a checkerboard manner. Since all the metallic components in this study are assumed to be a perfect electric conductor (PEC), blue and red colors just highlight the two sets of meshes with independent quasistatic potentials. According to the quasistatic potential analysis in [33], these 2D double-mesh MMs have one quasistatic mode lying at zone corner M.

We construct a circuit model for this structure to describe its Bloch modes at relatively low (but not zero) frequency. Note that the double-mesh MMs are a slab with finite thickness in the z-direction embedded in air background, and we only consider the EM modes with Bloch k near the M point. Since Bloch k lies below the light cone, the EM mode is guided in the slab and its leakage into air can be neglected. The EM field of Bloch mode would induce alternating currents on the metallic wires and charge accumulation. In Figure 2A, the total current passing through the top (bottom) wire along the u ̂ 1 ( u ̂ 2 ) direction is assumed to be I u ̂ 1 ( I u ̂ 2 ) . Taking into consideration the Bloch boundary condition, the currents at the unit cell corners should be I u ̂ 1 e i 2 k x + k y a and I u ̂ 1 e − i 2 k x + k y a ( I u ̂ 2 e i 2 k x − k y a and I u ̂ 2 e i 2 − k x + k y a ) along the u ̂ 1 ( u ̂ 2 ) direction. Due to these alternating currents, a time-varying charge Q should accumulate on the metallic block, which satisfies the following equation: i ω Q = I u ̂ 1 e i 2 k x + k y a − I u ̂ 1 e − i 2 k x + k y a + I u ̂ 2 e i 2 k x − k y a − I u ̂ 2 e i 2 − k x + k y a . (1)

Since we are dealing with the low-frequency Bloch modes at the k-point not far away from M, the charge on each block should be almost plus or minus alternatively in the square array. Also, as the two meshes serve as two conductors of a capacitor at quasistatic limit, the charge on one block should be proportional to the local potential difference between the two meshes and written as follows: Q = C Δ φ , (2) where Δφ is the local potential difference between a metallic block and its nearest four neighboring blocks. On the other hand, the alternating currents flowing on the wires would induce alternating magnetic fields, which are here specified by vector potential field A ⃗ r ⃗ . They satisfy the following relation: A u ̂ 1 A u ̂ 2 = L u ̂ 1 L u ̂ 2 I u ̂ 1 I u ̂ 2 , (3) where A u ̂ n (n = 1, 2) is a line integral of the vector potential on the path between two neighboring blocks (gray solid line in Figure 3B). The inductance matrix is a 2 × 2 matrix whose off-diagonal terms vanish due to the mirror symmetry of the structure. As we are considering the Bloch modes near the M point, I u ̂ 1 in Eq. (3) represents a series of currents (Figure 3A), rather than a single line current, in an infinite square lattice. Meanwhile, A u ̂ 1 represents the vector potential field (Figure 3B) induced by these collectively oscillating currents. Thus, the inductance L u ̂ 1 is actually the ratio between the current array and its induced vector potentials, rather than the self-inductance or mutual inductance of a single wire.

For a perfect conductor, the microscopic electric field inside should vanish, and we have the following equation: i ω A u ̂ 1 = Δ φ e i k x + k y a − Δ φ , i ω A u ̂ 2 = Δ φ e i − k x + k y a − Δ φ . , (4) where the local potential difference Δφ satisfies Bloch theorem, and the voltage drop between two adjacent unit cells can be evaluated by the difference between the local potential differences of the two unit cells.

Combining Eqs. 1–4, we obtain the following generalized eigenvalue equation: e i 2 k x + k y a − e − i 2 k x + k y a e i 2 k x − k y a − e i 2 − k x + k y a − i ω C − i ω L u ̂ 1 e i k x + k y a − 1 − i ω L u ̂ 2 e i − k x + k y a − 1 I u ̂ 1 I u ̂ 2 Δ φ = 0 0 0 . (5)

The equation has a nonzero eigenvector only when the determinant in Eq. 5 vanishes; then, the low-frequency dispersion can be derived as follows: ω 2 = 4 C − 1 L u ̂ 1 − 1 ⁡ sin 2 δ k x + δ k y a 2 + L u ̂ 2 − 1 ⁡ sin 2 δ k x − δ k y a 2 . (6)

The band dispersion is controlled by two equivalent inductances and one equivalent capacitance mentioned above, which can be calculated in magnetostatic and electrostatic limits. Also, this model applies well for low frequencies, as long as the Bloch modes near the quasistatic mode (at the M point in this case) are considered (δk _x and δk _y are sufficiently small compared to the reciprocal vector). Although the equivalent circuit model would deteriorate as the frequency increases, it still gives a qualitative prediction on the IFCs at midband frequency and can serve as a guideline for MM design.

From Eq. 6, the group velocities in u ̂ 1 and u ̂ 2 directions are written as v g , u ̂ 1 = a 2 C − 1 L u ̂ 1 − 1 , v g , u ̂ 2 = a 2 C − 1 L u ̂ 2 − 1 . (7) and can then be controlled by the equivalent capacitance and two equivalent inductances. Since the group velocities depend on the equivalent capacitance and inductances, the shape of the index ellipsoid can be controlled by adjusting the structural details of meshes, which determine the equivalent circuit parameters.

To investigate the dependence of group velocities on the structural parameters, we numerically calculate the equivalent inductances and capacitances of a series of double-mesh MMs. First, we alter the height of the metallic block centered at the unit cell while keeping other parameters unaltered (width w _p = 0.4a and w u ̂ 1 = w u ̂ 2 = 0.05 a ). Results from electrostatic and magnetostatic simulations are shown in Figure 4A. It is found that the change of block height would mostly alter the equivalent capacitance C between the two meshes, rather than inductances L u ̂ 1 and L u ̂ 2 . The reason is that the block height would increase the surface area of the meshes and, thus, the equivalent capacitance rises. On the other hand, the equivalent inductances are related to the width of the horizontal metallic wires and the distance between them [29]. The two inductances almost remain unchanged, as in Figure 4A. Hence, the group velocities at quasistatic limit along u ̂ 1 and u ̂ 2 directions should decrease as the block height changes. Figure 4B shows the results obtained from calculated band structures and the ones from equivalent parameters, which agreed well with each other.

In another case, we increase the u ̂ 1 wire width while keeping the other structural parameters unchanged (height d = 0.8a, width w _p = 0.4a, and w u ̂ 2 = 0.05 a ). The results are plotted in Figures 4C,D. In this case, the symmetry of the structure is reduced from S ₄ to C ₂, resulting in anisotropic inductances in u ̂ 1 and u ̂ 2 directions. The inductance L u ̂ 1 (purple open circle in Figure 4C) decreases, while L u ̂ 2 (blue solid circle in Figure 4C) remains almost unchanged. In the meantime, the equivalent capacitance slightly increases due to the surface area of the u ̂ 1 wire. To sum up, the group velocity v g , u ̂ 1 becomes faster as w u ̂ 2 increases while v g , u ̂ 2 becomes somewhat slower, see the symbols (model) and lines (simulation) in Figure 4D. Here, the index ellipsoids are oriented in u ̂ 1 and u ̂ 2 directions. In principle, the index ellipsoids can be oriented in other directions by properly designing the widths of four horizontal wires in the unit cell.

To further see the applicability of our model, Figures 5A,B plot the numerically simulated band structure and index ellipsoids for isotropic double-mesh MMs (height d = 0.8a, width w _p = 0.4a, and w u ̂ 1 = w u ̂ 2 = 0.05 a ) while Figure 5C gives the index ellipsoids predicted by our model. They agree well with each other and validate our theory. In both results, the index ellipsoids gradually evolve from a circle to a squircle shape as the frequency increases, due to the two sinusoidal terms in Eq. 6.

To see the influence for different equivalent inductances in u ̂ 1 and u ̂ 2 directions, we increase the bar width of the metallic wire in the u ̂ 1 direction in order to weaken its inductance L u ̂ 1 . Meanwhile, all the other parameters remain unchanged (height d = 0.8a, width w _p = 0.4a, w u ̂ 1 = 0.4 a , and w u ̂ 2 = 0.05 a ) (inset of Figure 5E). The group velocity in the u ̂ 1 direction of the first band becomes faster than that in the u ̂ 2 direction, leading to elliptical index ellipsoids simulated in Figure 5F. Also, as the frequency increases, the shape of index ellipsoids evolves gradually from an ellipse to a rounded rhombus. This is well predicted by our circuit model (Figure 5G).

According to our circuit model, we can achieve isotropic or anisotropic shape of index ellipsoids at quasistatic limit by adjusting the equivalent inductances in u ̂ 1 and u ̂ 2 directions. Then, the special-shaped IFCs can be achieved for midband frequencies of band 1, as proposed in Figure 1. For isotropic wire MM slabs in Figure 2A with appropriate structural parameters (height d = 0.42a, width w _p = 0.4a, and w u ̂ 1 = w u ̂ 2 = 0.2 a ), an index ellipsoid emerges from the M point due to the isotropic equivalent inductances. For midband frequencies (f = 0.2 (c/a) for example), the IFCs become cruciform shaped (Figure 6A). In another case, the index ellipsoid at M can be anisotropic when two equivalent inductance parameters are unequal. Considering the MM slabs with structural parameters (height d = 0.36a, width w _p = 0.4a, w u ̂ 1 = 0.3 a , and w u ̂ 2 = 0.05 a ), we have bowtie-shaped IFCs at around f = 0.2 (c/a) (Figure 6B). In addition, the frequency range of special-shaped IFCs can be adjusted by the mutual capacitance between double meshes.

In conclusion, we propose an equivalent circuit model for connection-type MMs to analyze its low-frequency dispersion away from quasistatic limit, taking into account the mesh connectivity. The model implies that the low-frequency IFCs can be engineered via the equivalent inductances and capacitance of the metallic meshes. For instance, we adjust the metallic block height and the wire width according to the model and demonstrate the feasibility of shaping the low-frequency IFCs. Coupled with the degree of freedom of the quasistatic mode at an arbitrary k-point, one can, in principle, realize the arbitrary topology of IFCs and manipulate the EM propagation properties using an artificial medium.