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 gross errors and approximations that are regularly used to evaluate the properties of this type of wave, we judge seriously important for successful interpretation and understanding of experiments to implement adapted formalism allowing to extract the relevant information. Through a comprehensive calculation supported by an analytical development, we establish a generalized formula for the propagation length which is different from what is usually employed in the literature. We also demonstrate that the Goos-H\"anchen shift becomes an extrinsic property that depends on the beam dimension with an asymptotic behavior limiting its value to that of the propagation length. The proposed theoretical scheme allows predicting some new and unforeseen results such as the effect due to a slight deviation of the angle of incidence or of the beam-waist position with respect to the structure. This formalism can be used to describe any polarization-dependent resonant structure illuminated by a polarized Gaussian beam.

1 Quantification of the BSW properties is very important to predict its effectiveness in being used in integrated or in surface optics. BSW s are electromagnetic surface modes used to design different configurations for applications ranging from sensing [1,2] to surface-optics [3][4][5][6][7][8][9][10] or micro-manipulation [11,12]. This kind of surface waves is of a great interest in integrated optics [13,14] due its very large propagation distance and the possibility to be excited in both T E and T M polarizations contrarily to surface plasmon. Similarly to the latter, BSW can either be excited in the Kretschmann configuration (total internal reflection) [15,16] or more simply by diffraction [17,18]. However, 3D BSW electromagnetic field distribution has never been theoretically reported except very recently by pure numerical methods (Finite Difference Time Domain [7] or Finite elements [19]). This is a prerequisite for evaluating the two most important properties of BSW , namely its propagation length P L and lateral or Goos-Hänchen shift L GH , which will be defined later on. This will be addressed through two different ways: (a) a rigorous method based on the Transfer Matrix Method T MM combined to the description of a 3D polarized Gaussian beam by an accurate Plane Wave Expansion P W E and, (b) an analytical calculation of the electromagnetic field associated to the BSW itself. As it will be demonstrated, the two methods converge to the same result that fails the commonly used formulas. For the P L, we establish a new 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 modes...). For the L GH , we demonstrate its value to be dependent on the incident beam dimension, which is completely innovative, compared to the accepted ideas for which this property is intrinsic to the structure itself.
As mentioned below, our findings are in great contradiction with commonly used formulas. On one hand, several studies [20][21][22][23] used theoretical formulas based on a development obtained for plane wave illumination [24,25]. For example, in ref. [21], the formula given in Eq. 1 of that paper is used to discuss the 2 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 then operated to evaluate the lateral shift of the reflected beam. Nevertheless, as in most theoretical studies [26][27][28], 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 in [21]) meaning that the incident beam is a 2D Gaussian beam (prismatic) instead of a realistic 3D beam. This certainly 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 as in [21], 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 and never its amplitude nor its 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. In this case, the relationship between the angular width of the reflectance dip and the BSW properties is no longer intuitive.

Proposed structure and plane wave analysis
First, we consider a typical configuration of 1D-PhC by optimizing its geometry using a plane wave illumination through a very simple algorithm based on T MM (see details in SI file) that links the electric incident and reflected field amplitudes to the transmitted and back reflected ones on the interface separating two different layers. The total transfer matrix, which is the product of all single matrices, allows determining the transmitted and reflected amplitudes over the entire multi-layered system as a function of the incident one (see Eqs E.4 and E.5 of the SI file) by taken into account all the geometrical and physical parameters of the structure (thicknesses and permittivities) and the 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 kind of calculation to adapt a multilayer design [20,29] that consists on N-periods of bi-layered stacks (see Fig. 1a) to operate at telecoms wavelength in T E polarization. All geometrical parameters are given in the caption the figure 1.
This exponential increasing of the Q-factor with N does not have a real physical meaning because, in practice, losses due to scattering by surface defects and by material absorption lead to a finite value of N (see Fig. 1d) for which the structure is optimized (maximum Q-factor) [30].
Where k BSW x is the tangential wave-vector component associated to the BSW excitation and t max is the value of the transmission coefficient for k x = k BSW x calculated through the T MM.

Modeling of the polarized 3D Gaussian beam
In the real experiment, a finite beam (commonly Gaussian spatial shape) is used to illuminate the multi-layered structure both in the Krestschmann configuration and either by diffraction. To model such a beam, the plane wave spectrum P W E (or angular spectrum) method is used by coherently summing the amplitude response of all the plane waves composing the Gaussian beam (see SI file for more details). This can be done over the entire structure even inside the layers. The angular spectrum of a 3D polarized Gaussian beams was described in ref. [31] and tested several times through comparison with experimental and/or results based on different methods [31][32][33]. An extended formalism from linear to elliptical or circular polarized beam is given by Eqs E.6 -to E.8 of the SI file.
The transmitted electric field distribution associated with the BSW is then E.14.
After fastidious algebra (see SI file), the transmitted electric field amplitude is analytically expressed as a function of the beam-waist W 0 and the F W HM (∆k) of the transmission coefficient through: Where erf is the error function defined by erf (x) = 2 √ π x 0 e −x 2 dx and ∆k = 2πn 1 cos θ m ∆θ T is the F W HM of the transmittance spectrum as defined above.
Equation 3 provides determining all the BSW properties (P L, L GH , maximum efficiency...) as it will be discussed below.

Results and discussion
Within T MM/P W E combination (Eq. E.13 of the SI file) one can calculate the electric field distribution over all the structure for any illumination direction, beamwaist or polarization. This versatile character is demonstrated through figures 2 that present the electric field intensity distribution in the mean plane of incidence (Oxz) across all the structure for a BSW excitation within a Gaussian beam of, W 0 = 10 µm in Fig. 1a, W 0 = 30 µm in Fig. 1b and W 0 = 1 mm in Fig. 1c. 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, the overlapping between the incident beam and the BSW excitation is weak and a small part 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 determination of the propagation length P L is very easy. When W 0 increases (Fig. 2b), the angular aperture of the beam decreases and the overlap grows resulting in a more efficient excitation of the BSW . 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 the 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 fixed differently 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 asymmetric beams when the incident beam waist is small enough due to the presence of out of BSW spectral (angular) components.
From such numerical results on can determine the BSW characteristics corresponding to experimental observed quantities that are recorded on the transmitted near-field, namely the lateral or Goos-Hänchen shift and the propagation length. Other properties could also be determined such as the Imbert or transverse shift [34], or the angular shift of the secondary reflected beam [35]. These two last quantities, deriving 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 [36]. Generally, they need circular or elliptical incident polarization to take place. Furthermore, two different definitions are still used for the L GH assuming it as, the displacement of the maximum of the intensity or, that of the intensity centroid [37]. Nonetheless, it is commonly agreed to consider the maximum intensity shift in cases where large propagation distances occur, such For the sake of clarity, the electric field was auto-normalized in three different zones: the incidence, the multilayer and the transmission zone. In addition, the scale in the incidence zone varies with the beam dimension as to show both incident and reflected beams.
as for surface plasmon resonance or BSW [23,31]. Consequently, we will restrict our calculation to this last definition as indicating on figure 1a. Note that Goos-Hänchen shift also exists for acoustic waves and was recently studied by analogy with optics [38]. Additional properties dealing with the reflected beam are also reachable as it will be discussed in the following. . c and d correspond to an angle of incidence θ m = θ BSW +1 o . Experimentally, this kind of distributions is measured by means of scanning near-field optical microscope to estimate both the GH shift (differential value between T E and T M) and the propagation length [23]. Note that the intensity maximum value is 70 times greater in T E than T M.
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 in figure 3b. For example, in figure 4b of reference [23] the cross-section made over the intensity map along the propagation direction of the BSW exhibits 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 o on the angle of incidence, numerical simulations allow reproducing an almost identical behavior as shown in figures 3c and d. From figure 3a or b, we determine both the spatial position of the intensity maximum that gives L GH = 649 µm and the P L = 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 Gaussian beam with W 0 = 5 µm. The first one (Fig. 4a) corresponds to a beam-waist 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 curved fringes similar to caustics resulting from the interference between the incident and the reflected angular-wide beams. All this demonstrates the difficulty of interpreting some experimental results but 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 [39]. Very recently, Petrova et al. [2] exploited the properties of the reflected beam for biosenseng applications. Several theoretical studies have been performed in this context [26][27][28] but all of them considered a 2D-Gaussian beam (prismatic beams) instead of a realistic 3D-beam. In those references, the authors studied the effect of the angular dispersion of the GH shift and they linked it to the fringe pattern that appear 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. Although, experimentally, the reflected beam is observed perpendicularly to its propagation direction as in Fig. 5c where the presented electric intensity distribution is evaluated through the T MM/P W E algorithm without any projection operation nor symmetry considerations. According to us, 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 similar pattern but with different oscillation features. At a first glance, the two results seem to be very similar. Unfortunately, even if the global shape is comparable, Fig. 6c (where the beam waist was fixed to W 0 = 6.87 µm for both simulations) disclaims it. Although, the oscillations are not at all concordant and their intensity level 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 T E, these out of incidenceplane waves exhibit T M 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 amplitude shape (see Fig. 2c  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 wavevector that is not located in the plane of incidence. This is true for both T E and T M polarized beams as shown by Eqs. E.11 and E.12 in the SI file where the field was also expressed in the T E, T M 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 in ref. [26].
In the latter, 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 fallacy is at the origin of the clear discrepancy, in terms of oscillation number and amplitude, between the experimental and theoretical results in Fig. 2 of ref. [2]. Therefore, a quantitative exploitation or comparison with experimental results must take into account the contribution of these components.

Goos-Hänchen shift and Propagation distance
The number of bi-layers is fixed to N = 7 in the following as in Fig. 2 from which one determined the L GH and P L for the three beam-waist values to be: c, Cross sections made over the two maps for W 0 = 6.87 µm are plotted showing a real discrepancy between the two oscillation patterns.
W 0 = 1 mm. Even if the propagation distance is almost constant, its value, in addition to the evolution of the L GH , is in clear contradiction with a simple theory based on plane wave analysis [24,25] estimating these two quantities to be: where ∆θ R is the F W HM of the dip resonance appearing in the reflectance spectrum and φ is the phase of the transmission coefficient through the whole structure.
Experimentally, the phase variation can hardly be measured. Nevertheless, as it is well-known, this phase value is equal to the half of the reflection coefficient one.
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 farfield detection of the reflected beam. Only direct measurement of the near-field allows access to this property. Theoretically, the variation of φ versus the angle of incidence is given in figure S3  To this end, we consider the same N = 7 bi-layer 1D-PhC and we calculate the evolution of the L GH as a function of the beam-waist value through the T MM/P W E 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 greatly varies with losses.
Second, the evolution of L GH depicted in figure 7a shows an asymptotic behavior limiting it to a maximum value independently from the beam dimension.
This behavior is in great contradiction with the results obtained by Konopsky in ref. [28] 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 wonderful 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 Q−factor of 1856 occurs.
To decide the issue, we make use of the analytical expression of Eq  For the propagation length (P L), 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 apparition of a dip in the calculated reflectance. Experimentally, this dip can bring all we need to determine the BSW propagation length P L due to the fact that ∆θ R ≈ ∆θ T even in the case of absorption. As determined from Fig. 2, we get a propagation length of P L ≈ 1.374 mm independently of the beam-waist value while Eq. 4 leads to an almost twice smaller value of P L = 769 µm. On notice that expression of LP given by Eq. 4 is commonly used to interpret or exploit experimental results [20,22,23]. Again, we are in front of a contradiction that needs to be clarified. For this purpose, 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 is the amplitude expression is e − ∆kx 2 (erf (∞) → 1) that corresponds to the electric field behavior far from its maximum. This allowed expressing the propagation length as: This is a 6 × 6 matrix linking the amplitudes in media m (incident ↑ and reflected ↓ on the interface separating the two m and m + 1 media) to the same field components in the m+1 media (see figure S1). In the general case, the matrix dimension is 6 × 6 and cannot be reduced to smaller dimension when we attempt to use it in the case of a Gaussian beam even if the later is T E or T M polarized.
This point was discussed in the paper. Nevertheless, it can be reduced to 2 × 2 matrix in the T E, T M frame (see below).
In the Cartesian frame, the M m,m+1 matrix is given by: Top layer x z y Figure S1 | Schematic of the stratified structure built upon N stacks of bilayers and terminated by the top layer where the BSW mode is supposed to take place. The total number of interfaces is then 2N + 2 separating 2N + 3 media including the substrate (media m = 1) and the superstrate or the transmission media. Each layer is characterized by its dielectric permittivity ε m and its thickness d m .
where w m is the orthogonal (here z) component of the wave-vector given by the dispersion relation w m = ε m − k 2 x − k 2 y , A and B correspond to "upward" and "downward" propagation operator given by A = e iwmdm , B = e −iwmdm and d m is the thickness of the layer m. Let us emphasize that w m can be real (propagating waves) or imaginary (evanescent waves), which implies that B is not the complex conjugate of A. The total transfer matrix T T is obtained by calculating the matrix product of all the individual matrices respecting the order imposed by the light propagation direction. In the case of 2N + 2 interfaces, its expression is : Thereby, the electric field in the transmission media is simply expressed as a function of the field in the substrate through: The sub-blocks T T ij are 3 × 3 upper triangular matrices meaning that they are invertible. Both reflected and transmitted fields are then expressed in terms of the incident one by: The knowledge of the three components of the incident field makes it easy to calculate the transmitted and reflected fields.
2 BSW and the other modes of the structure This property could be of interest for the enhancement of non-linear effects.

Plane Wave Expansion of a Gaussian beam
In addition to the T MM, the Plane Wave Expansion P W E (or angular spectrum) is used to take into account the finite size of the illumination (Gaussian beam) through its angular spectrum. Let us recall the plane wave expansion of a polarized Gaussian beam when expressed in the Cartesian frame Oxyz related to the structure (Oxy being the plane parallel to the interfaces as shown in figure  S1): ] is the amplitude of the plane wave characterized by its (u, v) transverse wave-vector components expressed in the proper frame of the Gaussian beam. k x , k y and k z are the same components in the Oxyz frame, I 0 is the electric intensity of the whole Gaussian beam and W 0 its beam-waist defined as the half width at 1/e of its maximum amplitude. α and β are given by: The polarization state of the whole Gaussian beam is defined by the couple (χ, a) as: • linear directed along the angle χ measured from the x-axis with a = 1 (χ = π/2 for T E and χ = 0 for T M) • circular with χ = π/4 and a = ± √ −1 • elliptical with major to minor axes ratio equal to |β/α| = tanχ and a = ± √ −1.
Let us emphasize the fact that the combination of the T MM and the P W E can be easily extended to any system of coordinates. Here, all numerical simulations are done in Cartesian coordinates but it is obviously possible to work in the T E, T M frame which is often the case in the literature. Nevertheless, the plane wave expansion of the incident beam must be adapted by electric field projection along these two axis having unit vectors: where e z is the unit vector along Oz and k // is the tangential (in xOy plane) component of the wave-vector. By projection of the electric field given in Eqs E.7 on these two axes, one gets: As expected, a T E-polarized Gaussian beam (α = 0, β = 1) exhibits both T E and T M components and likewise for T M-polarized beam (α = 1, β = 0).
In the same basis, the T MM is reduced to a diagonal 2 × 2 matrix meaning independent relations for T E and T M components. The transmission and reflection coefficients can be than derived from a matricial method [42] or by applying an algorithm based on the use of Einstein's addition law [43].

Analytical expression of the transmitted field
In a general case, the three electric field components are calculated by integrating over all the electric fields resulting (in transmission, reflection or inside any layer) 33 from each individual incident plane wave (k x , k y ) through the expression: For the transmitted electric field, E p (k x , k y ) is replaced by t p (k x , k y )E inc p (k x , k y ) in Eq. E.13. t p (k x , k y ) is the transmission coefficient of the entire structure calculated by the T MM for the incidence given by the two components k x , k y . In the case of a BSW excitation, the transmission coefficient spectrum exhibits an almost Lorentizan shape (similar to the transfer function of a Band-pass filter) leading to express it as: Where k BSW x is the x−component of the wave-vector associated to the BSW excitation (k BSW x = 2πn 1 λ sin θ m ). Figure S3a shows the amplitude and phase of the transmission coefficient calculated by the T MM method (solid lines) and by Eq.E.14 (stars) in the case of N = 7 structure. A very good agreement is obtained between T MM and Eq. E.14. We have verified that this analytical model still valid if we increase the N value and also in the case of BSW excited in T M polarization. Expression in Eq.E.14 can be extended to more than one Lorentzian resonance by summing the corresponding functions.
As clearly shown, the transmission coefficient is faithfully described by one complex Lorentizan function. This allows us simplifying the calculation of the integral in Eq. E.13 that becomes: Accordingly, the transmitted field in the direct space (Oxyz) is the convolution 35 product of two Fourier transform functions as follows: To go further, we need to express the FT of t T E that is given by: For the incident beam, the electric field in the direct space can be expressed through a simple variable change corresponding to a rotation operation (see Eq. 11 of ref. [44]) by: By injecting Eqs E.17 and E.18 into Eq. E.16 and after some tricky algebra, we obtain: Where erf is the error function defined by erf (x) = 2 √ π x 0 e −x 2 dx. where ∆θ T,R is the FWHM of BSW resonance peak/dip appearing in the square modulus of the transmission/reflection coefficient spectra.

L GH calculation
Meanwhile, the Goos-Hänchen shift corresponds to the value of x = L GH for which the derivative with respect to x of |E t (x, y, z)| 2 vanishes. This equation can be written as: