Hyperbolic metasurfaces exhibit unique dispersion and polarization properties, making them a promising platform for a plethora of photonic applications. At the same time, the ability to engineer the hyperbolicity via the predefined spectral positions of the metasurface resonances remains a notable challenge. Here, we analyze the dependencies of the spectral positions of the resonances corresponding to the limits of the hyperbolic regime for the metasurfaces based on square arrays of the rectangular nanopatches. We show that the spectral difference between the resonances increases linearly with stretching of the nanopatch, but this dependence becomes quadratic when the length of the stretched nanopatch exceeds 85% of the lattice constant, indicating the regime of extreme anisotropy. Finally, we demonstrate the characteristic feature of the engineered resonances by showing the canalization (divergenceless propagation) of the surface plasmon-polariton along the anisotropic nanopatch-based metasurface in the vicinity of the resonance. The results obtained may be used for the engineering of the anisotropic nanoparticle-based metasurfaces for a plethora of photonic applications.

## I. INTRODUCTION

Two-dimensional (2D) materials offer unique opportunities for photodetection, light emission, energy harvesting, and enhanced light–matter interactions.^{1,2} Even more interest brings the artificially engineered 2D micro- and nanostructures with on-demand properties, such as surface conductivity or surface polarizability, paving the way towards a plethora of specific applications and devices including lensing, holography, imaging, polarimetry, biosensing, etc.^{3–5} The rapidly developing use of 2D nanostructures poses new challenges for their proper engineering.

^{6}In other words, the components of the diagonalized 2D tensor have different signs leading to the topological transition from elliptical to hyperbolic shapes of the isofrequency contours.

^{7–9}The 2D tensor of the surface conductivity is usually studied for the plasmonic metasurfaces,

^{10,11}

^{12,13}In the near-infrared and visible spectrum, HMSs are represented by metallic gratings

^{14}and lattices of nanoparticles.

^{15–17}

While the HMS have the clear advantages for many applications, the selection of their specific designs requires the extensive numerical simulation, which can be time-consuming and require significant effort. Some specific applications requiring the proper engineering of HMS include, for instance, the high-directional excitation of surface waves,^{18,19} twisted hyperbolic optics,^{20–22} spin–orbit photonics,^{23–25} and mode engineering for high-Q resonators^{26,27} and sensors.^{28,29} For all these applications, it is highly important to set two well-separated resonances of the surface conductivity, namely, to adjust their wavelength and amplitude independently. These resonances may be characterized by the peaks in the reflectance spectra^{17,30} or the plasmon canalization regime, which is the collinear divergenceless propagation of the surface plasmon-polariton along the specific direction^{12,31} (Fig. 1). The plasmon canalization has been previously observed in hyperbolic metasurfaces,^{10,11,32} natural anisotropic materials,^{33–35} and twisted hyperbolic layers.^{20–22}

In this work, we analyze the spectral positions of the resonances and the bandwidth of the hyperbolic regime for dozens of plasmonic metasurface designs based on gold nanopatches. Starting with nanopatches of a square cross section, we gradually stretch them into gratings, introducing anisotropy. Thus, a nanopatch-based design allows one to study the full dynamics of the transition between isotropic and extreme anisotropic cases for the HMS resonances. We show that the resonant wavelength for the case of the electric field oriented along the stretching direction and the spectral bandwidth between two resonances are initially linearly proportional to the stretching factor, but at a certain degree of the stretching (anisotropy), their dependence becomes quadratic. Then, we demonstrate the plasmon canalization in the vicinity of one of the resonances highlighting the relevance of the metasurface engineering for the in-plane optical signal transferring. Our results bring the significant contribution to the engineering of the resonances for the anisotropic 2D materials and nanostructures.

## II. METHODS AND DEFINITIONS

We consider the metasurface based on the periodic array of gold nanopatches with a height of $h=20$ nm and a rectangular cross section $( a x\xd7 a y)$, packed in a square array with a lattice constant of $p=300$ nm at the interface of a quartz substrate ( $ n s=1.45$). The dispersion of gold is taken from the experimental data for the thin films.^{36} Initially, we consider a plasmonic metasurface based on the nanopatches with a square cross section with a side of $a=150$ nm [the unit cell is shown in Fig. 2(a 1)]. Then, we stretch the nanopatches along the $x$-direction transforming the square cross section into the rectangular one [Figs. 2(a 2)–2(a 4)]. The total area of the cross section is saved during the stretching; thus, the long and short sides of a rectangle correspond to $ a x=a\eta $ and $ a y=a/\eta $, respectively, where $\eta $ is the stretching factor characterizing the degree of the anisotropy for a metasurface. For further analysis, we introduce the parameter of a filling factor along the direction of the stretching defined as $ f x= a x/p$, i.e., the long side of a nanopatch normalized per period. The unit cells, shown in Figs. 2(a 1)–2(a 4), correspond to the stretching factors $\eta =1,1.4,1.8,2$ and filling factors along the $x$-axis $ f x=0.5,0.7,0.9,1$, respectively. The first ( $ a x= a y=a$) and last ( $ a x=p$) unit cells, shown in Figs. 2(a 1) and 2(a 4), correspond to the isotropic (within effective local medium approach) and extreme anisotropic cases, respectively. Namely, the nanopatch-based metasurface in extreme anisotropic case represents a grating with a period $p$ and a width $ a y$.

^{30}

^{,}

^{17,30}

^{37}The reflectance and transmittance spectra of a normally incident plane wave with electric field oriented along $x$- and $y$-axis were calculated by the finite-element method with COMSOL Multiphysics for a single unit cell with imposed periodic boundary conditions [Figs. 2(b 1)–2(b 4)]. Then, we retrieve the corresponding components of the surface conductivity tensor $ \sigma x$ and $ \sigma y$ by fitting the calculated reflectance spectra using the least-squares method according to Eq. (4); see Figs. 2(c 1)–2(c 4). The effect of the finite thickness of a metasurface can be compensated by properly correcting the reflection phase.

^{30}

The corresponding electric field distributions at the resonances for all the cases from Figs. 2(a 1)–2(a 4) are shown in Fig. 3. Figures 3(a 1) and 3(b 1) correspond to the isotropic case (square nanopatch) at the wavelength $ \lambda 0=790$ nm and, as a result, are the same with respect to the rotation by $ 90 \xb0$. Figures 3(a 2) and 3(a 3) and 3(b 2)–3(b 4) correspond to the wavelengths $ \lambda x$ (1000, and 1345 nm) and $ \lambda y$ (685, 600 and 580 nm), respectively. The field distribution for the non-resonant case, when the electric field is directed along the grating at $\lambda =1500$ nm, is shown in Fig. 3(a 4).

## III. RESULTS AND DISCUSSION

We analyze the spectral positions of the resonances and the SDBR during the stretching of the nanopatches up to the extreme near-grating case characterized by $ f x=0.99$. First, we set the value of the period $p=300$ nm and the stretching factor $\eta $ to 1 (isotropic case), 1.2, 1.5, and 2 and change the side of a square nanopatch $a$ from 100 to 297 nm, 125 to 247.5 nm, from 100 to 198 nm, and from 75 to 148.5 nm, respectively, so that $ f x$ changes from 0.5 to 0.99 for all three cases [Figs. 4(a) and 4(d)]. Second, we set the value of the period $p=300$ nm and the side of a square nanopatch $a$ to 100, 150, and 180 nm, that correspond to the resonant wavelength $ \lambda 0=705,790$, and $835$ nm, and change the stretching factor $\eta $ from 1.5 to 2.97, from 1 to 1.98, and from 1 to 1.65, respectively [Figs. 4(b) and 4(e)]. For the last case, $a=180$ nm, $p=300$ nm, the filling factor along the stretching direction $ f x$ changes from 0.6 to 0.99. Third, we fix the long side of a rectangular nanopatch as $ a x=225$ nm via three cases $a=180$ nm and $\eta =1.25$, $a=150$ nm and $\eta =1.5$, $a=125$ nm and $\eta =1.8$, and change the period $p$ from 227.3 to 450 nm, so that $ f x$ changes from 0.5 to 0.99 for all three cases [Figs. 4(c) and 4(f)]. The latter three cases are compared with the isotropic case of $a= a x=150$ nm with the period varying from 151.5 to 300 nm [gray line in Fig. 4(c)]. The corresponding values of the resonant wavelengths $ \lambda x$, $ \lambda y$ and the SDBR $\Delta \lambda $ are shown in Figs. 4(a)–4(c) and 4(d)–4(f), respectively.

The calculations with square nanopatches (isotropic case) are shown by the gray lines in Figs. 4(a) and 4(c). One can notice that resonant wavelength $ \lambda 0$ increases linearly up to $ f x\u22480.85$ and then quadratically for $ f x>0.85$ for both types of the induced anisotropy. The quadratic dependence of $ \lambda x( f x)$ in the range $0.85< f x<0.99$ is associated with the strong mutual interaction between the nanoparticles leading to the substantial modification of the eigenmodes dispersion.^{30,38} The wavelength of the smaller resonance $ \lambda y$ slightly linearly redshifts [Fig. 4(a)] and blueshifts [Figs. 4(b) and 4(c)] for the cases of the fixed $p$ and $\eta $, $p$ and $a$, and $ a x$, respectively, with the increase of $ f x$. The wavelength of the larger resonance $ \lambda x$ increases linearly up to $ f x\u22480.85$ and then quadratically for $ f x>0.85$ for all three cases [Figs. 4(a)–4(c)] as well as the resonant wavelength $ \lambda 0$ in the isotropic case. It leads to the significant change of the SDBR from 145, 270, and 415 nm at $ f x=0.5$ to 300, 520, and 720 nm at $ f x=0.85$, and then to 965, 1155, and 1620 nm at $ f x=0.99$ for the case of the fixed period $p=300$ nm and the stretching factors $\eta =1.2,1.5$, and $2$, respectively [Fig. 4(d)]. For the case of the fixed period $p=300$ nm and the side of a square nanopatch $a=150$ nm, the SDBR changes from 0 nm at $ f x=0.5$ ( $\eta =1$, isotropic case) to 610 nm at $ f x=0.85$ ( $\eta =1.7$, anisotropic case) and then to 1620 nm at $ f x=0.99$ ( $\eta =1.98$, extreme anisotropic case) as it is shown by the red line in Fig. 4(e). While for the case of the fixed long sides of a rectangular nanopatch $ a x=225$ nm and $ a y=100$ nm, the SDBR changes from 360 nm at $ f x=0.5$ ( $p=450$ nm) to 480 nm at $ f x=0.85$ ( $p=264.7$ nm) and then to 1420 nm at $ f x=0.99$ ( $p=227.3$ nm) as it is shown by the red line in Fig. 4(f). The complete picture can be represented by the color map in Fig. 5 demonstrating the dependence of the SDBR on the stretching factor $\eta $ and the side of a square nanopatch $a$ at the period of $p=300$ nm. The strong hyperbolicity is observed within the region $0.85< f x<0.99$ and corresponds to the extremely anisotropic case. The transition from linear to quadratic dependence of the resonant wavelength and SDBR on the anisotropy degree is counter-intuitive and important to engineer properly the ultra-anisotropic systems and the plasmon canalization phenomena, where both resonances must be tuned simultaneously in extreme anisotropic case.

Finally, we demonstrate the typical near-resonant spatial field distribution of the surface plasmon-polariton excited by the vertical ( $z$-oriented) dipole and propagating in the canalization regime along the plasmonic nanopatch-based metasurface (Fig. 6). We take the metasurface consisting of $15\xd715$ nanopatches taken from Fig. 2(a 2), i.e., $p=300$ nm, $a=150$ nm, $\eta =1.4$. The canalization regime is observed when one of the surface conductivity components is large, while the orthogonal one is near-zero. Therefore, we draw the spectral dependence of $Im( \sigma x)/Im( \sigma y)$ characterizing the efficiency of plasmon canalization in Fig. 6(a) by the blue line. One can notice the small peak in the vicinity of $ \lambda y=700$ nm and the pronounced peak in the vicinity of $ \lambda x=1000$ nm. For the small peak, both components of the surface conductivity are near-zero meaning inefficient propagation with a divergence, while for the large peak, $ \sigma y$ is near-zero and $ \sigma x$ is resonant leading to the divergenceless propagation of the surface wave. Nevertheless, the most efficient plasmon canalization takes place for the wavelength slightly shifted from the resonance that is explained by the high absorption losses in the vicinity of the resonance characterized by the large real part of the surface conductivity component along the canalization direction, as it is shown by the red dashed line in Fig. 6(a). As a consequence, we choose the wavelength $\lambda =1100$ nm for the optimum plasmon canalization for this specific metasurface design.

The spatial distributions of the normal component of the magnetic field $ H z$ are shown in Figs. 6(b) and 6(d) explicitly demonstrating the collinear divergenceless propagation of surface plasmon-polariton with a small decay along the $x$-axis. The canalization regime is characterized by the flat isofrequency contour, which also indicates the topological transition between the elliptical and hyperbolic shapes of the isofrequency contours. In order, to extract the isofrequency contour, we apply the two-dimensional Fourier transform to the calculated distribution of the $ H z$-field [Fig. 6(c)] following the standard procedure.^{32,39} On the other hand, we calculated the corresponding isofrequency contour analytically [blue dashed line in Fig. 6(c)] using the dispersion equation for hyperbolic plasmon-polaritons at metasurfaces^{10,11} with the retrieved surface conductivity tensor, shown in Fig. 2(c 2). We claim a good coincidence between the isofrequency contours calculated analytically with the effective surface conductivity and extracted from the full-wave numerical simulations for the full-scale metasurface consisting of $15\xd715$ nanopatches.

The experimental implementation of the surface waves excitation in the hyperbolic and canalization regimes includes the loop and probe near-field technique in the microwave range,^{32} the prism coupling,^{15,30} the defect-induced structures,^{14,40} the metallic^{41} and dielectric nanoantennas,^{42,43} the metal tip of the microscope,^{33,41,44} and the edge-oriented launching^{44–46} in the visible and infrared ranges.

## IV. SUMMARY

We have investigated a number of plasmonic metasurface designs based on gold nanopatches by varying the side of the nanoparticle, its stretching (anisotropy) degree and the period of a metasurface from the isotropic (square nanopatches) to extreme anisotropic (gratings) cases. We define the quadratic dependence of the spectral bandwidth for one of the resonances on the anisotropy degree when the electric field is oriented along the stretching direction. It leads to the broadband spectral bandwidth between the resonances characterizing the hyperbolic regime of a metasurface. As a possible application verifying the eligibility of the discussed results, we demonstrated the plasmon canalization in the vicinity of the engineered resonance of the plasmonic metasurface.

The results obtained may be extended to different meta-atom geometries, lattice types, and materials (e.g., silver, aluminum). We believe that these results may be used by optical researchers and engineers for the efficient designing, engineering, and optimization of HMS plasmon resonances for a wide range of applications.

## ACKNOWLEDGMENTS

The numerical calculations have been conducted with the support of the U.S. Civilian Research and Development Foundation (CRDF Global) (No. G-202401-71609). A.H. thanks the IEEE Microwave Theory and Technology Society (IEEE MTT-S) for the support within the Undergraduate/Pre-Graduate Scholarship. S.P. and O.D. acknowledge the EURIZON Fellowship Programme “Remote Research Grants EU#3003-EURIZON” (Agreement No. 871072). O.Y. acknowledges the support of the Alexander von Humboldt Foundation within the framework of the Humboldt Research Fellowship for postdoctoral researchers.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**A. Hrinchenko:** Data curation (lead); Formal analysis (lead); Investigation (lead); Software (lead); Validation (lead); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). **S. Polevoy:** Formal analysis (equal); Software (equal); Visualization (equal); Writing – review & editing (equal). **O. Demianyk:** Software (equal); Visualization (equal); Writing – review & editing (equal). **O. Yermakov:** Conceptualization (equal); Formal analysis (equal); Methodology (lead); Supervision (lead); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## REFERENCES

*Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light*