Nonlocal response of planar plasmonic layers

Nonlocal interactions of plasmonic nanostructures are currently being intensively investigated. It is generally accepted that nonlocal interactions are most pronounced on structures with nanometer unit sizes and affect the shape of spectral functions of characterizing quantities in the resonance region. Numerical analysis of nonlocal phenomena is very complicated and it is advantageous to use it for the analysis of geometrically more complex structures. For some simple structures it is possible to find analytical solutions which can then be used, among other purposes, to build semi-analytical approaches for the analysis of some more complex structures. This article is focused on the efficient description of nonlocal manifestations of the planar plasmonic layers using a hydrodynamic model. The originality is particulary focused on the extension from the single nonlocal layer towards the more general case of two adjacent nonlocal layers (bilayers). The results for both single layer and bilayer case are presented and discussed in detail. Also, one of the here presented important technical results is the novel presentation of the transfer matrices for the nonzero normal component of the current densities at the interfaces of the plasmonic layer.


Introduction
At present, the central interest of photonics and plasmonics is mainly via objects belonging to the nanodimensional scale.In this region, new effects, such as nonlocality, appear on the border of classical and quantum physics.The classical method of description based on the local approximation may therefore be inaccurate in certain circumstances.So far, new methods of description called nonlocal models (Wubs 2015;Xue et al. 2015;Ginzburg 1286 Page 2 of 33 and Zayats 2013) are still being developed.A number of possible approaches already exist, some of which are based mainly on quantum description, such as quantum nonlocal models (Ciracì and Della Sala 2016;Varas et al. 2016;Baghramyan et al. 2021;Yan 2015;Yan and Mortensen 2016).Due to the complexity of quantum models, here, we will use the standard hydrodynamic model (HD model or HDM) (Hiremath et al. 2012;Toscano et al. 2013;Raza et al. 2015) for our nonlocal thin film analysis.However, our model can be further modified to suit some HDM modifications.An example is a modification involving the effect of electron gas diffusion (Mortensen et al. 2014;Svendsen et al. 2020).The standard HD model has already been used to simulate the interaction of incident optical radiation with structures of basic geometric shapes such as a spherical particle, cylinder or metallic surface (Raza et al. 2011;Toscano et al. 2012;David and García de Abajo 2011;Moreau et al. 2013), respectively.It has also been proven that the so-called curl-free HDM approximation is inaccurate and therefore precise solutions need to be found (Raza et al. 2011).Our article is focused on the analysis of nonlocal properties of thin plasmonic layers.Compared to calculations based on the classical theory, solving the interaction described by the HD model is a much more difficult problem even for these elementary structures.
Simultaneously, considerable attention is being paid to metallo-dielectric layers, as it can be assumed that they may have unique properties in specific configurations and spectral ranges.As an example, we can mention hyperbolic metamaterials (Poddubny et al. 2013;Li and Khurgin 2016) which can theoretically transform evanescent waves into propagating waves, and thus allow subwavelength imaging or high-sensitivity sensors (Guo et al. 2020).Another theoretically possible realization that can be achieved via metallo-dielectric layers is an environment with a negative refractive index (Verhagen et al. 2010;Bang et al. 2019).Further, the development of sensitive sensors based on the specific properties of the surface plasmons (Roh et al. 2011;Nguyen et al. 2015;Uddin et al. 2021;Barchiesi et al. 2022) is very attractive, too.Clearly, the design of structures and their proper operation towards the mentioned applications always depend on the most accurate calculations.This corresponds to the great effort aimed at developing simulation models.Many nonlocal analytical or numerical models have recently been developed to analyze metallo-dielectric layers, each with its own usability or limitations.The most basic model is based on the modified Fresnel relations, i.e. on the nonlocal coefficients of reflection and transmission at the interface of two material half-spaces (Feis et al. 2018;David et al. 2013).Actually, modified Fresnel coefficients represent also the basis of many more complicated models.An example is a model analyzing the properties of a surface plasmon on a nonlocal metamaterial (Feis et al. 2018).Under certain metal layer thickness constraints, nonlocal properties of metallo-dielectric layers can be investigated (Pitelet et al. 2018), too.Modified Fresnel coefficients are also an important part of the surface HDM (Yang et al. 2018); this model takes into account the so-called spill-out phenomenon, i.e. the fact that a non-zero electron density can flow out of the metal interface (Skjølstrup et al. 2018;Li et al. 2015;Toscano et al. 2015;Taghizadeh and Pedersen 2019;Mortensen et al. 2021).Some other nonlocal models of metal layers have also been published, focusing mainly on the analysis of the waveguide properties of these layers and the surface plasmon propagation (Raza et al. 2013;Xiao et al. 2015).At the same time, nonlocal numerical methods based on rigorous formalism of scattering matrices (Benedicto et al. 2015) were developed which thus enabled the analysis of interactions of more complex structures.
In this contribution, our model, analysing plasmonic planar layers, is based on an analytical approach called the transfer matrix method.The specific novelty of our approach is the extension of the algorithm for plasmonic bilayers, i.e. two metallic films in direct contact.Therefore, the main benefit lies in the possibility of analysing such structures as IMMI (insulator-metal-metal-insulator) and other combinations of metallo-dielectric layers.Also, this model could be modifiable to include some selected quantum phenomena connected with a metal layer interface or to investigate nonlocal interactions in other materials, such as semiconductors.Recently, in fact, a HD model for semiconductor materials has been derived (Golestanizadeh et al. 2019), and thus our bilayer model based on the transfer matrix method would possibly enable the incorporation of nonlocality even for semiconductor layers.
This article is divided into eight main sections as follows.After this introduction, the Sect. 2 presents the fundamental relations related to the standard HDM, and the Sect. 3 is devoted to the revision of the interaction of the incident plane wave with a single plasmonic layer.The content of the Sect. 4 is the definition of the transfer matrix and its specific shape for such single layer.The correctness of the calculations of the analytical method presented here is also checked in this section.In the following Sect.5, the case of the metal bilayer, i.e. the IMMI structure, is analyzed in a similar way, and the obtained analytical solution of the transfer matrix is presented.Finally, the Sects.6 and 7 are focused on assessing the effect of nonlocal response on the reflectance of the gold and silver layers and the gold-silver bilayer, respectively.The paper is concluded in the Sect.8.The supplemental document to the article is devoted to some interesting cases of the predicted significant nonlocal response of gold and silver IMI (insulator-metal-insulator) and IMMI structures.

Fundamental relations
The nonlocal response of a material can be generally described as that the material properties at the considered point in space are also influenced by the physical quantities in its vicinity.Thus, the spatial response of a material is not a delta function which in the Fourier space of wave numbers means that the properties of the material are functions of these wave numbers.The HD model is based on the idea that it is possible to describe the dynamics of an electron gas similarly to a fluid with non-zero compressibility.As a result, HDM can be understood as a more sophisticated Drude-Lorentz model, in the sense that the electron gas moves on a positively charged background of metal ions, these ions form a dielectric portion of the metal permittivity and the remaining permittivity is generated by the excited current density response.The sum of these two permittivities gives both the permittivity of the metal for the transverse electromagnetic field but also the longitudinal permittivity for the longitudinal field.From the point of view of the classical theory, which takes into account only the transverse electromagnetic field, the existence of the longitudinal field and the longitudinal component of permittivity can be understood as a nonlocal response.The existence of a nonlocal response was expected mainly in localized nanostructures and therefore great attention was focused in this direction.Some approximations of the standard HD model have also been investigated to simplify the calculations, but these may lead to inaccurate predictions, such as new unexpected resonances below the plasma frequency.The hydrodynamic model in its linearized form Hiremath et al. (2012) is thus proven tool for describing nonlocal interaction in metallic media.It combines Maxwell's equations and, at the same time, the dynamics of the electron gas in an electric field.Simultaneously, the behavior of the electron gas is similar to a liquid, satisfies the continuity equation, and has non-zero compressibility.The standard form of the hydrodynamic model is defined by the following set of two equations (Raza et al. 2011;Hiremath et al. 2012) 1286 Page 4 of 33 The symbol b ( ) in HDM means the permittivity of a positive metal background.The Table 1 below lists the HDM parameters for gold and silver materials (Raza et al. 2015;Hiremath et al. 2012) that are used in the calculations in later sections of this article.
In general, one can assume that the total current density is the sum of the transverse and longitudinal components, given by the transverse E T (r, ) = ∇ × E (r, ) , J T (r, ) = ∇ × J (r, ) and longitudinal fields E L (r, ) = ∇ξ E (r, ) , J L (r, ) = ∇ξ J (r, ) .
Using this assumption, one can obtain, after simple algebraic manipulations, expressions for transversal and longitudinal wave numbers (David and García De Abajo 2011) In Eqs. ( 3) and ( 4) we introduced the following notation: where k T,y and k T,z denote the component of the transverse field wave vector in the directions ŷ and ẑ of axes of cartesian system, respectively.Similarly, we introduce the notation ‖k L ‖ ≡ k L and ‖k 0 ‖ ≡ k 0 .In Eq. (3), m ( ) corresponds to the permittivity of the material, so the transverse wave number has the same value as the classically defined wave number by the standard permittivity.The values of m ( ) can be determined for different materials by measurement, theoretically, using various models, or simply by a polynomial extrapolation of the tabulated values (Palik 1985).If we now define the conductivity separately, for the transverse and for the longitudinal field respectively, by the relations J(r, ) = (r, ) * E(r, ) , where symbol * means convolution through spatial variables, we can find from Eq. ( 2) the following relations in the Fourier space (3) Similarly, it is possible to define nonlocal formulas for permittivity.We can modify Eq. ( 1) by converting its right-hand side to the left and using the relations for conductivity, thus obtaining new relations for longitudinal and transverse permittivity (Raza et al. 2013) as It can also be shown that the value of L ( , k L ) is zero (Wubs 2015), this follows immedi- ately from (4).The origin of the nonlocal response is the existence of a longitudinal field which is evident from Eq. ( 6).The relation ( 7), is obvious at first sight, also be obtained immediately from (3).From the knowledge of the material dispersion of m ( ) , the permit- tivity of the positively charged background metal b ( ) can be easily determined from ( 7), which, in addition to (4), also appears in other important relations.

Revision of HD model for a single layer
A schematic of the interaction model is shown in Fig. 1.The plane of incidence of the transverse magnetic (TM) polarized plane wave is determined by the ẑ and ŷ axes.The first medium from above is a dielectric defined by the refractive index n 0 and the imped- ance z 0 = √ ∕ , where and denote relative permeability and permittivity of the first medium, respectively.The vacuum impedance defined by the vacuum values of permittivity and permeability is denoted by the symbol 0 , where 0 = √ 0 ∕ 0 .Another layer depicted with yellow color is a metal with a complex refractive index n 1 and impedance The last layer is again a dielectric with parameters n 2 , z 2 .Let us denote the norms of wave vectors of incident wave ‖k − 0 ‖ and the reflected wave ‖k + 0 ‖ as k 0 .Quantities with a lower index T will be used to describe transverse fields, similarly a lower index L will be bound to a longitudinal field.The incident wave in the first medium forms an angle 0 with the ẑ axis, the transmitted wave then makes an angle 1,T , with the ẑ axis in the sec- ond medium.
For the geometry of the planar layer considered here and the assumption of a steady state, the solution of Maxwell's equations in the form of a harmonic spatio-temporal evolution can be considered.So let us assume that we can define both transverse and longitudinal electric fields as follows: E T (r, ) ≡ T e i(k T r− t) , E L (r, ) ≡ E L e i(k L r− t) .The field amplitudes are key quantities for our further considerations, so let us introduce the following notation: ‖E T ‖ ≡ E T , ‖E L ‖ ≡ E L .From Maxwell's equations, we can easily derive the relations between the amplitudes of the magnetic and electric fields of a plane wave in the layer of the index a (a = 0, 1, 2) , in the form z a −1 E a = 0 H a .These equations together form the basics for our model.All the relevant further mentioned relations can be derived from the geometric relations between the electric and magnetic vectors of the waves shown in Fig. 1.However, the question of how to determine the 1286 Page 6 of 33 individual angles 1,T and 1,L of refraction has not yet been resolved.The answer is quite simple, since the electric and magnetic fields have the same frequency, it is enough to deal only with the spatial distribution of the fields on the interface.In other words, the projection of the wave vectors into the ŷ axis must be the same for all fields.It is actually a gener- alization of Snell's law of refraction, which can be written in the following form: The values of cos 1,T and cos 1,L can be determined directly from the trigonometric rela- tions and from the knowledge of the value of sin 1,T and sin 1,L , respectively.The formu- las can be further modified if we define: k T = ( ∕c)n 1,T , k L = ( ∕c)n 1,L .Thus we obtain relations If we were only concerned with a lossless dielectric layer with a refractive index n 1,T instead of a metal one, then for total reflection when (n 0 ∕n 1,T ) sin  0 > 1 applies, the cos 1,T will be imaginary and only the evanescent wave will reach the dielectric layer.In the case of a metal layer, the cos 1,T always has a complex value.Now we can determine the values of the tangential components of electrical and magnetic fields at the interfaces of the metal layer purely from geometric considerations.Due to Fig. 1 Scheme of the nonlocal interaction of an incident TM plane wave with a metal layer the continuity of the tangential components of the fields at the interfaces of the individual layers, all fields must have the same harmonic course along the interface.Therefore, we do not have to deal with this dependence, and it is sufficient to consider only the phase change between the total transmitted and reflected waves from both the first and second interfaces.
The transversal wave E s1,T (r, ) , propagating from the second interface of the metal layer to the first one, is the sum of all possible cases of multiple reflections from the first interface with complex reflection coefficient r 0,1 and from the second interface with the reflection coefficient r 1,2 .The relation for the amplitudes of total reflected transverse wave from the lower interface, that reached the upper interface, is expressed in the form of equation: In a similar way, we can determine the magnitudes of the total reflected and transmitted fields at the upper and lower interfaces for both the transverse and longitudinal fields.
The reflected field from the lower interface, that reaches the upper interface, differs from the originally reflected wave only by the phase factor.This phase change has a value of 1,T = k T,z d in the case of a transverse field, where k T,z has the meaning of the z-component of the wave vector.Using the propagation angles of individual fields from Eq. ( 9), we can determine this phase change 1,T for the transverse and 1,L for the longitudinal field in the form of From the values of phase changes, we can obtain important relations between individual norms of fields amplitudes, given with equations From simple geometric considerations, according to Fig. 1, and from the knowledge of the ratio between the magnitudes of electric and magnetic field intensity vectors for a plane wave, we can determine the norms of amplitudes of the total tangential, (in direction of ŷ and x axes), electric and magnetic fields E 0,1 = ‖E 0,1 ‖ , H 0,1 = ‖H 0,1 ‖ at the upper interface and similarly for electric and magnetic fields E 1,2 = ‖E 1,2 ‖ , H 1,2 = ‖H 1,2 ‖ at the lower interface of the metal layer, in the form of equations ( 12) 1286 Page 8 of 33 Similar to the above, H i0 , H r0 , E t1,T , E s1,T , E t1,L , E s1,L , E i1,T , E r1,T , E i1,L , E r1,L are the norms of the relevant vector amplitudes.
To proceed further, it is necessary to express the current density in terms of the transverse and longitudinal electric fields, which then gives the boundary conditions for the current density in a convenient form.From the equations proposed here, the nonlocal transfer matrix of the metal layer can then be determined, which will be discussed in the following sections.In general, the total electric field E tot (r, ) can be considered as the sum of the transversal and longitudinal field: E tot (r, ) = E T (r, ) + E L (r, ) .As mentioned in the Sect. 2 of the article, the transversal and longitudinal electric fields can be defined as E T (r, ) = ∇ × E (r, ) and E L (r, ) = ∇ξ E (r, ) .Using the known operator identities ∇ ⋅ (∇×) = 0 and ∇ × (∇) = 0 , it immediately follows: ∇ ⋅ E T (r, ) = 0 and ∇ × E L (r, ) = 0 .Substituting the total electric field E tot (r, ) into the first Eq.( 2) of the HDM, after simple modifications we obtain: If we further consider that for a transverse field which takes the form of a plane wave, T E T (r, ) , then we can easily obtain from ( 21) the following general relation where the coefficients F ≡ F( ) and S ≡ S( ) are given by the expressions: Let k T,y and k L,y denote the y-components of the transverse and longitudinal wavenum- ber vectors, respectively.From (9) it automatically follows: k T,y = k L,y .In accordance with the already established convention, let us define the current densities at the upper and lower interfaces of the metal layer oriented perpendicularly to the interfaces in the positive direction of the ẑ axis as J u (y, ) ≡ u e i(k T,y y− t) and J b (y, ) ≡ b e i(k T,y y− t) , respectively.Next, let us introduce the norms (numerical values) of normal components of vector amplitudes of current densities with the following symbols: ‖J u ‖ ≡ J u and ‖J b ‖ ≡ J b in a similar way as for the previous quantities.We can now use Eq. ( 22) to define the boundary conditions for the normal component of the current density at the interfaces of the metal layer in the form The first boundary condition ( 23) determines the value of the normal component of the current density directed outside the metal layer at its upper interface.Analogously, the value of the current density flowing out of the lower interface of the metal layer is defined with (24).If the metal layer is surrounded by dielectric materials, as shown in Fig. 1, J u and J b are zero.On the other hand, if one of the surrounding environments is metal, then current may flow at this metal-metal interface.Such generally defined boundary conditions allow us to study another very interesting case where two metal layers (or bilayers) are in contact to each other allowing nonlocal and even quantum effects (such as a tunneling current) to occur.First, however, it is necessary to find a general mathematical formula for the transfer matrix between the analyzed layers.We can perform this step without specific knowledge of norms of normal current amplitudes to the interface J u and J b .

Transfer matrix of a metal layer
The transfer matrix method is a very useful method for determining the optical properties of isotropic layered structures (Macleod 2001;Born et al. 2019;Acquaroli 2018).This method assumes a harmonic time and space dependence of all fields.A transfer matrix can generally be defined as a matrix that gives the relation between the magnitudes of tangential electric and magnetic field at the interfaces of a given layer.The transfer matrix (Born et al. 2019) for the jth layer and TM polarization has the following form: In (25), the term j is equal to the phase change for the transversal field according to (13), and the symbol j is a substitution for z j −1 cos −1 j,T .A very important characteristic property is the transmission t and reflection r coefficient of the given layer, these quantities can be defined for TM polarization by relations Consider again the case of one layer surrounded at the top by an environment of index 0 (superstrate) and at the bottom by an environment of index 2 (substrate).In order to be able to determine the transmission and reflection coefficients, we have to state additional relations for the field amplitudes in the tangential direction to the interface, which are as follows: If we use the general expression of the transfer matrix and the supplementary relations ( 26)-( 29), then after simple algebraic modifications, we obtain the important relation: Equation ( 30) actually represents a matrix notation of a system of two equations with two unknowns r and t.By solving this system, we obtain the required relations for the transmission and reflection coefficients (Chilwell and Hodgkinson 1984) of the considered plane layer for TM polarization in the general form: (25) .
1286 Page 10 of 33 In the case of the classical local model, we can substitute the individual terms of the transfer matrix ( 25) into ( 31) and ( 32).Since the relations for the reflection and transmission coefficients were derived in a very general way, they can also be used for the nonlocal HD model.It is only necessary to know the specific shape of the transfer matrix of the nonlocal medium.
Let us now consider the problem of how to determine the transfer matrix for the HD model.This matrix can be derived from Eqs. ( 17) to ( 20) and the boundary conditions for the normal component of the current density ( 23) and ( 24).The mentioned equations can be seen as a system of six equations for six unknowns: E r1,T ,E i1,T , E r1,L , E i1,L , E 0,1 , 0 H 0,1 .The remaining unknown amplitudes are related to the unknowns defined via relations ( 14) and ( 15), the quantities E 1,2 , 0 H 1,2 , J u and J b are appearing in the equations as parameters.Our goal is to identify all unknowns using these parameters only.In the simplest case, the metal layer is surrounded on both sides by a dielectric, electrons from the metal cannot enter the dielectric under normal conditions, and thus the magnitude of the normal component of the current density at both interfaces is zero, J u = J b = 0 .Even in this simplest case, the derivation of the transfer matrix is relatively laborious, so we will only present the result here.The individual elements of the transfer matrix are represented with rather complicated expressions, so we will present them in abbreviated form together with substitutions allowing a simplified notation: giving finally the matrix elements as Thus, relations (37)-(40) determine the transfer matrix for the metal layer and the incident TM plane wave.
Page 11 of 33 1286 The publication (Pitelet et al. 2018) represents a valuable opportunity to compare the calculations obtained by our model with a nonlocal model focused on the nonlocal response of metallo-dielectric layer structures where the nonlocal deviation of reflectance is analyzed on the MIM structure, schematically shown in Fig. 2.
The proposed comparison of local and nonlocal reflectance, as a function of the incident angle i of TM plane wave, of the structure shown in Fig. 2 is presented in Fig. 3.For a more convenient visual comparison, we have taken Fig. 3a, c from Pitelet et al. (2018), our calculations represent Fig. 3b, d.
As can be seen, our model assumes a slightly smaller nonlocal deviation.The difference in results is most likely due to the fact that our model does not neglect the influence of the transverse current density on the overall nonlocal response.
It should be noted that the method in Benedicto et al. (2015) alternatively uses the effective polarization quantity P f (r, ) to describe the interaction instead of the current density J(r, ) which can be introduced by Maxwell's equation in the form: Based on the commonly stated form of Maxwell's equation and ( 41), the constitutive relation can be derived: i P f (r, ) = J(r, ) .As can be read in Benedicto et al. (2015), the assumption P f (r, ) = 0 is introduced for the case of the transverse field (for ∇ ⋅ E(r, ) = 0 ) which is analogous to the approximation J(r, ) with only the longitudinal field, neglecting the transverse component J T (r, ) of the current density.
The transfer matrix discussed here represents an important basis for further considerations such as the analysis of two coupled nonlocal metal layers.This case will be discussed in the next chapter.

Transfer matrix of two connected metal layers
In the classical view of linear spatially local response (LRA-local response approximation), all properties of a metallic material are included in the permittivity of this material which can be understood as a function of frequency (Maier 2010;Mortensen 2013).From the point of view of the HD model, the overall response of a metal is a combination of the dielectric properties of the positively charged background and the response of the electron gas in the form of induced current densities.The nonlocal deviation is due to the interplay of two differently behaving metal subparts.For these reasons, it can be assumed that the boundary conditions of the current density should be significantly reflected in the overall behavior of the induced current density and thus of the metal layer as a whole.The simplest structure on which the influence of boundary conditions of current densities can be investigated is a metal bilayer surrounded by dielectric media, thus the investigation of this structure is logically our next goal.Before we start analyzing the problem, it is necessary to introduce proper index notation.Let us denote the first metal layer by index 1, the metal layer below it by index 2, according to Fig. 4, and the quantities belonging to the given layer by the same index as the given layer.As a first step in solving the above problem, we have to solve the general case of a metal layer with nonzero boundary conditions ( 23) and ( 24).Nonzero boundary conditions change the relation between the amplitudes of the tangent fields at the edges of the metal layer.In the case of the first layer, we can write it as follows: The matrix M 1 in (42) has the same form as in the already discussed case of isolated metal layers.Therefore, the individual elements of M 1 are defined by the expressions (37) to ( 40), but the relation between the tangent components of the fields is now supplemented by the vector N 1 ensuring the fulfillment of the new boundary conditions for the current density.In the general form, the vector N 1 is defined by the expressions: (43) S sin 1,L .Such general relations ( 43) and ( 44) can be used to describe the middle layer in a threelayer metal structure.However, in the case of a metal bilayer, the first metal layer is surrounded from above by a dielectric environment and therefore J 1u = 0.
For later consideration, it is advantageous to introduce a matrix expression of the vector N 1 , where matrix elements Γ 1 ,Π 1 ,Θ 1 , and Ω 1 are given by Eqs. ( 43) and ( 44), as follows: Now, the problem remains how to determine the normal component J 1b of the current densities to the interfaces.The answer will give us a simple reasoning.Consider a metal bilayer as in Fig. 4 and assume in both materials a different magnitude of the current density in the tangential direction, but the same in the normal direction at a common interface.In both metals, the current density satisfies the continuity equation ∇ ⋅ J(r, t) + t (r, t) = 0 which should be met independently of the choice of axes of the Cartesian system.
Calculating the divergence in the Cartesian system defined by the axes ŷ and ẑ , where the axis ŷ is parallel to the interface, will give a finite number if the calculation is per- formed separately in the half-space above the interface and the half-space below the interface.The problem arises if we want to determine the divergence of the current density ∇ � ⋅ J � (r � , t) in a Cartesian system with axes ŷ′ and ẑ′ which is rotated as compared to the original one.For this purpose, however, it is necessary to unambiguously determine the value of the current density J � (r � , t) in the rotated system at the common interface of a bilayer, otherwise inconsistencies cannot be avoided.Moreover, in accordance with thermodynamic principles, it can be assumed that a metallic bilayer in the state of thermodynamic equilibrium exhibites an equilibrium value of the free electron gas density at the mutual interface of the metals and a smooth functional course even in the vicinity of the interface.The continuous evolution of the density necessarily corresponds to a continuous evolution of the values of the effective mass m * of electrons.Classically, the current density is defined as J(r, t) = env(r, t) , where n is the electron density, e is the electron charge and v(r, t) is the electron velocity.From the classical point of view, the density element has the momentum p(r, t) = (m * ∕e)J(r, t) .Therefore, the assumption of the discontinuity J(r, t) at the interface implies automatically the discontinuity of the momentum field, which contradicts the law of conservation of a momentum.
The only way to avoid this problem is to introduce additional boundary conditions defining the continuity of normal and tangent components of current densities at the common interface of metal layers.The mentioned boundary conditions for the current density are actually additional relations which make it possible to determine the magnitude of the normal current density.
Next, we will focus on the specific determination of boundary conditions for current density, the continuity of normal components can be defined using Eqs.( 23) and ( 24).The upper metal layer has a nonzero normal current density at the lower boundary, so Eq. ( 24) with the corresponding parameters with index 1 can be applied.In contrast, the lower metal layer has a nonzero value of J u only, so we use Eq. ( 23) with the corresponding index 2. Introducing this index 2, it will be appropriately applied for all relevant quantities within this model (electric fields, current densities, and transfer matrices).Thus we obtain the boundary condition of continuity of normal components in the form of: Page 15 of 33 1286 The values of the sine functions for the boundary condition are obtained by a simple modification of Eq. ( 9), specifically in the form: The amplitude of tangential current density at the common interface of the metal layers can be derived by substituting the transverse and longitudinal components of E 1,2 from Eq. ( 19) into ( 22).The same value must be obtained by substituting the transverse and longitudinal components of E 1,2 , but determined by a relation analogous to (17) for the second layer, into (22).By subsequent comparison, we get: The overall model of a metal bilayer thus consists of two sets of equations which separately describe the first and second metal layers.These two sets of equations are bound together by boundary conditions of continuity of current densities at a common interface of the bilayer ( 47) and ( 49).The first mentioned set of equations consists of a system of equations ( 16) to (20).The second set is obtained simply by shifting the indices of all quantities in the first set, but with the exception of the vacuum impedance quantity 0 .The considered sets of equations together with the boundary condition for the current density already allow us to determine the normal component of the current density at the interface as linear combination of tangential amplitudes H 0,1 , H 2,3 , E 1,2 , and H 1,2 .Deriving this important relation is quite complicated but straightforward, and therefore we will only present here the following result: Formula ( 50) is of course also valid for J 2u .Due to the complexity of the relation itself, the substitutions below were used: (47) (50) 1286 Page 16 of 33 In order to further facilitate and clarify the derivation of the transfer matrix of the bilayer itself, it is necessary to rewrite the relation (50) via substitutions of complex expressions into a simplified form: Since we already know the relation for the normal component of the current density J 1b and J 2u , we can substitute it into (46) and thus determine the individual modification vector for the upper layer.The upper metal layer is surrounded on top by a dielectric medium, so this layer has a nonzero component J 1b only.For the modification vector of the first layer, we obtain the relation: By straightforward mathematical manipulations in (59), we obtained defining relations for matrices N 1b , N 2b , and N 3b .In the case of the second layer, we use the appropriate equivalent of ( 46).Unlike the first layer, the second layer is surrounded at the bottom by a dielectric medium and thus has a nonzero component J 2u only.In an analogous way, we obtain similar matrices N 1u , N 2u , and N 3u for the second layer: We can now return to Eq. ( 42) or an equivalent of this equation for the lower metal layer.Since the values of the tangential fields at the upper interface of the layer can be (54) (58) . determined using the transfer matrix from the knowledge of the values of the tangential fields at the lower interface, it is necessary to first deal with the lower layer.By the above procedure and after a simple modification, we obtain the relation ( 61), where the expression I + N 2u −1 has the meaning of the inverse matrix to the matrix I + N 2u : We can now apply the same procedure as above in the case of the first layer, substituting (59) into (42), we get then the relation: Now it remains to substitute relation (61) into Eq.( 62).This gives the relation between the tangential fields at the lower and upper interfaces of the bilayer.After simple mathematical adjustments, the required expression for the transfer matrix of the bilayer can be obtained.
In order to simplify the desired relation, it is appropriate to introduce the substitution: It is now possible to establish expression for the transfer matrix M bil of a metallic bilayer in the following simplified form: Hence, we are now able to determine the transfer matrix for the two adjacent metal layers.
Clearly, if one wants to calculate the transfer matrix or individual characteristics or both for a metal bilayer made of the same metal in this way, one should always obtain the same results as if one used the relations for transfer matrix of one metal layer with the appropriate thickness ].The following identities of the transfer matrices and reflection coefficients must therefore apply: In Eq. ( 65), M one represents the transfer matrix of one metal layer on both sides surrounded by dielectric materials.The selected metal layer of thickness d can be divided into two metal layers of thicknesses d 1 and d 2 in any way, i.e. d = d 1 + d 2 .This criterion can also be seen as one way to verify the correctness of the relation (64).Due to the complexity of expression (64), we verified the required criterion (66) numerically, and found that the deviation from identity (66) was practically negligible, at the level of numerical rounding errors.
In Fig. 5, we present a comparison of the reflectance curves as a function of wavelength.The green curve belongs to the classical local model, the blue and red curve to the nonlocal calculation of reflectance of 8 nm thick gold layer considered as a gold bilayer with a thickness of the upper and lower layer 2.5 and 5.5 nm , respectively.The blue curve repre- sents the correct way of calculations by the nonlocal transfer matrix method for bilayers 1286 Page 18 of 33 presented here by relations (64).The red curve is the result for case of the incorrect (hard wall) boundary conditions for the current density.The above comparison (red and blue curves) clearly shows the importance of using the correct boundary conditions.

Plasmonic single layers: results and discussion
In this part, we will firstly focus on the comparison of calculations according to classical local theory and according to nonlocal theory based on the HD model.Given that the strongest nonlocal deviation from the classical local theory can be expected for the cases of grazing angles of incidence of a TM wave on a metal layer, the experimental arrangement must correspond to this fact.An optical prism must be used to prevent total reflection when the plane wave passes from the air to the dielectric with a larger refractive index.To avoid undesirable reflections from the dielectric interface with the external environment (air), the edge of the optical prism should be at right angles to the reflected beam.Furthermore, undesired reflections of the reflected wave at the edge of the optical prism should be avoided.For this reason, the edge of the optical prism should be at a perpendicular angle to the reflected beam.Similarly, it is possible to grind the surface of the lower dielectric so that the transmitted wave propagates perpendicularly to the interface of the dielectric and the external environment.Our theoretical calculations of nonlocal responses could therefore be confronted in the future with a hypothetical experimental measurement setup, according to Fig. 6.Now we will return to the calculations of the nonlocal deviations from the local model.For this purpose, it is ideal to choose an relatively easily measurable energy quantity-the reflectance defined as R = rr * .Our goal is to find, within certain limits, the range of the main parameters of our model, where the largest deviations in the calculations of reflectance between classical local and nonlocal theory will occur.Since this is a multi-parametric problem, the question of setting up the simulations appropriately is not entirely trivial.Therefore, we proceeded as follows.Based on many calculations, we gradually (by the iterative way) determined the values of some parameters and subsequently used them for further calculations.After several repetitions of this procedure, we arrived at the optimal values of the parameters of the calculations presented below.In this paper, we focused on the simplest cases, namely the case shown in Fig. 7.As indicated in Fig. 7, the key parameters that influence the reflectance are the angle of incidence 0 measured from the perpendicular direction, the refractive indices n 0 , n 2 of the environments surrounding the metal layer, and the thickness d of the metal layer itself.Our aim was to perform the calculations for realistic material values of refractive indices and for gold or silver or both as metals commonly used in plasmonics.For this reason, we have chosen a range of refractive index values between 1 and 2 for the calculations, in this range there is a number of materials which imaginary component of the refractive index can be neglected.
Figure 8a, b represent an illustrative guide to the issue of choosing the values of the refractive indices of the external environment.These calculations show that the n 0 and n 2 cannot be too different for the nonlocal interaction to be strong.Therefore, for many calculations we chose n 0 = 1.7 and n 2 = 1.6 as a compromise value for both analyzed met- als.Manifestations of nonlocal effects occur mainly for nanostructures where at least one dimension is in the order of a few nanometers.For our calculations, we chose d = 8 nm as the acceptable thickness so that the nonlocal interactions were strong enough, and at the same time, we could neglect quantum phenomena.The silver and gold metal layers have, at the same parameter values ( n 0 , n 2 , d, 0 ), the maximum nonlocal deviation at different wavelengths, see Fig. 9a, b, therefore it is necessary to define appropriate wavelength for both cases of metals specifically.With this, we retrospectively justified the determination of the wavelength for the calculations in Fig. 8a, b.For the analysis of the gold layer, we chose = 532 nm which corresponds to the second harmonic frequency of the Nd:YAG laser.In the case of the silver layer, we chose = 450 nm , i.e. radiation that emits, for example, blue laser diode.As Fig. 10a, b indicate, nonlocal effects are strongest in the case where the plane wave propagates to the metal interface at a grazing angle.Therefore, we set the value of the plane wave incident angle 0 to 0.95 ⋅ ∕2 = 85.5 • , as a compromise value for both metals.
The first simulations presented here in Fig. 8a, b show the magnitude of the nonlocal deviation of reflectance as a function of the refractive index values n 0 of the medium above (superstrate) and below the metal layer n 2 (substrate).Both Fig. 8a, b indicate a similar character of the increasing intensity of the nonlocal deviation with increasing refractive indices n 0 and n 2 , where n 0 is slightly larger than n 2 .From large number of our analyses, we have observed that the nature of this dependence remains very similar even when the wavelength shift is in the range of ±50 nm around chosen wavelengths (532 nm and 450 nm in the case of gold and silver layer, respectively).
One of the other possibilities of the analysis is to display the magnitude of the reflectance deviation as a function depending on the layer thickness and the wavelength of the incident plane wave.Figure 9a for gold and Fig. 9b for silver layer show this case.As in previous calculations, we set 0 = 85.5 • .For both calculations, we set the refractive indices to n 0 = 1.7 and n 2 = 1.6 which corresponds to the state with a significant nonlocal response according to both Fig. 8a, b.In both Fig. 9a, b, the range of parameters at which the reflectance deviates significantly from the classical predictions is clearly evident.In the case of a gold layer, i.e. Figure 9a, the largest deviation occurs at a layer thickness of d = 9.5 nm and at the wavelength of 523 nm .Position of the maximum of the nonlocal response, see Fig. 9b in the case of a silver layer, is determined by the parameters d = 3 nm and = 339 nm .The silver layer shows a stronger nonlocal response, even for the 8 nm thick layer.The deviation at this thickness and wavelength of 409 nm reaches a value of 0.11.In both cases of metal layers, the most significant deviation shifts with increasing layer thickness to longer wavelengths.
From both Fig. 9a, b, it can be concluded that a non-negligibly strong nonlocal response can be expected even for layers about 40 nm thick.In the case of the gold layer with the mentioned thickness, the maximum deviation has the value 0.049 and is determined by the wavelength = 761 nm , for the silver layer the maximum is given with the following data: (R HD − R) = 0.061 and = 728 nm .Another parameter which depends on the strength of the nonlocal interaction is the angle of incidence of the plane wave on the metal layer.Figure 10a for the gold and Fig. 10b for the silver layer represent this possibility of analysis.The wavelength was selected as the second functional parameter, for this purpose.We left the values of the previously used parameters at the same values for a more objective comparison, namely n 0 = 1.7 , n 2 = 1.6 , d = 8 nm .Figure 10a, b clearly show the range of values of the angle of incidence at which the nonlocal interaction occurs.According to our calculations, the gold layer shows the strongest nonlocal behavior around the wavelength of 500 nm and for angles of incidence greater than 70 degrees in the case of the silver layer, it is around the wavelength of 400 nm and the angular range of 80 • to 90 • .As can be seen in Fig. 10a, in the vicinity of an incidence angle of 70 degrees, specifically at 0 = 70.247• , colour transition is visible.This striking line of rapid change in the reflectance deviation corresponds to the critical angle.At higher values of the incident angle, the incident TM plane wave is coupled to the plasmon-polariton mode, resulting in a significant decrease in the reflectance (Solis-Tinoco et al. 2021).Unlike the classical model, the nonlocal one predicts a lower (more realistic) decrease in the reflectance when the critical angle value is exceeded.This appears in Fig. 10a, the region with the highest nonlocal reflectance deviation.The colour transition line given by the critical angle value is also slightly visible in the case of the silver layer in Fig. 10b, again in the vicinity of the incidence angle value of 70 degrees.

Plasmonic bilayers: results and discussion
So far we have dealt with the single metal layer, now, we will move to the analysis of the nonlocal behavior of the metal bilayer.Clearly, one can compare the above calculations for a single metallic layer with their analogy regarding the metal bilayer.As before, we will rely on the above presented theory, namely the relations for the transfer matrix of the metal bilayer ], together with the formula for the reflectivity (32).In order to build on the previous analysis of one metal layer, we will keep the same values of refractive indices of external environments and also focus our attention only on the metals as gold or silver or both.Clearly, the tuning of the bilayer thickness parameter will be considered.In accordance with the previous analysis of the influence of the single layer thickness on the strength of the nonlocal behavior, it seems most appropriate to keep the total thickness of the bilayer at 8 nm .Before we start commenting on the selected calculations for the metal bilayer, it is necessary to mention the phenomenon related to the attenuation correction that we have implemented.
During the validation of our model, we were interested in a specific situation, where the reflectance takes the maximum value ( R = 1 ).At the same time, we mainly focused on the case of a significant nonlocal deviation.As was previously shown, a noticeable nonlocal response is predicted for such settings, where the refractive indices of the external medium are close (see Fig. 8a, b).From our analyses, it can be concluded that this condition is broadly valid for different wavelengths.
According to classical calculations, both the gold and silver single layers, even the 8 nm thin, show a reflectance close to 1 for the parameters of the refractive indices and the angle of incidence of the plane wave chosen as n 0 = 2.1 , n 3 = 2 , 0 = 81 • .If we use this param- eter setting to calculate the reflectance using this nonlocal model for the metal bilayer, we find that according to our model, the reflectance is higher than one for some cases.This is illustrated in Fig. 13a, where d 1 denotes the thickness of the upper gold layer, d 2 denotes the thickness of the lower silver layer.At first glance, it may seem that our model must be wrong, because it predicts non-physical results (approximate range of parameters: ).However, we came to the conclusion that these non-physical results are caused by the inadequate setting of the electron gas attenuation in both metal layers.The next Fig.11 shows a hypothetical notion of electron gas attenuation.The incident plane wave initiates the motion of the electron gas in a metal.However, randomly moving electrons cause a damping of this excited motion of the current density.This process is schematically illustrated in Fig. 11, where the incident wave-induced current density is illustrated by a series of electrons in the form of green spheres connected by a red zigzag line.This red zigzag line shows the force coupling between the regularly moving electrons of the current density wave.The electric force bond between the individual electrons of the current density can be compared to a kind of spring.If the distance between two electrons is less than the equilibrium distance, then they are pushed apart.If, on the other hand, these electrons are located at a greater distance than their mutual equilibrium distance, the total system of electrons is also out of equilibrium which results in the pressure of the outer electrons pushing the two electrons together.This process of electron gas interaction means that the movement of individual electrons cannot be considered as isolated objects and it is relevant to consider the dynamics of the electron gas or rather the hydrodynamics of the electron fluid.
Electromagnetic wave propagation in a metal is inextricably linked with electron density waves, and if one wants to avoid inconsistencies leading to infinite density values or infinitely large time changes, as already explained above, this electron density wave must be continuous at the common boundary of both metals.Therefore, we can immediately rule out the case that waves of electron densities from individual metal layers do not interact with each other.In other words, the movement of the electron density at the boundary of two metal layers and in its vicinity must necessarily respond to the movement of electron densities which exhibit different behavior, i.e. have different densities, attenuation, and other parameters.Therefore, the electron current density near the boundary of the two metals cannot be attributed to the attenuation value of either the first or the second metal, but some average value must be taken.Since we have focused in this article only on the two metallic materials of gold and silver, see Table 1, they have almost identical values of the main parameters of the HD model, except for the attenuation.For this reason, we will focus on the idea of averaging only the attenuation constant, although in the case of another pair of metals, we can also think about averaging the values of other parameters.From a classical point of view, the orderly movement of current density electrons is disturbed by other processes, so that the electron liquid tends to thermodynamic equilibrium (Jakoby 2009).This manifests itself in the attenuation of the ordered motion of the electron liquid.According to Drude (Yang et al. 2015) and other classical models of metal permittivity, the attenuation constant is defined as a quantity inversely proportional to the relaxation time, i.e. in our simplified case as the collision frequency = 1∕t c , where t c is the average time interval between two electron inelastic collisions.Suppose one can follow a thought construct, see Fig. 11.A part of the electrons of the electron density is subject to an ordered motion, caused by a propagating plane wave in the metal, and is damped by collisions of randomly moving electrons.The electrons, that are currently part of the electron gas wave due to force bonds, are linked together and form a hypothetical coherent chain of electrons.Such a hypothetical chain of electrons passing through the boundary of two different metal layers is exposed to a different attenuation in these two metals.Due to a certain coherence and non-zero inertia of the individual electrons of the mentioned chain, energy is supplied from electrons with lower attenuation from one metal to electrons with higher attenuation in the second metal.The described phenomenon applies analogously in the opposite direction.The electrons involved in the orderly motion in the metal with a lower attenuation thus lose their kinetic energy faster and show greater attenuation.Conversely, electrons in a higher attenuation metal lose their energy more slowly and exhibit lower attenuation due to inertial forces from electrons in a lower attenuation metal.This process should lead to the averaging of the attenuation constants of both metals near their common boundary.For further considerations, it is convenient to focus the attention on a single electron.An important argument for introducing an average (effective) value of the attenuation constant is the following.Commonly stated values of either the relaxation time or the electron free path (Mortensen 2021;Ketterson 2016) (up to 100 nm for pure noble metals) mean that the electron can with great probability freely overcome a path greater than the thickness of at least one layer of the metal bilayer, which we analyze below.Therefore, it is reasonable to relate the collision frequency (attenuation constant) at least to the total thickness of the bilayer.Let us assume that the trajectory of the electron, at a distance smaller than the free path of the electron, forms an angle e with the ŷ axis.Due to the boundary conditions of continuity for the current density, this angle is the same in both layers of the metallic bilayer.The distance traveled by the electron in the first metal layer is therefore 1 + sin 2 e and s 2 = d 2 √ 1 + sin 2 e in the second layer, according to Fig. 12.If the velocity of the electron is v and the collision frequency of the electrons in the first and second metal is 1 = 1∕t c1 and 2 = 1∕t c2 , then the number of collisions related to the con- sidered electron is: p = (s 1 ∕v)∕t c1 + (s 2 ∕v)∕t c2 .The average collision frequency (attenuation constant) of electrons in the bilayer is therefore av = p∕t m , where t m = (s 1 + s 2 )∕v is the time of electron movement.Sufficient for our purposes, refined attenuation constant to the first order is: The modification of the attenuation constant of the HD model, of course, has an effect on the wave numbers of the longitudinal and transverse fields in both metal layers.If we denote by the index i the metal layer (for upper i = 1 , lower i = 2 ), we can then derive the expression of the transverse and longitudinal wave number as: As in (3), the term b i ( ) means here the permittivity of a positively charged background and m i is the permittivity of i-th metal.The quantities plasma frequency p i and attenuation constant i represent the original unchanged HDM parameters for a given metal.It only remains to determine the permittivity of the positively charged background of both metals which can be simply derived in the form: Thus, it is obvious that, in contrast to the classical local theory and the HD model for the IMI structure, the transverse wave numbers are modified, i.e. k T i ≠ ( ∕c) √ m i ( ) .In addition, the wave numbers are now functionally dependent on the thicknesses of the metal layers.
After this explanation, one can focus on the nonlocal response of a metal bilayer.The structure of the metal bilayer and the main parameters on which the interaction depends are indicated in Fig. 12.In order to assess the behavior of the metal bilayer, we chose the same reflectance deviation calculations as in the case of one metal layer.
We can now focus on comparing the calculation of the deviation (R HD − 1) , i.e. the anal- ogy of the calculations shown in Fig. 13a, if the averaging of the attenuation constant is used.As can be seen in Fig. 13b, the value of the reflectance deviation is negative for all the thicknesses of the upper gold and lower silver layer of this bilayer.We have verified the assumption that the reflectance does not exceed the maximum value of 1 in a wider range of wavelengths, and it clearly shows that our proposed averaging technique of the attenuation constant is a correct method leading to realistic results.

Nonlocal symmetric bilayers
To illustrate the influence of the nonlocal response of a metal bilayer, we present here several figures of calculations for the case of a bilayer which consists of gold and silver layer, (67 .  of refractive index values, manifested by increased reflectance, occupies roughly the same area in both Fig. 8a, b for both the silver and gold layers at a given angle of incidence, regardless of the wavelength.In the case of a metal bilayer, see Fig. 14a, the range of refractive index values with the expected increased reflectance is slightly shifted to higher values of the refractive index n 0 of the upper environment, but n 3 > n 0 still applies.
Another of the nonlocal response characteristics discussed above is the dependence of the reflectance deviation on the angle of incidence and the wavelength of the plane wave.The graphical representation of the mentioned reflectance deviation for the gold-silver bilayer is presented in Figs.15a, b, the calculation was again chosen for the symmetric case of 4 nm thick both gold and silver layers.As we mentioned, our analyses of results show that in the case of very thin bilayers, the order of the individual layers has very little effect on the reflectance deviations examined here.An example is a silver-gold bilayer composed of 4 nm thick layers, the nonlocal response of this bilayer as a deviation of the nonlocal reflectance as compared to the local one, with respect to the angle of incidence and wavelength is shown in Fig. 15a.If we compare Fig. 15a with the analogous result of the reflectivities in Fig. 15b, for the gold-silver bilayer, i.e. the case of the opposite order of layers, we find that the differences between the two responses are very small.The influence of the gold layer on the overall nonlocal response of the analyzed gold-silver bilayer is manifested to a large extent in the form of a pattern defining the maximum values of the deviation, this is evident from the comparison of Figs.15b and 10a.On the other hand, the silver layer contributes to the total nonlocal deviation of the bilayer in a different way.The pattern in the Fig. 15b, which is formed by the most pronounced values of the reflectance deviation in the range of wavelengths 330 nm to 400 nm , can be attributed to the influence of the silver layer.Compared to its counterpart in Fig. 10b, the pattern already has a more significantly changed form, smaller range of values, and shifted position to shorter wavelengths.Figure 16 shows the nonlocal deviation of the reflectance of the silvergold bilayer which is composed of a silver and a gold layer, both of the same thickness.
The situation of the asymmetric case, when one of the metal layers of the metal bilayer is thicker than the second one, is briefly discussed in the following Sect.7.2 of this paper.

Nonlocal asymmetric bilayers
In order to follow the analyses of the reflectivities presented in the previous Sect.7.1, we will focus mainly on the case of 8 nm thick bilayers, but composed of layers with thick- nesses of 1 nm and 7 nm .The other two Fig. 17a, b again show the evolution of the reflec- tance deviation as a function of the angle of incidence and the wavelength of the plane wave.While Fig. 17a relates to a bilayer formed of an upper gold layer 7 nm thick and a lower silver layer 1 nm thick, Fig. 17b again relates to a gold-silver bilayer but with thicknesses of 1 nm and 7 nm , respectively.Both of the above-mentioned Fig. 17a, b show the response of the thicker layer in a given bilayer, but at specific wavelengths a very strong nonlocal interaction of the thinner layer prevails.The other two Fig. 18a, b show the dependence of the reflectance deviation of the gold-silver bilayer on the values of the refractive index n 0 of the medium above the gold layer and similarly n 3 below the silver layer.In the case of Fig. 18a, the parameters were set to these values, i.e. the wavelength at 532 nm , angle of incidence 0 at 85.5 • , thickness of gold layer at 7 nm and silver layer at 1 nm .The above-mentioned Fig. 18b shows an unexpectedly strong response, which concerns the interaction of a plane wave with the parameters = 480 nm , 0 = 85.5 • with a bilayer of a gold and a silver layer 1 nm and 7 nm thick, respectively.It is now proposed to compare the nonlocal responses of a single gold layer (in previous Sect.6) Fig. 8a, and a gold-silver bilayer Fig. 18a of the same thickness under the same conditions.As can be seen from the comparison, the single layer shows a slightly stronger nonlocal response, but not fundamentally.On the other hand, as already indicated, under specific conditions, the nonlocal response can be very strong.Similarly, we can compare the response in Fig. 18c of a 8 nm thick silver layer interacting with a plane wave of = 480 nm and 0 = 85.5 • with response in Fig. 18b.Thus, a significant difference between Fig. 18b, c is due to the nonlocal behavior of the 1 nm thick gold layer.The maximum reflectance deviation, see Fig. 18b, reaches 0.242, so the nonlocal model predicts, under certain circumstances, 39.6% higher reflectance compared to the classical calculation ( R ≈ 0.6).
Finally, the following Fig. 19 shows the reflectance deviation as a function depending on the wavelength of the incident plane wave and the total thickness of the bilayer.This bilayer always consists of an upper gold layer with a thickness of 1 nm , the rest of the total thickness of the bilayer is occupied by a silver layer.The main parameters of the calculation belonging to Fig. 19 are the angle of incidence 0 = 85.5 • and the refractive index of the medium above and below the bilayer n 0 = 2 and n 3 = 1.9 , respectively.As can be seen from Fig. 19, our nonlocal model predicts a very significant reflectance deviation near the wavelength of 500 nm and the total bilayer thickness of 8 nm .Note that according to Fig. 19, the nonlocal model at the parameters = 495 nm and d = 8 nm indicates a reflec- tance value of 0.632.In other words, this means that compared to the classical reflectance, which is approximately 0.33, the nonlocal model gives a 91% higher value.
As can be seen from the reflection deviation graphs, the total nonlocal response of a metal bilayer is a combination of the responses of the individual layers, but it is not just a simple average of these individual responses.Therefore, for a more detailed description of the nonlocal processes, especially for two nonlocal materials interconnected by the flow of current densities, it is necessary to apply an appropriate model that more accurately captures the physical reality, exactly as was presented here.

Conclusion
The aim of this paper was to present an analytical nonlocal model for planar plasmonic layers which could analyze the nonlocal response of both an isolated metal layer and two connected metal layers.As was shown, our model predicts a relatively strong nonlocal response for the specific cases of TM plane wave interaction with both IMI and IMMI structures.For the case of plasmonic bilayers, since our model does not neglect the effect of the flow of current densities from one metal layer to another, it brings much more reliable data for such bilayers, to our knowledge for the first time in the literature.Specifically, we discussed both symmetric and asymmetric bilayers with respect to their thickness.Other preferences of the model are its simple implementation and its modifiability taking into account quantum phenomena as the corrections at the interface, such as the quantum tunneling or the spill out effect.

Fig. 2
Fig. 2 Schematic representation of the MIM structure chosen to compare nonlocal models

Fig. 3
Fig. 3 Reflectivity of the structure shown in Fig. 2 as a function of the incident angle i for two values of d 2 .The wavelength of the incident TM plane wave is = 600 nm .The medium above the top gold layer is TiO 2 with the refraction index n 0 = √ 6.78 and the medium below is air.The permittivity of gold has value m = −8.44 + 1.41i and nonlocal parameter = 1.35 × 10 6 m s −1 .The upper metallic layer thickness is d 1 = 18 nm .a, c are taken from Pitelet et al. (2018)

Fig. 4
Fig.4Scheme of the nonlocal interaction of metal bilayer and incident TM plane wave

Fig. 5
Fig.5Comparison of reflectance as a function of wavelength for a gold layer with a thickness of 8 nm .Refractive indices of the environment above and below the layer have the value n 0 = 1.7 , n 2 = 1.6 .The angle of incidence 0 of the TM plane wave is 85.5 •

Fig. 6 Fig. 7
Fig. 6 Possible experimental setup for measuring the nonlocal reflectance deviation

Fig. 8
Fig. 8 Dependence of nonlocal deviation of reflectance of 8 nm thick metal layer on refractive indices of external environments n 0 and n 2 .The angle of incidence of the TM plane wave is 85.5 • , wavelengths are selected in the case of a gold layer: 532 nm and b silver layer: 450 nm

Fig. 10
Fig. 10 Dependence of nonlocal deviation of reflectance of 8 nm thick a gold and b silver layer on wavelength and angle of incidence of the TM plane.The refractive indices of external environments are set to n 0 = 1.7 and n 2 = 1.6

Fig. 11
Fig. 11 Schematic representation of electron gas attenuation at the interface of two different metals.The interface of the two metal layers is indicated by a dashed line.The green circles represent the electrons, the blue arrows the direction of their random movement, and the red zigzag line indicates the force bond between electrons of the induced current density driven by the EM field of the interacting incident plane wave.The purple arrows indicate the direction and magnitude of the damping force which is directed against the direction of movement of the current density.(Color figure online)

Fig. 12
Fig. 12Illustration of the bilayer structure and main parameters of interaction.The TM plane wave propagates through a dielectric medium with refractive index n 0 and impacts on the metal bilayer at the angle of incidence 0 .The upper layer of the metal bilayer has a thickness of d 1 , the second lower one of d 2 .The medium under the metal bilayer-substrate-is dielectric with refractive index n 3

1286
Page 26 of 33 both of the same thickness.(The reader will find additional detailed simulations for some other types of asymmetric bilayers in the Sect.7.2).Figures 14a and 15b are analogous to Figs. 8a and 10a, for a single metallic layer, the simulations were performed with the same settings of the main parameters, i.e. refractive indices, wavelength, and angle of incidence of a plane wave.As can already be seen from Figs. 8a, b relating to a single metal layer, there is a certain range of values of the refractive indices of the external environment in which the nonlocal response causes a significantly increased reflectance.This range

Fig. 13
Fig. 13 The of (R HD − 1) the HD model with a modified attenuation constant as a function of the thicknesses of the individual metal layers d 1 and d 2 of the gold-silver bilayer.The original setting of the attenuation constants gives inaccurate results in the region around d 1 = 20 nm and d 2 = 7 nm.The angle of incidence and wavelength of the TM plane wave are 81 • and 900 nm , respectively.The refractive indices of external environments are set to n 0 = 2.1 and n 3 = 2

Fig. 15
Fig. 15 Dependence of the nonlocal deviation of the reflectivity on the angle of incidence 0 and the wavelength .Case a represents the response of a silver-gold bilayer, case b and its 3D representation c a goldsilver bilayer.The bilayer has total thickness of 8 nm and is composed of a gold and silver layer of the same thickness 4 nm .Refractive indices of external environments have values of n 0 = 1.7 and n 3 = 1.6 .A modified HD model with an averaged value of the attenuation constant was used

Fig. 16
Fig. 16 Dependence of nonlocal deviation of reflectance on total thickness of bilayer d and wavelength .The angle of incidence of TM wave is 85 • .The refractive indices of external environments have values n 0 = 1.7 , n 3 = 1.6 .The HD model with a modifying attenuation constant was used

Fig.
Fig. Dependence of nonlocal deviation of reflectance on refractive indices of external environments at the wavelength of a 532 nm , b, c 480 nm , and angle of incidence 85 • .The metal bilayer has a total thickness of 8 nm , gold layer in bilayer has thickness of a 7 nm and b 1 nm .The thickness of a single silver layer c is also 8 nm .In the case of bilayers (a and b), the HD model with a modifying attenuation constant was used