Introduction

Quantification 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 sensing1,2,3 to surface-optics4,5,6,7,8,9,10,11 or micro-manipulation12,13. BSW is of a great interest in integrated optics14,15 due to its very large (>3 mm16) 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 diffraction19,20. However, 3D BSW electromagnetic-field distribution has never been theoretically reported except very recently by pure numerical methods (Finite Difference Time Domain8 or Finite elements21). 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 (LGH), 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 LGH, 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 studies16,22,23,24 used theoretical formulas based on a development obtained for plane wave illumination25,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 studies27,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 studies22,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 denominator31. 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 Lorentzian-shape 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 design22,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.

Fig. 1: Schematic of the studied 1D-photonic crystal structure and transmission properties.
figure 1

a The incident beam illuminates the structure from a glass substrate at an angle θm. It is linearly polarized along the y-direction (Transverse Electric polarization). The bilayer stack is composed of a layer of high-index media (nh = 2.23) with a thickness dh =  294 nm deposited on a second layer made in low-index material (nl = 1.75) of thickness dl = 240 nm. The total number of stacks is named N and the structure is terminated by a top layer of high-index material of dtop = 5503 thickness. b Transmitted electric intensity in logarithmic scale versus the number of bilayers N. Additional modes H1, H2, and H3 occur when N increases corresponding to smaller values of the angle of incidence. c Variations of θBSW as a function of the number of bilayers (N) in blue solid line and its Full Width to Half Maximum (FWHM) variation Δθ (in log-scale) in red dashed line. d Variations of the quality factor Q of the Bloch Surface Wave excitation with the bilayer number N for different values of the optical refractive index imaginary part n″.

Full size image

Figure 1b shows the square modulus of the transmitted electric-field 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 H1 to H3 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:

$$Q=\frac{{\theta }_{{\rm{BSW}}}}{\Delta {\theta }_{{\rm{T}}}}$$
(1)

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 increases33 (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 = e0.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}_{{\rm{eff}}}={n}_{{\rm{1}}}\sin {\theta }_{{\rm{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:

$$t({k}_{x},{k}_{y})=\frac{{t}_{\max }}{1+\frac{2i}{\Delta k}({k}_{x}-{k}_{x}^{{\rm{BSW}}})}\cdot {1}_{{k}_{y}}$$
(2)

where \({k}_{x}^{{\rm{BSW}}}\) 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}_{x}^{{\rm{BSW}}}\) calculated through the TMM. Δk is the FWMH of the square modulus of the reflection or the transmission coefficient given by: \(\Delta k=\frac{2\pi }{\lambda }{n}_{{\rm{1}}}\cos {\theta }_{{\rm{m}}}\Delta {\theta }_{{\rm{T/R}}}\).

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 199934 and verified by comparison with experimental and/or results based on different methods34,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 as \(\tilde{t}({k}_{x},{k}_{y})=-T{T}_{{\rm{21}}}^{-1}\times T{T}_{{\rm{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 W0 and the FWHM (Δk) of the transmission coefficient through:

$${E}_{{\rm{t}}}(x,y,z=0) = \,\, \frac{\sqrt{{I}_{{\rm{0}}}}{t}_{\max }\Delta k}{4\cos {\theta }_{{\rm{m}}}}{e}^{-\frac{8\Delta k\cos {\theta }_{{\rm{m}}}^{2}x-{W}_{{\rm{0}}}^{2}{(\Delta k)}^{2}}{16\cos {\theta }_{{\rm{m}}}^{2}}}\\ \times \left[{\rm{erf}}\left(\frac{4\cos {\theta }_{{\rm{m}}}^{2}\ x-{W}_{{\rm{0}}}^{2}\Delta k}{4{W}_{{\rm{0}}}\cos {\theta }_{{\rm{m}}}}\right) +1\right]{e}^{-\frac{{y}^{2}}{{W}_{{\rm{0}}}^{2}}}{e}^{-i{k}_{x}^{{\rm{BSW}}}x}$$
(3)

Where erf is the error function defined by \(erf(x)=\frac{2}{\sqrt{\pi }}\mathop{\int}\nolimits_{0}^{x}{e}^{-{x}^{2}}dx\) and Δk is the FWHM of the transmittance spectrum defined before.

Equation  (3) provides all the BSW properties (PL, LGH, 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, W0 = 10 μm in Fig. 2a, W0 = 30 μm in Fig. 2b and W0 = 1 mm in Fig. 2c. The spatial shape of the excited BSW greatly depends on the value of the incident beam-waist W0. 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 W0 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.

Fig. 2: Electric-field amplitude distributions for three different beam-waist values.
figure 2

In the three cases, the distributions are calculated in a vertical plane passing by the center of the incident beam in the case of a structure composed of N = 7 bilayer stacks. The beam-waist of the incident beam is set as W0 = 10 μm in a, b, c, W0 = 30 μm in d, e, f, and W0 = 1 mm in g, h, i and it is Transverse Electric (TE)-polarized. Maps in (a, d, g) (b, e, h) and (c, f, i) correspond to the electric-field amplitude in the superstrate, the bilayer stacks and the substrate, respectively. For the sake of clarity, the electric field was auto-normalized in each zone so that the maximum value of the scale bar varies from one to another map allowing the visualization of both incident and reflected beams.

Full size image

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 shift37, or the angular shift of the secondary reflected beam38. 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 coefficient39. Generally, they occur only with circular or elliptical incident polarized beams. Furthermore, two different definitions are still used for the LGH, which can be considered as the displacement of the maximum of the intensity or, that of the intensity centroid40. 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 BSW16,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 optics41. 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.

Fig. 3: 3D map of the electric-field intensity distribution.
figure 3

In all figures, this quantity is calculated at the top interface (z = 0 nm) for θ = θBSW, N = 7 and W0 = 300 μm. a is a 3D map of the electric-intensity distribution of the Bloch Surface Wave (BSW). b Top-view maps of the electric intensity for the two polarizations (Transverse Electric TE on top and Transverse Magnetic TM on bottom). c and d correspond to an angle of incidence θm = θBSW + 1∘. Experimentally, this kind of distributions is measured by means of scanning near-field optical microscope to estimate both the Goos–Hänchen shift LGH (difference between TE and TM) and the Propagation Length PL16. Note that the intensity maximum value is 70 times greater in TE than in TM.

Full size image

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 LGH = 649 μm and the PL = 1.37 mm.

Fig. 4: Intensity map distributions for defocused Gaussian beams.
figure 4

a The waist is at 300 μm under the first interface and 100 μm above this interface in b. The beam-waist is set to W0 = 5 μm and the mean angle of incidence corresponds to the Bloch Surface Wave (BSW) excitation. The dashed blue rectangle in a corresponds to a zoom-in over the region where incident and reflected beam overlap. The green dashed dotted-dotted line on Fig. c corresponds to the electric-field intensity along the Ox direction in the case of a waist perfectly centered on the first interface. The red dotted line corresponds to 20× the same quantity in Transverse Magnetic (TM) polarization (no excited BSW). The blue dashed and the black solid lines correspond to a and b cases.

Full size image

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 (W0 = 5 μm) Gaussian beam. 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. 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 excitation42. Very recently, Petrova et al.3 exploited the properties of the reflected beam for biosensing applications. Several theoretical studies have been performed in this context27,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 W0 = 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 electric-intensity 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.

Fig. 5: Electric-field amplitude distributions.
figure 5

a In Oxz plane at y = 0, b, in xy plane at z = −1 mm, and c, in ξy plane that is perpendicular to the propagation direction of the reflected beam (see the ξ axis depicted in a). We recall the geometric parameters: N = 7, a beam-waist W0 = 5 μm and a Transverse Electric (TE) polarized incident beam.

Full size image

 Figure 6a, b shows a transverse cross-section (along the Ox axis) 1 mm under the first interface (substrate-PhC) for a beam-waist varying from W0 = 5 μm to W0 = 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 W0 = 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 given by Supplementary Eq. S.6 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.

Fig. 6: Reflected beam behavior with regard to the calculation method.
figure 6

Intensity map of the reflected beam calculated along the Ox axis at 1 mm under the Photonic Crystal (PhC) as a function of the beam waist W0 when simulations are done for a realistic 3D Gaussian beam in a or a prismatic 2D-Gaussian beam in b. c Cross sections made over the two maps for W0 = 6.87 μm showing a real discrepancy between the two oscillation patterns.

Full size image

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 kx 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.

Goos–Hänchen shift and PL

The number of bilayers is fixed to N = 7 in the following as in Fig. 2 from which one determined the LGH and PL for the three beam-waist values to be: {LGH = 49.85 μm, PL = 1.3736 mm} for W0 = 10 μm, {LGH = 124.34 μm, PL =  1.3740 mm} for W0 = 30 μm and {LGH = 1.07 mm, PL = 1.3745 mm} for W0 = 1 mm. Even if the propagation length is almost constant, its value, in addition to the evolution of the LGH, is in clear contradiction with a simple theory based on plane wave analysis25,26 estimating these two quantities to be:

$${\mathrm{PL}}=\frac{\lambda }{\pi \Delta {\theta }_{{\rm{R}}}},\quad {L}_{{\rm{GH}}}=\frac{-\lambda }{2\pi }\frac{\partial \phi }{\partial \theta }$$
(4)

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 LGH 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 (LGH = 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 LGH as a function of the beam-waist value through the TMM/PWE algorithm. Figure 7a shows that LGH significantly varies with W0 as long as the angular width of the beam (\(\frac{\lambda }{\pi {W}_{{\rm{0}}}}\)) is 22 times larger than the angular width of the BSW (here ΔθT = 0.642 mrad) corresponding to W0 ≃ 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 LGH. As expected, the latter varies dramatically with losses. Second, the evolution of LGH 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 LGH 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 Q-factor of 1856 occurs.

Fig. 7: Variations of the Goos–Hänchen lateral shift.
figure 7

a, in the case of N = 7 as a function of the beam-waist for three different values of the imaginary part of the refractive index of the media (all the layers except the glass substrate and air obviously). The blue curve corresponds to the loss-less case while the red and green ones correspond to n″ = 10−4 and n″ = 10−3 respectively. The yellow, ocher and purple vertical thick lines correspond to three different values of the beam-waist W0. All the three curves in a correspond to values calculated within the Transfer Matrix Method (TMM) combined to the Plane Wave Expansion (PWE) algorithm while the circles are obtained from Eq. (5). b, Comparison between TMM/PWE results (solid line) and values calculated from Eq. (5) (green circles) of the LGH as a function of n″ in the case of a fixed value of the beam-waist (W0 = 300 μm).

Full size image

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 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 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.

$$ \frac{{W}_{{\rm{0}}}\sqrt{\pi }}{2{\mathrm{PL}}\cos {\theta }_{{\rm{m}}}}\left[{\mathrm{er}}f\left(\frac{\cos {\theta }_{{\rm{m}}}{L}_{{\rm{GH}}}}{{W}_{{\rm{0}}}}-\frac{{W}_{{\rm{0}}}}{2{\mathrm{PL}}\cos {\theta }_{{\rm{m}}}}\right)+1\right]\\ \quad-{e}^{-{\left(\frac{\cos {\theta }_{{\rm{m}}}{L}_{{\rm{GH}}}}{{W}_{{\rm{0}}}}-\frac{{W}_{{\rm{0}}}}{2{\mathrm{PL}}\cos {\theta }_{{\rm{m}}}}\right)}^{2}}=0$$
(5)

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 beam-waist 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 results16,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 ≫ W0), the predominant term in the amplitude expression is \({e}^{-\frac{\Delta kx}{2}}\) (erf(∞) → 1) that corresponds to the electric-field behavior far from its maximum. This allowed expressing the propagation length as:

$${\mathrm{PL}}=\frac{2}{\Delta k}=\frac{\lambda }{\pi {n}_{{\rm{1}}}\cos {\theta }_{{\rm{m}}}\Delta {\theta }_{{\rm{T,R}}}}$$
(6)

Replacing θm = 1.189 rad, n1 = 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 : \(\Delta {\theta }_{{\rm{R}}}^{{\rm{air}}}={n}_{{\rm{1}}}\Delta {\theta }_{{\rm{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}_{{\rm{PL}}}^{S}=\frac{1}{{n}_{{\rm{1}}}\cos {\theta }_{{\rm{m}}}}\) in the substrate or \({C}_{{\rm{PL}}}^{A}=\frac{1}{\cos {\theta }_{{\rm{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}_{{\rm{PL}}}^{S}\approx =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 \({\mathrm{PL}}_{\exp }=470\ \mu 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 \({\mathrm{PL}}=448.5/(1.501\cos 50.3{2}^{\!\deg})=468\ {\mathrm{\mu 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 CPL = 1.64, which would correspond to an error of 64%.

For the Goos–Hänchen shift LGH, the error depends on the beam-waist value (W0). In our case, this error can reach 86% for highly focused beams (W0 = 10 μm) and is only twice smaller (43%) in the case of wide beams as considered in the study by Soboleva23.

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 pattern5,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.