Scattering Amplitude of Surface Plasmon Polariton Excited by a Finite Grating

Unusual optical properties of laser-ablated metal surfaces arise from the excitation of local plasmon resonances in nano- and microstructures produced by laser-processing and from the mutual interaction of those structures through surface plasmon polariton (SPP) waves. This interaction provides a synergistic effect, which can make the optical properties of the composite nanostructure drastically different from the properties of its elements. At the same time, the prediction and analysis of these properties are hampered by the complexity of the analytical solution to the problem of SPP excitation by surface objects of arbitrary configuration. Such a problem can be reduced to a simpler one if one considers the geometry of a structured surface as a superposition of harmonic Fourier components. Therefore, the analytical solution to the problem of surface plasmon polariton excitation through the scattering of light by a sinusoidally perturbed plasmonic metal/vacuum boundary becomes very important. In this work, we show that this problem can be solved using a well-known method for calculating guided-mode amplitudes in the presence of current sources, which is used widely in the waveguide theory. The calculations are carried out for the simplest 2D cases of (1) a sinusoidal current of finite length and (2) a finite-length sinusoidal corrugation on a plasmonic metal surface illuminated by a normally incident plane wave. The analytical solution is compared with the results of numerical simulations. It is shown that, in the first case, the analytical and numerical solutions agree almost perfectly. In the second case, the analytical solution correctly predicts the optimum height of the corrugation xopt, providing the maximum SPP excitation efficiency. At the same time, the analytical and numerical values of the SPP amplitude agree very well when the corrugation height x turns out to be x≪xopt or x≫xopt (at least up to 3xopt); at x=xopt, the mismatch of those does not exceed 25%. The limitations of the analytical model leading to such a mismatch are discussed. We believe that the presented approach is useful for modeling various phenomena associated with SPP excitation.


Introduction
Nanomaterials require a specialized toolkit to construct on-chip electronics and develop highly efficient devices capable of concentrating light at the nanoscale. To overcome the limitations of traditional dielectric waveguide nanotechnology, one promising approach involves the utilization of metallic structures that support plasmon modes. The field of plasmonics offers a compelling solution by harnessing the collective oscillations of electrons at the surface of these metallic nanostructures. This allows for intense light-matter interactions and enables the precise manipulation and control of light [1].
Nowadays, the field of plasmonics is a rapidly growing area of research that has gained significant attention in recent years due to its potential applications in various fields, including biosensing [2][3][4], fiber sensing [5,6], imaging [7][8][9], data storage [10], and energy conversion [11,12]. The field of plasmonics focuses on the interaction between electromagnetic radiation and free electrons in metal structures, known as surface plasmon polaritons (SPPs), which can lead to unusual optical properties and strong light-matter interaction.
The history of plasmonics dates back to the early 20th century, when SPPs were first discovered as anomalies in the reflectance spectra [13]. However, it was not until the 1950s that researchers started to study the interaction between light and metal surfaces and realized the potential of surface plasmons for enhancing the absorption and scattering of light. Since then, the field of plasmonics has evolved rapidly, with numerous studies exploring the properties and applications of SPPs in various structures, including metallic nanoparticles [14], thin films [15], and nanostructures created by laser ablation [16,17].
Optical phenomena involving SPPs can occur at subwavelength scales and are not limited by the diffraction limit. Therefore, one of the most promising applications of plasmonics is in the field of nanophotonics, which involves the manipulation of light at the nanoscale. This field has been driven by the need to develop smaller and more efficient optical components for a wide range of applications, from telecommunications [18] to medical imaging and drug delivery [19]. Plasmonic coupling between nanoantennas over a conductive surface is the basis for the design of highly sensitive sensors [20] and thermoplasmonic devices [21].
However, the design and optimization of plasmonic nanostructures for specific applications can be challenging, as it requires a deep understanding of the complex interactions between light and matter at the nanoscale. The excitation of SPPs in these structures is particularly important, as it can enhance the local electromagnetic field and lead to increased light-matter interactions. Therefore, there is a need for accurate and efficient computational tools for modeling the behavior of SPPs in plasmonic nanostructures, which can aid in the design and optimization of these structures for various applications.
The study and analysis of nanostructured plasmonic surfaces raise even broader questions. In addition to the excitation of SPPs in individual structures, the mutual interactions between neighboring structures can also lead to novel optical effects. This interaction arises from the coupling of SPPs between adjacent structures, which can result in the formation of new modes with different spectral and spatial properties. The resulting composite structures exhibit a synergistic effect that is not present in the individual structures, leading to significant changes in the optical properties.
In this work, we extend the approach initially suggested by us in [22,23] for a description of SPP excitation by confined nanoantennas to the case of SPP excitation by nanostructured substrates. While modern manufacturing techniques, including laser ablation, allow the creation of arbitrarily shaped profiles, we investigate the simplest case of a sinusoidally perturbed surface because the general scattering problem can be decomposed into the set of Fourier components. We pay extra attention to properly accounting for SPP radiation losses and the impact of the finite length of the profiled area. All the analytical results are confirmed by the numerical calculations.

Excitation of SPP within Born Approximation
First, we are going to determine the efficiency of the induced SPP wave based on the Lorentz reciprocity theorem [24]. Suppose we have a finite-length grating on the metal-air interface under the plane wave illumination (Figure 1a). To calculate the amplitude of the excited plasmon polaritons, we use the expression following the unconjugated form of the Lorentz reciprocity theorem [24] where e SPP = e SPP x n x + e SPP z n z is the SPP mode's electric field with e SPP z (x) = e 0 e −γ(x)|x| and e SPP ; n x , n y are the unit vectors in the X and Y directions, respectively; J is the excitation current density distribution; and k 0 = 2π/λ is the vacuum wavenumber.
Let us assume that this grating represents a corrugation with periodically varying height x g = x 0 cos k g z and total length z 0 (Figure 1a), and the height of the corrugation x 0 is small compared to its period Λ = 2π/k g . Suppose that a plane electromagnetic wave with a wavelength λ not too different from the period Λ falls normally on this corrugation from the vacuum (see Figure 1a). The density of the equivalent current, due to the presence of the corrugation, can be described as in [24]: where E is the electric field amplitude, the dependencies ε(x, z) and ε s (x) set the modulation of the dielectric permittivity of the layer of the considered plasmonic metal with corrugated and initial flat surfaces, respectively. The field E is a superposition of the incident wave and that reflected from the metal surface with some additional contribution associated with light scattering.
While the E itself depends on the SPP field and the latter is described by the a SPP amplitude, (1) defines a self-consistent problem. In order to simplify the problem, we consider using the first Born approximation [25]. For this purpose, we will ignore all the scattering terms in (2), and assume that where E 0_ins is the incident wave field amplitude, r = 1 − √ ε Me / 1 + √ ε Me the reflection coefficient from the flat metal surface, and γ 0 (1) can be evaluated analytically, provided that x 0 λ by a power decomposition of a small parameter k 0 x 0 . After taking into account the linear terms only, one can obtain the following expression for the Z component of the electric field right after the end of the grating: where One can find a derivation of (4) for the case when the total length of the grating possesses the form z 0 = n + 1 2 Λ (where n is an integer) from the Lorentz reciprocity theorem in Appendix A.
To verify the analytical results obtained, we performed numerical simulations using the commercial software COMSOL Multiphysics. In this subsection, we create a twodimensional numerical model simulating half of a grating with a total length of (2n + 1/2)Λ (with the perfect electric conductor boundary condition to ensure mirror symmetry; n is considered an integer) surrounded by a perfectly matched layer. The solution to the Maxwell equation was obtained using the finite-element method in the frequency domain (FEFD). The SPP contribution was determined from Equation (1) with the calculated electric field distribution substituted into (2). Integration in Equation (1) was performed along a vertical line perpendicular to the interface, positioned after the sinusoidal perturbation. We chose the simulation area size and mesh size to ensure that the results did not change significantly under small variations of the model's parameters. Figure 2 compares the analytical calculations according to (4) and the results of fullwave numerical calculations for different grating periods. Hereinafter, we consider gold as a plasmonic metal (dielectric permittivity of gold is taken from [26]). The figure demonstrates that, overall, (4) agrees with the results of the simulations for relatively low gratings (namely, for grating heights up to 30 nm). As the length of the grating increases, the distinction of the dependencies grows stronger. Moreover, the degree of divergence depends on the grating period. When the difference between the resonant period and the actual one (red color in Figure 2; the period equals 778 nm) suppresses the SPP generation, the analytical expression describes the behavior of the numerical results even for longer grating lengths of up to 300 µm. However, while the grating period approaches the resonance, the divergence becomes more pronounced. One can explain it as follows. The Born approximation assumes that the background field largely determines the total electric field value. If this condition is met, the additional terms related to the "self-action" of the SPP field are negligible In the case of resonant excitation, the SPP field is compatible in strength with the background field; therefore, the approximation is less accurate.  Another factor limiting the accuracy of the described approach is the grating amplitude. Indeed, the background field decays as it penetrates the metal. By contrast, the total electric field is essentially non-zero near the grating, even for the larger values of x 0 . Thus, the background field is no longer the dominant contribution to the total field. Figure 3 supports this conclusion. Indeed, (4) adequately describes the numerical experiment for small grating heights (up to 5 nm). A further increase in the grating amplitude leads to a divergence of the numerical results from the analytically predicted ones. Importantly, the difference is retained even when the next orders of magnitude of the parameter k 0 x 0 are taken into account, despite the corresponding expressions becoming more complicated.  Therefore, there is a need for another approach, which will better describe the SPP generation for larger grating amplitudes, as well as the case of resonant SPP excitation.

Excitation of SPP by a Flat Harmonic Current
Now, we will explore a slightly different approach. Firstly, let us consider the problem of plasmon excitation by a harmonic current distributed on a flat metal-dielectric interface.
The current density possesses the form J surf = J surf · n z , where n z is a unit vector in the Z-axis direction, δ(x) is the Dirac delta function, k g = 2π Λ , I surf is the current amplitude, which implicitly includes the time dependence exp(iωt), Λ is the current modulation period along the Z-axis (Figure 1b). The current is not infinite in this direction, and its range is limited by the length z 0 .
Assuming in (1) that J = J surf , we note that of the two terms, exp(ik g z) and exp(−ik g z), which are in the Euler representation for the function cos k g z, defining the spatial modulation of the excitation current density in relation (5), only the former one makes a significant contribution to the amplitude of the forward-propagating (i.e., in the Z axis direction) SPP, while the contribution of the second one is negligibly small. The situation is the opposite for the back-propagating plasmon. Taking into account only the main contribution to the amplitude of the forward-propagating SPP, the Z component of its electric field and Me −1 . The results of calculations for the amplitude of this component, as carried out by expression (6) for various excitation current parameters when the plasmonic metal is gold, are shown in Figure 4; the current amplitude I surf is assumed constant within its range and equal to 1 A/m. Under these conditions, Figure 4a illustrates the dependence of the SPP amplitude E 0z on the spatial modulation period when the wavelength of SPP remains constant. The vacuum wavelength is set to be λ = 0.8 µm, which yields the corresponding plasmon wavelength λ SPP = 0.783 µm with its typical propagation length l SPP = 1/k I = 92 µm. Figure 4a demonstrates that when the values of λ SPP and Λ are inconsistent, the plasmon polariton amplitude oscillates as z 0 increases. These oscillations decay when z 0 l SPP , then the amplitude reaches the stationary value. If λ SPP = Λ, the stationary value is the largest and reached without oscillations.
The behavior of the curves can be interpreted as a result of excitation in the current zone of the "free" and "forced" SPP waves with the wave numbers k SPP and k g , respectively. These two waves have opposite-signed amplitudes, so at the beginning of the excitation current zone, at z = −z 0 (Figure 1b), the amplitude of the total right-propagating wave is zero. At the end of this zone, at z = 0, the total wave is the sum of the beats of the "free" and "forced" waves, the first of which gradually decays due to the non-zero imaginary part of the plasmon polariton propagation constant k SPP . It is possible to show that, in this representation, the total field is proportional to 1 − e −z 0 k I +i∆kz 0 , i.e., as it is assumed in expression (6). It gives the initial amplitude for an SPP wave that is free from the external action excited outside the current range at z > 0. If λ SPP = Λ, fluctuations in the amplitude of the excited SPP are consequences of the beats mentioned above. In resonance, there is no beating due to the equality of the propagation constants of "free" and "forced" SPP waves. Figure 4b shows the results of calculations of the dependence of the SPP amplitude on the current modulation period. The current amplitude and plasmon parameters are the same as in the previous example. The curves obtained for z 0 < l SPP have pronounced traces of "free" and "forced" waves' beats. As the length z 0 increases, the beats disappear, and the frequency response narrows and acquires the Lorentzian profile with a half-width equal to λ 2 SPP (πl SPP ) −1 . To illustrate the case when the current modulation period can go beyond the 770-800 µm range, as assumed in Figure 4a  To extract the initial amplitude of the excited SPP wave from the entire electromagnetic field generated by the source current, we used the fact that sufficiently far from the excitation current, the longitudinal electric field at the interface decays exponentially along the Z-axis with a damping constant corresponding to k SPP , which indicates the dominant contribution of the SPP wave to the total field. Considering this, the initial amplitude of the SPP wave was taken from the extrapolation from the distant point to z = 0 using the known exponential law of SPP attenuation. This extraction method provides the same results compared to the one used in the first subsection despite being based on a slightly different approach.
One can see that the numerical and analytical results agree well with each other. Therefore, we argue that the considered approach describes well the process of SPP generation by a surface current. Thus, the problem of plasmon polariton excitation due to the sinusoidal perturbation of the metal boundary (i.e., grating) reduces effectively to the appropriate choice of the current density. Of course, one can later introduce additional adjustments for better agreement with the experimental data since the current distribution may not be purely surface.

SPP Excitation by a Plane Wave Illumination
We are now ready to describe the SPP excitation with the approach developed in the previous section. We will start with the same corrugation as in Section 2.1, but introduce another approximation for E in (2), namely E = (1 + r)E 0_ins n z . Like in the previously investigated case, E 0_ins is the incident wave field amplitude, and r = 1 − √ ε Me / 1 + √ ε Me is the reflection coefficient from the flat metal surface. While the grating's amplitude is chosen to be much less than the vacuum wavelength, it suggests to be a more accurate approximation.
Note that the current density J g will be zero everywhere except in the corrugated regions highlighted in Figure 1a in orange and blue, where it is J g0 n z and −J g0 n z , respectively, and This gives the corresponding limits of integration over the variable x in (5), when substituted into it J g as the excitation current J. After performing the substitution, the expression for the SPP amplitude will include integrals of the form where the half-wave rectified cosine functions cos + x and cos − x are given by cos + x = cos x, cos x 0 0, cos x < 0 and cos − x = 0, cos x 0 cos x, cos x < 0 , so that the half-wave rectified cosine functions keep their sign for all x. The analytical calculation of these integrals becomes possible by decomposing exponents containing half-wave rectified cosine functions into a power series. If we limit the accuracy of such an expansion to the quadratic terms in x 0 , after integration, the amplitude Z component of the SPP field turns out to be where Let us add that we could use a less accurate way to account for the excitation current in the grating, which consists of first summing this current along the vertical axis (the X axis in Figure 1a) and then representing it in a purely surface form (5). Later, we will refer to this approach as the "flat current" approximation. It is not difficult to show that the amplitude of such a current will be (6) gives almost the same result as (7), except for the absence of a quadratic term in the parameter x 0 in (8). It is not difficult to understand that this term disappears in the flat current approximation because one does not consider the plasmonic mode's finite depth of penetration into metal and air when integrating in (6). The power expansion of the exponents in the transverse SPP mode accounts for this feature. However, the necessity of introducing the terms of the third or higher degree of x 0 to improve the accuracy of the SPP amplitude calculation, as compared to the flat current approximation, is still obscure. We will address these questions when comparing the analytical and numerical calculations.
So far, we have not discussed the impact of radiation losses introduced by the scattering of the SPP wave on the grating. We will account for them in the following considerations. Let us differentiate expression (7) with respect to z 0 . We obtain that dE 0z dz 0 = −γ E E 0_ins − (k I − i∆k) E 0z . The last expression is a differential equation for the plasmon amplitude, in which the energy dissipation is given by the terms k I E 0z . It is reasonable to assume that the radiative losses, as well as thermal losses, are proportional to the SPP amplitude. Adding the corresponding term transforms the considered differential equation to the following form: One should also make a corresponding substitution into (7): which, consequently, becomes an exact solution of Equation (10). Note that the introduction of radiation losses increases which is the typical grating length at which the SPP amplitude reaches saturation, which now turns out to be smaller than l SPP . Hereafter, we will refer to the case where the grating length is much greater than z sat as the long grating case. The value of z sat also specifies the half-width of the resonance contour formed by the dependence of the SPP intensity on the period of the long grating, according to the expression

Evaluation of γ rad
Despite the fact that (11) takes proper account of the impact of losses on the SPP amplitude, the value of γ rad remains unknown. Note that the differential Equation (11) turns out to be very similar to the one considered in [27], although the latter was obtained from different considerations. Using the approach proposed in the mentioned paper, let us add to the differential equation one more relation for the field intensity of the reflected light: where is the electric field intensity of the wave reflected directly from the flat metal surface, and is the electric field intensity of the wave arising as a result of plasmon energy leakage with κ as a coefficient characterizing this leakage. When analyzing the resulting system of Equations (10) and (14), one should ensure that the grating length is sufficient to consider the incident and reflected waves as plane waves. In this case, they may not be separated into separate spatial Fourier components, as was done in [27]. Otherwise, the approach to solving the system of equations remains the same as in the above work. Using the transformations based on energy conservation, time-reversal, and geometry mirror symmetry [27], and the results for the lossless case (k I = 0), we find that the parameters γ E , γ rad , and κ are connected by the following relation: where Me . This relation differs from the similar one given in [27] mainly due to the presence of the normalization factor N 1 , which appears due to the change in the approach to the description of the spatial frequency spectrum of radiation. Assuming that this relationship remains valid in the lossy case, we then obtain the expression for the electric field of reflected light E ref = ρE 0_ins , where ρ is the light reflection coefficient from the grating, which, if z 0 > z sat , is equal to where φ is the phase shift in the reflection from a flat metallic surface. Finally, the expression for the SPP amplitude is obtained by substituting the value of the parameter γ rad found from the relation (17) into expression (11).
It is important to note that application of the "lossless" relation (17) to the system of Equations (10) and (14) in the lossy case is an approximation that is justified only at sufficiently small values of parameter k I , and which, otherwise, can lead to errors in calculating the SPP amplitude, which will be discussed later. Hereafter, we will refer to this approximation as the "k I → 0 approximation".
Note that within this approximation, provided the dissipative and radiative losses of SPP are equal, the output of the long resonant (∆k = 0, Λ = λ SPP ) grating achieves the maximum possible SPP amplitude E (SUP) 0z . The reflection coefficient in this case, as (18) suggests, becomes zero. The height of the grating profile, ensuring that γ rad = k I , as follows from (8) and (17), can be calculated as where for some wavelengths of the optical range are presented in Table 1. is the maximum resonant amplitude of the longitudinal component of the SPP electric field at z 0 > z sat and x 0 = x opt .
The results of the calculations of the SPP amplitude dependence on the parameters of the gold grating, carried out following expressions (11) and (17) for the case of a normal incidence of a plane wave with λ = 0.8 µm, are shown in Figure 5. Figure 5 also illustrates the dependence of the SPP amplitude on the grating length z 0 when the height of its profile is x 0 = 5 nm. Note that the magnitude of the resonant grating period coincides with the SPP wavelength equal to λ SPP = 0.783 µm. As one would expect, the curve corresponding to this period in Figure 5a exhibits monotonic growth up to the saturation value. The oscillations in non-resonant curves are similar to those in Figure 4a for the SPP current excitation case. Figure 5b shows the dependence of the resonant (Λ = λ SPP ) SPP amplitude E (rns) 0z on the grating length z 0 . Note that the saturation length z sat for these dependencies is no longer constant. It is nearly constant as long as the grating profile height remains much smaller than x opt = 17.7 nm (solid curves at x 0 = 1, 2, 4 nm). Otherwise, it decreases sharply due to an increase in radiation damping (solid curves at x 0 = 25 and 75 nm), as follows from expression (12).  0z of the SPP excited by the long grating on its period. As expected, these curves have a Lorentzian profile similar to that for SPP excitation by a surface current. The resonance period, equal to 0.783 µm, is constant for all dependencies because the used model does not allow variation in the SPP wavelength as the height of the corrugation changes.
In Figure 5d, which shows the dependence of the amplitude E The nature of the dependence of the SPP amplitude on the grating parameters does not change for other wavelengths of the optical range compared to the one in Figure 5; therefore, these dependencies are not given. Of course, the values x opt and E (SUP) 0z for other wavelengths are different, as supported by Table 1.
Although this paper is devoted to the problem of SPP excitation by the grating and not to the reflective properties of the latter, some intriguing features of SPP excitation, as we will see later, are also apparent in the analysis of the grating reflection coefficient. In this regard, Figure 6 presents the results of analytical calculations of the reflectance R = |ρ| 2 , depending on the parameters of the long grating obtained per expression (18) at λ = 0.8 µm.
In this case, the solid curves in Figure 6a illustrate the dependence of the coefficient R on the grating period, and Figure 6b shows the dependence of the resonance value of this coefficient on the height of the grating. As the amplitude x 0 increases, the presented dependencies show first an increase in the depth of the dip in the grating reflection, up to the case R = 0, which occurs when the condition x 0 = x opt (γ rad = k I ) is satisfied. After that, the dip depth decreases, and the dependence expands due to the rapid growth of radiation attenuation. The results of numerical calculations of the SPP amplitude at λ = 0.8 µm performed for the geometry shown in Figure 1a, with the same technique of extracting the plasmon polariton contribution from the scattered field, as in the case of SPP wave excitation by a surface current, are presented by dotted and dashed curves in Figure 5. As can be seen, they agree fairly well with the analytical ones at sufficiently small (no more than 5 nm) values of the grating profile height. With a further increase in the profile height, one can observe a mismatch between the numerical and analytical results. In particular, in Figure 5c, the resonant wavelength shifts toward a shorter grating period, while for the analytical dependencies, the value of the resonant period is constant. This difference appears to be due to the effect of changes in the SPP mode propagation constant in the perturbed waveguide [24], the consideration of which is beyond the scope of this paper. However, we emphasize that we take the numerical values of the SPP resonance amplitude as the maximum values of the dashed dependencies in Figure 5c. Figure 5b shows that these values agree well with the analytical one, not only when x 0 x opt (curves at x 0 = 1, 2, 4 nm), but also at much higher than x opt heights (see the curve at x 0 = 75 nm). The differences, however, are mainly observed at x 0 ≈ x opt (curve at x 0 = 25 nm). This is also clearly noticeable when comparing the numerical (curve 2) and analytical (curve 1) results of calculations of the dependence of the SPP resonance amplitude on the height of the long grating profile x 0 in Figure 5d. This figure shows that the magnitude of the discrepancy between the analytical and numerical results reaches approximately 20% at x 0 = x opt (for other x 0 up to 3x opt , the discrepancy does not exceed 25%). Despite this, the positions of the maxima of the numerical and analytical curves visually coincide at x opt ≈ 18 nm. The agreement between the analytical and numerical data at a high grating height continues, at least up to the value of 100 nm (x 0 ≈ 6x opt ). The same picture is observed for other calculated wavelengths, as illustrated by the results given in Table 1.
It is also interesting to compare the widths of the analytical and numerical resonance dependencies presented in Figure 5c. When the height of the grating profile is low, these curves virtually coincide in width (curves at x 0 = 4 nm). However, as the height x 0 increases, the numerical resonance curves are narrower than the analytical ones. This feature is illustrated by curves 4 and 5 in Figure 5d, the first of which represents the results of analytical calculations of the FWHM value, according to (13), while the second is derived from numerical simulations.
Considering the reasons for the discrepancy between numerical and analytical results, we note that, in expression (8), if we reject the quadratic on the parameter x 0 term and pass to the "flat current" approximation, it leads to a much more pronounced difference between analytical curve 3, obtained in the approximation, and numerical dependence 2 in Figure 5d. Such a mismatch occurs not only at x 0 ≈ x opt , as it was when comparing curves 1 and 2, but also if x 0 > x opt , when it can reach 100% or more. In addition, the maximum of curve 3 ("flat current" approximation) shifts to the right, and now its position is remarkably different from the numerical value of x opt . The additional quadratic term in (8) leads to a substantially better agreement between the corresponding analytic curve 1 and the numerical curve 2. Although one would expect even better agreement between the analytical and numerical results after introducing the cubic terms in (8), direct analytical calculations show that taking this term into account practically does not change the form of the analytical curve 1 in Figure 5d for heights in the range of 0-100 nm. Therefore, one should look for another source of discrepancy between the analytical and numerical calculations. Such a difference may be related to the limitations of the Born and k I → 0 approximations.
A comparison of the results of numerical calculations of the grating reflectance, which are presented as dashed curves, and in Figure 6 for the λ = 0.8 µm case, with solid analytical curves, suggests the possible origins of this issue. The numerical results presented in the figure were obtained in the approximation of an infinitely long grating illuminated by an unbounded normal incident plane wave. Thanks to the periodicity, we performed all calculations within one grating period with periodic boundary conditions along the Z axis and a perfectly matched layer (PML) to absorb reflected light at the upper boundary of the calculation domain. Figure 6 demonstrates that the numerical values of the R coefficient tend to the value of the reflectance from a flat gold surface, equal to 0.988, in two cases: at x 0 → 0 and away from resonance. The analytical value for R under the same conditions tends to 1. This difference is an obvious consequence of the k I → 0 approximation since this approximation ignores dissipative losses, without which a flat surface of any metal has R = 1.
Another consequence of using this approximation seems to be more significant. A comparison of the analytical and numerical curves in Figure 6a exhibits a marked difference in the dip depth even at x 0 x opt (the curve at x 0 = 4 nm in Figure 6a). At the same time, the results of numerical and analytical calculations of the SPP amplitude at this grating's height presented in Figure 5 agree well with each other, which advocates the hypothesis that there is an error in the analytical calculation of the γ rad value due to the k I → 0 approximation.
Indeed, it follows from the relation (11) that the amplitude of the SPP excited by the long resonant grating is limited by the sum of dissipative and radiative losses according to the expression E (sat+rsn) 0z x opt , the radiative losses are negligible compared to the dissipative losses. Therefore, the amplitude E (rsn) 0z weakly depends on γ rad and, thus, on the error of its calculation for small values of x 0 . It is not the case, however, for the reflected wave. The first term in (14) for the reflected wave amplitude does not depend on the radiation losses, and the second term can be written as E lkg = − √ 2N 1 γ rad E 0z , as in (17). Therefore, the inaccuracy of the calculation of the parameter γ rad will directly affect the value of this summand and, as a consequence, the intensity of the reflected wave. Thus, the error in the calculation of the value of the coefficient R will appear regardless of the fulfillment of the condition γ rad k I , or equivalently, whether or not the grating's height is relatively small.
Further discussing the source of the discrepancy between the analytical and numerical calculations, we note that if we fix the height of the grating profile at 4 nm but increase the value of γ rad by a factor of 1.87, the corresponding analytical dependence for R(Λ) in Figure 6a (grey dashed curve) will almost completely coincide with the numerical one. With this modification of radiation attenuation, the analytical dependence E (rsn) 0z (z 0 ) for the resonant plasmon amplitude will shift slightly down in Figure 5b and coincide completely with the dashed curve for x 0 = 4 nm (not shown in this figure). One might expect that for other grating heights, such a technique would help to eliminate the discrepancy between the analytical and numerical dependencies for the SPP amplitude in Figure 5b. However, the magnitude of this correction differs for different x 0 . Thus, in the case where x 0 = x opt , the coincidence of the analytical and numerical values of E (rsn) 0z (z 0 ) is achieved when γ rad is 1.52 times larger, and when x 0 > x opt no additional correction is required. One can interpret these features as follows. As long as x 0 x opt , dissipative SPP losses dominate and cannot be assumed to be small, thus reducing the accuracy of the γ rad parameter calculation within the k I → 0 approximation and, as a result, requiring the maximum value of the correction factor. As the height of the grating profile increases, the contribution of dissipative losses to the SPP attenuation gradually decreases and eventually becomes negligibly small compared to radiation losses when x 0 x opt . In this case, the accuracy of the approximation k I → 0 increases, leading to the complete coincidence of analytical results with numerical ones. Figure 6a shows that as the profile height continues to increase, the resonance position shifts toward lower values of Λ (just as it does for the SPP resonance amplitude in Figure 5c). Therefore, we take the minimum value of the numerical dependence R(Λ) as the numerical value of the resonant reflection coefficient R rns .
In Figure 6b, curve 2 shows the results of numerical calculations of R rns as functions of the profile height x 0 of the long grating at λ = 0.8 µm. The minimum of this curve is slightly shifted to the left, relative to the minimum of analytical dependence 1, and is observed at the grating profile height x R=0 = 14.5 nm, which differs from x opt = 17.7 nm. Note that the numerically calculated height x R=0 for other wavelengths of the optical range also turn out to be smaller than x opt , which is illustrated by the data in the results of the calculations of this height in Table 1.
The limitations of the k I → 0 approximation can partially explain the difference in the position of the minimum of the analytical and numerical curves for R rns . Indeed, if one accepts the validity of the conclusion drawn above (that the maximum of the SPP resonance amplitude observed when the condition x 0 = x opt (γ rad = k I ) is satisfied), then the reflectance R rns cannot be zero at any k I = 0 for the real lossy metal. Indeed, it follows from expressions (14)-(16) that γ rad = k I provides E lkg = |E 0_ins |, |E r | = |rE 0_ins |. Within the k I → 0 approximation, |r| = 1, so the waves leaking out and reflected from the interface have equal amplitudes E lkg = |E r |. As a result, the destructive interferences of E lkg and E r waves nullify the reflection coefficient. However, due to the dissipative losses in a real metal, the value of |r| is always slightly less than 1, so E lkg > |E r |, and there is no complete mutual damping of the reflected and outgoing waves in the considered case. Therefore, the reflectance from the lossy metal becomes zero at a slightly lower grating height than x opt , when the amplitude of the outgoing wave E lkg equals E r , attenuated due to the reflection losses. The value x R=0 calculated from these considerations at λ = 0.8 µm is 17.5 nm, which is still noticeably (by 3 nm) greater than the numerical value x R=0 for this wavelength. This mismatch can seemingly no longer be attributed to the k I → 0 approximation and its source is in the limitations of the Born approximation. Indeed, the accuracy of the latter should decrease in the presence of a strong SPP wave arising at x 0 ≈ x opt . Therefore, at x 0 x opt , when such a wave attenuates, the accuracy of the k I → 0 approximation increases (as noted above), as well as the accuracy of the Born approximation. As a consequence, the numerical and analytical dependencies for R rns (x 0 ) at large values of the grating profile height converge to each other, as in Figure 6b (at least while this height remains smaller than 6x opt ).
Despite the difference in the dip depths, the widths of the analytical and numerical dependencies R(Λ) are the same for the small grating heights (curves for x 0 = 4 nm in Figure 6a). However, as x 0 increases, the numerical dependencies appear slightly wider than the analytical dependencies. This feature is illustrated by curves 6 and 4 in Figure 5d, the first of which represents the results of numerical calculations of the width of the resonance at half its depth (FWHD), and the second one represents the results of the corresponding analytical predictions, which, as expected, coincide with the same calculations for the half-width of analytical dependencies E (sat) 0z (Λ). Surprisingly, the numerical curve 6 for FWHD is shifted in the opposite direction from analytical curve 4 compared to numerical curve 5 for FWHM. This bidirectional shift most likely cannot be explained by the limitations of the k I → 0 approximation due to the use of the Born approximation.
Thus, the analytical model considered in this paper describes generally well the qualitative dependence of the amplitude of the excited SPPs on the parameters of the sinusoidal corrugation on the plasmonic metal surface. The quantitative coincidence of the analytical and numerical calculation results is observed for small (x 0 x opt ) and large (x 0 > x opt ) heights of the corrugation. However, for grating amplitudes in between, the numerical value of the resonance period of the grating turns out to be somewhat smaller than the analytical one, which requires further specification of the value of the SPP propagation constant for this case. Nevertheless, the model under consideration fairly accurately predicts the grating height x opt , at which it provides the maximal SPP amplitude. However, the value of the amplitude differs in the analytical and numerical calculations by less than 25%, which appears to be the result of the constraints of k I → 0 and Born approximations.

Conclusions
We presented an analytical description of SPP excitation by arbitrary configurations of nanostructured surfaces using a well-known method for calculating guided-mode amplitudes in waveguide theory. In this work, we focused on two cases: the case of a finite-length sinusoidal flat current and the case of a finite-length sinusoidal corrugation on a plasmonic metal surface under the normal incidence of a plane wave. For the first case, the described approach perfectly agrees with the numerical experiment. For the second case, our model correctly predicts the optimum amplitude of the corrugation to maximize the SPP excitation efficiency. We compared our analytical solution with numerical simulations and found excellent agreement between the two for non-resonant excitation and when the corrugation height was much smaller or larger than the optimum grating height. When the grating amplitude approaches the optimal value, the mismatch between the analytical and numerical values of the SPP amplitude was no more than 25%. We discussed the limitations of our analytical model that led to this mismatch. We believe that our approach is useful for modeling various phenomena associated with SPP excitation in metal nanostructures fabricated by laser processing or other methods.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
We start with the Lorentz reciprocity theorem in its unconjugated form [24], which is useful for describing lossy systems: where E i and H i are the solutions of Maxwell's equations for a system with sub-index i, defined by the spatial distribution of permittivity ε i (r) and magnetic susceptibility µ i (r); J i is the current density; ε 0 and µ 0 are the vacuum permeability and vacuum susceptibility, respectively.
Let us choose a gold-vacuum interface in a current-free planar configuration, where the back-propagating SPP is a system that corresponds to sub-index 1, so that H ← SPP = n y e ik SPP z e −κ 1 x , x 0 e κ 2 x , x < 0 (A4) where k SPP = k 0 ε 1 ε 2 /(ε 1 + ε 2 ), κ 1,2 = ε 1,2 k 2 0 − k 2 SPP , k 0 = 2π/λ. For symmetrical reasons, the origin is chosen at the center of the grating (see Figure A1), so that it would be useful for numerical simulations. Later, the origin can be moved to the end of the grating to follow the main text's thread.
x z � 1 � � Figure A1. Geometry of the problem. As the second system, we choose the same system as with sub-index 1 (with no grating), but with the presence of the incident plane wave and all the reflected and scattered waves. The grating is to be modeled with the non-zero current density J 2 , which eventually will be chosen to be the Born-like current density.
where E → SPP and H → SPP differ from (A3) and (A4) by the sign of k SPP , A describes the efficiency of the SPP excitation, andẼ 2 andH 2 combine all the terms that do not include forwardpropagating SPP.
With all the considerations, the right-hand side of (A1) simplifies down to the only term −E 1 J 2 . After integrating (A1) over the dashed region and exploiting the orthogonality relations [24], one can express A in terms of J 2 : where is the normalization constant. Importantly, (A6) allows one to calculate the amplitude of the excited SPP for any given current density distribution J 2 . Thus, the only issue here is to determine an appropriate function J 2 (r).
As part of the straightforward approach, let us define J 2 to be equal where E 0 is the incident field E 0 = e z e −ik 1 x + re ik 1 x , x 0 (1 + r)e −ik 2 x , x < 0 (A9) with k 1,2 = √ ε 1,2 k 0 , r = ( √ ε 1 − √ ε 2 )/( √ ε 1 + √ ε 2 ), and ∆ε = (ε 2 − ε 1 )1 −L/2<z<L/2      1, 0 < x < x 0 cos k g z −1, x 0 cos k g z < x < 0 0, otherwise where L is the total length of the grating. During the evaluation, there is a need to calculate the integrals of the following form: While the inner integration can be conducted rather easily, it results in the appearance of "double exponents", such as exp(a cos(bz)), with a and b being some combination of the problem's parameters. In order to acquire the compact and easy-to-use form of the analytical expression, exp(a cos(bz)) should be decomposed in series by a parameter a: exp(a cos(bz)) ≈ 1 + a cos(bz) + a 2 cos 2 (bz)/2 + . . . . Taking the first two summands in the series corresponds to the linear on the x 0 order.
We will now focus on the specific case of L = (n + 1/2)Λ, provided that n is an integer. After performing the corresponding calculations and substituting the found value of A into (A5), one obtains the expression for the electric field component parallel to the interface (1 + r)x 0 (1 + e ik SPP (n+1/2)Λ) ).
While it is possible to derive a similar expression for the general case of z 0 , even the considered case allows one to understand the limitations of this method. Examples of such limitations are given in the main text of this article.