Influence of the electron spill-out and nonlocality on gap-plasmons in the limit of vanishing gaps

We study the effect of electron spill-out and of nonlocality on the propagation of light inside a gap between two semi-infinite metallic regions. We compare the predictions of a local response model taking into account only the spill-out, to the predictions of a quantum hydrodynamic model able to take both phenomena into account. We show that only the latter is able to correctly retrieve the correct limit when the gap closes, while the local model suffers from undesirable features (divergence of the fields, overestimation of the losses). Finally, we show that, to a certain extent, the correct results can be retrieved using a simple local approach without spill-out or conventional Thomas-Fermi approximation, but considering an effective gap width.


I. INTRODUCTION
The possibility to design extremely miniaturized devices, even optical circuits, has fostered research in plasmonics for decades. The most recent advances in fabrication processes have led to metallic nanostructures where extremely narrow gaps play a central role [1]. Metalinsulator-metal (MIM) structures support a mode with a particularly large effective index called the gap-plasmon [2,3]. The insertion of MIM sections in finite systems can thus lead to resonances in exceptionally small volumes [4,5], while the quality factor of the resonance [6] and its absorption cross-section [7] are essentially preserved. The insulator gaps can be so small (typically smaller than 1 nm [8]) that many questions arise concerning the accuracy of the electromagnetic description at such scales. Indeed, at these scales the classical local response approximation (LRA), based on the Drude model, reaches its limits [9], as the spatial dispersion due to the repulsion between electrons inside the metal has a direct impact on the optical response of the system. Electron-electron interactions can be taken into account in a simple way in the framework of the hydrodynamic theory (HT) with hard-wall boundary conditions (i.e., no electron spill-out is allowed), by adding the Thomas-Fermi (TF) electron pressure term [10,11].
When dealing with the near-atomic subnanometer length scales, it becomes inevitable to consider the quantum nature of electrons, and the quantum mechanical effects, such as electron spill-out and quantum tunneling, become important. The TFHT cannot address these phenomena properly and, therefore, more sophisticated * muhammad.khalid@iit.it † cristian.ciraci@iit.it theories are required for an accurate description of electron dynamics [12][13][14]. Ab-initio quantum mechanical approaches, such as time-dependent density functional theory (TD-DFT), provide an appealing solution to such problems although they are often computationally demanding, especially for large plasmonic structures [15][16][17].
Within the LRA framework it is still possible to account for the tunneling by interposing between the metal regions a fictitious material described by an effective conductivity [18] or by using a finite-current boundary condition at the metal-insulator interface [19][20][21]. However, when the LRA is used to account for spill-out, nonphysical features (including a divergence of the associated electric field) arise [22,23]. Due to this artificial divergence of the field, the imaginary part of the propagation constant of the gap-plasmon can be largely overestimated, which consequently results in very high modepropagation losses.
In the past few years, a hybrid approach called quantum hydrodynamic theory (QHT)-which adds a ∇ndependent correction (n being the charge density) to the TF kinetic energy-has been developed to study nonlocal and quantum effects (electron spill-out/tunneling) in plasmonic systems more accurately [24][25][26][27][28]. The QHT method intrinsically account for the electron spill-out and tunneling, but also spatial dispersion (nonlocality) which is relevant in the unidirectional plasmonic waveguides [29] and nonlinear optical dynamics [30,31]. This approach is based on orbital-free description of quantum electronic systems and its accuracy relies on the exactness of non-interacting kinetic energy functional. Ideally, if the exact expression of such functional were known, QHT would be as accurate as TD-DFT [32].
Here we apply the QHT approach to the gap-plasmon in the vanishing gap limit, and compare to the LRA model in the presence of spill-out. We show that the unphysical divergences predicted by the LRA lead to a large and artificial overestimation of the absorption when the gap-plasmon propagates, but that including spatial dispersion in the description of the electron gas makes such features disappear. In fact, the gap-mode propagation constant is influenced by the presence of longitudinal modes (due to nonlocality) inside the metallic region, which are determined by the hydrodynamic equations and, in particular, by the quantum pressure term. The correct limit, as the gap closes, is thus retrieved only when both the spill-out and nonlocality are taken into account.
In section II, we detail the theoretical background of our study, first by deriving simplified equations corresponding to the LRA approximation in the presence of spill-out and then by detailing the differences with the QHT framework. In section III, we discuss our numerical results by comparing the behavior of the effective index for the different approaches of interest (LRA, TFHT, LRA with spill-out and QHT). Finally, in section IV, we discuss the implications of our work regarding the definition of the actual optical edges of a metal. We show that it makes sense to define an effective gap for the gapplasmon, allowing to understand why in so many works LRA without spill-out has been so successfully applied.

II. THEORETICAL FRAMEWORK
The propagation of the electric and magnetic fields, E and H, of the gap-plasmon is described by Maxwell's equations: where µ 0 and ε 0 are the vacuum permeability and permittivity, respectively; ε ∞ is the relative permittivity accounting for any dielectric local response. J is the current inside the electron gas. This current is given by J = −enu where u is the velocity of the electron gas, considered as a fluid, e is the electron charge in absolute value, and n is the electron density. In this framework the velocity satisfies the following equations [33,34]: where m is the electron mass and the energy functional G[n] contains all the internal energy of the electronic system (here the operator δ δn indicates the functional derivative with respect to the density n). The last term in the second equation gives rise to a nonlocal contribution to the current J. Using the previous equations, we obtain the following evolution equation for the current: Neglecting higher order terms and taking n(r, t) = n 0 (r) + n 1 (r, t), and E(r, t) = E 0 (r) + E 1 (r, t), with n 0 being the equilibrium charge density and E 0 the electrostatic field, we first obtain an equation for the dynamic part, are the zero-th and first order terms of the potential, respectively. Second, we get an equation linking the static electron density to the static electric field, where E 0 and n 0 can be found using Gauss's law: n + being the positive background density. Finally, combining Eqs. (1) and assuming a harmonic field oscillating at ω we obtain: where k 0 = ω/c. Combining Eqs. (5) and (6), we obtain: Note also that using Eq. (5) we can simplify Eq. (8) as: A. LRA in the presence of spill-out We first derive the equations corresponding to the LRA [18,22,23]. They have the merit to allow a better grasp of the physics at play. The main assumption here is to neglect the pressure term, i.e. G[n] 0, making the description purely local.
The impact of nonlocality alone in the case of vanishing gaps (and hard-wall boundaries) has already been investigated in Ref. [35], where it has been shown that the effective index of the gap-plasmon increases although staying finite (contrarily to the diverging behavior of the LRA case without spill-out [36,37]). Yet, it does not tend to the expected value of the bulk metal refractive index either.
Starting from Eq. (10) we get The current can always be incorporated into Eqs. (1) as an effective polarization P f , by simply writing ∂ t P f = J. This finally allows to introduce a local susceptibility by writing, just as for a Drude model, that except that the susceptibility is proportional to the electron density and thus depends on the position here. Now if we consider the propagation of a guided mode in the z direction, and assume that all fields are invariant in the x direction, Maxwell's equations in the p-polarization reduce to: where ε(y) = ε ∞ + χ f (y) is the local permittivity.
Combining the previous equations to get rid of H x we obtain, Since we are looking for a guided mode we take ∂ z = ik z so that we obtain and finally manage to keep only the E z field to write the wave equation under the form We underline that such an approach is not selfconsistent, as it requires to choose arbitrarily the electron density. The choice of electron density put aside, our local model with spill-out is identical to the model used in [22] and yields perfectly similar results. In Ref. 22, the authors compute n 0 using a quantum mechanical approach whereas we use a simple model density to approximate the electron spill-out (details are given in the next section).
We have then converted this equation into a numerical problem using a finite difference scheme, in order to look for values of k z in the complex plan which allows to make the determinant of the finite difference matrix vanish. As such a problem is extremely unstable numerically, obviously because of the different scales involved (the gap being much smaller than the skin depth), we have used the exponentially decaying analytic solution inside the metal. This solution is known when ε(y) is uniform and corresponds to the classical solution of Maxwell's equations considering Drude's model for the metal. Connecting the finite difference problem with the analytic solution allows to reduce the size of the problem and to improve the stability, finally allowing to compute the wave vector of the gap-plasmon, even for vanishing gaps. We compared it with a finite element based approach, and both are in almost perfect agreement (see Appendix). B. Self-consistent quantum hydrodynamic theory Equation (9) with Eqs. (7) and (10) constitute the basis of the QHT [24,26] linear response and allow to calculate the detailed electron gas dynamic at the metal surface. In this work, we employ the following approximation: where T TF is the TF kinetic energy, T vW is the gradientdependent correction, namely the von Weizsäcker term, to the TF functional, which allows to take into account effects depending on the gradient of the electron density n and thus should not be neglected when spatially dependent densities are considered. E XC is the exchangecorrelation energy functional. The parameter λ vW in front of the von Weizsäcker term can be related to the decay of the charge density from a metal surface, and hence controls the amount of the electron spill-out and tunneling. It has been shown in Ref. 26 that for QHT λ vW is inversely proportional to the decay of the charge density and can be related to the "equilibrium" and "induced" electron spill-out. In general, λ vW is a free parameter of the kinetic functional that can be tuned depending on what the figure of merit is (i.e., plasmon energy, equilibrium density asymptotic decay, tunneling, etc). Most often, in the literature a value in the [1/9, 1]-range is used for λ vW ; in which λ vW = 1/9 corresponds to a lower electron spill-out and a faster decay of n 0 from a metal surface, whereas λ vW = 1 gives a higher spill-out and a slower decay of the charge density. In general, using λ vW = 1/9 approximates well the Lindhard function for small k-vectors while λ vW = 1 gives a good estimation for large k-vectors [? ]. In fact, the von Weizsäcker correction in Eq. (20) is the first-order term in the expansion of the kinetic energy [? ] and to construct a more generally valid functional one would need to consider higherorder terms (i.e., Laplacian dependence), which would introduce more free parameters [42]. The explicit expressions of the energy functional G[n] and its functional derivative, i.e. δG δn 1 , can be found in Ref. 26. We further note that Landau damping is another important factor to consider and it can also be incorporated in Eq. (10), as was done in Ref. 27 in the context of QHT. We expect that the Landau damping will result in increased damping inversely proportional to the gap size. Since the main objective of the present work is to analyze the impact of nonlocality and electron spill-out on the gap plasmons, we neglect Landau damping. It is informative to note that considering a spatially-constant n 0 = n b , n b being the charge density of the bulk metal, and G[n] = 0 in Eq. (20) leads to the usual Drude-type relation (LRA) for the polarization, and using T vW = E xc = 0 returns G[n] = T TF , the standard hydrodynamic theory in the Thomas-Fermi approximation (TFHT).
In contrast to the conventional theories in plasmonics, QHT can efficiently take into account the nonlocal, electron spill-out and tunneling effects, and can provide the details of microscopic distribution of the fields which usually cannot be described by conventional approaches. The QHT is highly dependent on the spatially varying equilibrium charge density n 0 and, consequently, approximation of the optical response of a nanoplasmonic system relies strongly on the description of n 0 . Over the past few years, QHT has been applied to probe linear plasmonic properties of a variety of individual nanoparticles as well as nearly-touching nanostructures (i.e., in the tunneling regime) [26,[38][39][40] and a fairly good agreement between the QHT and TD-DFT calculations has been reported [38]. Very recently, electron spill-out effects in singular metasurfaces have been analyzed using QHT, showing that the spill-out effectively blunts the sharp singularities [41]. QHT has also been applied to study nonlinear optical properties of thick metal films and it has been shown that there exist spill-out induced resonances which can be exploited to achieve very large second-harmonic generation efficiencies [31]. These resonances cannot be excited in the LRA and TFHT, as electron spill-out is overlooked in these approximations. Therefore, it becomes important to take into account the nonlocal and quantum effects when dealing with plasmonic structure with sub-nanometer geometrical features. Despite the fact that QHT gives a very good prediction of the plasmon resonances by adding a ∇n−dependent correction to the kinetic energy, however, it also suffers some limitations. For example, it predicts unexpected modes between the surface and bulk plasmon frequencies lying in the far tail region of the exponentially decaying charge density. Very recently, a development in the QHT theory based on the Laplacian-level kinetic energy functionals has been presented [42]. The introduction of Laplaciandependent charge density results in more robust numeri-cal solutions but its implementation becomes more complex. Since in the present study we only consider very small gaps, the tail of the charge density never extends to the problematic region. Therefore, we expect that the QHT approach in the limit of the von Weizsäcker approximation used herein gives us results with the same level of accuracy and with the advantage of avoiding additional numerical complexity, arising from the Laplacian-level functionals. We adopt the jellium approximation [43] in implementation of the QHT, which assumes that the electrons in a metal are confined by a constant positive background charge n + = ( 4 3 πr 3 s ) −1 , where r s indicates the Wigner-Seitz radius (r s = 4 for Na and r s = 3 for Ag).
In order to compute the gap-plasmon mode and its dispersion relation we first compute the space-dependent equilibrium charge density using Eq. (9) and then solve Eqs. (7) and (10) assuming a solution of the type E(r) = E(y) e ikgpz (with similar expressions for J), where k gp is the propagation constant of the gap plasmon mode. For a given angular frequency ω, we solve this system of equations using the finite-element method to compute the associated propagation constant k gp (ω), as well as the modeẼ(y).

III. RESULTS
Let us consider a MIM configuration as schematically shown in Fig. 1, that is, a dielectric gap g sandwiched between two semi-infinite metallic regions. For simplicity we use air as the dielectric. We examine the impact of nonlocality and the electron spill-out on the gap-plasmon, the unique mode guided in p-polarization in the structure and propagating along the z direction. More precisely, we compute the effective index n eff = kz k0 of the gap-plasmon according to the different models described above and compare the results, in particular in the gap size range (< 1 nm) where the charge densities from the two surfaces significantly overlap, as schematically shown in Fig. 1. We begin with the LRA with spill-out for simple Drude metals, i.e., sodium (Na), which we compare with the LRA and TFHT approximations without spillout as well as with the QHT. Finally, we consider the case of silver (Ag) in the QHT framework and compare to the commonly used LRA and TFHT without spill-out.

A. LRA with spill-out
We begin with the simplified case presented above in which the equilibrium electron density n 0 (r) can be described by an exponentially decaying profile with the following analytical function, where the parameter κ defines the spill-out or the decay of the charge density from the metal surface. In what follows, we set κ (0.7a 0 ) −1 , with a 0 being the Bohr radius, which we get by fitting of the analytical function with the decay of the self-consistent charge density (see Eq. (9)), considering λ vW = 1/9. We take r s = 4, ε ∞ = 1 andhγ = 0.16 eV, values that correspond to Na. Although it is impossible to work with Na in practical applications, due to its high reactivity, it presents an ideal Drude response as the impact of interband transitions can be omitted, thus offering a convenient theoretical study platform.
In Fig. 2 we compare the results of our simplified LRA with spill-out against the commonly used LRA without any spill-out. The input equilibrium charge density profiles used in the both approaches are plotted in Fig. 2(a) for g = 1 nm. In Fig. 2(b) and (c), we compare the behavior of the real and imaginary parts of the effective mode index. As expected, for large gaps (g > 1 nm) spillout can be neglected, the LRA with and without electron spill-out agree very well. The spill-out has a noticeable impact on the real part of the effective index for gaps that are typically around 1 nm or less: it increases the effective index, signaling a plasmonic slow down of the mode which is larger than without any form of spill-out. This is understandable, as the Poynting vector inside a plasma is opposite to the propagation direction, leading to a plasmonic drag [44], to a lower group velocity and finally to a larger wave vector.
When the gap is as small as 0.5 nm it begins to "close", electromagnetically speaking: the electron density inside the gap becomes so large that the propagation of the gapplasmon is hindered. The real part of the effective index decreases quickly, while its imaginary part is much higher when the spill-out is taken into account. It is important to note that in the limit of vanishing gaps, where the spill-out is expected to virtually close the gap, both real and imaginary parts of n eff should ideally converge to the index of the bulk metal. In Fig. 2(c), we notice that when g < 0.35 nm, i.e., the gap distances where the gap is expected to vanish completely, the imaginary part of n eff tends to bend back. However, we have found to be very hard to obtain a proper solution in this region. Indeed the dotted line as shown in Figs. 2(b) and (c) represents a region where we could not find a solution due to numerical artifacts. When g < 0.28 nm, we find n eff = n Bulk , where n Bulk being the refractive index of the bulk metal. These results are in agreement with the ones reported in [22]. We attribute the large imaginary part of the effective index to peaks in the electric field that can be seen when the electron density reaches a value so that χ f 0, meaning a vanishing local permittivity ε(y). The components of the electric field in the transverse and longitudinal directions with respect to mode propagation are shown in Figs. 2(d) and (e), respectively. The extremely narrow peaks the electric field exhibits have already been seen in a similar context [22] with a different electron gas density. This means that they are fundamentally linked to the assumption that a local permittivity is able to de-scribe the spill-out accurately. A clue that nonlocality may play an important role here is that ε = 0 is the condition for which volume plasmons are expected to be supported inside the electron gas.
The peaks in the electric field drastically increase the losses experienced by the gap plasmon, since the local absorption by the electron gas is proportional to E 1 · J * , where * denotes the complex conjugate [45]. This explains the larger imaginary part. We underline that, on the contrary, the magnetic field remains smooth (see Fig. 2(f)).

B. QHT for simple metals
In the framework of QHT, we compute the selfconsistent ground-state density of a symmetric MIM waveguide from Eq. (9). First we consider that the guiding structure is composed of an air gap sandwiched between two semi-infinite Na metals, to better grasp the difference with the simplified spill-out model.
The self-consistent ground-state density n 0 (r) for a Naair-Na configuration is shown in Fig. 3 (a) considering different gaps. The metal-dielectric interfaces for each gap size are indicated by the dashed lines. It can be seen that when the inter-distance of the metal surfaces is sufficiently large, there is no overlap of the charge densities. In this scenario, the system can be well described using a classical formulation. However, when the gap spacing between the metal boundaries decreases, the charge density profiles start to overlap and the electrons from the two surfaces start to tunnel across the gap. For sufficiently small gaps the charge transfer between the metals is large enough to close the gap optically speaking, despite the fact that the metal ion surfaces do not touch physically. Figure 3 (b) presents a comparison of n 0 between the LRA, TFHT and QHT for gap size g = 1 nm. In the LRA and TFHT approximations, we assume a constant n 0 = n b in the metal and zero outside. In the QHT case, the space-dependent n 0 (r) is plotted for both λ vW = 1/9 and 1. The parameter λ vW = 1/9 represents a lower spill-out and a faster decay of the charge density, whereas λ vW = 1 gives a higher spill-out and a slower decay, as shown in Fig. 3 (c).
In Fig. 4, we plot the effective refractive index as a function of gap distance g. The results computed within the QHT, considering the self-consistent space-dependent charge density, are compared against the LRA and TFHT without spill-out. It is interesting to note that for extremely small gap separation, LRA predicts very large values of both real and imaginary parts of n eff , whereas the TFHT gives relatively smaller values. Nonetheless, the trend is quite similar, i.e. n eff increases monotonically as the gap decreases. In contrast, in the QHT approximation, the n eff first increases with the decreasing gap size and then rapidly decrease after a certain gap distance, getting a value of the refractive index of the bulk metal. This is due to fact that at a certain gap size, a rea- sonable charge transfer through the gap occurs between the metal surfaces due to the overlapping of the charge densities and a virtual fade-out of the physical gap appears when the spacing between the metal boundaries is further decreased. For sufficiently small gaps, the mode virtually vanishes due to the electron tunneling and both the real and imaginary parts of the n eff converges to the real and imaginary parts of the refractive index of the bulk metal, as shown in the inset of Fig. 4(a) and (b), respectively. The expected limit for the effective index of the gap-plasmon is in that case actually retrieved. In the QHT case, results for both λ vW = 1/9 and λ vW = 1 are presented. Since λ vW = 1/9 entails the lower spill-out, the gap virtually vanishes when g < 0.3 nm, however, the gap disappears even at higher separation (g < 0.45 nm) for λ vW = 1, as it implies higher electron spill-out. Now, to show the influence of the electron spill-out on the mode profile, we plot here the electric |E y | and magnetic |H x | fields in the orthogonal direction to the metal-air interfaces under different approximations for g = 1 nm and E = 1 eV, as shown in Fig. 5. Both in the LRA and TFHT, |E y | is constant throughout the gap and then exponentially decays in the metal, as can be seen in Figs. 5 (a) and (c). QHT, on the other hand, predicts that |E y | has a maximum value at the center of the gap and it decreases as moving away from the center. The effect of the electron pressure in the metal in the TFHT can also be observed on the field profile. The |H x |, as plotted in Fig. 5 (b), remains quite similar in  the gap in all approximations expect near the interfaces where QHT predict a bit lower value. The decay of fields is given in Fig. 5 (c) and (d), which show that the fields decay faster in QHT and slower in TFHT with respect to the LRA fields, essentially because the wavevector along z and thus the decay along the y axis is slightly different. Now, we analyze the impact of electron spill-out and tunneling on the induced charge density. We, in fact, consider the spill-out and tunneling in both the equilibrium charge density and in the linear response, as the modification of the equilibrium density gives a more enhanced overlap of the induced charge densities. The induced charge density indeed is the equivalent of the transition density in the TD-DFT, and, in the limit of the kinetic energy functional approximation, represents the sum of all possible electron transitions, including those describing the tunneling. As an example, we plot in Fig. 6 the induced charge density n 1 along with the equilibrium charge density n 0 in the same scenario as considered in Fig. 5. The peak of the induced charge density from each interface lies outside the jellium edge and the tails of the densities go more into the dielectric medium and consequently showing more overlap.
Finally, it is interesting to compare the QHT and the simplified LRA with spill-out, as shown in Fig. 7. The real part of n eff in both cases accord well until the gap closes ( Fig. 7 (a)) whereas, as described previously, the imaginary part is much higher in the simplified spill-out model ( Fig. 7 (b)). This really means that nonlocality avoids the existence of the peaks in the electric field when the electron density reaches a particular value, making the imaginary part much closer to the prediction of the  LRA without spill-out.

C. QHT for noble metals
Noble metals, such as Ag, show in general a more complex behavior compared to Drude-type metals due the contribution of interband transitions and core electrons.
In what follows, we use r s = 3 a.u. andhγ = 0.03 eV for Ag, and account for interband effects with a constant background dielectric constant ε ∞ = 5.9 [31], which adds a local contribution to the free-electron metal permittivity. The effective refractive index for the Ag-air-Ag configuration as a function of the gap size computed using different models are plotted in Fig. 8. The real and imaginary parts of n eff are plotted in Fig. 8 (a) and (b), respectively, showing quite similar behavior as seen for the Na-Air-Na case in the previous section. LRA and TFHT clearly show unphysical increasing character of effective index when sub-nanometer gaps are in question. QHT, on the other hand, predicts that in the tunneling regime, for sufficiently small metal separation, the gap virtually fades away completely and the effective index converges to the refractive index of the bulk Ag. Note that in Fig. 8, the QHT gives a different threshold for g at which the gap virtually vanishes, depending on the amount of electron spill-out considered from the metal surface, i.e., λ vW = 1/9 or λ vW = 1. So far, we have analyzed the effective mode index of the MIM waveguide as function of air gap between the metals. It is also important to examine the behavior of n eff as a function of energy. Figure 9 (a) and (b) present the real and imaginary parts of n eff , respectively, as a function of energy (E in eV) for a Ag-air-Ag waveguide, considering a gap size g = 1 nm -a width which is experimentally rather easy to obtain [9].
This underlines the differences between the predictions of the different approaches regarding the effective index of the mode. The LRA predicts a bend back for the real part and a much larger imaginary part in this region. The TFHT on the contrary predicts a much lower imagi- nary part at any frequency. The QHT approach predicts a more reasonable behavior of the gap-plasmon with no bend back in the effective index and values of its imaginary part that can be much lower than for the LRA at higher frequencies.
The lower the imaginary part of the effective index, the larger the typical propagation length of the gap-plasmon, and finally the larger the quality factor of an eventual cavity for the gap-plasmon [46]. This underlines the importance of an accurate prediction in this regime.

IV. EFFECTIVE GAP
To summarize qualitatively the above results, the spillout seems to be able to close the gap electromagnetically speaking before it is actually spatially or mechanically closed. In addition, we have seen that when the spill-out is taken into account, whether in the simplified model or within the QHT framework, the effective index of the gap-plasmon is slightly larger than when the spill-out is neglected, however it presents a very similar trend. All of these results suggest that everything occurs as if the gap were narrower when the spill-out is taken into account, which seems physically sound. The fact that, for the simplified model, peaks seem to appear when the electron density reaches a certain value suggests that it should be possible to define the electromagnetic edges of the metal, as well as an effective value for the gap width. We define this value as the width that should be considered in the LRA model or in the TFHT approach for the gap-plasmon to fit the effective index of the gap-plasmon provided by the QHT. The effective gap is thus given by where g is the actual gap between the edges of the metal. Said otherwise, in order to retrieve the results of the QHT, one can use the LRA (or the TFHT) with hard wall boundaries, but simply use a diminished size of the gap to take the spill-out into account. For silver, shifting the LRA curve of δ g = 0.15 nm or the TFHT curve of δ g = 0.22 nm allows to retrieve the effective index of the gap-plasmon as computed with the QHT extremely accurately, as long as the gap does not start closing (which occurs around g 0.6 nm). The results for the QHT reported in Fig. 10 are for λ vW = 1/9, however, the QHT result for λ vW = 1 can also be nicely reproduced from the LRA and TF approximations considering δ g 0.25 nm and δ g 0.32 nm, respectively. The idea that the gap-plasmon is sensitive to the effective gap between the metals is thus well founded, as it yields quantitative results. This can be understood because the spill out makes the response of the electron gas outside the metal similar to the response inside the metal, as soon as the electron density is high enough. However, we point out that we have not found another way to define the effective gap based on the profile of the electron density or of the fields. The values we give are simply the ones that yield the best fit and there does not seem to be any simple way to predict such values.
This discussion finally underlines that at such scales, for the gap-plasmon, the tiny difference (a rule of thumb would be, according to our results, to consider a 0.1 nm 3 . 0  Figure 10. The effective mode index for Ag-air-Ag waveguide (a) real and (b) imaginary part as a function of gap distance g between the metal edges plotted at E = 1 eV. The LRA and TFHT curves are shifted by δg to cope with the QHT results, which allows to introduce the concept of an effective gap that can be defined as ge = g − δg.
The ge seems to appear smaller than the mechanical gap g by an amount of δg. In particular, using δg = 0.15 nm in the LRA and δg = 0.22 nm in the TF approximation gives good fit with QHT. The insets plot n eff as a function of actual gap.
difference) that exists between the actual edge of the metal and the "electromagnetic edge" has its importance and must be considered in the limit of vanishing gaps.

V. CONCLUSIONS
In this article, we analyzed the impact of both electron spill-out and nonlocality on the propagation of gapplasmons in MIM waveguides using the full self-consistent QHT, which can efficiently take into account both spatial dispersion inside the electron gas and the spatial variation of the electron density outside of the metal. We also compare our results with the commonly used approaches (LRA and TFHT) without any spill-out and with a popular LRA approach for the spill out [18,22,23].
Peaks in the electric field profile and a large imaginary part for the effective index of the gap plasmon are predicted by the latter approach. Such features are however absent in the framework of the QHT, showing that nonlocality is key to an accurate description of the gap plasmon. We underline that only in the QHT framework does the effective mode index converge to the refractive index of the bulk metal, as expected from physical considerations. Hence QHT, without being affected by the artifacts present in the local description, provides a good measure of the effective mode index and of the local field amplitudes, which are particularly relevant for nonlinear optical processes where response of the system is strongly affected by the local intensities. These points reinforce the idea that QHT is the right tool for plasmonics at the subnanometer level while using a LRA model for the spill-out may lead, at least in the case of gap-plasmon resonators, to an overestimation of the losses and an underestimation of the quality factor of the resonance entirely due to the fact that nonlocality has been neglected.
Our study thus shows that the spill-out indeed imposes an upper limit on the effective index of the gap-plasmon, confirming previous studies [22]. A gap-plasmon is very unlikely to present an effective index larger than 15 for noble metals in the visible, which puts a theoretical limit to the miniaturization of gap-plasmon resonators in the visible. Importantly, whatever model is chosen to describe the spill-out, it seems possible to introduce an effective gap width between metals, smaller than the gap size defined as the distance between the atoms of the metal across the gap. This effective gap width allows to take simply the effect of the spill-out into account, as long as the gap can be considered open to the propagation of the mode. In many experiments, the definition of the gap between metals is of crucial importance [7,9]. Our reasoning suggests using optical characterization methods (like ellipsometry typically) to determine the thickness of a ma-terial deposited at the surface of a metal is particularly relevant in the context of gap-plasmon resonators, as the optical thickness which is retrieved probably takes the spill-out partially into account. The gap width deduced from such measurements is in that way probably closer to the effective gap width we have defined here. This may thus explain why hard-wall approaches have been so far surprisingly successful at predicting the experimental behavior of gap-plasmons in the limit of vanishing gaps.

ACKNOWLEDGMENTS
Antoine Moreau is an Academy CAP 20-25 chair holder. He acknowledges the support received from the Agence Nationale de la Recherche of the French government through the program "Investissements d'Avenir" (16-IDEX-0001 CAP 20-25). a second order central scheme where E i = E z (y i ) represents the field at the discrete position y i = (i − 1)∆ and v ± i .
= v(y i±1 ) + v(y i ). At the boundary, where we use that for y > L the field has the exponential form E z = E z (L)e −ikz(y−L) . We obtain a homogeneous linear system for the discrete field E(y i ) parameterized by the wave number k z . The plasmon resonances are associated to the k z for which the system has a non trivial kernel. We use a iterative procedure based on the Nelder-Mead simplex algorithm [47] where we vary k z in the complex plane in order to minimize the singular values of the matrix associated to the linear system. The results are depicted in Fig