Correcting the formalism governing Bloch Surface Waves excited by 3D Gaussian beams

Due to the growing number of publications and applications based on the exploitation of Bloch Surface Waves and the numerous errors and approximations that are used to evaluate their properties, we judge important for the successful interpretation and understanding of experiments to implement an adapted formalism allowing to extract the relevant information. Through comprehensive calculations supported by an analytical development, we establish generalized formula for the propagation length and the Goos-Hänchen shift, which are different from what is usually employed in the literature. The relative errors in the estimation of these two quantities are evaluated to vary between 50% and 200%. The effect due to a slight deviation of the angle of incidence or of the beam-waist position with respect to the structure are studied showing high effects on the Bloch Surface Waves properties. This formalism is adapted to any polarization-dependent Lorentzian-shape resonant structures illuminated by a polarized Gaussian beam. Theoretical models capable of accurately capturing the behaviour of plasmonic, Bloch surfaces waves are vital for the interpretation of experimental results. Here, the authors demonstrate the importance of extrinsic factors in determining the Goos-Hänchen shift in a generalised model of the propagation length.

Q uantification of the Bloch Surface Wave (BSW) properties is very important to predict its effectiveness in being used in integrated or in surface optics. BSWs are electromagnetic surface modes used to design different configurations for applications ranging from sensing 1-3 to surface-optics [4][5][6][7][8][9][10][11] or micro-manipulation 12,13 . BSW is of a great interest in integrated optics 14,15 due to its very large (>3 mm 16 ) Propagation Length (PL) and the possibility to be excited in both Transverse Electric (TE) and Transverse Magnetic (TM) polarizations contrarily to surface plasmon. Similarly to the latter, BSW can either be excited in the Kretschmann configuration (total internal reflection) 17,18 or more simply by diffraction 19,20 . However, 3D BSW electromagnetic-field distribution has never been theoretically reported except very recently by pure numerical methods (Finite Difference Time Domain 8 or Finite elements 21 ). This is a prerequisite for evaluating the two most important properties of BSW, namely its propagation length (PL) and lateral or Goos-Hänchen shift (L GH ), which will be defined later. This will be addressed through two different ways: (a) a rigorous method based on the Transfer Matrix Method (TMM) combined to the description of a 3D polarized Gaussian beam by an accurate Plane Wave Expansion (PWE) and, (b) an analytical calculation of the electromagnetic field associated to the BSW itself. As it will be demonstrated, both methods converge to the same result that fails the commonly used formulas in the literature. For the PL, we establish an equation that is widely valid for any surface wave excited within high quality-factor resonance having a Lorentzian shape (surface plasmon, Fano, membrane mode, symmetry protected modes, Bounded In the Continuum (BIC) modes...). For the L GH , we demonstrate its value to be dependent on the incident beam dimension while it is currently considered as intrinsic to the structure itself.
On the one hand, several studies 16,[22][23][24] used theoretical formulas based on a development obtained for plane wave illumination 25,26 . For example, Soboleva et al. 23 used the formula given in Eq. 1 of that paper to discuss the occurrence of a giant Goos-Hänchen shift on the reflected beam issued from the excitation of a BSW. In that paper, the measured reflectance angular spectrum is used to estimate the Fano profile of the resonance and subsequently operated to evaluate the lateral shift of the reflected beam. Nevertheless, as in most theoretical studies [27][28][29] , the approach used to assess the reflected beam distribution is based on the consideration of one-dimensional angular distribution for the beam (see Eq. 3 of ref. 23 ) assuming that the incident beam is a 2D-Gaussian beam (prismatic) instead of a realistic 3D beam. Such approximation leads to less reliable physical properties of the studied phenomenon as it will be discussed in more details below. On the other hand, in diverse studies 22,23,30 , the use of the reflectance spectrum to estimate the BSW properties is somewhat questionable. In fact, the BSW corresponds to a surface mode that is excited in the total internal reflection condition meaning that the reflection coefficient is equal to 100% in amplitude for purely dielectric flat layers. Consequently, the signature of the BSW excitation on the reflection coefficient only involves its phase but neither the amplitude nor the intensity that is usually experimentally measured. When the reflectance spectrum exhibits a dip resonance, this gives directly the effective index of the BSW (through the tangential wave-vector component) but means above all that losses occur by scattering or by absorption as explicitly studied by Michelotti et al. 30 . In this case, the relationship between the angular width of the reflectance dip and the BSW properties is no longer intuitive. More, precisely, when absorption (imaginary part of the dielectric constant) occurs, the dip in reflectance and the maximum in transmission coefficient do not correspond to the same value of the angle of incidence (as usually encountered in the case of surface plasmon resonance). In fact, the mathematical solution corresponding to the extrema of reflectance and transmission leads to two different values of the angle of incidence while the analytical expression of these two quantities have the same denominator 31 . Assuming that the angle of incidence is a real value causes a small shift between their angular positions.
Here we provide theoretical models capable of accurately determine the two main properties of surface wave that are vital for the interpretation of experimental results. We also demonstrate the importance of extrinsic factors (illumination conditions) in determining the Goos-Hänchen shift and the propagation length occurring within the excitation of such surface waves through a generalised formalism applicable to Lorentzianshape resonances.

Results
Proposed structure and plane wave analysis. For illustration purposes, we consider a typical configuration of one-dimensional Photonic Crystal (1D-PhC) by optimizing its geometry using a plane wave illumination through a very simple algorithm based on TMM (see details in Supplementary Notes 1 and 2) that links the electric incident and reflected field amplitudes to the transmitted and back reflected ones on the interface of two different layers. The total transfer matrix, which is the product of all single matrices, determines the transmitted and reflected amplitudes over the entire multilayered system (see Supplementary Eqs. S.4 and S.5) taking into account all the geometrical and physical parameters of the structure (thicknesses and permittivities) and the incident plane wave properties (polarization, wavelength, angle of incidence). The eigenvalues of this total matrix are the eigenmodes of the structure that can be simply calculated through basic inverse matrix algorithm.
We use this TMM method to adapt a multilayer design 22,32 that consists of N-periods of bilayered stacks (see Fig. 1a) to operate at telecom wavelength in TE polarization. All geometrical parameters are given in the caption of Fig. 1. Figure 1b shows the square modulus of the transmitted electricfield amplitude (in logarithmic scale) at the upper interface as a function of the bilayer number (N) at λ = 1550 nm. Note that additional Bloch modes (indicated by H 1 to H 3 on Fig. 1b) occur inside the PhC structure as discussed in Supplementary Note 3 and presented on Supplementary Fig. 2. The angle of incidence and the natural logarithm of the Full Width at Half Maximum FWHM) Δθ T of the BSW resonance are given on Fig. 1c. As expected, the angular position θ BSW converges asymptotically but promptly (N ≃ 7) to the value corresponding to the pole of the transmission coefficient of the infinite structure. Δθ T varies exponentially with N (see red with stars line on Fig. 1c) meaning rather similar variations for the BSW resonance quality factor defined by: This exponential increase of the Q-factor with N will be very hard to be experimentally verified because, in practice, losses due to scattering by surface defects and by material absorption lead to a finite value of the Q-factor even if the number of stacks increases 33 (see Fig. 1d). Although such loss mechanism is quite hard to be quantified, most of the authors agreed to model it by adding a small imaginary part n″ to the optical refractive index. Unfortunately, when introducing such absorption (One notes that absorption is fundamentally proportional to the imaginary part of the permittivity and not to that of the optical index meaning that θ BSW will be also affected by absorption (Kramers-Kronig relations)) through n″ for all media (except glass substrate), the BSW angular position remains the same while the BSW efficiency becomes weaker. In our case, we estimate that BSW excitation is negligible for n″ > 10 −3 . Figure 1d shows the quality factor Q variations with the number of b-layers N for different values of the imaginary part n″. For loss-less materials, Q tends to infinity as Q = e 0.8623N+0.0586 while asymptotic behaviors occur for n″ ≠ 0. We have verified that the constants in the last relation only depend on the effective index associated with the BSW excitation (here n eff ¼ n 1 sin θ BSW ¼ 1:3928).
To go further through formal calculation, both the Gaussian shape of the illumination beam and the transmission coefficient of the structure must be analytically expressed (see Supplementary Notes 4 and 5, respectively). Fortunately, in the case of a BSW excitation, the transmission coefficient spectrum can be realistically approached by a Lorentzian function (see more details in the Supplementary Note 5) leading to express it as: where k BSW x is the tangential wave-vector component associated with the BSW excitation and t max is the value of the transmission coefficient for k x ¼ k BSW x calculated through the TMM. Δk is the FWMH of the square modulus of the reflection or the transmission coefficient given by: Modeling of the polarized 3D Gaussian beam. In a real experiment, a finite beam (commonly Gaussian spatial shape) is used to illuminate the multilayered structure either by the Kretschmann configuration or by diffraction. To model such a beam, the plane wave spectrum PWE (or angular spectrum) method is used by coherently summing the amplitude response of all the plane waves composing the Gaussian beam (see Supplementary Note 4 for more details). This can be done over the entire structure even inside the layers. The angular spectrum of a 3D polarized Gaussian beam has been described since 1999 34 and verified by comparison with experimental and/or results based on different methods [34][35][36] . An extended formalism from linear to elliptical or circular polarized beam is given by Supplementary Eqs. S.6-S.8. The transmitted electric-field distribution associated with the BSW is then calculated in the direct space Oxyz through the Supplementary Eq. S.15. The latter involves the transmission Jones matrix of the structure that is basically given by the TMM astðk x ; k y Þ ¼ ÀTT À1

21
TT 22 (see the Supplementary Eq. S.4). All results calculated through this integral are obtained without any approximation meaning that the vectorial character of both the incident field and the transmission coefficient is taken into account. Nonetheless, due to the resonant character of the transmission, one can reduce the calculation to a scalar equation by considering only the resonant term of the transmission (for instance the TE term in our case). Replacing the transmission coefficient through its expression given by Eq. (2) and after fastidious algebra (see Supplementary Note 5), the transmitted electric-field amplitude is analytically expressed as a function of the beam-waist W 0 and the FWHM (Δk) of the transmission coefficient through: Where erf is the error function defined by erf ðxÞ ¼ 2 ffiffi π p R x 0 e Àx 2 dx and Δk is the FWHM of the transmittance spectrum defined before.
Equation (3) provides all the BSW properties (PL, L GH , maximum efficiency...) as it will be discussed below.
Numerical simulations. Based on combination of TMM and PWE methods introduced in the previous section (see Supplementary Eq. S.13), one can calculate the electric-field distribution over all the structure for any illumination direction, beam-waist or polarization. This versatile character is demonstrated in Fig. 2 that shows the electric-field intensity distribution in the mean plane of incidence (Oxz) across the whole structure for a BSW excitation with a Gaussian beam of, W 0 = 10 μm in Fig. 2a, Fig. 2b and W 0 = 1 mm in Fig. 2c. The spatial shape of the excited BSW greatly depends on the value of the incident beam-waist W 0 . When the beam waist is small, only a portion of the incident energy is coupled to the BSW giving rise to a comet shape for the intensity distribution of the BSW at the top interface. In this case, the calculation of the PL is simple. When W 0 increases (Fig. 2b), the angular aperture of the beam decreases and the overlap with the BSW grows resulting in a more efficient excitation of the latter. Nevertheless, the comet shape becomes less evident due to the competition between the propagation length and the beam width itself. When the beam-waist is very large (Fig. 2c), the comet shape completely disappears in the face of the Gaussian shape. In all three cases, we can clearly see that large electric-field confinement occurs in the top layer. For the sake of clarity, the vertical scale in the substrate zone is chosen to be large enough to see both the incident and reflected beams. The latter is greatly affected by the BSW excitation and appears to be split into two asymmetrical beams when the incident beam waist is small enough due to the presence of large out-of-BSW spectral (angular) components.
From the numerical results on can determine the BSW characteristics corresponding to experimentally observed quantities that are recorded within the transmitted near-field, namely the lateral or Goos-Hänchen shift and the PL. Other properties could also be explained such as the Imbert or transverse shift 37 , or the angular shift of the secondary reflected beam 38 . The latter phenomena originated from the spin-orbit coupling between light and a flat interface, occur on the reflected beam and are mediated by the angular dispersion of the reflection coefficient 39 . Generally, they occur only with circular or elliptical incident polarized beams. Furthermore, two different definitions are still used for the L GH , which can be considered as the displacement of the maximum of the intensity or, that of the intensity centroid 40 . Nonetheless, it is commonly agreed to consider the maximum intensity shift in cases where large propagation lengths occur such as for surface plasmon resonance or BSW 16,34 . Consequently, we will restrict our definition to this last one as indicated on Fig. 1a. Note that the Goos-Hänchen shift also exists for acoustic waves and was recently studied by analogy with optics 41 . Additional properties dealing with the reflected beam are also reachable as it will be discussed in the following. Figure 3a shows the 3D map of the BSW electric near-field intensity distribution at the top interface in a xOy plane as it can be measured by means of Scanning Near-field Optical Microscope (SNOM). We can clearly see the surface wave character through the intensity decay that occurs along the propagation direction (Ox here). Top-view distributions are given in Fig. 3b showing the excitation of the BSW (comet shape) only in TE polarization and highlighting the lateral shift that accompanies this excitation.
In some recent experimental studies, it was sometimes found that the near-field images of the BSW present a different behavior compared to what is expected (a pure comet shape) as pointed out in Fig. 3b. For example, Dubey et al. 16 show (see Fig. 4b of that paper) the cross-section made over the intensity map along the propagation direction of the BSW to exhibit a depletion next to the maximum. At a first glance, this effect can be attributed to a surface irregularity of the top interface. In fact, by introducing an angular mismatch less than 1 ∘ on the angle of incidence, numerical simulations allow reproducing an almost identical behavior as shown in Fig 3c, d. From Fig. 3a or 3b, we determine both the spatial position of the intensity maximum that gives L GH = 649 μm and the PL = 1.37 mm. Nonetheless, there is another parameter, which is difficult to experimentally estimate, and which could also affect the excitation of the BSW, namely the incident beam defocusing. In fact, in all numerical simulations the beam-waist is supposed to be centered on the top of the substrate. Figure 4 shows two different cases of defocusing. Both of them correspond to the N = 7-structure illuminated by a beam-waist (W 0 = 5 μm) Gaussian beam. The first one (Fig. 4a) corresponds to a beamwaist located 300 μm under the PhC structure while it is supposed to be 100 μm above the substrate-PhC interface in the second (Fig. 4b). As in Fig. 2, the calculated amplitudes of the electric field are mapped in the Oxz plane with different spatial scales. In the first case, oscillations affect the BSW itself especially near its intensity maximum (see the blue dashed line at the top of Fig. 4c) while an additional lateral shift of this maximum occurs in the second case (solid black line). This demonstrates how the BSW shape can be significantly affected by a small focusing default of the incident beam. In addition, another effect arises on the interference pattern appearing in the reflected beam due to the spatial broadening of the beam falling the s. In fact, the total lateral shift at reflection becomes greater and leads to increase the spatial separation between the different angular components of the incident beam. The region where the beam impinges the first interface is emphasized in the blue rectangle in the bottom of Fig. 4a. One can see the occurrence of curved fringes similar to caustics resulting from the interference between two highly focused incident and reflected beams. This demonstrates the difficulty of interpreting some experimental results but it also shows the way to have an effective excitation of the BSW.
The reflected beam. Experimentally, the excitation of BSW is controlled by exploiting the reflected beam (presence of a dip in the reflectance). Consequently, the properties of the latter deserve to be understood to extract information about the BSW excitation. In particular, the oscillation pattern appearing on the reflected beam in the case of strongly focused beams is often highlighted as a signature of the BSW excitation 42 . Very recently, Petrova et al. 3 exploited the properties of the reflected beam for biosensing applications. Several theoretical studies have been performed in this context [27][28][29] but all of them considered a 2D-Gaussian beam (prismatic beam) instead of a realistic 3D beam. In those references, the authors studied the effect of the angular dispersion of the Goos-Hänchen shift and they linked it to the fringe pattern that appears on reflection. To point out this phenomenon, which occurs also in Surface Plasmon excitation within the Kretschmann configuration, we consider an incident beam with W 0 = 5 μm illuminating the 1D-PhC in the case of N = 7 and we calculate the electric-field distribution in three different planes. Figure 5a shows the electric-field amplitude in the Oxz plane as in Fig. 2. The fringe pattern is clearly apparent on the reflected beam. A zoom-in over the reflected beam cross-section in the Oxy plane in the substrate, at z = 1 mm below the first interface, is shown in Fig. 5b. The spatial oscillations of the electric-field intensity are perfectly visible. However, experimentally, the reflected beam is observed perpendicularly to its propagation direction as in Fig. 5c where the presented electricintensity distribution is evaluated through the TMM/PWE algorithm without any projection operation nor symmetry considerations. To the best of our knowledge, this is the first time that such images are calculated in the case of a real 3D Gaussian beam. In fact, the 2D calculations lead to a similar pattern but with different oscillation features. Figure 6a, b shows a transverse cross-section (along the Ox axis) 1 mm under the first interface (substrate-PhC) for a beamwaist varying from W 0 = 5 μm to W 0 = 50 μm in the case of 3D and 2D-Gaussian beams respectively. At a first glance, the two results seem to be very similar. However, even if the global shapes are comparable, the cross sections shown on Fig. 6c (where the beam waist was fixed to W 0 = 6.87 μm for both simulations) show a quantitatively different result. The oscillations are not at all concordant and their intensity levels are clearly different. This is directly due to the contribution of the plane waves that are out of the incidence plane. Indeed, even if the global polarization of the beam is TE, these out-of-incidence plane waves exhibit TM components whose weight increases as their propagation direction falls out from the plane of incidence. Nonetheless, the Gaussian envelop of the beam amplitude produces a two-lobes shape as shown on Fig. 2c by Bouhelier et al. 36 . This is explicitly where the x − and z − components of the electric field are not zero even in TE polarization (χ = π/2 → α = 0 and β = 1). Although, the contribution of these TM components to the BSW (in transmission) is negligible because only TE components resonate, this does not prevent their contribution to the reflected beam. This interpretation is derived from the Maxwell-Gauss equation (DivE = 0), which must be fulfilled for each incident plane wave of the field expansion in the Fourier space. This implies a depolarization term that appears for all waves with wave-vector that is not located in the plane of incidence. This is true for both TE and TM polarized beams as shown by Supplementary Eqs. S.11 and S.12 where the field was also expressed in the (TE,TM) basis. Consequently, the response of a realistic Gaussian beam cannot be calculated by limiting the plane wave expansion over only one spatial frequency component (the k x one) as it is done by Andaloro et al. 27 . In the latter paper, authors claimed that the y-dependence can be suppressed because it does not affect the beam-interface interaction, which is rigorously false especially if we deal with the reflected beam. According to us, the same argument is at the origin of the clear discrepancy, in terms of number of oscillations and amplitude, between the experimental and theoretical results as shown on Fig. 2 of the paper by Petrova et al. 3 . Therefore, a quantitative exploitation or comparison with experimental results must take into account the contribution of these components.   25,26 estimating these two quantities to be: where Δθ R is the FWHM of the dip resonance appearing in the reflectance spectrum and ϕ is the phase of the transmission coefficient through the whole structure. Experimentally, the transmission phase variation can hardly be measured. Nevertheless, as it is well-known, this phase is equal to the half of the reflection coefficient phase. Consequently, assuming an interferometric detection (heterodyne), one can reach the reflection phase value. Unfortunately, this proportionality between the two phases of transmission and reflection coefficients is no longer valid when dealing with absorption. However, the L GH of the BSW cannot be obtained by any far-field detection of the reflected beam. Only direct measurement of the near-field allows access to this property. Theoretically, the variation of ϕ with the angle of incidence is given in Supplementary Fig. 3. From this figure, and according to Eq. (4), we estimate the theoretical values of the Goos-Hänchen shift to be constant (L GH = 770 μm), which is inconsistent with the calculated values from Fig. 2 that depend on the beam dimension. This discrepancy needs to be elucidated.
To this end, we consider at first the same N = 7 bilayer 1D-PhC and we calculate the evolution of the L GH as a function of the beam-waist value through the TMM/PWE algorithm. Figure 7a shows that L GH significantly varies with W 0 as long as the angular width of the beam ( λ πW 0 ) is 22 times larger than the angular width of the BSW (here Δθ T = 0.642 mrad) corresponding to W 0 ≃ 1 cm as shown by the solid blue line on Fig. 7a. The effect of absorption is also studied on the same Fig. 7a where we consider the two cases of n″ = 10 −4 (red dashed line) and n″ = 10 −3 (dashed dotted green line) and study their impact on the L GH . As expected, the latter varies dramatically with losses. Second, the evolution of L GH depicted in Fig. 7a shows an asymptotic behavior limiting it to a maximum value to the propagation length PL (see Supplementary Note 6) independently from the beam dimension. This obviously contradicts the results obtained by Konopsky et al. 29 where formula 2 of that paper states that L GH is proportional to the square root of the beam diameter. Nonetheless, all the demonstrations made in that paper are formulated for Gaussian beams illuminating an interface near the critical angle in total internal reflection configuration, which is different from our case where a sharp resonance with a large Qfactor of 1856 occurs.
To clarify the issue, we make use of the analytical expression of Eq. (3). The spatial position of the transmitted intensity maximum corresponds to the value of x for which the xderivative of the square modulus of the electric-field amplitude given by Eq. (3) vanishes. This condition leads to Eq. (5) that was numerically solved (see the Supplementary Note 6) for the three cases of Fig. 7a. They correspond to the colorful thick vertical lines depicted on the same figure showing a perfect agreement with the TMM/PWE (see large circles corresponding to intersection of these thick vertical lines with the solid blue line). We have verified the relative error to be less than 5 × 10 −3 .
For the PL, Eq. (4) cannot be exploited if we assume purely dielectric structure without any absorption because, as mentioned above, the reflectance is equal to 100% and does not exhibit any dip. Nevertheless, introducing a small absorption allows the presence of a dip in the calculated reflectance. Experimentally, this dip can bring all one needs to determine the BSW propagation length PL due to the fact that Δθ R ≈ Δθ T even if absorption occurs. As determined from Fig. 2, we obtain a propagation length of PL ≈ 1.374 mm independently of the beamwaist value while Eq. (4) leads to an almost twice smaller value of PL = 769 μm. We notice that expression of PL given by Eq. (4) is commonly used to interpret or to exploit experimental results 16,22,24 . Again, we are in front of a contradiction that needs to be clarified.
For this purpose (see Supplementary Note 7 for more detail), we will still consider the analytical expression of the transmitted field given by Eq. (3) where we can clearly see that for x → ∞(x ≫ W 0 ), the predominant term in the amplitude expression is e À Δkx 2 (erf(∞) → 1) that corresponds to the electricfield behavior far from its maximum. This allowed expressing the propagation length as: Replacing θ m = 1.189 rad, n 1 = 1.501, and Δθ T = 0.642 mrad into Eq. (6) leads to PL = 1.376 mm, which perfectly agrees the value (1.374 mm) graphically determined through the results of the TMM/PWE algorithm. We have verified the good agreement between the TMM/PWE results (solid line) and this analytical formula (green circles) on Fig. 7b when the imaginary part of the refractive optical index vary from 10 −6 to 10 −1 . This perfect agreement between a rigorous numerical method and the mathematical formulation of the transmitted field is an indisputable proof of the accuracy of the two methods.
In a real experiment, the detection is made in air not in the substrate medium. Assuming that a hemispherical lens is used as a substrate (normal incidence between air and lens), one can express the FWHM of the reflectance in air by : Δθ air R ¼ n 1 Δθ R . Unfortunately, in most experimental studies, it is very difficult to know if the presented reflectance spectrum corresponds to the angle measured in air or in the substrate.

Discussion
From Eqs. (4) and (6), we can determine the correction term for the propagation length to be C S PL ¼ 1 n 1 cos θ m in the substrate or C A PL ¼ 1 cos θ m in air. The latter becomes significantly important for larger value of θ m (larger effective index BSWs). For the studied structure in this paper, this correction is equal to C S PL %¼ 1:8 meaning a relative error between the two PL values of 80%.
To validate our formalism through experimental measurements, let us consider the configuration studied by Descrovi et al. 22 : from near-field detection, authors measured an experimental value of PL exp ¼ 470 μm, while the use of Eq. (4) leads to a lower value of LP = 448 μm. By considering the correcting term given by our formalism, and taken into account the presence of the substrate, we find PL ¼ 448:5=ð1:501 cos 50:32 Þ ¼ 468 μm that agrees perfectly the experimental measured value. Generally, BSW configurations are designed to operate at large angle of incidence θ m compared to the critical angle to avoid a direct transmission of a part of the incident light especially if we work with highly focused beams. This causes the relative error to increase drastically reaching 580% at θ m ≈ 80°. For the configuration studied by Michelotti et al. 30 , the correcting term is equal to C PL = 1.64, which would correspond to an error of 64%.
For the Goos-Hänchen shift L GH , the error depends on the beam-waist value (W 0 ). In our case, this error can reach 86% for highly focused beams (W 0 = 10 μm) and is only twice smaller (43%) in the case of wide beams as considered in the study by Soboleva 23 .
In summary, the combination of the PWE with the TMM, and the use of an accurate angular spectrum expansion of a Gaussian beam, turn out to be a powerful tool for simulating and conceiving 1D-PhC structures dedicated to surface wave excitation. The use of the PWE can be extended to integrate any other method (Rigorous Coupled Wave Analysis (RCWA) for instance) able to take into account diffraction by grating (periodic) or by individual pattern 5,14,43,44 . The examples discussed in this paper demonstrate the versatility of this tool that allows highlighting and estimating the unexpected effects of some external parameters (alignment error, focusing default, presence of adhesion layer on the top surface) on the excitation of the surface wave. This combination is essential for the calculation of the reflected beam through the consideration of a realistic incident 3D beam. The major result of this paper is obtained through rigorous analytical mathematical development that leads to a significant correction of the two important properties of the BSW, namely the lateral shift (Eq. (3)) and the propagation length (Eq. (6)), for which inaccurate formulas are so far commonly used in the literature.

Data availability
Materials and data that support the findings of this research are available within the paper. All data are available from the corresponding author upon request.

Code availability
Numerical codes that support the analytical demonstration related to this research are available from the corresponding author upon request.