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

TL;DR

This paper develops a new formalism for accurately analyzing Bloch surface waves excited by 3D Gaussian beams, correcting previous approximations and enabling better interpretation of experiments.

Contribution

It introduces a generalized formula for propagation length and shows the Goos-H"anchen shift depends on beam dimensions, enhancing understanding of Bloch surface wave behavior.

Findings

New formula for propagation length differs from literature

Goos-H"anchen shift depends on beam size and asymptotically limits

Predicts effects of slight deviations in incidence angle or beam position

Abstract

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.

Equations

$$Q=\frac{{\theta }_{BSW}}{\Delta {\theta }_T}$$

$$t(k_x,k_y)=\frac{t_{max}}{1+\frac{2i}{\Delta k_x}(k_x-k_x^{BSW})}\cdot 1_{k_y}$$

$$E_t(x,y,z=0)$$

$$PL=\frac{\lambda }{\pi \Delta {\theta }_R}\text{vspace}0.cm\text{,}{L}_{GH}=\frac{-\lambda }{2\pi }\frac{\partial \varphi }{\partial \theta }$$

$$\frac{W_0\sqrt{\pi }\Delta k}{4\cos {\theta }_m}\left[erf\left(\frac{\cos {\theta }_m{L}_{GH}}{W_0}-\frac{W_0\Delta k}{4\cos {\theta }_m}\right)+1\right]-e^{-{\left(\frac{\cos {\theta }_m{L}_{GH}}{W_0}-\frac{W_0\Delta k}{4\cos {\theta }_m}\right)}^2}=0$$

$$LP=\frac{2}{\Delta k_x}=\frac{\lambda }{\pi n_{1}cos{\theta }_m\Delta {\theta }_{T,R}}$$

$$\centering\left\{\begin{array}[]{c}E_{m+1~{}x}^{\uparrow}\\ E_{m+1~{}y}^{\uparrow}\\ E_{m+1~{}z}^{\uparrow}\\ E_{m+1~{}x}^{\downarrow}\\ E_{m+1~{}y}^{\downarrow}\\ E_{m+1~{}z}^{\downarrow}\\ \end{array}\right\}=M_{m,m+1}\times\left\{\begin{array}[]{c}E_{m~{}x}^{\uparrow}\\ E_{m~{}y}^{\uparrow}\\ E_{m~{}z}^{\uparrow}\\ E_{m~{}x}^{\downarrow}\\ E_{m~{}y}^{\downarrow}\\ E_{m~{}z}^{\downarrow}\\ \end{array}\right\}\@add@centering$$

$$\left(\begin{array}[]{cccccc}\cfrac{w_{m+1}+w_{m}}{2w_{m+1}}A&0&\cfrac{-k_{x}(\varepsilon_{m+1}-\varepsilon_{m})}{2w_{m+1}\varepsilon_{m+1}}A&\cfrac{w_{m+1}-w_{m}}{2w_{m+1}}B&0&\cfrac{-k_{x}(\varepsilon_{m+1}-\varepsilon_{m})}{2w_{m+1}\varepsilon_{m+1}}B\\ 0&\cfrac{w_{m+1}+w_{m}}{2w_{m+1}}A&\cfrac{-k_{y}(\varepsilon_{m+1}-\varepsilon_{m})}{2w_{m+1}\varepsilon_{m+1}}A&0&\cfrac{w_{m+1}-w_{m}}{2w_{m+1}}B&\cfrac{-k_{y}(\varepsilon_{m+1}-\varepsilon_{m})}{2w_{m+1}\varepsilon_{m+1}}B\\ 0&0&\cfrac{\varepsilon_{m}w_{m+1}+\varepsilon_{m+1}w_{m}}{2w_{m+1}\varepsilon_{m+1}}A&0&0&\cfrac{\varepsilon_{m}w_{m+1}-\varepsilon_{m+1}w_{m}}{2w_{m+1}\varepsilon_{m+1}}B\\ \cfrac{w_{m+1}-w_{m}}{2w_{m+1}}A&0&\cfrac{k_{x}(\varepsilon_{m+1}-\varepsilon_{m})}{2w_{m+1}\varepsilon_{m+1}}A&\cfrac{w_{m+1}+w_{m}}{2w_{m+1}}B&0&\cfrac{k_{x}(\varepsilon_{m+1}-\varepsilon_{m})}{2w_{m+1}\varepsilon_{m+1}}B\\ 0&\cfrac{w_{m+1}-w_{m}}{2w_{m+1}}A&\cfrac{k_{y}(\varepsilon_{m+1}-\varepsilon_{m})}{2w_{m+1}\varepsilon_{m+1}}A&0&\cfrac{w_{m+1}+w_{m}}{2w_{m+1}}B&\cfrac{k_{y}(\varepsilon_{m+1}-\varepsilon_{m})}{2w_{m+1}\varepsilon_{m+1}}B\\ 0&0&\cfrac{\varepsilon_{m}w_{m+1}-\varepsilon_{m+1}w_{m}}{2w_{m+1}\varepsilon_{m+1}}A&0&0&\cfrac{\varepsilon_{m}w_{m+1}+\varepsilon_{m+1}w_{m}}{2w_{m+1}\varepsilon_{m+1}}B\\ \end{array}\right)$$

$$TT=\prod _{m=2N+1}^1{M}_{m,m+1}=M_{2N+1,2N+2}\times M_{2N,2N+1}\cdots \times M_{2,3}\times M_{1,2}$$

$$\centering\left\{\begin{array}[]{c}E_{t~{}x}\\ E_{t~{}y}\\ E_{t~{}z}\\ 0\\ 0\\ 0\\ \end{array}\right\}=TT\times\left\{\begin{array}[]{c}E_{i~{}x}\\ E_{i~{}y}\\ E_{i~{}z}\\ E_{r~{}x}\\ E_{r~{}y}\\ E_{r~{}z}\\ \end{array}\right\}=\left\{\begin{array}[]{c c}TT_{11}&TT_{12}\\ TT_{21}&TT_{22}\\ \end{array}\right\}\times\left\{\begin{array}[]{c}E_{i~{}x}\\ E_{i~{}y}\\ E_{i~{}z}\\ E_{r~{}x}\\ E_{r~{}y}\\ E_{r~{}z}\\ \end{array}\right\}\@add@centering$$

$${\vec{E}}_r(k_x,k_y)$$

$${\vec{E}}_t(k_x,k_y)$$

$$E_x(k_x,k_y)$$

$$E_y(k_x,k_y)$$

$$E_z(k_x,k_y)$$

$$\alpha$$

$$\beta$$

$${\vec{e}}_{TE}$$

$${\vec{e}}_{TM}$$

$$E_{TE}(k_x,k_y)$$

$$E_{TM}(k_x,k_y)$$

$$E_p(x,y,z)=\iint E_p^R(k_x,k_y)e^{i(k_{x}x+k_{y}y+\sqrt{{\epsilon }_m\frac{{\omega }^2}{c^2}-k_x^2-k_y^2}z)}d{k}_{x}d{k}_y$$

$$t(k_x,k_y)=\frac{t_{max}}{1+\frac{2i}{\Delta k_x}(k_x-k_x^{BSW})}$$

$${\vec{E}}_t(x,y,z)={\iint }_{-\infty }^{\infty }\tilde{t}(k_x,k_y)\cdot {\vec{E}}_{inc}(k_x,k_y)e^{i(k_{x}x+k_{y}y+\sqrt{{\epsilon }_m\frac{{\omega }^2}{c^2}-k_x^2-k_y^2}d{k}_{x}d{k}_y}$$

$${\vec{E}}_t(x,y,z)=FT\{\tilde{t}(k_x,k_y)\}\otimes FT\{{\vec{E}}_{inc}(k_x,k_y)\}$$

$$t_{TE}(x,y,z=z_t)=FT\{t_{TE}\}=\frac{2}{\Delta k_x}H(x)e^{-\frac{\Delta k_{x}x}{2}}{e}^{i{k}_x^{BSW}x}$$

$$E_{TE}(x,y,z=0)=\frac{1}{\sqrt{\pi }{W}_0}{e}^{-\frac{x^2{\cos }^2{\theta }_m}{W_0^2}+\frac{y^2}{W_0^2}+i{k}_x^{BSW}x}$$

$$E_t(x,y,z=0)$$

$$LP=\frac{2}{\Delta k_x}=\frac{\lambda }{\pi n_{1}cos{\theta }_m\Delta {\theta }_{T,R}}$$

$$\frac{W_0\sqrt{\pi }}{2LP\cos {\theta }_m}\left[erf\left(\frac{\cos {\theta }_m{L}_{GH}}{W_0}-\frac{W_0}{2LP\cos {\theta }_m}\right)+1\right]-e^{-{\left(\frac{\cos {\theta }_m{L}_{GH}}{W_0}-\frac{W_0}{2LP\cos {\theta }_m}\right)}^2}=0$$

Full text

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 $TE$ and $TM$ 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 $PL$ 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 $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, the two methods converge to the same result that fails the commonly used formulas. For the $PL$, 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 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 $TMM$ (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 $TE$ polarization. All geometrical parameters are given in the caption the figure 1.

Figure 1b shows the square modulus of the transmitted electric field amplitude (in logarithmic scale) at the upper interface as a function of the bi-layer number ($N$) at $\lambda=1550~{}nm$. The angle of incidence and the natural logarithm of the $FWHM$ (Full Width at Half Maximum) $\Delta\theta_{T}$ of the $BSW$ resonance are given on figure 1c. As expected, the angular position $\theta_{BSW}$ converges asymptotically but promptly ($N\simeq 7$) to the value corresponding to the pole of the transmission coefficient of the infinite structure. $\Delta\theta_{T}$ varies exponentially with $N$ (see red with stars line on figure 1c) meaning rather similar variations for the $BSW$ resonance quality factor defined by:

[TABLE]

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]. Even if these losses are quite hard to be quantified, most of authors agreed introducing them in calculations by adding a small imaginary part to the optical index.

Unfortunately, when introducing such absorption 111One notes that absorption is fundamentally proportional to the imaginary part of the permittivity and not to that of the optical index meaning that $\theta_{SBW}$ will be affected. through $n^{\prime\prime}$ for all media (except glass substrate), the $BSW$ angular position remains unchanged while the $BSW$ efficiency becomes weaker. In our case, we estimate that $BSW$ excitation is completely canceled for $n^{\prime\prime}>10^{-3}$. Figure 1d shows the quality factor $Q$ variations versus the number of bi-layers $N$ for different values of the imaginary part $n^{\prime\prime}$. For loss-less materials, $Q$ tends to infinity as $Q=e^{0.8623N+0.0586}$ while asymptotic behaviors occur for $n^{\prime\prime}\neq 0$. We have verified that the numerical values appearing in the last relation only depends on the effective index associated to the $BSW$ excitation (here $n_{eff}=n_{1}\sin{\theta_{BSW}}=1.3928$). More importantly, to go further through analytical calculation, it is essential that the transmission coefficient be expressed explicitly as a function of the wave-vector components. Fortunately, in the case of a $BSW$ excitation, the transmission coefficient spectrum can be approached very realistically by a Lorentizan function (see more details in the SI file) leading to express it as:

[TABLE]

Where $k_{x}^{BSW}$ 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_{x}^{BSW}$ calculated through the $TMM$.

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 $PWE$ (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 calculated in the direct space $Oxyz$ through the equation E.15 of the SI file. The latter involves the transmission Jones matrix of the structure that is basically given by the $TMM$ as $\tilde{t}(k_{x},k_{y})=-TT_{21}^{-1}\times TT_{22}$ (Eq. E.4 of the SI file). 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 only considering the resonant term of the transmission (for instance the $TE$ term in our case) and replacing the transmission coefficient by its expression given by Eq. 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 $FWHM$ ($\Delta k$) of the transmission coefficient through:

[TABLE]

Where $erf$ is the error function defined by $erf(x)=\frac{2}{\sqrt{\pi}}\displaystyle\int^{x}_{0}{e^{-x^{2}}}dx$ and $\Delta k=2\pi n_{1}\cos\theta_{m}\Delta\theta_{T}$ is the $FWHM$ of the transmittance spectrum as defined above.

Equation 3 provides determining all the $BSW$ properties ($PL$, $L_{GH}$, maximum efficiency…) as it will be discussed below.

Results and discussion

Within $TMM/PWE$ combination (Eq. E.13 of the SI file) one can calculate the electric field distribution over all the structure for any illumination direction, beam-waist 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~{}\mu m$ in Fig. 1a, $W_{0}=30~{}\mu 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 $PL$ 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 as for surface plasmon resonance or $BSW$ [31, 23]. 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.

Figure 3a shows the 3D map of the $BSW$ electric near-field intensity distribution at $z=0$ from 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 distance ($Ox$ here). Top-view distributions are given in Fig. 3b allowing identifying the excitation of the $BSW$ (comet shape) in only $TE$ polarization and highlighting the lateral shift that accompanies the excitation of the $BSW$.

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~{}\mu 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 Gaussian beam with $W_{0}=5~{}\mu m$. The first one (Fig. 4a) corresponds to a beam-waist located $300~{}\mu m$ under the PhC structure while it is supposed to be $100~{}\mu m$ above the substrate-PhC interface in the second (Fig. 4b). As in Fig. 2, the calculated amplitude of the electric field in Fig. 4 is 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 the figure (b)) while additional lateral shift of this maximum occurs in the second case (solid black line). This demonstrates how the $BSW$ shape can be affected by a lack of beam focus. In addition, another effect arises on the interference pattern appearing in the reflected beam due to the spatial broadening of the beam falling the PhC. 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 hurts 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 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~{}\mu 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 $TMM/PWE$ 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.

Figures 6a and 6b show a transverse cross-section (along the $Ox$ axis) made $1~{}mm$ under the first interface (substrate-PhC) when the beam-waist varies from $W_{0}=5~{}\mu m$ to $W_{0}=50~{}\mu m$ for a 3D and 2D Gaussian beams respectively. 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~{}\mu 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 $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 amplitude shape (see Fig. 2c of ref. [33]) for these components. This is explicitly written in Eqs E.6 of the SI file where the $x-$ and $z-$components of the electric field are not zero even in $TE$ polarization ($\chi=\pi/2\rightarrow\alpha=0$ and $\beta=1$). The contribution of these $TM$ components to the $BSW$ (in transmission) is negligible because only $TE$ components resonate but this does not prevent their contribution to the reflected beam.

This interpretation is derived from the Maxwell-Gauss equation ($Div{\vec{E}}=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 wavevector that is not located in the plane of incidence. This is true for both $TE$ and $TM$ polarized beams as shown by Eqs. E.11 and E.12 in the SI file 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 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 $PL$ for the three beam-waist values to be: $\{L_{GH}=49.85~{}\mu m,LP=1.3736~{}mm\}$ for $W_{0}=10~{}\mu m$, $\{L_{GH}=124.34~{}\mu m,LP=1.3740~{}mm\}$ for $W_{0}=30~{}\mu m$ and $\{L_{GH}=1.07~{}mm,LP=1.3745~{}mm\}$ for $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 [25, 24] estimating these two quantities to be:

[TABLE]

where $\Delta\theta_{R}$ is the $FWHM$ of the dip resonance appearing in the reflectance spectrum and $\phi$ 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 far-field detection of the reflected beam. Only direct measurement of the near-field allows access to this property. Theoretically, the variation of $\phi$ versus the angle of incidence is given in figure S3 of the SI file. 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~{}\mu m$), which cannot be consistent with the calculated values from Fig. 2 that depend on the beam dimension. This discrepancy needs to be elucidated.

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 $TMM/PWE$ algorithm. Figure 7a shows that $L_{GH}$ significantly varies with $W_{0}$ as long as the angular width of the beam ($\frac{\lambda}{\pi W_{0}}$) is $22$ times larger than the angular width of the $BSW$ (here $\Delta\theta_{T}=0.642~{}mrad$) corresponding to $W_{0}\simeq 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. 3. The spatial position of the transmitted intensity maximum corresponds to the value of $x$ for which the $x-$derivative of the square modulus of the electric field amplitude given by Eq. 3 vanishes. This condition leads to Eq. 5 that was numerically revolved (see the SI file) for the three cases of Fig. 7a. The obtained values are depicted by circles on the same figure showing a perfect agreement with the $TMM/PWE$. We have verified the absolute error to be less than $5\times 10^{-3}$.

[TABLE]

For the propagation length ($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 apparition of a dip in the calculated reflectance. Experimentally, this dip can bring all we need to determine the $BSW$ propagation length $PL$ due to the fact that $\Delta\theta_{R}\approx\Delta\theta_{T}$ even in the case of absorption. As determined from Fig. 2, we get a propagation length of $PL\approx 1.374~{}mm$ independently of the beam-waist value while Eq. 4 leads to an almost twice smaller value of $PL=769~{}\mu 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\rightarrow\infty(x>>W_{0})$, the predominant term is the amplitude expression is $e^{-\frac{\Delta kx}{2}}$ ($erf(\infty)\rightarrow 1$) that corresponds to the electric field behavior far from its maximum. This allowed expressing the propagation length as:

[TABLE]

Replacing $\theta_{m}=1.189~{}rad$, $n_{1}=1.501$ and $\Delta\theta_{T}=0.642~{}mrad$ into Eq.6 leads to $LP=1.376~{}mm$ which perfectly agrees the estimated value ($1.374~{}mm$) by $TMM/PWE$ algorithm. We have verified the good agreement between the $TMM/PWE$ results (solid line) and this analytical formula (green circles) on figure 7b. 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. Conclusion

Combining the $PWE$ with the $TMM$, and using an accurate angular spectrum expansion of a Gaussian beam, turns 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 ($RCWA$ for instance) able to take into account diffraction by grating (periodic) or by individual pattern such as in [4, 13, 40, 41]. 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, or focusing, presence of adhesion layer on the top surface, …) on the excitation of the surface wave. The major result of this paper is obtained through analytical 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.

Acknowledgements

Several computations have been performed on the supercomputer facilities of the "Mésocentre de calcul de Franche-Comté".

Supplementary Information file

All the numerical simulations were done using codes developed under Matlab environment. These codes combine the $TMM$ and the $PWE$ in Cartesian coordinates. They allows simulations with an arbitrary number of planar membrane stack. The $TMM$ was used alone in the first part of the paper to correctly design a $BSW$ excitation at $\lambda=1550~{}nm$ as shown on figure 2 of the paper.

1 The $TMM$ method

Let us first recall the principle of the $TMM$. For a given stratified structure of $N$ stacks of bilayers (see schematic of figure S1), let $M_{m,m+1}$ the transfer matrix linking the electric field components in the medium $m+1$ to those of the medium $m$ in Cartesian coordinates. $M$ is then given by:

[TABLE]

This is a $6\times 6$ matrix linking the amplitudes in media $m$ (incident $\uparrow$ and reflected $\downarrow$ 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\times 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 $TE$ or $TM$ polarized. This point was discussed in the paper. Nevertheless, it can be reduced to $2\times 2$ matrix in the $TE,TM$ frame (see below).

In the Cartesian frame, the $M_{m,m+1}$ matrix is given by:

$M_{m,m+1}(k_{x},k_{y})=$

[TABLE]

where $w_{m}$ is the orthogonal (here $z$) component of the wave-vector given by the dispersion relation $w_{m}=\sqrt{\varepsilon_{m}-k_{x}^{2}-k_{y}^{2}}$, $A$ and $B$ correspond to "upward" and "downward" propagation operator given by $A=e^{iw_{m}d_{m}}$, $B=e^{-iw_{m}d_{m}}$ 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 $TT$ 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 :

[TABLE]

Thereby, the electric field in the transmission media is simply expressed as a function of the field in the substrate through:

[TABLE]

The sub-blocks $TT_{ij}$ are $3\times 3$ upper triangular matrices meaning that they are invertible. Both reflected and transmitted fields are then expressed in terms of the incident one by:

[TABLE]

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

As mentioned in the paper and according to figure 1b, supplementary Bloch modes exist at smaller values of the angle of incidence (smaller effective index). They correspond to the cavity modes of the stratified finite structure. Peaks $H_{1}$ to $H_{3}$ of figure 1b are the three first harmonics of the total cavity formed by the structure. In addition to the $BSW$ mode presented on figure S2a, we show the normalized (to the incident) electric field intensity distribution inside the structure for the three angles of incidence (b-d) in the case of $N=20$-structure illuminated by a plane wave in $TE$ polarization. Light confinement occurs in all cases. As expected, the electric field is localized in the upper layer and leads to large electric field confinement at the upper interface (evanescent wave) when the $BSW$ is excited. This property could be of interest for the enhancement of non-linear effects.

3 Plane Wave Expansion of a Gaussian beam

In addition to the $TMM$, the Plane Wave Expansion $PWE$ (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):

[TABLE]

Where $E_{0}(u,v)=\sqrt{\frac{I_{0}}{2\pi}}W_{0}\exp\left[-W_{0}^{2}(u^{2}+v^{2})/4\right]$ 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. $\alpha$ and $\beta$ are given by:

[TABLE]

The polarization state of the whole Gaussian beam is defined by the couple ($\chi$, $a$) as:

• •

linear directed along the angle $\chi$ measured from the $x$-axis with $a=1$ ($\chi=\pi/2$ for $TE$ and $\chi=0$ for $TM$)

• •

circular with $\chi=\pi/4$ and $a=\pm\sqrt{-1}$

• •

elliptical with major to minor axes ratio equal to $|\beta/\alpha|=tan\chi$ and $a=\pm\sqrt{-1}$.

Let us emphasize the fact that the combination of the $TMM$ and the $PWE$ 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 $TE,TM$ 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:

[TABLE]

where $\vec{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 3 on these two axes, one gets:

[TABLE]

As expected, a $TE$-polarized Gaussian beam ($\alpha=0,\beta=1$) exhibits both $TE$ and $TM$ components and likewise for $TM$-polarized beam ($\alpha=1,\beta=0$).

In the same basis, the $TMM$ is reduced to a diagonal $2\times 2$ matrix meaning independent relations for $TE$ and $TM$ 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].

4 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) from each individual incident plane wave ($k_{x},k_{y}$) through the expression:

[TABLE]

with $p=\{x,y,z\}$ or $\{TE,TM\}$, $\varepsilon_{m}$ is the dielectric permittivity of media $m$ where the field is calculated, $\omega$ is the angular frequency (pulsation) of the incident light and $c$ is the light velocity in vacuum. The component $E_{p}^{R}(k_{x},k_{y})$ contains the response of the structure to the illumination by the spectral component $(k_{x},k_{y})$ calculated through the $TMM$ as described previously. It could correspond to the transmitted, reflected or any field inside the structure in medium $m$.

Rigorously, the integral in equation E.13 ranges from $-\infty$ to $+\infty$. Practically, a truncation is numerically introduced by considering only $N_{x}$ values of $k_{x}$ and $N_{y}$ values of $k_{y}$ in such a way as to correctly describe both the Gaussian shape of the beam and the structure response (resonant modes for instance). Therefore, on the one hand, we consider all plane waves with angular deviation $\delta\theta=\theta_{BSW}-\theta_{i}\in[-\cfrac{5\lambda}{\pi W_{0}};+\cfrac{5\lambda}{\pi W_{0}}]$ which allows taking into account all the angular components of the incident beam with amplitude larger than $\sqrt{I_{0}}\times 10^{-2}$. On the other hand, the values of $N_{x}$ and $N_{y}$ must be as large as possible to prevent any aliasing effect especially if large spatial window are considered. For example, values of $(N_{x},N_{y})=(6000,350)$ are usually used to calculated the spatial distribution of the electric field over a window of $\Delta x\times\Delta y=150W_{0}\times 10W_{0}$. In all cases, convergence tests are systematically applied.

For the transmitted electric field, $E_{p}(k_{x},k_{y})$ is replaced by $t_{p}(k_{x},k_{y})E_{p}^{inc}(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 $TMM$ 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:

[TABLE]

Where $k_{x}^{BSW}$ is the $x-$component of the wave-vector associated to the $BSW$ excitation ($k_{x}^{BSW}=\frac{2\pi n_{1}}{\lambda}\sin\theta_{m}$). Figure S3a shows the amplitude and phase of the transmission coefficient calculated by the $TMM$ method (solid lines) and by Eq.E.14 (stars) in the case of $N=7$ structure. A very good agreement is obtained between $TMM$ 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 $TM$ 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:

[TABLE]

Where $\tilde{t}(k_{x},k_{y})$ is the $3\times 3$ transmission Jones matrix expressed in the $Oxyz$ frame linking the transmission electric field to the incident one. In general, it is given by Eq.E.5 (term between brackets) but in the case of $BSW$excitation, Eq. E.15 turns into a scalar equation only involving the resonant electric field component (here $TE$ component) and the transmission coefficient given by Eq. E.14. Let us emphasize that this assumption is only valid for the calculation of the transmitted field and cannot be used for the reflected beam because all the three (in Cartesian) or the two (in $TE,TM$ basis) components of the incident electric field undergo the same reflection coefficient in amplitude (in the case of loss-less materials).

Accordingly, the transmitted field in the direct space ($Oxyz$) is the convolution product of two Fourier transform functions as follows:

[TABLE]

To go further, we need to express the FT of $t_{TE}$ that is given by:

[TABLE]

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:

[TABLE]

By injecting Eqs E.17 and E.18 into Eq. E.16 and after some tricky algebra, we obtain:

[TABLE]

Where $erf$ is the error function defined by $erf(x)=\frac{2}{\sqrt{\pi}}\displaystyle\int^{x}_{0}{e^{-x^{2}}}dx$.

5 $LP$ analytical expression

From the last equation, we can clearly see that for $x\rightarrow\infty$, the predominant term is $e^{-\cfrac{\Delta kx}{2}}$ ($erf(\infty)\rightarrow 1$) that corresponds to the electric field behavior far from its maximum. This allowed expressing the propagation length as:

[TABLE]

where $\Delta\theta_{T,R}$ is the FWHM of $BSW$ resonance peak/dip appearing in the square modulus of the transmission/reflection coefficient spectra.

6 $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:

[TABLE]

Unfortunately, this equation cannot be solved analytically. A Root-finding algorithm based on Newton’s method is then developed and used to calculate the solution of Eq. E.21. Nevertheless, we verified that the asymptotic value of $L_{GH}$ corresponds to the $LP$ value. This still valid for any configuration with or without absorption as shown on figure S3b. In the case of $N=7$, the $L_{GH}$ variations versus the incident beam-waist are presented on figure S3b (solid red line) in comparison with the values obtained through the combined $TMM/PWE$ method (blue circles). We can see a very good agreement between the two curves confirming the extrinsic character of this shift.

All figures (curves and maps) are obtained within Matlab and presentations are then improved with CorelDRAW. Maple 16.00 version was also used to verify Eq.E.19.