Influence of spatial dispersion on surface plasmons, nanoparticles and grating couplers

TL;DR

This paper demonstrates that spatial dispersion significantly influences surface plasmons, nanoparticles, and grating couplers, especially in high permittivity dielectrics, highlighting the need to consider nonlocal effects in plasmonic modeling.

Contribution

It shows that nonlocality affects surface plasmons and nanoparticle responses, emphasizing the importance of using advanced models like the hydrodynamic model in high permittivity environments.

Findings

Spatial dispersion impacts surface plasmon propagation.

Nanoparticles larger than 20 nm are sensitive to spatial dispersion.

High permittivity dielectrics require nonlocality considerations in plasmonic structures.

Abstract

Recent experiments have shown that spatial dispersion may have a conspicuous impact on the response of plasmonic structures. This suggests that in some cases the Drude model should be replaced by more advanced descriptions that take spatial dispersion into account, like the hydrodynamic model. Here we show that nonlocality in the metallic response affects surface plasmons propagating at the interface between a metal and a dielectric with high permittivity. As a direct consequence, any nanoparticle with a radius larger than 20 nm can be expected to be sensitive to spatial dispersion whatever its size. The same behavior is expected for a simple metallic grating allowing the excitation of surface plasmons, just as in Wood's famous experiments. Importantly, our work suggests that for any plasmonic structure in a high permittivity dielectric, nonlocality should be taken into account.

Tables

Table 1. Table 1: Uncertain parameters extracted from telemetry

Parameter	$\overline{z}$	${\delta }_z$	Units
$h_{\text{G}}$	$68.1$	$0.1$	(nm)
$d_{\text{G}}$	$499.2$	$1.0$	(nm)
$a_{\text{G}}$	$165.4$	$1.5$	(nm)
$h_{\text{D}}$	$84.8$	$0.3$	(nm)

Equations

$$n_{\mathrm{eff}}=\sqrt{\frac{{\epsilon }_{\mathrm{d}}{\epsilon }_{\mathrm{m}}}{{\epsilon }_{\mathrm{d}}+{\epsilon }_{\mathrm{m}}}},$$

$$-{\beta }^2\nabla (\nabla \cdot \mathbf{J})+\ddot{\mathbf{J}}+\gamma \dot{\mathbf{J}}={\epsilon }_0{\omega }_{\mathrm{P}}^2\dot{\mathbf{E}},$$

$$\frac{k_{\mathrm{m}}}{{\epsilon }_{\mathrm{m}}}+\frac{k_{\mathrm{d}}}{{\epsilon }_{\mathrm{d}}}-i\Omega =0.$$

$$\Omega =\frac{k_{\mathrm{SP}}^2}{{\kappa }_{\mathrm{l}}}\left(\frac{1}{{\epsilon }_{\mathrm{m}}}-\frac{1}{1+{\chi }_{\mathrm{b}}}\right),$$

$$2\pi R=m\frac{{\lambda }_0}{n_{\mathrm{eff}}}$$

$${\mathbf{P}}_{\mathrm{f}}.\mathbf{n}=\left(\mathbf{E}-(1+{\chi }_b)\frac{{\beta }^2}{{\omega }_p^2}\nabla (\nabla .\mathbf{E})\right).\mathbf{n}=0$$

$$a_n=\frac{\sqrt{{ϵ}_d}{J}_n(k_{d}R)(J_n^{\prime }(k_{m}R)+{\alpha }_n)-\sqrt{{ϵ}_m}{J}_n(k_{m}R)J_n^{\prime }(k_{d}R)}{\sqrt{{ϵ}_m}{H}_n^{\prime }(k_{d}R)J_n(k_{m}R)-\sqrt{{ϵ}_d}{H}_n(k_{d}R)(J_n^{\prime }(k_{m}R)+{\alpha }_n)}$$

$${\alpha }_n=\frac{\frac{n^2}{k_m{R}^2}{J}_n(k_{L}R)J_n(k_{m}R)}{-k_L{J}^{\prime }n(k_{L}R)+(1+{\chi }_b)\frac{{\beta }^2}{{\omega }_p^2}\left[k_L^2{J}_n^{\prime \prime \prime }(k_{L}R)+\frac{k_L^2}{R}{J}_n^{\prime \prime }(k_{L}R)-\frac{k_L}{R^2}{J}_n^{\prime }(k_{L}R)+\frac{2n^2}{R^3}{J}_n(k_{L}R)\right]}$$

$$k_L=\sqrt{-\frac{{\omega }_p^2}{{\beta }^2}\left[\frac{1}{1+{\chi }_b}+\frac{1}{{\chi }_f}\right]}.$$

$$C_e=-\frac{2}{k_{d}R}\sum _{n=-\infty }^{\infty }\Re (a_n).$$

$$k_{\mathrm{SP}}=m\frac{\lambda }{d_{\mathrm{G}}},$$

$${\lambda }_{\mathrm{c},m}=\frac{\Re (n_{\mathrm{eff}})}{m}{d}_{\mathrm{G}},$$

$${\left|\left|R_{\mathrm{meas}}(\lambda )-R(\lambda )\right|\right|}_{L^2({\lambda }_1,{\lambda }_2)}^2={\int }_{{\lambda }_1}^{{\lambda }_2}{\left(R_{\mathrm{meas}}(\lambda )-R(\lambda )\right)}^2\text{mbox}d\lambda .$$

$$f\left(\mathbf{z}\right)\approx \sum _{m=1}^{M}f\left({\mathbf{z}}^{(m)}\right){\Psi }_m\left(\mathbf{z}\right),$$

$$\mathbb{E}\left[R\right]$$

$$\mathbb{V}\left[R\right]$$

$$\mathbb{E}\left[R\right]$$

$$\mathbb{V}\left[R\right]$$

Full text

††preprint: APS/123-QED

For more than a century now, Drude’s modelDrude (1900), coupled to Maxwell’s equationsMaxwell (1881), has been able to describe very accurately the optical response of metals, even for quite extreme geometriesMoreau et al. (2012); Chen et al. (2014); Akselrod et al. (2014); Powell et al. (2016); Ayata et al. (2017). Many advanced theories describing metal-vacuum interfacesFuchs and Kliewer (1971); Feibelman (1975) have been developed during the second half of the twentieth century, especially with the development of Electron Energy Loss Spectroscopy (EELS) which provided experimental data to better ground the theoretical discussionsRocca (1995); Rocca et al. (1995); Park and Palmer (2009). These new approaches were able to take into account complex phenomena, like spatial dispersion or electron spill-out, and allowed to better understand the success of Drude’s model. It turns out in fact that the Drude’s theory is the zero-th order approximation of all more advanced descriptions introduced later on. All these studies however, seemed to conclude that spatial dispersion and spill-out have a limited impactLiebsch (1993a, b); Feibelman (1994); Liebsch (1994) on surface plasmons (SPs) such that optical experiments were not likely to show any difference from Drude’s predictionsBoardman (1982); Frostmann and Gerhardts (1986); Chang et al. (2003). Moreover, because most resonances in metallic structures can actually be explained as cavity resonances of some sort for SPs, it has become widely accepted that plasmonic resonances could be very accurately described by Drude’s model. For decades, then, there has not been any urge to adopt advanced descriptions of the response of metals in plasmonics. Only in the case of metallic clusters, due to the extremely small size, spatial dispersion and spill-out were expected to play a significant role, requiring the most advanced descriptionsLiebsch (1993a); Scholl et al. (2012); Raza et al. (2015a); Esteban et al. (2012).

However, a recent experiment with film-coupled nanoparticles showed that the frequency of the resonance of modes that are localized in small volumes (of the order of 1 nm3) is simply not correctly predicted by Drude’s model, whereas the linearized hydrodynamic modelCiracì et al. (2013a), in its simplest formulation (the Thomas-Fermi approximation) seemed to be accurate enough. This is the case even with an optical excitation and relatively large metallic particles. This can be linked to the fact that the small gaps between the nanoparticle and the metal support a gap-plasmon – a guided mode which is particularly sensitive to spatial dispersion because it has a very large wavevectorMoreau et al. (2013); Raza et al. (2013a); David et al. (2013); Wiener et al. (2013); Toscano et al. (2013) whatever the frequency. This allows to better understand why small gaps, which are more and more common in plasmonicsAkselrod et al. (2014); Moreau et al. (2012); Lassiter et al. (2014); Haffner et al. (2015); Nielsen et al. (2017), may require more advanced descriptions of the metallic response. Furthermore, this explains why the hydrodynamic model, despite its well documented deficienciesHaberland (2013), is probably a good replacement for Drude’s model in plasmonicsToscano et al. (2015); Ciracì and Della Sala (2016): high wavevector plasmonic modes enhance spatial dispersion effectsWiener et al. (2013); Toscano et al. (2013); David et al. (2013), but not the impact of the spill-out, especially since noble metals present a high extraction work. The spatial dispersion which arises from the repulsion between free electrons inside the metal is taken into account accurately by a hydrodynamic model, which presents in addition the advantage of being relatively easy to implement numericallyToscano et al. (2012); Benedicto et al. (2015); Schmitt et al. (2016, 2018).

Here we show that the impact of spatial dispersion on the SP propagating at the interface between a metal and a dielectric is enhanced when the refractive index of the dielectric is large. We show, as a direct consequence, that even large nanoparticles can be expected to be sensitive to nonlocality and that, using a grating coupler, it is theoretically possible to estimate the main parameter of the linearized hydrodynamic model, in a situation where other more complex phenomena can be ruled out. We underline that such a well controlled setup differs strongly from previous experiments which all involved chemically synthesized nanoparticlesCiracì et al. (2012); Scholl et al. (2012); Ciracì et al. (2014); Raza et al. (2013b) whose geometry cannot always be fully controlled.

I Surface Plasmon

In this first part we study the influence of spatial dispersion on a SP propagating along a metal-dielectric interface (especially for dielectrics with high refractive indices). A SP can be seen as current loops propagating beneath the surface of a metal. Such a phenomenon is accompanied by an electromagnetic field in the metal (with relative permittivity $\varepsilon_{\mathrm{m}}$) and in the dielectric (with relative permittivity $\varepsilon_{\mathrm{d}}$) which is transversely evanescent in both media due to the fact that its effective index $n_{\mathrm{eff}}=\frac{k_{\mathrm{SP}}}{k_{0}}$ is always larger than the refractive index of the dielectric medium.

First, neglecting losses, the dispersion relation can be written as

[TABLE]

where we assumed a simple Drude model $\varepsilon_{\mathrm{m}}=1-\frac{\omega_{\mathrm{P}}^{2}}{\omega^{2}}$. In this case, the curve of a SP has a horizontal asymptote at $\omega_{\mathrm{SP}}=\frac{\omega_{\mathrm{P}}}{\sqrt{1+\varepsilon_{\mathrm{d}}}}$. It is however unrealistic to neglect the losses inside the metal because the frequency $\omega_{\mathrm{SP}}$ is usually in a wavelength range, the UV, where the interband transitions make the metal highly lossy. As a consequence, a bend-back can be seen on the dispersion curve of the SP, which thus never reaches very high wavevectors (see local dispersion curves in Fig. 1).

However, when the permittivity of the dielectric increases, the frequency $\omega_{\text{SP}}$ decreases. At the same time, losses due to the interband transitions can be expected to be low enough as to enable a support of high wavevector SPs. Such modes are more likely to be sensitive to spatial dispersionMoreau et al. (2013).

In order to take nonlocality into account we rely on the linearized hydrodynamic model for the free electrons, already introduced in previous works Raza et al. (2011); Chapuis et al. (2008); Toscano et al. (2012); Ruppin (2005); Ruppin and Kempa (2005); Wiener et al. (2012); Ciracì et al. (2013b). The electric current $\mathbf{J}$ inside the metal is linked to the electric field $\mathbf{E}$ by

[TABLE]

where $\omega_{\mathrm{P}}$ is the plasma frequency, $\varepsilon_{0}$ the vacuum permittivity, $\gamma$ the damping factor and $\beta$ the nonlocal parameter. The $\beta$ factor represents the increase of the internal pressure in the electron gas due to exchange interaction and Coulomb repulsion. There are actually several theoretical expressions for this parameterCrouseilles et al. (2008). We rely on the experimental data availableCiracì et al. (2012, 2014), which consistently point to a value of $\beta=1.35\times 10^{6}\ \nicefrac{{\mathrm{m}}}{{\mathrm{s}}}$. Finally, the electric current inside the metal can always be considered as an effective polarization $\mathbf{P}_{\mathrm{f}}$ due to the free electrons and is then given by $\dot{\mathbf{P}}_{\mathrm{f}}=\mathbf{J}$. In that framework, the metal can then be described as a nonlocally polarizable medium.

We use accurate material parametersRakic et al. (1998) that allow a distinction between the response of the free electrons, that is subject to spatial dispersion, and the response of bound electrons, that can be considered to be purely localMaradudin and Mills (1973). In the following the metal is always assumed to be silver, which is favorable since silver is less lossy than gold. The total metal polarization reads then $\mathbf{P}=\mathbf{P}_{\mathrm{f}}+\mathbf{P}_{\mathrm{b}}$, where $\mathbf{P}_{\mathrm{b}}=\varepsilon_{0}\chi_{\mathrm{b}}\mathbf{E}$ with $\chi_{\mathrm{b}}$ being the local susceptibility associated to the bound electrons. In order do obtain $\chi_{\mathrm{b}}$, we fit the experimental data of silver permittivity with a generalized dispersion model based on Padé seriesLanteri et al. (2017) (which is sufficient to perform nonlocal time-domain simulations) and subtract the local Drude contribution $\chi_{\mathrm{f}}$. Throughout this work, we rely on the DIOGENeS 111DIOGENeS, https://diogenes.inria.fr/ Discontinuous Galerkin Time Domain (DGTD) suiteSchmitt et al. (2016). Finally we consider additional boundary conditions that are both the most natural (vanishing normal component of the polarization current $\mathbf{J}$ at the metal boundary, thus forbidding free electrons to escape the metal) and the most conservative (they reduce the impact of nonlocal effects)Moreau et al. (2013).

We first study the influence of nonlocality by considering the dispersion relation of a SP propagating along a perfectly plane dielectric-metal interface. Assuming time dependence of the form $e^{-i\omega t}$ , solving relation (2) for $\mathbf{J}$, injecting the result in Maxwell’s equations and considering both classical and additional boundary conditions, the nonlocal dispersion reads

[TABLE]

Here, $k_{j}=\sqrt{\varepsilon_{j}k_{0}^{2}-k_{\mathrm{SP}}^{2}}$, $j=\mathrm{m,d}$, are the vertical components of the SP wavevector and $\Omega$ is the parameter including $\beta$ and thus taking nonlocality into account. Presuming $\Omega=0$ allows to retrieve the usual dispersion relation for SPs, the $k_{j}=k_{j}^{\prime}+ik_{j}^{\prime\prime}$ being here essentially imaginary ($k_{j}^{\prime\prime}>>k_{j}^{\prime}$). The expression of this parameter isMoreau et al. (2013)

[TABLE]

where \kappa_{\mathrm{l}}^{2}=k_{\mathrm{SP}}^{2}+\Big{(}\frac{\omega_{\mathrm{P}}^{2}}{\beta^{2}}\Big{)}\Big{(}\frac{1}{\chi_{\mathrm{f}}}+\frac{1}{1+\chi_{\mathrm{b}}}\Big{)} is the vertical component of the wavevector associated with the longitudinal part of the SP appearing only when considering non-localityMoreau et al. (2013).

As can be seen from (4), the parameter $\Omega$ is roughly proportional to the square of $k_{\mathrm{SP}}$, which clearly indicates that the higher the wavevector, the higher the impact of nonlocality. Fig. 1 shows local and nonlocal SP dispersion curves for different values of the dielectric permittivity $\varepsilon_{\mathrm{d}}$ ranging from 1 to 10. We can easily see that, for a fixed frequency, an increasing $\varepsilon_{\mathrm{d}}$ pushes the SP towards higher $n_{\mathrm{eff}}$, and thus causes the SP to be more sensitive to nonlocality. The leftmost curve, obtained using air as dielectric ($\varepsilon_{\mathrm{d}}=1$) clearly shows that the effective index $n_{\mathrm{eff}}$ is too small to show any impact of nonlocality. However, after a certain $\varepsilon_{\mathrm{d}}$ value, the characteristic bend-back occurring at the SP frequency $\omega_{\mathrm{SP}}$ just disappears when spatial dispersion is taken into account.

This corresponds exactly to a recent theoretical studyRaza et al. (2013a) which shows that an artificial decrease of the metallic losses can induce exactly the same behavior on the dispersion curves of plasmonic guided modes. This occurs when the impact of nonlocality overcomes the influence of the metallic losses.

The impact of the dielectric’s permittivity overcomes the one of losses for two reasons. First, as explained above, increasing $\varepsilon_{\mathrm{d}}$ lowers $\omega_{\mathrm{SP}}$ and thus takes the frequency away from the interband transitions. Additionally, a higher dielectric permittivity directly gives to the SP a higher $n_{\mathrm{eff}}$ and thus enhances the influence of spatial dispersion, leading to a large impact of nonlocality even well below $\omega_{\mathrm{SP}}$.

We underline that, despite extensive studies on nanoparticles, such a behavior of the SP has seemingly not been reported previously. Straightforwardly, this suggests that nonlocality will have an impact on (i) the resonances of relatively large nanoparticles of noble metals (with a diameter well above 20 nm, as they can be considered as resonant cavities for the surface plasmon) and (ii) SP grating couplers very similar to the canonical experiment of WoodWood (1902) provided the grating is buried in high index dielectrics.

II Nanoparticles

The resonance of large metallic nanoparticles can be interpreted as cavity resonances for the surface mode with a resonance conditionKreibig and Vollmer (2013) which can be written as

[TABLE]

or simply as $k_{\mathrm{SP}}=\frac{m}{R}$, where $R$ is the radius of the particle. Such a condition is strictly valid only (i) for a cylinder instead of a sphere and (ii) if the curvature of the particle can be neglected, which is almost never the case. However, this condition being roughly valid even for spherical nanoparticles instead of cylindersKreibig and Vollmer (2013); Klar et al. (1998), it allows to understand that if the wavevector of the SP is influenced by nonlocality, the resonance frequencies of a nanoparticle should be influenced as well, irrespective of its size.

For decades now, the community has actually focused on nanoparticles with a diameter much smaller than 20 nmScholl et al. (2012); Raza et al. (2013b), hoping that enhanced nonlocal effects would take place - since spatial dispersion is linked to supplementary pressure terms in the description of the electron gasCrouseilles et al. (2008); Scalora et al. (2010). However, with very small nanoparticles (i) other effects like the spill-outScholl et al. (2013); Esteban et al. (2012); Teperik et al. (2013); Toscano et al. (2015) kick in and (ii) the geometry of the nanoparticles is not well controlledScholl et al. (2012); Raza et al. (2013b). Given the relatively poor agreement between experiments and the prediction of the hydrodynamic modelRaza et al. (2013b, 2015b), it seems difficult to consider the hydrodynamic model to be sufficient at such small scales. It may even be a little bit early to introduce further improvements of the hydrodynamic modelMortensen et al. (2014) based on these results. The simple analysis above suggests that larger nanoparticles buried in a high permittivity medium could actually be a better setup to test the hydrodynamic model.

In order to further strengthen this analysis, we used Mie theoryKlar et al. (1998); Kreibig and Vollmer (2013) and adapted the formalism proposed by RuppinRuppin (2001) for metallic cylinders, for which formula 5 is the most relevant. We use the most realistic material parameters possibleRakic et al. (1998) and a supplementary boundary condition which can be writtenMoreau et al. (2013)

[TABLE]

where $\mathbf{n}$ is the unitary vector normal to the interface. When the field is decomposed on the cylindrical harmonics, the Mie coefficientsRuppin (2001) are

[TABLE]

where

[TABLE]

with $k_{m}=\sqrt{\epsilon_{m}}\,k_{0}$, $k_{d}=\sqrt{\epsilon_{d}}\,k_{0}$ and

[TABLE]

We then compute the absorption cross-section as

[TABLE]

We underline that when $\alpha_{n}=0$, the expression (7) reduces to the regular Mie expression.

We have computed the local and nonlocal response of nanoparticles with different sizes when they are in water (as is common) or in TiO2 (see Fig. 2). The latter is a good choice to enhance the influence of spatial dispersion because of its high refractive index. Its permittivity $\varepsilon_{\mathrm{TiO}_{2}}$ is described by a generalized dispersion model fitted to experimental data corresponding to thin films of $\mathrm{TiO}_{2}$ grown by atomic layer depositionSiefke et al. (2016). This results in a real part of the refractive index of $\mathrm{TiO}_{2}$ comprised between 2.25 and 2.5 over a wavelength spectrum ranging from 2000 nm to 400 nm. The extinction coefficient is of the order of $10^{-7}$.

In water, nonlocality has a noticeable impact only for a radius approaching 10 nm, whereas in TiO2 nonlocality tends to blueshift all the resonances of more than 4.8 nm in wavelength even for a radius of 100 nm. We stress here that a high index dielectric is able to sufficiently enhance the magnitude of the nonlocal effects to make it observable on the response of particles/cylinders 5 times larger than the ones usually considered by the community. This should be enough to rule out any other effect like the spill-out and with such a large size, the geometry of the nanoparticles are better controlled - or the nanoparticles could even be probed individuallyKlar et al. (1998); Raza et al. (2015b).

III Grating coupler

We now discuss the structure shown on Fig. 3(a), which is a simple 1D metallic grating buried in a high index dielectric with infinite extent in the lateral directions. Normal incident illumination is assumed from the top and we recover the zero-th reflected order for a broadband spectrum of the incident wavelength. In such a simple configuration many diffraction orders, including evanescent ones, are excited.

For the $m$-th order, and since we only consider normal incidence here, the coupling condition to the SP can be simply written as

[TABLE]

where $\lambda$ is the wavelength, $m$ the diffraction order and $d_{\mathrm{G}}$ the grating pitch. Such a relation is valid only for a very shallow grating when the surface mode can be considered undisturbed. It is equivalent to assuming that the spatial periodicity of the grating is a multiple of the periodicity of the SP. Such a condition is proper to the periodicity of the structureMaystre (2012) and thus is expected to be valid whether spatial dispersion (which will modify $k_{\mathrm{SP}}$) is taken into account or not. This allows to find the wavelength $\lambda_{\mathrm{c},m}$ for which the grating is able to excite the SP

[TABLE]

where $n_{\mathrm{eff}}=\frac{k_{\mathrm{SP}}}{k_{0}}$ is the SP effective index.

Each time this condition is satisfied, a dip due to the excitation of a SP will appear in the reflectivity. Spatial dispersion should cause blueshifts of the resonances with respect to a fully local approach because the effective index is always smaller when nonlocality is taken into account. We define this blueshift as the positive quantity $\Delta\lambda_{\mathrm{c},m}=\lambda_{\mathrm{c},m}^{\mathrm{local}}-\lambda_{\mathrm{c},m}^{\mathrm{nonlocal}}\propto n_{\mathrm{eff}}^{\mathrm{local}}-n_{\mathrm{eff}}^{\mathrm{nonlocal}}$. As the different resonances correspond to different orders being coupled to the SP, the corresponding wavevectors will be different and thus the impact of spatial dispersion will change from one resonance to another. To better understand the impact of nonlocality on the resonances, it is crucial to be able to identify them, hence the interest in the relation (12). We have thus taken the finite thickness of the dielectric layer into account by computing the properties of the guided mode of the non-corrugated structure (see Fig. 3(b)) using an open-access numerical toolBenedicto et al. (2015); Defrance et al. (2016). Aiming for a proper interpretation of the resonances, we have then carefully increased the depth of the grating, in order to allow a thorough physical discussion (see Fig. 5). We use Rigorous Coupled Wave Analysis (RCWA)Granet and Guizal (1996); Lalanne and Morris (1996) for these simulations.

III.1 Parameters of the grating coupler

III.1.1 Materials

Although we restrict ourselves to silver in the scope of this work, we underline that gold or any other metal and even semi-conductorsMaack et al. (2017) (claiming an equivalent electronic mean free pathGall (2016)) could have been used and would have led to the same conclusions. Again, in order to observe the largest possible effects, we consider TiO2 as a high permittivity dielectric.

III.1.2 Thickness of the dielectric and metallic layer

We have taken $h_{\mathrm{m}}=150\ \mathrm{nm}$ for the thickness of the metallic layer, thus being several times thicker than the skin depth and ensuring that the substrate does not play any significant role.

The thickness of the dielectric $h_{\mathrm{D}}$ layer is a distinctly more crucial parameter. In order to excite a SP with the highest possible wavevector, a dielectric of the highest possible thickness would be desirable. However, since the dielectric layer is finitely thick, increasing $h_{\mathrm{D}}$ leads to a higher number of classical guided modes. That hinders a clear interpretation of the resonances or renders it even impossible. Luckily, since we are considering the excitation of high wavevector guided modes that show a fairly low vertical extension, a thickness of $h_{\mathrm{D}}=85.0$ nm turns out to be a good trade-off. No other guided mode than the SP exist for this choice of $h_{\mathrm{D}}$ in most of the spectrum.

Fig. 4(a) shows the dispersion curve for a mode propagating at an air-$\mathrm{TiO}_{2}$-$\mathrm{Ag}$ interface for $h_{\mathrm{D}}=85$ nm with and without spatial dispersion. Obviously, the difference with the dispersion curve of the SP at a $\mathrm{TiO}_{2}$-$\mathrm{Ag}$ interface is small. The impact of nonlocality is clearly the same and the bend-back disappears.

III.1.3 Grating pitch

Using the coupling condition (12) and the local and nonlocal dispersion curves as shown in Fig. 4(a), a raw estimation of the blueshift due to spatial dispersion can be made. Fig. 4(b) shows the expected blueshift for different diffraction orders and sweeps over the grating periodicity $d_{\mathrm{G}}$. The higher the diffraction order, the higher the wavevector of the excited SP – which leads to an increased impact of nonlocality. We have finally chosen a pitch length $d_{\mathrm{G}}=500\ \mathrm{nm}$ according to a maximum of the predicted blueshift of about 5.5 nm for $m=5$. Keeping the pitch below the shortest working wavelength guarantees that only the zero-th order propagates, even if the evanescent orders of diffraction are coupled to the SP. In other words, all of the non-reflected light must be absorbed by the guided modes along the structure.

III.1.4 Groove width and grating thickness

The determination of the two remaining geometrical parameters, i.e. the groove $a_{\mathrm{G}}$ and the grating thickness $h_{\mathrm{G}}$, is less straightforward. We have to establish a trade-off, such that the excitation of the SP for different orders of diffraction can be done efficiently without perturbing the guided mode too much. In order to avoid gap plasmons to build up in the slitsDechaux et al. (2016), a large enough $a_{\mathrm{G}}$ value is needed. We choose $a_{\mathrm{G}}=\frac{1}{3}d_{\mathrm{G}}$ here. Starting from $h_{\mathrm{G}}=2$ nm we have increased the grating depth until an efficient coupling to the SP mode was found - we relied on a pure RCWA methodGranet and Guizal (1996); Lalanne and Morris (1996) to adjust this parameter.

It is easy to associate a diffraction order $m$ for most of the resonances supported by the grating described above (see Fig. 5). However, starting from almost zero and progressively increasing the height $h_{\mathrm{G}}$, the resonance that we first attributed to be of order $m=1$ exhibits a splitting (see Fig. 5).

The field maps on Fig. 5 clearly show that one of the resonances corresponds to a cavity-like resonance, which is entirely located in the grooves of the grating, i.e it is reflected back and forth horizontally. The other resonances are cavity resonances of SPs that propagate on top of the grating and are reflected by the edges of the metal. We label those resonances with $1\mathrm{b}$ (where ’b’ stands for ’bottom’) and $1\mathrm{t}$ (where ’t’ stands for ’top’), respectively. Since the thickness of the dielectric is different for both kinds of plasmons, their wavevector cannot be the same, therefore the splitting. We finally chose to take $h_{\mathrm{G}}=68$ nm.

Were the depth of the grating much larger, the 1t mode would enter a spoof SP regimePendry et al. (2004); Khanikaev et al. (2010); Rusina et al. (2010); Kats et al. (2011); Pors et al. (2012); Ng et al. (2013); Shen and Cui (2014) where the grating could be considered as an effective medium. In the present case, the conditions to be in this regimeKats et al. (2011) are however not fulfilled, as the pitch is roughly one third of the wavelength at most – for the 1t mode. In addition, the grating is too shallow to allow for spoof SPs to be excited, making 1b and 1t modes more similar to the hybrid SPs evoked in the literaturePendry et al. (2004). As the following will show, the 1t mode can be considered as insensitive to nonlocality, suggesting that spoof SPs may be equally insensitive.

III.2 Influence of nonlocality on the grating reflectivity

Now that the grating has been designed and its physics is well understood, we use a software taking spatial dispersion into account in order to assess its influence on the reflectivity of the grating. This is necessary, since analytical solutions do not exist anymore and we have to rely on a numerical algorithm. Here, we use a DGTD method for the simulations of the grating. The result is shown in Fig. 6 for wavelengths ranging from $\lambda=500$ nm to $\lambda=2000$ nm. The small discrepancy between local DGTD and RCWA is due to a known problem of the latter in very peculiar conditionsLi and Granet (2011) and occurs only far off the resonances. We can clearly identify 5 dips due to SP excitation corresponding to 4 diffraction orders, the first order being split as discussed above. The results show a considerable influence of the spatial dispersion, which is significantly stronger than what has been predicted theoretically. This can be attributed to the grating itself and to the shift of the resonances towards large wavelengths (and larger wavevectors for the SP). The resonances linked to the diffraction orders $m=4$ and $m=2$ (see inset on Fig. 6) experience blueshifts of respectively 1.8 nm and 6.5 nm compared to the local prediction. Especially for $m=2$ this is clearly higher than the expected shift of 3 nm (see Fig. 4(b)). For $m=4$ the shift is slightly smaller than expected. Finally, we can identify two wavelength bands of interest. The first one for $\lambda=[550,800]$ nm (see inset) shows the highest difference between local and nonlocal simulations. It comprises two blueshifted diffraction orders, and the whole response between the two is also clearly affected by nonlocality - making this regime a good choice to estimate the parameter $\beta$. The second region of interest is the rightmost part of Fig. 6 for $\lambda=[800,2000]$ nm. It comprises the two suborders $1t$ and $1b$, which are less sensitive to nonlocality. For this reason, we believe that this region is not useful to probe nonlocality, but well-suited for a geometrical parameter characterization. We have to keep in mind that any estimation of $\beta$ relies merely on a comparison between material models. Such comparisons are very sensitive to the geometrical parameters. Since we are trying to measure discrepancies of the order of $1\%$ of the wavelength, we must ensure that nonlocality will not be concealed by uncertainties on local parameters.

IV Telemetry and parameter estimation

In the previous section we have theoretically predicted that the metallic grating that we have designed is sensitive to nonlocality. The goal of the present section is to present the challenges which would have to be faced by experimentalists and to propose a methodology relying on the solution of inverse problems and Uncertainty Quantification (UQ), which could be used to assess realistically the presence of nonlocality and to estimate the crucial $\beta$ parameter.

In order to take into account the fact that, whatever the geometry and the imperfections of the grating, the spectra will be noisy, we have generated an artificial but realistic reflectance spectrum. We use nonlocal simulations and then added a noise whose characteristics are similar to actual experimental spectra provided by experimentalistsJamadi et al. (2018). The noise has been chosen with a correlation length of $0.001$ nm and a maximum difference of 0.04 with the original, unperturbed spectrum. We have then tested whether the parameters of the grating and of the model could be retrieved despite this realistic level of noise.

IV.1 Post-fabrication telemetry

IV.1.1 Grating parameters

The grating parameters would have to be determined before any dielectrics is deposited. A natural way of determining these geometrical parameters would be to rely on a Scanning Electron Microscope (SEM) or an Atomic Force Microscope (AFM) which would both allow to directly measure the precise parameters for each of the grooves of the grating. However, for a structure without any dielectrics, the impact of nonlocality is negligible - which means that the geometrical parameters can be determined using telemetry, without any assumption on the non local parameter, i.e with a local model. The optical response of a grating actually depends more on average geometrical parametersKaliteevski et al. (2002), as the grooves may be all slightly different. Telemetry has the advantage of allowing the determination of these average parameters (especially the period of the grating), which will constitute our geometric model.

We generate an artificial measurement spectrum, that we denote $R_{\mathrm{meas}}(\lambda)$. Furthermore, it is possible, for given $d_{\mathrm{G}}$, $a_{\mathrm{G}}$ and $h_{\mathrm{G}}$ to compute a theoretical spectrum $R(\lambda)$. We define the distance between the two spectra as

[TABLE]

The integration interval $[\lambda_{1},\lambda_{2}]$ will be chosen differently in the following, depending on which parameters have to be retrieved. We performed multiple optimization runs with different algorithms (particle swarm, pattern search and a derivative-based optimizer combined with a kriging-based meta-model) in order to find the parameters ($d_{G}$, $a_{G}$ and $h_{G}$) which would produce the minimum distance between the theoretical spectrum and the artificially generated one. We underline that we constrained the optimization to look for geometrical parameters in intervals which would be in accordance with the precision of the etching process. The fabrication tolerances of state-of-the-art nano-processing still lead to a priori tolerances of about $\pm 5\ \mathrm{nm}$ for the etching process Langereis et al. (2009) and about $\pm 11\%$ Chen (2015) for the dielectric deposition.

Considering spectra between 400 nm and 1200 nm to retrieve the parameters, we have found that the derivative-based optimization method performed best, followed by pattern search and particle swarm, the latter seemingly being less appropriate for this type of optimization problem. The geometrical parameters could be retrieved with an excellent accuracy despite the noise (see Fig.7).

IV.1.2 Dielectric thickness

The dielectric thickness has to be determined by telemetry. While the previous step can be performed with purely local simulations, here nonlocality clearly plays a role. This time, our artificially measured data are generated using the nonlocal spectrum corresponding to the right part of Fig. 6 i.e $[\lambda_{1},\lambda_{2}]=[800,2000]$ nm. Since the resonances 1t and 1b are almost insensitive to nonlocality, trying to retrieve $h_{\mathrm{D}}$ by minimizing the distance between the measured non local spectrum and a local in this wavelength range makes sense and actually yields results that are very close to the real value (the retrieved value was $h_{\mathrm{D}}=84.8\ \mathrm{nm}$, with the real value being $h_{\mathrm{D,init}}=85\ \mathrm{nm}$).

IV.2 Geometric uncertainty vs. nonlocality

The geometric telemetry has led to the fabricated geometries in Table 1, where, for each parameter $z$ with initial values $z_{\mathrm{init}}$ used to create $R_{\mathrm{meas}}$, $\overline{z}=z_{\mathrm{opt}}$ is the mean value and $\delta_{z}={\left|\overline{z}-z_{\mathrm{init}}\right|}$ the maximum deviation.

Given the uncertainties due to the retrieval process (with other means of measuring the geometric parameters these uncertainties would likely be of the same magnitude), we need to be sure that we will be able to distinguish the impact of nonlocality from an unavoidable small error in the determination of the geometric parameters. We use Uncertainty Quantification (UQ) to provide us with answers and, in order to estimate the impact of nonlocality we use, again, a $\beta$ value from the literatureCiracì et al. (2012).

We recast the geometrical parameters $h_{\text{G}}$, $d_{\text{G}}$, $a_{\text{G}}$ and $h_{\text{D}}$ as random variables (RVs) following uniform distributions $\mathcal{U}\left[\overline{z}-\delta_{z},\overline{z}+\delta_{z}\right]$ (see 1). We then perform an UQ Smith (2014) study, in order to estimate the expectation value and variance of the reflectivity in two resonance regions corresponding to the diffraction orders 2 and 4. In particular, we consider $25$ wavelengths in the range $[550,600]\ \mathrm{nm}$ ($m=4$) and $50$ wavelengths in the range $[700,800]\ \mathrm{nm}$ ($m=2$).

Since the underlying computational model is a complex one, we rely on black-box UQ methods, i.e. the model and its numerical solvers are used without any modifications. In the context of the present work, we employ a spectral method Le Maitre and Knio (2010); Xiu (2010), in particular the stochastic collocation method Babuska et al. (2010); Nobile et al. (2008); Xiu and Hesthaven (2005), taking advantage of a number of factors. First of all, we deal with a small number of RVs, therefore the costs of the method remain affordable. Moreover, numerical tests indicate that the dependence of the reflectivity upon the RVs is smooth, which is a prerequisite for fast convergence. Finally, we assume that the RVs are mutually independent, which greatly simplifies the method’s implementation.

The first step is to choose $M$ different set of values for the geometrical parameters $\mathbf{z}=\left(h_{\text{G}},d_{\text{G}},a_{\text{G}},h_{\text{D}}\right)$, called the collocation points. For each wavelength mentioned above, the reflectance $R(\lambda)$ can be considered as a function $f(z)$ which is approximated by

[TABLE]

where $\mathbf{z}^{(m)}$ are realizations of the random vector (the collocation points) and $\Psi_{m}$ are multivariate Lagrange polynomials. The collocation points are based on univariate Clenshaw-Curtis quadrature nodes and are produced by Smolyak sparse grid rules Bungartz and Griebel (2004); Smolyak (1963). The multivariate polynomials are formatted as products of univariate Lagrange polynomials, defined by the corresponding univariate Clenshaw-Curtis nodes.

The reflectivity’s expectation value $\mathbb{E}\left[R\right]$ and variance $\mathbb{V}\left[R\right]$ can now be estimated by post-processing the approximation terms of (14). Starting with the definitions

[TABLE]

where $\varrho(\mathbf{z})$ is the joint probability density function, we approximate the corresponding integrals with the multi-dimensional Gauss quadrature formulas

[TABLE]

where $w_{m}$ denote the corresponding quadrature weights. We use the UQ study results in order to estimate $\pm 2\sigma$ intervals around the optimized local reflectivity curve, where $\sigma=\sqrt{\mathbb{V}[R]}$ refers to the standard deviation. The results corresponding to each resonance area are presented in Fig. 8 and Fig. 9 (in blue), respectively. Since the nonlocal DGTD simulations are too expensive for an UQ of the same kind as we have performed for the local model, i.e. the evaluation of (17) and (18), we rely on a min-max study. Here, min-max represents solver calls for all interval bound combinations (in orange) as depicted in the corresponding figures.

According to Fig. 8, a clear measurement of the resonance $m=4$ is almost impossible due to the small difference between the local and nonlocal curves. Nevertheless, the second resonance, i.e. $m=2$ (see Fig. 9) is significantly more sensitive to nonlocality and stays distinguishable. For the sake of robustness, we have also performed simulations (data not shown) with the lowest theoretically acceptable value for $\beta=0.85{\times 10^{6}}$ m.s*-1*Hiremath et al. (2012), which still guarantees a blueshift of $5\ \mathrm{nm}$ for the second resonance, i.e. stronger than the geometric uncertainty.

IV.3 Model calibration

Knowing the geometry and estimating the impact of nonlocality to be greater than geometric uncertainties, we now want to extract the nonlocal material parameter $\beta$. We underline that different theoretical expressions can be found in the literature for this constantCrouseilles et al. (2008); Hiremath et al. (2012); Raza et al. (2011), so that the theoretically acceptable values for $\beta$ may lie between $\beta_{\mathrm{min}}\approx 0.85{\times 10^{6}}$ m.s*-1* Hiremath et al. (2012) and $\beta_{\mathrm{max}}\approx 1.4{\times 10^{6}}$ m.s*-1*. Raza et al. (2011). However, the experimental results available so farCiracì et al. (2012, 2014) point consistently towards a value close to the upper estimation of $\beta=1.35{\times 10^{6}}$ m.s*-1*. This underlines how important the determination of $\beta$ can be and explains why we have considered this value so far.

In order to estimate how precise our estimation of $\beta$ could be with the grating setup, we proceed in the same fashion as for the geometric telemetry but we use a wavelength range of $[550,800]$ nm. Using DIOGENeS and DGTD Schmitt et al. (2018), we find the $\beta$ value which minimizes the distance between $R$ and $R_{\mathrm{meas}}$. The geometric size of the structure, in combination with the very small effective wavelengths and the short interaction range of nonlocal effects, which is in the range of several nm at the metallo-dielectric interface, result in a computationally expensive procedure. The solution of the inverse problem can be significantly accelerated by meta-model-based optimization algorithms. We have used the kriging (Gaussian process) meta-model in combination with a derivative-based optimization implementation of the FAMOSADuvigneau optimization toolbox. We find $\beta_{\mathrm{inverse}}=1.385{\times 10^{6}}$ m.s*-1*, which is reasonably close to the $\beta_{\mathrm{init}}=1.35{\times 10^{6}}$ m.s*-1* (the value used to generate $R_{\mathrm{meas}}$), indicating that the value of $\beta$ can be retrieved with an error smaller than 10%.

V Conclusion

We have first shown that, in the framework of the hydrodynamic model, SPs can be sensitive enough to spatial dispersion – provided that the dielectric considered has a sufficiently high permittivity, like $\mathrm{TiO}_{2}$. Such a conclusion is in contrast with previous works that suggested that the impact of spatial dispersion could be too difficult to measure optically – which is only true for an interface between metal and air.

Since there is a link between such guided modes and the localized resonances of metallic nanoparticles, this lead us to expect an impact of nonlocality on essentially any metallic nanoparticle with a radius much larger than 20 nm, for which the geometry is more likely to be well controlled, buried in a high index medium. Using Mie theory, we estimate the blueshift brought by nonlocality to be at least of 5 nm in wavelength in TiO2, an effect that could potentially be observed experimentally.

Since high wavevector guided modes cannot be excited using prism couplers, we have then studied how the SPs can be excited using a grating coupler. We have shown, using state of the art numerical tools, that such a structure would allow the observation of spatial dispersion by means of blueshifted resonances up to almost 7 nm (around 1% of the wavelength). Using uncertainty quantification and inverse problem solving, we have identified which resonance precisely could be used to estimate the main parameter of the hydrodynamic model and shown how such an estimation could be made. We underline that such a procedure could well be applied to nanoparticles as well.

As already evoked in several earlier worksRaza et al. (2015b); Pitelet et al. (2017, 2018) and clearly demonstrated in the present work, any plasmonic structure surrounded by a high refractive index medium like TiO2 will only be accurately described if spatial dispersion is taken into account.

We believe that, by proposing a structure with realistic parameters and a procedure to carefully estimate the impact of spatial dispersion, this work will pave the way for future experiments that shall give reliable answers to the community on the limits of Drude’s model for plasmonics and its potential replacements.