Electron Spill-Out Effect in Singular Metasurfaces

: The electron spill-out effect is considered in a singular metasurface. Using the hydrodynamic model, we found that electron spill-out effectively smears the sharp singularity. The introduction of the electron spill-out effect also signiﬁcantly changes the reﬂection spectrum, charge distribution, ﬁeld proﬁle for a singular metasurface. Therefore, this spill-out contribution is crucial and cannot be ignored for a realistic description of optical response in a singular system.


Introduction
Plasmonic metasurfaces [1] have been a hot topic in the past decade and have been used in many applications, such as wave-front control [2,3], wave absorbing [4,5], imaging [6,7], and so forth. As a subset of plasmonic metasurfaces, singular metasurfaces [8,9] have one additional hidden dimension. In contrast, to conventional metasurfaces which require specific incident wave-vectors for surface plasmon excitation, singular metasurfaces give rise to excitations at arbitrary incident wave vectors. Because of this unique property, singular metasurfaces can be used as broadband electromagnetic absorbers, an example being a black gold surface [10].
A singular metasurface has sharp points or narrow gaps, which are called "singularities". These singularities are so tiny that a macroscopic permittivity cannot accurately describe the local optical properties because of the quantum nature of electrons. Thus, the quantum effect of electrons should be considered near a singularity. These quantum effects mainly include nonlocality [11] and electron spill-out [12]. Nonlocality is attributed to the wave nature of electrons combined with the Pauli exclusion principle, and point like charges are surrounded in the electron gas not by a point like screening charge but by a finite cloud of electrons whose extent is determined by the Thomas-Fermi screening length (typically less than one angstrom) [13]. This limitation sets the bound for the electric field enhancement. Besides, the electron spill-out effect means that electrons can leak out of the metal surface, where the length scale of spill-out is characterized by Feibelman's d parameters (typically a few angstroms) [14]. Such a small length scale makes the spill-out effect observable only in a nanometer-size metallic particle [15] or a sub-nanometer gap [16].
To fully account for these quantum effects in plasmonic systems, a time-dependent density functional theory (TD-DFT) [17][18][19] approach is required. However, this DFT approach is quite time-consuming, especially for a large-number electron system. Luckily, alternative methods, such as the hydrodynamic model, can handle these quantum effects in a semi-classical way, avoiding solving many-body Schrödinger equations. The simplest model is the hard-wall hydrodynamic model (HW-HD) [20,21], which incorporates the nonlocal contribution from electrons. More recently, a more advanced model that considers both nonlocality and spill-out is proposed, named as the self-consistent hydrodynamic model (SC-HD) [22][23][24][25].
In our previous work on singular plasmonic systems, we have only incorporated the contribution of nonlocality, which shows that the nonlocal effect strongly changes the spectrum of the system. A continuous spectrum in the local case becomes discrete once the nonlocal contribution of the electron is taken into account [26][27][28]. However, the electron spill-out effect was ignored in these calculations, and its role in the singular system is still unknown. In this paper, we utilize two models (SC-HD and HW-HD) in the calculation of the optical response of singular metasurfaces and compare them with classical local response approximation (LRA). For simplicity, a simple metal without inter-band transitions is studied to clarify the role of electron spill-out from the singular metasurface.

Methods
The schematic for this paper is shown in Figure 1. We consider both groove (upper one in the left panel) and wedge (lower one in the left panel) singular metasurfaces, which are parameterized with the period T and vertex angle θ [9]. On the right panel of Figure 1, we conceptually illustrate the electron density distribution near the metal interface ∂Ω. The total electron density n can be decomposed into the ground state n 0 (r) and the excited state n 1 (r), where n 0 corresponds to equilibrium electron density in the absence of field excitation, while n 1 to electron density induced by an incident field [29]. For SC-HD, the normalized electron density of the ground state (n 0 normalized by ion density n ion ) continuously transitions from 1 to 0 across the boundary ∂Ω, where the electron density exhibits some oscillations termed Friedel oscillations [29]. Likewise, the induced electron density, n 1 , spills out of the metal surface and strongly oscillates near the interface. In comparison, the equilibrium electron density in both HW-HD and LRA becomes a step function [22], which means n 0 is discontinuous at ∂Ω. Moreover, the induced density for HW-HD only oscillates inside the metal, while for LRA n 1 becomes a delta function at ∂Ω. In the LRA, there is no constraint on the electron density so that it becomes infinite. Once imposing a constraint, such as that in SC-HD and HW-HD, the induced electron density becomes finite.
For the electrons in the metal, the equation of the motion can be expressed as [22][23][24][25] m e n( ∂v ∂t together with the continuity equation for electron density, n ∂n ∂t where m e , −e, and v are the mass, the charge, and the velocity of the electron, respectively. Note that the Γ in the Equation (1) is the damping term. The functional G[n(r, t)] represents the internal energy of the electron gas, which includes internal kinetic energy and the exchange-correlation energy [18]. The detailed expression of the functional derivative of G[n] is given in Yan's work [23]. The electron density in the above system of equations can be written as n(r) = n 0 (r) + n 1 (r), where the induced electron density n 1 is treated as a perturbation. The electron at equilibrium is also assumed to be motionless so that v = v 1 , and hence the magnetic field appears when the system is perturbed, expressed as B = B 1 . However, the static electrons support an electric field E 0 , so the electric field is formulated as E = E 0 + E 1 . Likewise, we also calculate the perturbation for the derivative of the energy functional G[n], written as δG δn = δG δn 0 + δG δn 1 [22].
Using perturbation theory, the equation system Equations (1) and (2) can be divided into two sets of equations: zero-order ground state and first-order excited state. The zeroorder equation system is given as [22][23][24][25] which collectively determine the ground state charge density of plasmonic system. Note that in the ground state calculation, the jellium model is used for the ion, in which the ion density is uniform inside the interface ∂Ω but 0 outside [29]. In equation system above, n + is the jellium background density such that n + = n ion , µ is the chemical potential, and φ 0 is the electrostatic potential defined by E 0 = −∇φ 0 . In contrast, the first-order system of equations is written as [22][23][24][25] which couples with Maxwell's equation where the induced current J 1 = −en 0 v 1 and the induced charge density ρ 1 = −en 1 . E 1 is the total electric field including incident and scattered fields. Equations (6)-(8) together determine the excited state.
The difference between SC-HD, HW-HD and LRA mainly comes from the functional derivative δG/δn. For SC-HD, this functional derivative includes contributions from the Thomas-Fermi kinetic energy, the von Weizsäcker(VW) energy and the exchange-correlation energy [30,31]. When only Thomas-Fermi kinetic energy is considered, we arrive at HW-HD. Thomas-Fermi kinetic energy only accounts for the homogeneous electron gas. To account for the inhomogeneity of electron gas, we need to consider high-order terms, where von Weizsäcker(VW) energy is the second order correction. Note that von Weizsäcker(VW) energy contains derivative terms like ∇n, which becomes important and indispensable near the interface where the electron density, n, oscillates strongly. For the LRA, the functional derivative is completely ignored, so that everything becomes local.
The ground state and excited state equations can be implemented with the finiteelement method in the commercial software COMSOL Multiphysics. We first numerically solve Equations (3)-(5) to obtain the ground-state electron distribution n 0 . After that, a second calculation solving Equations (6)- (8) gives the excited state, including reflection spectrum, induced charge distribution, electric field, and so forth. The detailed implementation in Comsol can be found in Ding's work [25].

Results
In this section, we employ the SC-HD, HW-HD and LRA models to calculate the optical response of two singular metasurfaces in Figure 1. In the calculation, we consider the case of a simple metal, like sodium. The corresponding ion density is n + = 2.5173 × 10 28 m −3 , from which we obtain the plasma frequency ω p = e 2 n ion /m e ε 0 ≈ 5.89 eV/h [32]. Besides, the damping parameter is Γ = 0.17 eV/h [22]. The geometric setup for the singular metasurface is T = 10 nm and θ = 0.6π rad. Excitation is by a plane wave incident normally on the singular metasurface.
The ground state calculation is shown in Figure 2. The normalized equilibrium electron densities n 0 /n ion of groove and wedge singular metasurfaces are presented in Figure 2a,b, respectively. It is obvious that the electron has some distribution outside the metal surface (the black curve), which is the spill-out effect we are concerned about. Because of this electron spill-out, the region immediately outside ∂Ω has a negative net charge, while the region immediately inside has a positive net charge, which contributes to a dipole layer in the ground state. Despite geometric singularities in the jellium background n + , the ground state electron smears out this sharp point, making it as far as the electron are concerned an effectively blunt singular metasurface [9]. We have also marked the length scale in the figure, which shows this blunt tip has an around 2 Å diameter. In contrast, the ground states for both HW-HD and LRA possess a step-function distribution rather than a continuous transition in the case of SC-HD. With ground-state electron distribution determined, we can next calculate the excited state. We compare the reflection spectrum for SC-HD, HW-HD, and LRA cases in Figure 3, where Figure 3a,b correspond to groove and wedge singular metasurfaces, respectively. For LRA, we have previously demonstrated that the reflection is a continuous spectrum [9]. For a groove singular metasurface, the lower band Lg (ω c 1 < ω < ω sp ) is exited while the upper band (ω sp < ω < ω c 2 ) is a dark mode under normal incidence, where the notation "Lg" means the LRA mode for groove singular metasurface ("L" is short for LRA and "g" stands for groove), ω sp = ω p / √ 2 is the surface plasmon frequency and ω c 1,2 is the cut-off frequency. On the contrary, the wedge singular metasurface only supports upper-band mode Lw, leaving the lower band as a dark mode. However, when the Thomas-Fermi kinetic energy is considered, the reflection spectrum ceases to be continuous but instead has a discrete spectrum [28]. In both Figure 3a,b, a broadband continuous spectrum in LRA becomes a few discrete peaks in HW-HD, that is, Hg1, Hw1 and Hw2. The notation rule is the same as that in LRA except for the additional number which defines the discrete mode index. Now, we take the electron spill-out effect into account, which requires a full consideration of all energy contributions in functional derivative δG/δn, as shown in SC-HD case in Figure 3. The reflection spectra of SC-HD exhibit distinct features comparing with other two cases. For the groove singular metasurface, we observe that the Sg1 and Sg2 peaks tend to redshift relative to the Lg band, while the Hg1 peak relatively blueshifts. Similar peak shifting can be observed in the wedge singular metasurface.
As a free parameter, θ characterizes the shape of the singular metasurface. In Figure 3c,d, we implement the same calculation as in Figure 3a,b except that the vertex angles for the corresponding groove and wedge singular metasurfaces are θ = 0.8π rad. For LRA, the bandwidth of both groove and wedge singular metasurfaces become narrower when the vertex angle is increased from 0.6π to 0.8π. From the reflection spectrum of the LRA, we observe that the shifting of the cut-off frequency effectively blueshifts the central frequency of band Lg for the groove singular metasurface, while redshifting the central frequency of band Lw for the wedge singular metasurface. In Figure 3e,f, we further increase the vertex angle to 0.9π, making the singular surface close to a flat surface (θ = π). The spectra at θ = 0.9π for SC-HD, HW-HD and LRA all show weak resonances and behave like a flat metal surface.
For SC-HD and HW-HD, when increasing the vertex angle the discrete peaks (Sg1, Sg2, Sg3 and Hg1) for the groove singular metasurface blueshift, while the peaks (Sw1, Sw2, Hw1 and Hw2) for the wedge singular metasurface redshift. Specifically, for the groove singular metasurface, the blue-shift of peak Hg1 from Figure 3a-c is around 0.04ω p , while for peak Sg1 is approximately 0.03ω p . Therefore, changing the vertex angle only overall moves the spectrum. However, the electron spill-out contribution not only shifts the peak position but also strongly changes the spectrum features, by creating some additional peaks (Sg3 or Sw2).
From reflection spectra, we notice an isolated peak (Sg3 or Sw2) in the case of SC-HD. These additional peaks are far from the mode discretized from the continuous mode in LRA. To understand the physics of the additional peak in the reflection spectrum of SC-HD case, we have compared the induced charge density distribution ρ 1 = −en 1 near the singular surface in Figure 4. In this figure, we give the induced charge distribution for all the marked peaks in Figure 3. Since the induced charge density for LRA is an infinite delta function, we therefore only show the distribution for SC-HD and HW-HD. For S-type peaks in both groove and wedge singular metasurfaces, the induced charge leaks out of the metal surface but with different patterns. Let us take the groove singular metasurface as an example, that is, Figure 4a, the Sg1 peak corresponds to the first-order surface plasmon mode, while the Sg2 is the second-order surface plasmon mode. These modes come from the discretization of the continuous mode Lg. We can tell that for Sg1 and Sg2, the induced charge near the interface is either primarily positive or primarily negative. These charge density profiles are similar to that of peak Hg1 of HW-HD. Although the induced charge is confined inside the metal surface for Hg1 mode, it is the same kind of surface plasmon mode as Sg1 and Sg2. However, the modes Sg1 and Sg2 differ from Hg1 in the centroid of induced charge density. From Figure 4a, we observe that for Sg1 and Sg2 the centroid of their induced charge density ρ 1 lies outside the jellium edge ∂Ω, while the absence of spill-out in Hg1 makes the centroid always inside the jellium edge. It is these spill-out and "spill-in" effects that lead to the different peak shiftings in the reflection spectrum. Likewise, in Figure 4b the spill-out of Sw1 mode and the spill-in of Hw1 and Hw2 modes makes the wedge singular metasurface exhibit a similar peak shifting behavior.
Special attention should be paid to the Sg3 mode in Figure 4a. From the charge profile, we can see that nearly the same amount of positive and negative charges on the two sides of the interface, which differs from that of Sg1 and Sg2 whose charge is dominated by either positive or negative charge. This additional mode is called Bennett mode, or multipole surface plasmon mode [29,33]. For traditional surface plasmon mode such as Sg1 and Sg2, the surface charge forms an effective monopole normal to the interface. However, for a multipole surface plasmon, an effective dipole forms normal to the interface. Therefore, the formation of a Bennett mode requires the spill-out of the electron at the metal interface, which explains why this multipole mode is absent in the HW-HD and LRA. Moreover, a similar Bennett mode is observed in Sw2 peak in wedge singular metasurface in Figure 4b. Similar to Sg3, induced electron for Sw2 mode also forms a dipole momentum normal to the interface. Finally, we illustrate the electric field enhancement of a singular metasurface for SC-HD, HW-HD and LRA in Figure 5. In our previous work, we show that the electric field enhancement of singular metasurface in the LRA case diverges as long as the vertex angle θ is larger than the critical angle θ c [9,34]. For a LRA mode (ω = 0.65ω p ) in Figure 5a, the corresponding critical angle θ c = Im[ln( ε−1 ε+1 )] = 0.237 rad < θ, making the electric field diverge at the terminal point. For a LRA mode (ω = 0.75ω p ) in Figure 5b, the critical angle θ c = −Im[ln( 1−ε ε+1 )] = 0.333 rad < θ, leading to a divergence of electric for a wedge singular metasurface. Figure 5. Electric field enhancement of singular metasurface (θ = 0.6π rad) with SC-HD, HW-HD and LRA, which corresponds to the peak marked in Figure 3. (a) Groove case; (b) Wedge case. Note that the color function has been fixed as the same for all of contour plot so as to compare field enhancement properties. Besides, the electric field of LRA mode diverges at the singularity.
Comparing LRA field with that of SC-HD and HW-HD in Figure 5, we conclude that the electric field in LRA is more localized at the sharp point. This delocalization is mainly attributed to the Pauli exclusion principle, that is, nonlocality, which hinders the localization of electrons to a very large density. This constraint leads the field enhancement of SC-HD and HW-HD to be finite and smaller than that of LRA case. Moreover, another comparison can be made between traditional surface plasmon mode and multipole surface plasmon mode. For Bennett mode in Figure 5 (Sg3 in the left panel or Sw2 in the right panel), we can observe that this multipole mode is highly localized at the interface, which means it shows a stronger decay normal to the interface compared with other traditional monopole surface plasmon modes. Note that in Figure 5 some modes exhibit maximum field enhancement at the cusp, while some are not. This is mainly because of different reflection phase at the cusp experienced by these surface plasmon modes [9,35].

Conclusions
In this paper, we consider the electron spill-out effect in the singular metasurface. The ground state calculation shows that our sharp singular point is effectively blunt. Besides, for the excited state, the introduction of the spill-out effect strongly changes the reflection spectrum, charge distribution and near-field profile. Therefore, it is vital to include this spill-out effect in the singular system, especially for a simple metal where the inter-band transition is absent. Finally, the electron spill-out calculation is still numerical in this paper, and an analytical approach of spill-out effect in the singular plasmonic system will be left to future work.
Author Contributions: F.Y. designed and performed research; F.Y. and K.D. analyzed data; F.Y. and K.D. wrote the paper. All authors have read and agreed to the published version of the manuscript.
Funding: This research was funded by Gordon and Betty Moore Foundation.

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