Nonlocal effects in singular plasmonic metasurfaces

TL;DR

This paper demonstrates that nonlocal effects significantly modify the optical response of singular plasmonic metasurfaces, changing their absorption spectrum from continuous to discrete and reducing energy concentration at singularities.

Contribution

It reveals the impact of nonlocality on singular plasmonic metasurfaces, providing a more accurate model of their optical behavior compared to local theories.

Findings

Spectrum becomes discrete due to nonlocal effects

Energy concentration at singularities is reduced

Nonlocality alters absorption characteristics of metasurfaces

Abstract

A local model of the dielectric response of a metal predicts that singular surfaces, such as sharp-edged structures, have a continuous absorption spectrum and extreme concentration of energy at the singularity. Here we show that nonlocality drastically alters this picture: the spectrum is now discrete and energy concentration, though still substantial, is greatly reduced.

Figures

Equations

$$\begin{array}{rl}{\epsilon }_L^z(\omega ,\mathbf{k},z)& ={\epsilon }_L^{z^{{}^{\prime }}}(\omega ,\left|\frac{d{z}^{{}^{\prime }}}{dz}\right|\mathbf{k})\\ & ={\epsilon }_{\infty }-\frac{{\omega }_p^2}{\omega (\omega +i\Gamma )-\frac{{\beta }^2{d}^2|\mathbf{k}{|}^2}{T^2}{\left|\sinh \left(2\pi \frac{z}{d}\right)\right|}^2}\\ & ={\epsilon }_{\infty }-\frac{{\omega }_p^2}{\omega (\omega +i\Gamma )-{\beta }_{eff}^2(z)|\mathbf{k}{|}^2}\end{array}$$

$$\begin{array}{r}(\epsilon -1)\sqrt{k_x^2}(e^{\sqrt{k_x^2}(d_1+d_2)}+1)(e^{\sqrt{k_x^2}{d}_3}+1)\frac{e^{\kappa d_3}-1}{\kappa (e^{\kappa d_3}+1)}+((\epsilon -1)(e^{\sqrt{k_x^2}(d_1+d_2)}-e^{\sqrt{k_x^2}{d}_3})+(\epsilon +1)(e^{\sqrt{k_x^2}d}-1))=0\end{array}$$

$$\begin{array}{r}(\epsilon -1)\sqrt{k_x^2}(e^{\sqrt{k_x^2}(d_1+d_2)}-1)(e^{\sqrt{k_x^2}{d}_3}-1)\frac{e^{\kappa d_3}+1}{\kappa (e^{\kappa d_3}-1)}-((\epsilon -1)(e^{\sqrt{k_x^2}(d_1+d_2)}-e^{\sqrt{k_x^2}{d}_3})-(\epsilon +1)(e^{\sqrt{k_x^2}d}-1))=0\end{array}$$

$$\begin{array}{r}{\int }_{-L/2}^{L/2}{k}_{px}(x)dx+\varphi =2n\pi \end{array}$$

$$H_{z}(k_{x},y)=\left\{{\begin{array}[]{lr}(1-r)a_{a}\frac{e^{-|k_{x}||y|}}{|k_{x}|}+(1+r)a_{s}\text{sgn}(k_{x})\frac{e^{-|k_{x}||y|}}{\text{sgn}(y)|k_{x}|}+{b_{+}}e^{-|k_{x}|y}+{b_{-}}e^{|k_{x}|y},&{-d_{2}<y<d_{1}}\\ -t\frac{k_{0x}^{{}^{\prime}}}{k_{0x}}a_{a}\frac{e^{-|k_{x}||y+\frac{d}{2}|}}{|k_{x}|}-ta_{s}\text{sgn}(k_{x})\frac{e^{-|k_{x}||y+\frac{d}{2}|}}{\text{sgn}(y+\frac{d}{2})|k_{x}|}+{c_{+}}e^{-|k_{x}|y}+{c_{-}}e^{|k_{x}|y},&{-(d_{2}+d_{3})<y<-d_{2}}\end{array}}\right.$$

$$\varphi(k_{x},y)=\left\{{\begin{array}[]{lr}0,&{-d_{2}<y<d_{1}}\\ \frac{\mathrm{sgn}(k_{x})}{\omega\varepsilon_{0}\varepsilon}({d_{+}}e^{-\kappa y}-{d_{-}}e^{\kappa y}),&{-(d_{2}+d_{3})<y<-d_{2}}\end{array}}\right.$$

$$\begin{array}{rl}{E}_x(k_x,y)& =\frac{i}{\omega \epsilon }\frac{\partial H_z}{\partial y}-\frac{\partial \phi }{\partial x}\\ E_y(k_x,y)& =-\frac{i}{\omega \epsilon }\frac{\partial H_z}{\partial x}-\frac{\partial \phi }{\partial y}\end{array}.$$

$$H_{z}^{a}(x,y)=\left\{{\begin{array}[]{lr}\begin{split}&i2\pi a_{a}({\Gamma_{a+}}e^{-\sqrt{k_{px}^{2}}y}+{\Gamma_{a-}}e^{\sqrt{k_{px}^{2}}y})\\ &\times(e^{i\int_{0}^{|x|}k_{px}(x^{\prime})dx^{\prime}}+e^{-i\int_{0}^{|x|}k_{px}(x^{\prime})dx^{\prime}+i(\phi_{0}+\phi)})\frac{1}{1-e^{i(\phi_{0}+\phi)}}\end{split},&{-d_{2}<y<d_{1}}\\ \begin{split}&i2\pi a_{a}({\Lambda_{a+}}e^{-\sqrt{k_{px}^{2}}y}+{\Lambda_{a-}}e^{\sqrt{k_{px}^{2}}y})\\ &\times(e^{i\int_{0}^{|x|}k_{px}(x^{\prime})dx^{\prime}}+e^{-i\int_{0}^{|x|}k_{px}(x^{\prime})dx^{\prime}+i(\phi_{0}+\phi)})\frac{1}{1-e^{i(\phi_{0}+\phi)}}\end{split},&{-(d_{2}+d_{3})<y<-d_{2}}\end{array}}\right.$$

$$\varphi^{a}(x,y)=\left\{{\begin{array}[]{lr}0,&{-d_{2}<y<d_{1}}\\ \begin{split}&\frac{i2\pi a_{a}\mathrm{sgn}(k_{px})\mathrm{sgn}(x)}{\omega\varepsilon_{0}\varepsilon}({\Omega_{a+}}e^{-\kappa_{p}y}-{\Omega_{a-}}e^{\kappa_{p}y})\\ &\times(e^{i\int_{0}^{|x|}k_{px}(x^{\prime})dx^{\prime}}-e^{-i\int_{0}^{|x|}k_{px}(x^{\prime})dx^{\prime}+i(\phi_{0}+\phi)})\frac{1}{1-e^{i(\phi_{0}+\phi)}}\end{split},&{-(d_{2}+d_{3})<y<-d_{2}}\end{array}}\right.$$

$$P_{abs}=2{\int }_{-(d_2+d_3)}^{d_1}{S}_{x}dy$$

$$\begin{array}{r}r=\frac{\sqrt{ϵ}-1+\frac{{\sigma }_e}{{\sigma }_{e0}}}{\sqrt{ϵ}+1+\frac{{\sigma }_e}{{\sigma }_{e0}}}\end{array}.$$

$$\begin{array}{rl}& \nabla \times \nabla \times \mathbf{E}=k_0^2\mathbf{E}+i\omega {\mu }_0\mathbf{J}\\ & {\beta }^2\nabla (\nabla \cdot \mathbf{J})+\omega (\omega +i\Gamma )\mathbf{J}=i\omega {\omega }_p^2{\epsilon }_0\mathbf{E}\end{array}$$

Full text

Imperial College London] The Blackett Laboratory, Department of Physics, Imperial College London, London SW7 2AZ, United Kingdom

Imperial College London] The Blackett Laboratory, Department of Physics, Imperial College London, London SW7 2AZ, United Kingdom

Imperial College London] The Blackett Laboratory, Department of Physics, Imperial College London, London SW7 2AZ, United Kingdom

Imperial College London] The Blackett Laboratory, Department of Physics, Imperial College London, London SW7 2AZ, United Kingdom

\abbreviations

1 Introduction

In the past decade, rapid development in nanofabrication has prompted great interest in plasmonic systems which collect and concentrate light into subwavelength volumes 1. Among these systems, singular plasmonic structures, such as narrow gaps and sharp edges 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, can even concentrate the field down to the Coulomb screening length where electron-electron interaction should be considered 14, 15, 16, 17, 18. The hydrodynamical model 19, 20, 21 takes this interaction into account and successfully explains nonlocal effects in the gap of a dimer 22, 23, 24.

In our recent work25, 26, we proposed and studied singular metasurfaces within a local description where the optical response of metals is described by a $\omega$-dependent permittivity such as the Drude model, $\varepsilon(\omega)=\varepsilon_{\infty}-\frac{\omega_{p}^{2}}{\omega(\omega+i\Gamma)}$. However, the local model assumes that the electron gas is a continuum taking no account of the Fermi surface and the finite density of electrons. In our singular metasurface, the SPPs will propagate to the singularities, where electrons accumulate and reach infinite density, which leads to divergences of the electric field. However, the infinite density is not physical because of the discrete nature of the electron gas, which determines the screening length ($\delta_{C}\sim 0.1$ nm, for noble metals) 27 and prevents the electron density from blowing up. In this work, we take this nonlocal effect into account for singular metasurfaces where the scale of the singularity could go below the screening length. Our results show that the singular metasurface is very sensitive to the nonlocal effects of the electron gas and that nonlocality in the metal has to be considered for an accurate description of our singular system. We discuss how singular metasurfaces could be used to reveal nonlocal effects.

2 Results and discussion

Fig. 1(a) is a schematic of our singular metasurface: a periodic array of grooves on a metal surface with sharp edges illuminated by a plane wave. In the presence of nonlocality, we have both transverse and longitudinal modes in the metal. For the transverse mode, the dielectric function is modeled with classical Drude permittivity with $\varepsilon_{\infty}=1$, $\omega_{p}=8.95$ eV/$\hbar$ and $\Gamma=65.8$ meV/$\hbar$ 28, while for the longitudinal mode, the dielectric function is described within the hydrodynamic model as $\varepsilon_{L}(\omega,\mathbf{k})=\varepsilon_{\infty}-\frac{\omega_{p}^{2}}{\omega(\omega+i\Gamma)-\beta^{2}|\mathbf{k}|^{2}}$ 19, where $\beta$ is the nonlocal parameter, measuring the degree of nonlocality. The decay length of the longitudinal mode is $\delta=\frac{\beta}{\sqrt{\omega_{p}^{2}-\omega^{2}}}\approx\frac{\beta}{\omega_{p}}$ when the frequency is well below the plasmon frequency. In the presence of nonlocality, the surface charge becomes a volume charge with a decay length $\delta$ inside the metal surface, which is depicted as an electron density layer (purple layer) on top of the metasurface with period $T$ shown in Fig. 1(a). The existence of this layer smooths the singularity and has remarkable effects on the far field spectrum.

The complex boundary of a singular metasurface complicates the calculation, so we shall map this singular surface into a slab array with period $d$ (see Fig. 1(b)), where the thickness of the slabs is $d_{3}$ and the thickness of the dielectric region is $d_{1}+d_{2}$ 26. In the slab frame, the transverse electric permittivity remains the same since it is $k$-independent and conformal mapping conserves the permittivity in the $x-y$ plane. However, the dependence of the longitudinal permittivity on the $k$ vector complicates the problem because $k$ is not invariant under conformal mapping. Therefore, the longitudinal permittivity in the slab frame becomes coordinate dependent, and is written as

[TABLE]

where the primed complex coordinate $z^{{}^{\prime}}=x^{{}^{\prime}}+iy^{{}^{\prime}}$ is for the metasurface frame, while the unprimed coordinate $z=x+iy$ is for the slab frame. The effective nonlocal parameter $\beta_{eff}$ in the slab frame, $\beta_{eff}(z)=\frac{\beta d}{T}\left|\sinh\left(2\pi\frac{z}{d}\right)\right|$, indicates that the decay length is smaller near the origin and increases when $x\rightarrow\pm\infty$ (see Fig. 1(b)), i.e. the permittivity is more and more nonlocal along the slab.

We start by calculating the SPP dispersion in the slab frame. Plane waves (incident, reflected and transmitted waves) in the metasurface frame are modeled as monopole sources in the slab frame and a non-retarded approximation is taken 26. These three sorts of waves excite the SPPs along the slab. To determine the coefficients of this SPP mode, the tangential component of the magnetic and electric fields ($H_{z}$ and $E_{x}$ in our coordinate system) should be matched at the interface between metal and air. However, imposing the continuity of $H_{z}$ and $E_{x}$ is not sufficient as there exists an additional longitudinal mode. To solve this indeterminacy, an additional boundary condition, the continuity of $E_{y}$, is imposed at the interface between metal and air, where $E_{y}$ is the electric field component normal to the interface. This continuity comes from the assumption of no surface charge in the hydrodynamic model. By imposing these boundary conditions, we can calculate the mode in the slab frame. Since $\beta_{eff}$ is a function of the slab frame coordinate, $k_{x}$ varies along the slab. To make an analytic solution possible for this inhomogeneity, the WKB approximation is introduced. The WKB approximation, which has been successfully used in plasmonic systems22, 23, 29, applies when the phase of the SPP wave changes more quickly than its amplitude. With the above additional boundary condition and the WKB approximation, the SPP mode and the corresponding dispersion relation are determined. The dispersion relation for anti-symmetric mode is

[TABLE]

while the dispersion relation for the symmetric mode is given by

[TABLE]

where $\kappa=\sqrt{k_{x}^{2}+\frac{\omega_{p}^{2}}{\frac{\beta^{2}d^{2}}{T^{2}}\left|\sinh\left(2\pi\frac{z}{d}\right)\right|^{2}}\frac{\varepsilon}{\varepsilon-1}}$ and $k_{px}$ stands for the plasmon pole. The symmetry of the SPP mode is defined by the parity of $E_{x}(y)$ in the slab frame: The odd(even) function corresponds to the anti-symmetric(symmetric) mode. The dispersion relation at the center of the slab $k_{px}(\omega,x=0)$ is presented in Fig. 1(c), together with the local result for comparison 26. For the local case, both anti-symmetric and symmetric modes asymptotically approach $\omega=\omega_{sp}$. In contrast, within the nonlocal description, this SPP mode asymptotically approaches the longitudinal bulk mode $\omega=\sqrt{\beta^{2}k^{2}+\omega_{p}^{2}}$. In the large $k$ limit, $\omega\approx\beta k$, resulting in a linear dispersion relation, see Fig. 1(c). In the following, we will focus on the anti-symmetric band which can be excited by a normally incident plane wave. For this anti-symmetric mode, the profile of $k_{px}^{2}$ on the interface of the slab at a representative frequency within the band $\omega=0.6\omega_{p}$ is shown in Fig. 1(d), where the point when $k_{px}^{2}=0$ is called a "turning point" 30. At the turning point, the SPP wave is reflected and the two turning points ($x=-L/2$, $x=L/2$) confine SPP waves to form a Fabry-Perot resonant cavity with width $L$.

With the mode profile in the slab frame (see Methods section for detailed calculations), we can calculate the power flow at $x=0$ and model this power flow as the energy consumption by an effective surface conductivity, which has been previously introduced as flat surface model 26. This surface conductivity gives us the reflection spectrum. In all the following calculations, the parameter $d_{3}=0.9d$ is used, fixing the shape of the metasurface.

In Fig. 2 we show how nonlocality greatly affects the optical response of singular metasurfaces. Here we present the reflection spectrum of a groove metasurface with period $T=10$ nm for different nonlocal parameters $\beta$. In Fig. 2(a), a realistic value of $\beta=1.27\times 10^{6}$ m/s is used 21. The theoretical calculation (blue solid line) agrees very well with the simulation results (red dashed line) from the commercial finite element method solver COMSOL MULTIPHYSICS, which are also compared with a local calculation (gray line). It is clear that the spectrum becomes discrete when nonlocality is introduced. To further clarify the effect of nonlocality on the spectrum, we conduct a series of calculations with different $\beta$. When $\beta$ is 10 times smaller ($\beta=1.27\times 10^{5}$ m/s), more resonances appear in the spectrum, shown in Fig. 2(b), in which the analytic result also agrees well with numerics. We can reduce $\beta$ further to $1.27\times 10^{4}$ m/s (Fig. 2(c)), and even to $1.27$ m/s(Fig. 2(d)) where Comsol simulations do not converge. The nonlocal spectrum supports more and more peaks as $\beta$ is reduced. In the limit $\beta\rightarrow 0$, the spectrum will become continuous. The origin of the spectrum discretization is that nonlocality blunts the singularity. The resonance condition for different peaks is

[TABLE]

where $n$ is an integer and $\phi$ is the phase change at the turning point which is discussed in Methods. The finite range of integration yields a finite number of resonances. As $\beta$ decreases, the turning point is further from the origin and the $k$ integration is larger, thereby allowing for more peaks and appreciably changing the spectrum. In a nutshell, this microscopic feature of nonlocality greatly affects the optical response of the singular metasurface.

Next, we study the effect of the other length scale of the problem, the metasurface period, while keeping a realistic nonlocal parameter ($\beta=1.27\times 10^{6}$ m/s). In Fig. 3(a) we reproduce the reflection spectrum for 10 nm period where we have already seen an excellent agreement with full electrodynamics simulations as in this regime the quasi-static approximation is well satisfied. When the period is increased to 50 nm (Fig. 3(b)), more resonances are excited and at the same time the spectrum moves to lower frequencies. The reason for this is that when the period is increased the length scale where nonlocality is relevant is effectively smaller. Hence, increasing the period has a similar effect to reducing $\beta$ since these are the two length scales which compete in this problem.

Then we further increase the period to 100 nm, where the quasi-static approximation is not expected to be accurate. However, there is still a good agreement between our theory and simulations. In particular, the linewidth is not broadened in comparison with $T=10$ nm and $T=50$ nm, except for the low order resonance peaks whose linewidth gets broadened slightly. This is a remarkable difference from localized plasmonic structures of sizes comparable to the period considered here, where radiative broadening is important and a theory beyond quasi-statics is needed to accurately describe these systems22. Different from localized plasmonic structures, the periodic singular metasurface does not have a strong radiative broadening because by increasing the period we have less periods within overall length. The first resonance peak which has the highest broadening forms a stronger dipole moment, so its broadening is larger than the higher order peaks whose electron distribution is more homogeneous.

After studying the far field spectrum, we turn to the near field profile in the presence of nonlocality. Fig. 4 illustrates the electric field profile of the metasurface, in which the first and the second column depict the electric field ($E_{x}$ and $E_{y}$) distribution in one period ($-T/2<y<T/2$). Also, the electric field $E_{n}$ (the component normal to the interface, which dominates for the anti-symmetric mode) along the surface is presented in the last two columns with different scales in the $y$-axis. Fig. 4(a) shows the electric field profile for the second resonance peak in the reflection spectrum of Fig. 2(a) (period $T=10$ nm, $\beta=1.27\times 10^{6}$ m/s), at a frequency of $\omega=0.7\omega_{p}$. From this field profile, we see that the electric field is continuous across the interface because the hydrodynamic model assumes no surface charge. On the contrary, the existence of a surface charge in the local description makes the normal component of electric field discontinuous. On the other hand, the field near the singularity is not as localized as the field in the local case. This is because the smearing of the singularity introduced by nonlocality results in a weaker compression of the plasmon wavelength as it travels towards the singularity. This effect causes a strong reduction of the field enhancement in the vicinity of singularity. In fact, while the electric field diverges at the singularity in the local case, it has a value of approximately 70 when nonlocality is taken into account, as shown by the plot of the field at the interface. In addition, the agreement between analytical and numerical results is excellent except near the singularity. This discrepancy between theory and Comsol is due to the WKB approximation, which does not work well when the $k$-vector is small.

To further unveil the effect of nonlocality on the field enhancement, in Fig. 4(b) we consider a metasurface with the same period as in Fig. 4(a) but with a smaller nonlocal parameter ($\beta=1.27\times 10^{5}$ m/s). We plot the fields of the fifth resonance peak in Fig. 2(b), at the same frequency $\omega=0.7\omega_{p}$ as Fig. 4(a). Comparing to the case with a realistic description of nonlocality (panel (a)), the field in this case shows more oscillations and presents a sharper change at the metal surface, looking more alike to the local situation. In addition, the plot of the field along the interface demonstrates how the field enhancement strongly increases as compared to panel (a), which is consistent with a divergent field in the limit of $\beta\rightarrow 0$.

Finally, we also consider the near field profiles for a metasurface with a long period ($T=100$ nm) and nonlocal metal (realistic $\beta$). The field profile for the fifth-order resonance in the metasurface (the fifth peak in spectrum in Fig. 3(c), $\omega=0.7\omega_{p}$), is shown in Fig. 4(c). Similar to decreasing the value of $\beta$ for a fixed period, keeping $\beta$ and increasing the period effectively reduces the length scale of nonlocality. Hence, the fields have more oscillations and do not smooth out over the metal/dielectric interface but present jumps. Besides, the electric field is larger than that for the 10 nm case in Fig. 4(a), since the effectively reduced nonlocal length scale implies that the plasmons have more time to travel towards the singularity while their wavelength is compressed, and confinement increases. Hence, higher field enhancements are to be found in metasurfaces of larger periods, which are surprisingly unaffected by radiative broadening.

3 Conclusions

In this paper we have explored the consequences of nonlocality for surfaces in the form of metasurfaces containing sharp edge singularities. Local theory predicts that external radiation will excite a continuous spectrum of modes with infinite energy density at the singularity. However, nonlocality forbids infinite concentration of charge at the surface with dramatic consequences for the spectra which are now discrete, and for the energy density which is no longer singular, although it takes a large value at the structural singularity. In effect nonlocality blunts sharp edges. This result suggests a way of measuring nonlocality, particularly by observing the mode spacing in the discrete spectrum it implies. When the metasurface period is very short compared to the free space wavelength of light, magnetism plays little part in the modes which are almost entirely electrostatic. In a local theory the spectrum is scale invariant hence the period of the metasurface does not affect the spectrum. This ceases to be true when nonlocality introduces a length scale in the form of the screening length. As a result, there is a balance between the nonlocal parameter, $\beta$, and the metasurface period: by increasing the period we can compensate for nonlocality at least until the period approaches the free space wavelength. Another point we noted was the insensitivity of radiative damping of the modes to the metasurface period. The charge currents generated by a mode couple weakly to external radiation and increasing the period might be expected to increase this coupling. For an isolated resonant particle this is indeed the case and radiative damping increases dramatically once the particle is much bigger than about 100 nm. In contrast, for a given length of metasurface, increasing the period means including less periods in the overall length hence the small effect.

4 Methods

4.1 Field in the slab frame

Since the conformal mapping conserves the spectrum of the system, we carry out our analytic calculation in the slab frame where the geometry is simpler. In the slab frame, we can write the field in each region by taking the quasi-static approximation as

[TABLE]

[TABLE]

where $H_{z}$ contributes to the transverse mode, and $\varphi$ to the longitudinal mode. Besides, the anti-symmetric source magnitude is $a_{a}=-i\frac{k_{0x}T}{4\pi}H_{0}$, while the symmetric one is $a_{s}=\frac{k_{0y}T}{4\pi}H_{0}$. The electric components can be obtained by

[TABLE]

Then by using the continuity of $H_{z}$, $E_{x}$ and $E_{y}$ at interface between metal and air, we can calculate the coefficients ($b_{\pm}$, $c_{\pm}$ and $d_{\pm}$) for the SPP mode in $k$-space. The anti-symmetric excitation gives us a mode with anti-symmetry, the dispersion relation for which is given in Eq. 2. In contrast, the symmetric excitation gives a symmetric mode, whose dispersion relation is given by Eq. 3. For the anti-symmetric mode studied in this paper, the field distribution in real space ($H_{z}(x,y)$ and $\varphi(x,y)$) can be obtained by taking a Fourier transformation of the field in k-space($H_{z}(k_{x},y)$ and $\varphi(k_{x},y)$). The field distribution in real space can be written as

[TABLE]

[TABLE]

where $\phi_{0}=\int_{-L/2}^{L/2}k_{px}(x^{\prime})dx^{\prime}$ and $\phi$ is the phase change at the singularity which will be discussed in the next subsection.

4.2 Phase change at the singularity

The WKB approximation fails near the turning point as $k$ is small30. Therefore, we cannot use $k$ from the dispersion relation based on WKB to calculate the phase change $\phi$. Instead, the phase change at the tip is calculated in the metasurface frame. The result in Fig. 5(a) shows that the $k$-vector near the tip saturates as the width of the cavity $\delta$ decreases. The phase change $\phi$ of this cavity array can be obtained by matching the field at the terminus26, 31, 32. The calculated results show that both $k$-vector and phase change saturate near the singularity. Therefore, we use this value of the phase $\phi$ in the calculation.

4.3 Reflection spectrum

In order to calculate the reflection, we need to calculate the absorption in the slab frame. Our trick is to evaluate the power flow at the excitation point ($x=0$) by

[TABLE]

where $S_{x}=\frac{1}{2}\mathrm{Re}[E_{y}^{*}H_{z}]$. The absorption by the slab is modeled as a real surface conductivity $\sigma_{er}$. Then using Kramers-Kronig relation 26, 33, 34, 27, a causal complex surface conductivity $\sigma_{e}=\sigma_{er}+i\sigma_{ei}$ gives the reflection spectrum of singular metasurface by

[TABLE]

4.4 Comsol modelling

All of our numerical simulations are based on the RF and PDE modules in Comsol35, 20, where the hydrodynamic system of equations in the metal are implemented as

[TABLE]

where the $\mathbf{J}$ is the current density for the electron. When $\beta\rightarrow 0$, we arrive at classical Drude local-response function $\varepsilon(\omega)=1-\omega_{p}^{2}/(\omega(\omega+i\Gamma))$. By solving the above coupled equations, the response of our singular metasurface can be obtained. {acknowledgement} The authors thank A. I. Fernandez-Dominguez, Y. Luo and C. Ciraci for fruitful discussion. F.Y. acknowledges a Lee Family Scholarship for financial support. Y.T.W. acknowledges funding from the Leverhulme Trust. P.A.H and J.B.P. acknowledge funding from the Gordon and Betty Moore Foundation.