Domino plasmons for subwavelength terahertz circuitry

A new approach for the spatial and temporal modulation of electromagnetic fields at terahertz frequencies is presented. The waveguiding elements are based on plasmonic and metamaterial notions and consist of an easy-to-manufacture periodic chain of metallic box-shaped elements protruding out of a metallic surface. It is shown that the dispersion relation of the corresponding electromagnetic modes is rather insensitive to the waveguide width, preserving tight confinement and reasonable absorption loss even when the waveguide transverse dimensions are well in the subwavelength regime. This property enables the simple implementation of key devices, such as tapers and power dividers. Additionally, directional couplers, waveguide bends, and ring resonators are characterized, demonstrating the flexibility of the proposed concept and the prospects for terahertz applications requiring high integration density.


I. INTRODUCTION
Electromagnetic radiation in the terahertz (THz) regime is a central resource for many scientific fields such as spectroscopy, sensing, imaging, and communications. We are currently witnessing the take off of THz technologies 1,2,3 with applications as diverse as astronomy 4 , medicine 5 , or security 6 . Within this general endeavour, the building of compact THz circuits would stand out as an important accomplishment. This requires the design of THz waveguides carrying tightly confined electromagnetic (EM) modes, preferably with subwavelength transverse cross section. Besides circuit integration and compact device design, sub-λ localization may be advantageous for waveguide THz time-domain spectroscopy 7 and non-diffraction-limited imaging 8 . Although a number of structures has been put forward, none of them passes all the following requirements. First, structures should be easily manufactured and, if possible, planar and monolithic. Second, as above motivated, subwavelength transverse size is also needed. Finally, absorption and bend losses should be small, and the waveguides sufficiently versatile for the design of key functional devices. Particularly important are in/out-couplers since they work as the interface to external waves. In this context, compact tapers able to laterally compress the modes down to the sub-λ level seem essential.
With these requirements in mind it is illuminating to compare optical waveguides with those proposed for the THz range. Optical fiber modes exhibit a cut-off diameter such that fibers cannot be of subwavelength size and, moreover, modal size grows as this cut-off is approached. In the THz regime dielectric wires 9 meet the same limitation. The photonic crystal fiber approach 10 shares the inability of sub-λ confinement, and low-index discontinuity waveguides 11 perform better regarding modal size but are neither planar nor monolithic. Metals are an alternative when operating at optical frequencies because surface plasmon polaritons (SPPs) propagating on flat metaldielectric interfaces provide subwavelength confinement in the direction perpendicular to the surface. This metallic route has been tried for the THz regime, but the corresponding EM modes, known as Zenneck or Sommerfeld waves, loose their confined character at these frequencies 12,13 .
Geometrically-induced surface plasmons 14 supported by periodically corrugated metallic surfaces overcome the low localization of Zenneck waves. A number of waveguiding structures based on this concept has been already suggested 15,16,17,18 , but they have failed to meet simultaneously all the above mentioned requirements. In this work we present structures based on this metamaterial approach that, thank to their superior properties, fulfill the requisite easy fabrication, subwavelength confinement, low loss, and device flexibility.

II. RESULTS
The basic structure consists of a periodic arrangement of metallic parallelepipeds standing on top of a metallic surface and resembling a chain of domino pieces, see inset of Fig. 1(a). In contrast with corrugated wires 16 , this is a planar and monolithic system and should not pose significant manufacturing problems, its geometry being much simpler than corrugated V-grooves 17 or wedges 18 . The properties of its guided modes, hereafter referred to as domino plasmons (DPs), are mainly controlled by the geometric parameters defining the dominoes, i.e., periodicity d, parallelepiped height h, lateral width L, and inter-domino spacing a.
First, the basic properties of DPs are described starting with the dispersion relation. To gain physical insight we model the metal as a perfect electric conductor (PEC). Inset: diagram of the domino structure and geometric parameters (the arrow depicts the mode propagation direction). b and c, Modal shape of DPs: transverse (xy) electric field (arrows) and horizontal (xz) electric field (color shading) for DPs of height h = 1.5d and widths L = 0.5d (b) and L = 8d (c). d and e, Same as in panels b and c, but now for a 1D array of free-standing metallic rods (h = 3d). The designations even and odd label the symmetries of the modes with respect to the corresponding white lines. The fields in b-d are plotted for d/λ = 0.125, marked with a red open dot in a, the white bar in c being the wavelength (valid for panels b-d).
The field in e is computed for the same k as in panels b-d, corresponding now to d/λ = 0.056. In panels a-e metals are modelled as PECs. f, DP modal effective index as a function of lateral dimension L in units of wavelength. Various operating frequency regimes are considered: λ = 1.6 mm (red), λ = 0.16 mm (green), λ = 0.016 mm (blue), and λ = 1.5 µm (magenta). To compute panel f, a realistic description of the metals is used. As described in the main text, the periodicity d is different for the various operating frequencies, and h = 1.5d, a = 0.5d, L = 0.5d, . . . , 24d.
Later, realistic metals will be considered, letting us to discern what is the frequency regime where DPs display the most promising behaviour. The numerical technique is described in the Methods section. Due to the simple scaling properties of PECs, we choose the periodicity d as the length unit. The value of a is not critical for the properties of DPs and is set as a = 0.5d in this work. DP bands present a typical plasmonic character, i.e., they approach the light line for low frequencies and reach a horizontal frequency limit at the edge (k edge = π/d) of the first Brillouin zone ( Fig. 1(a), notice that only fundamental modes are plotted). In all cases DP modes lie outside the light cone, a fact that explains their non-radiative character. While the limit frequency of SPPs for large k is related to the plasma frequency, the corresponding value for DPs is controlled by the geometry. The influence of the height h is clear: the band frequency rises for short dominoes (h = 0.75d, blue line) as compared with that of taller ones (h = 1.5d, black line). This can be understood as follows. When L = ∞, dominoes become a one-dimensional (1D) array of grooves which supports a geometrically-induced plasmon mode. Its dispersion relation is represented, for h = 1.5d, with a dashed line in Fig. 1(a). An approximation for the dispersion relation in this limit, neglecting diffraction effects and for λ >> d, where k is the modal wave vector, k 0 = 2π/λ, and the y component of the wave vector inside the grooves, q y , is related with the corresponding x and z components by q 2 x + q 2 y + q 2 z = k 2 0 . The EM fields inside the grooves are independent of x and z for the chosen longitudinal polarization and thus q x = q z = 0, q y = k 0 . This means that Eq. 1 presents an asymptote for k 0 h = π/2 or, in other words, for a normalized frequency d/4h. Such dependence with h explains the behaviour observed in Fig. 1

(a).
For h = 1.5d, the value chosen in the rest of the paper, this estimation for the limit normalized frequency provides a figure of about 0.17, not far from the numerically computed value of 0.14.
The most striking characteristic of DPs is their behaviour when the lateral width L is changed. All bands in the range L = 0.5d, . . . , 3d lie almost on top of each other ( Fig. 1(a), grey curves). In other words, the modal effective index, n eff = k/k 0 , is rather insensitive to lateral width. If we now look at the numerically computed dispersion relations of Fig. 1(a) in the light of Eq. 1, it seems that q y ≈ k 0 for a wide range of L. Remarkably, the bands remain almost unchanged even for L = 0.5d, whose modal size is well inside the subwavelength regime, as observed in the field distribution plotted in panel (b) for a narrow structure (L = 0.5d). Notice that the white bar representing λ in panel (c) is valid for panels (b)-(d), and the modal size, defined as the FWHM of the modulus of the Poynting field distribution, is 0.21 λ for panel (b). The described behaviour is to be contrasted with that of conventional plasmonic modes in the optical regime for which sub-λ lateral confinement is not a trivial issue 20 . When the lateral size of the structure supporting such modes becomes subwavelength, either the modal size grows, or the effective index increases. This is best illustrated for a free standing metallic stripe 21 . A stripe of large lateral width essentially supports an SPP mode but, as this width is decreased and becomes sub-λ, the plasmonic modes that are not cut-off become either highly confined but quite lossy, or otherwise highly delocalized.
The insensitivity of DPs to lateral size, which constitutes the foundation of some of the devices presented later, can be linked in an intuitive way to the structure geometry. The key role played by the metallic substrate can be explained by comparison of DPs and the modes of a Yagi-Uda structure, i.e., a periodic array of free standing metallic rods 22,23,24 . For small L the dispersion of DPs of height h is identical to that of the fundamental mode of a Yagi-Uda structure of height 2h, and the field of the DP is simply the upper half of the fundamental mode of a rod array of height 2h, see Figs. 1(b) and (d). For large L this correspondence does not hold anymore, as demonstrated by comparison of panels (c) and (e) which display completely different symmetries. For dominoes, a mode with the same symmetry as the Yagi-Uda mode shown in panel (e) is forbidden because the horizontal electric field has to be zero at the metallic ground. This mode suppression is critical because the frequency of the fundamental mode of a Yagi-Uda structure, which is essentially controlled by the length of the long axis of the rods, is very sensitive to the lateral width L, being very different for short vertical rods and long horizontal rods. It is similarly important that the air gaps between the dominoes are laterally open. Let us imagine that, instead of being open, they were laterally closed by PEC walls. Since the electric field inside the gap points longitudinally from one domino to the next one, it should then vanish at both PEC lateral walls, giving a cosine-like x dependence. It is then clear that the mode frequency should increase when the lateral size L is decreased, arriving to the conclusion that the modes of such an structure have L-sensitive bands. In contrast, for DPs the boundary conditions at the laterally open gaps are more akin to perfect magnetic walls, and the corresponding field is approximately x-independent inside the air gaps even when L → 0.
Let us scrutinize in more detail the role played by the lateral dimension L, considering now realistic metals and paying attention to the spectral regime. The periodicity d is chosen to set the operating wavelength in the desired region of the EM spectrum, and L is varied in the range L = 0.5d, . . . , 24d, while the remaining parameters are kept constant (a = 0.5d, h = 1.5d). Aluminum is selected for low frequencies, where metals behave almost like PECs. In order to work at λ = 1.6 mm, providing an operating angular frequency of the order of 1 THz, we first consider d = 200 µm. The evolution of the modal effective index as a function of the lateral dimen-sion normalized to the wavelength is plotted in Fig. 1(f). The curves corresponding to a PEC (black line) and aluminum at λ = 1.6 mm (red line) are, as expected, almost identical. We can now quantify the sensitivity of the effective index to L, its variation being only about 12% even when L goes from L = ∞ to L = 0.5d = λ/16, well inside the sub-λ regime. The rapid variation of n eff when L << d can be understood with the help of Eq. 1. In this limit the EM fields start to be z-dependent inside the air gaps, as verified by inspection of the fields (not shown here), letting us to conclude that q z = 0 and consequently q y < k 0 . Then, if we assume that Eq. 1 is still valid in this limit, it predicts that the modal effective index is reduced when L → 0. To investigate the performance of DPs at higher frequencies, the structures have been scaled down by factors 1/10 and 1/100. The fact that the curves corresponding to λ = 0.16 mm (green line) and λ = 0.016 mm (blue line) do not lie on top of the previous ones is a signature of the departure of aluminum from the PEC behaviour. Nevertheless, even at λ = 0.016 mm, the variation of the effective index is still smaller than 15%, a property that can be exploited for taper design in the far infrared, albeit with slightly reduced performance as compared with the THz case. When the operating frequency moves to the telecom regime (λ = 1.5 µm, magenta line) the variation of the effective index is much larger (about 38%). Let us mention in passing that, for this wavelength, we have used gold because it has lower loss than aluminum, and that the geometric scaling had to be different (the periodicity is now d = 0.13 µm) because with the same simple scaling as before no DP mode is supported due to the value of the permittivity of gold. The important message of panel (f) is that, although a variation of n eff begins to be noticeable when the lateral dimension L goes below λ/2, DP bands are fairly insensitive to L in the range L = 0.5d, . . . , 24d when operating at low frequencies. Such property has, to our knowledge, never been reported for geometrically-induced plasmon modes in corrugated wedges, channels, or wires. Now that we have identified the THz range as the most promising operation regime, we will characterize the corresponding absorption (ohmic) and bending losses. . The only source of damping in the domino structure is absorption by the metal as no radiation losses occur in such straight waveguide. Propagation length is very short for small λ, when the Brillouin zone edge is reached, but it rises for increasing wavelength as the dispersion curve approaches the light line. Several lateral widths L are considered and it is observed that DPs in structures of smaller L travel longer distances before being damped. The radiation losses occurring in a waveguide of L = 0.5d with a 90 degree bend are analyzed in panel (b). In this case the metal is modelled as a PEC for easier evaluation of radiation losses. Data for three radii of curvature of the bend, R = λ/2, λ, 2λ, are presented. In all three cases, losses are low for small λ but they grow as the wavelength is increased and the modal wave vector approaches k 0 . The dependence with the radius of curvature of the bend is not surprising either, radiation losses being higher for smaller radius of curvature. Consideration of absorption and bending losses suggests that λ = 1.6 mm can be an appropriate operating wavelength. For λ = 1.6 mm and L = 0.5d = 100 µm, the modal size is 0.34 mm, i.e., about λ/5. The propagation length is about 20λ, much larger than the size of the devices presented later, and the bend loss for R = λ is reasonably low, about 10%. Notice that for corrugated V-grooves, the bend loss for R = 2λ is around 50% 17 .
In the rest of the paper we present a number of devices enabling the spatial and temporal modulation of EM fields. The operating wavelength is λ = 1.6 mm and d = 200 µm. In all ensuing computations the metal is considered as a PEC, to easily quantify the device radiation loss. Waveguide tapering 16,18,25,26 is interesting for field concentration and amplification, both important requirements for imaging applications 27,28 and circuit integration 29 . Tapering is expected to be easy in domino structures due to the insensitivity to lateral width L of the dispersion relation of DPs. This hypothesis is confirmed in Fig. 3(a), which shows a top view of the propagation of power along a tapered domino structure. The lateral widths of the waveguide are L in = 16d at the input port (bottom) and L out = 0.5d at the output port (top). The length of the taper is 16d, i.e., two wavelengths at the operating frequency, corresponding to a taper semiangle of θ = 31 • . Remarkably, reflection is smaller than 2% and only 5% of the incoming power is lost as radiation due to the shrinkage of the lateral dimension along the taper. Panels (b)-(e) are cross sections at different positions along the taper, vividly showing the process of field subwavelength concentration and enhancement. As mentioned above, the same design does not work in the telecom regime due to the sensitivity of DP bands to the lateral dimension L at those frequencies. A more adiabatic taper, i.e., with smaller change of n eff per unit length, would not circumvent the problem because absorption loss is high in the telecom regime (for instance, for λ = 1.5 µm, d = 0.13 µm, and L = 16d, the propagation length is = 1.2 µm).
We now discuss how to construct power dividers and directional couplers 30,31 . Figure 4 plots the dispersion relation of the modes supported by one and two parallel domino chains. Black solid line corresponds to the DP of L = 0.5d. For two parallel waveguides (both with L = 0.5d), interaction gives rise to symmetric and antisymmetric supermodes which, for a separation s = 0.5d, are shown by magenta curves with the solid line depicting the even supermode whereas the dashed one represents the odd supermode. Blue dashed-dotted line corresponds to a wider structure of L = 1.5d and, importantly, this band stays very close to that of the s = 0.5d even supermode. Due to this good overlapping, a simple linking of a wide domino structure with L = 1.5d to two parallel domino chains (widths L = 0.5d and separation s = 0.5d) should not result in an important impedance mismatch. This is verified in Fig. 5(a), displaying a top view of a power divider based on the previous idea. The input port (L = 1.5d) at the bottom feeds the wide domino which, after 5 periods, is linked to two parallel domino structures (L = 0.5d, s = 0.5d). Due to symmetry only the even supermode is excited. The lateral separation s between the individual waveguides is then increased, followed by 90 • waveguide bends. Simulations show that reflection is also negligible in this power splitter. Total radiation loss in the device is equal to that happening in the bends (previously discussed in Fig. 2(b)), each being 5% of the incoming power, which is small taking into account that the bend radius is equal to the operating wavelength.
Coupling between parallel domino structures can also be easily understood with the help of Fig. 4, displaying with green lines the even (solid) and odd (dashed) supermodes for parallel domino chains separated a distance s = 1.25d. The insets show the modal shape of both modes. A larger s is chosen in this simulation to prevent the cut-off of the odd supermode. The coupling length needed to exchange the carried power between the waveguides is estimated as π/|k even − k odd | = 2.05 mm. This value agrees with the result found in the simulation of a directional coupler. Figure 5(b) displays a top view of two parallel domino chains (L = 0.5d) separated a distance s = 1.25d. A mode is launched in the left bottom waveguide and characteristic beatings with the predicted coupling length are observed when a second waveguide is located close to it.
As a final example a waveguide ring resonator is demonstrated, able to route various frequencies to different output ports 32 . To show once more the robustness of our proposal, this device has been designed with posts of circular cross section as seen from above (radius r = 0.28d), instead of square cross section parallelepipeds. We have checked that similar results are obtained for the square geometry. The top view of this device is displayed in Fig. 5(c). For the chosen geometric parameters (ring radius R = 2 mm, minimum separation in the coupling sections s = 1.25d), it is possible to drop a resonant wavelength λ = 1.54 mm at the output port. Non-negligible radiation losses only occur at the bends in the waveguide ring and they are about 25% of the input power, reasonably small considering that no optimization was attempted.
In summary, we have presented a new class of surface electromagnetic modes, termed domino plasmons, which feature an extraordinary characteristic: their insensibility to the lateral dimension of the structure supporting them. We have shown how the guiding properties of these modes enable the planar routing of THz radiation at a deep subwavelength level. Finally, the flexibility and versatility of waveguides based on domino plasmons have been demonstrated through the implementation of a variety of functional devices, which may thus prove very useful as a new concept for subwavelength THz circuitry.

III. METHODS
All results in this paper have been obtained by means of numerical simulations performed with the Finite Element Method (FEM) using commercial (COMSOL Multiphysics) software. A summary of the main modelling details follows. The technique is a volume FEM operating in frequency domain. Due to the various length scales involved, highly non-uniform meshes are used. In order to ensure an adequate representation of the electromagnetic fields, the size of the elements is a fraction of the skin depth inside the metallic parallelepipeds, a fraction of the characteristic geometric dimension (e.g., a) in the neighborhood of the structure, and a fraction of the operating wavelength at the boundaries of the simulation domain. The mesh size is refined until the results are stable with an accuracy of 1%. For instance, in the case of band structure computations, the convergence criterium is based on the value of both the real and imaginary parts of the wave vector. The final mesh is different for each simulation, but typical values of the tetrahedra sizes in different regions of the simulation domain are of the order of 1/5 of the skin depth and char- acteristic geometric dimension, and 1/10 of the wavelength. The typical number of degrees of freedom lies between 3 × 10 5 and 3 × 10 6 , a range in which the matrix equations are inverted with direct solvers. For the computation of band structures, the corresponding eigenvalue problem is posed in a single unit cell where Bloch boundary conditions together with Scattering boundary conditions are used. For device modelling, the excitation of the corresponding scattering problem is provided by Port boundary conditions, the ports being fed by the solutions of the corresponding eigenvalue problem. Open space is mimicked with Perfectly Matched Layers or/and Scattering boundary conditions. For the simulation of ideal metals, Perfect Electric Conductor (PEC) boundary conditions are employed. For realistic metals the dielectric permittivities are taken from Ref. [33] for Al ( Al = −3.39 × 10 4 + i3.5 × 10 6 at λ = 1.6 mm), and Ref. [34] for Au ( Au = −103 + i8.7 at λ = 1.5 µm). At low frequencies it was often sufficient to simulate the domains outside the conductors together with Surface Impedance Boundary Conditions, as verified by comparison with simulations fully accounting for the EM fields inside the metals.