Simulation and analysis of exotic non-specular phenomena

F. Krayzel Clermont Université, Université Blaise Pascal, LASMEA, BP 10448, F-63000 Clermont-Ferrand CNRS, UMR6602, LASMEA, F-63177 Aubière R. Pollès Clermont Université, Université Blaise Pascal, LASMEA, BP 10448, F-63000 Clermont-Ferrand CNRS, UMR6602, LASMEA, F-63177 Aubière A. Moreau antoine.moreau@univ-bpclermont.fr Clermont Université, Université Blaise Pascal, LASMEA, BP 10448, F-63000 Clermont-Ferrand CNRS, UMR6602, LASMEA, F-63177 Aubière M. Mihailovic Clermont Université, Université Blaise Pascal, LASMEA, BP 10448, F-63000 Clermont-Ferrand CNRS, UMR6602, LASMEA, F-63177 Aubière G. Granet Clermont Université, Université Blaise Pascal, LASMEA, BP 10448, F-63000 Clermont-Ferrand CNRS, UMR6602, LASMEA, F-63177 Aubière


INTRODUCTION
Non-specular phenomena in the reflection of beams have been studied since the middle of the nineteenth century (see [1] for a review) from a theoretical point of view and since the middle of the twentieth century with the work of Goos and Hänchen [2] for the experimental part.The phenomenon initially considered is the lateral displacement of a beam in the case of total internal reflexion.A formula giving the asymptotic shift of the beam was soon proposed by Artmann [3] to describe the Goos-Hänchen effect, but it finally proved to be more general since it could account for all lateral shifts, even the giant lateral shifts due to leaky mode excitation [4].This is probably the reason why the displacement of beams is often called "giant Goos-Hänchen effect" though it has little to do with what happens in the case of mere total reflection.At that time it was underlined that these leaky modes could even be responsible for negative lateral shifts [4,5].Later, negative lateral shifts have been predicted for metals [6] and lossy structures [7] -a completely different phenomenon.As a more complete description of the reflected beam's profile was furnished [4,7], lateral shifts appeared as a particular non-specular phenomenon among others [8].By that time it seemed that everything had been said and done about nonspecular phenomena and lateral shifts.
When left-handed materials (LHM) [9] were shown to be feasible [10], it appeared that a large part of the physics of lamellar structures had been missed out [11]- [13].The study of lat-eral shifts thus attracted renewed interest: the Goos-Hänchen effect was shown to be negative for total reflection on LHM [14,15], giant lateral shifts were predicted for structures containing LHM [16]- [20].An analog of the Goos-Hänchen effect was proposed in the case of a beam reflecting on photonic crystals [21] and structures including photonic crystals and producing giant positive or negative lateral shifts were proposed [22,23].Negative shifts were predicted for various structures such as a simple slab when all fields are propagative [24] or for lossy structures [25,26].It should be stressed that lateral shifts are not just a theoretical subject but that they have been measured in various cases [27,28] and have even found applications [29].
All these non-specular phenomena in the reflection of beams are linked to the behaviour of the underlying multilayered structure : they are a signature of a phenomenon.But Artmann's formula or Tamir's approach, since they are tied to the reflection coefficient's properties, do not give any insight into the physics of the structure.The only way to analyze the behaviour of the structure is to use numerical simulation.
The purpose of this paper is to study multilayered structures presenting exotic non-specular phenomena and to interpret previous works, using reliable numerical tools that are given with this article.
The paper is organized as follows.In a first part we recall the

CLASSICAL APPROACHES OF NON-SPECULAR PHENOMENA
The two classical approaches of non-specular phenomena are (i) the study of the asymptotic regime for which simple analytical results like Artmann's formula are available (ii) Tamir's description of the reflected beam's profile when a resonance of the structure is excited.

Asymptotic regime
When the width of the incident beam is large enough, the lateral shift of the reflected beam is given by Artmann's formula where θ is the angle of incidence of the beam in a medium characterized by an optical index n, k 0 = ω/c being the wavevector's modulus in the vacuum, and φ the phase of the reflection coefficient of the structure.
The reflected beam's enlargment in the asymptotic regime was given for the first time by Tamir [8] under the form where α = n k 0 sin θ, w r (resp.w i ) is the waist of the reflected (resp.incident) beam and r the reflection coefficient.
A straightforward calculation shows that this result can also be written where ρ = |r| and ρ = ∂ρ ∂α .
The asymptotic formulae are valid only when the incident beam is so large (and thus spectrally so narrow) that φ can be considered as linear -at least for Artmann's formula.For the enlargement's formulae, which has been largely ignored by the community, the conditions to fulfill are unclear and would probably require a more detailed analysis which is beyond the scope of the present work.Anyway, for the genuine Goos-Hänchen effect the lateral shift predicted by Artmann's formula cannot be reached [30], and very large beams may be required to reach the asymptotic lateral shift when a leaky wave is excited.In general, the lateral shift is thus smaller than Artmann's formulae prediction [31] -which nevertheless gives strong indications that something is happening inside the considered structure.

Description of the reflected beam's profile
The analytical approach of Tamir and co-authors [4,7] relies on the assumption that the reflexion coefficient can be approximated (in the angular domain covered by the beam, centered on α 0 ) by a simple expression based on the presence of a pole and a zero associated to a particular resonance of the structure : where α p (resp.α z ) is the location of the pole (resp.the zero) in the complex plane.Provided the above assumption is founded, this approach is then valid for any width of the incident beam, which makes it both more precise and general.But it cannot describe the original Goos-Hänchen effect, since the variation of the phase is due to the presence of a branch point instead of a pole.Some rather complex calculations [4] lead to the conclusion that the reflected field presents two parts (i) a gaussian part due to the reflection of the incident beam on the first interface and (ii) a part due to the leakage of the resonance excited in the structure for which an analytical approximation can be found.
With this approach, it is possible to conclude that the lateral shift due to the excitation of a leaky mode is much larger than the Goos-Hänchen shift (but that the reflected beam is distorted) and that it can be negative either when the leaky mode is backward [4] or when the structure is lossy [7].In the latter case, the zero is closer to the real axis than the pole which gives rise to a negative lateral shift.
This particular type of negative lateral shift can be predicted using Artmann's formula but in this case the modulus of the reflection coefficient is close to zero and the variation of the phase occurs on so narrow an angular spectrum that the asymptotic regime is extremely difficult to reach.Analytical [7] and numerical [31] studies of such cases have shown that the only effect that can be seen is a splitting of the beam in two and almost no negative lateral shift.In these conditions, it is much more relevant to use a precise description of the reflected beam's profile than Artmann's formula, which may be misleading.

SIMULATION OF MULTILAYERED STRUCTURES WITH UNCONVENTIONAL MEDIA
As already mentioned above, the study of reflected beams profile has recently been extended to the case of structures with a high number of conventional layers [21,22] or containing left-handed materials [12].Multilayered structures may behave as a left-handed material [32] and meta-materials can be considered as multi-layered structures [33].Reliable and efficient simulation tools are of course necessary to study all these cases.

Description of the problem
The problem is to compute the field in a multilayered structure illuminated by a gaussian beam coming from above as shown Figure 1.The spectral amplitude of the beam is given by where x 0 is the position of the beam's center, given by x|A| 2 dx |A| 2 dx, and α 0 = √ µ k 0 sin θ 0 , θ 0 being the angle of incidence of the beam.
For the E (resp.H ) polarization and for a given α the field E y (resp.H y ) inside the j-th layer can be written where , z j is the position of the top of the jth layer and A + j and B + j are computed for an incident plane wave of amplitude 1.
Above (respectively under) the structure, the expression of the field is simply ) by definition of r and t.
The system which has to be solved to find r and t (and eventually all the A + j and B + j coefficients) is then, in the case of the x FIG. 1 The structure presents N − 1 layers.Each layer is characterized by its permittivity ε j and its permeability µ j .The problem is invariant under translation in the y direction.

E polarization and ∀
For the H polarization, µ j only has to be replaced by ε j .

Choice of the numerical method
Numerically, Eq. ( 6) is computed using a regular discretization of the spectral domain so that the field is pseudoperiodical in the x direction.Only a few discrete values of α are thus considered.Using a double sum method as described in [34] could probably avoid the periodization of the problem but it is complicated and costly so that this technique has not been used here.
The main source of numerical instabilities is usually the cancellation phenomenon : since the number of digits retained in any floating number is finite, the addition of two floating numbers is rarely exact.In order to avoid numerical problems, it is recommended to manipulate floating numbers of the same magnitude.A good test for the numerical stability of a method is thus to consider the case of a total frustrated reflexion : due to the evanescent waves, the method will have to deal with both very large and very small numbers.If only numbers of the same magnitude are combined, the method is stable.
We have therefore used all these methods to compute the transmission coefficient in the case of a simple frustrated total reflexion (a slab of air in a medium with a permittivity = 4 and with an incidence angle of π/4).The T matrix and Abèles methods become unstable for a thickness of air greater than 5.5λ.The S matrix and the Dirichlet to Neumann methods remain perfectly stable and provide exactly the same results.The T matrix method is very often used because of its simplicity, but it is unstable even in the case of a simple Bragg mirror when the field is propagative everywhere as shown Figure 2.
The Dirichlet to Neumann (DtN) and scattering matrix method may still present numerical problems of another type : the division by zero.This may happen for particular values of a layer thickness for the DtN method.Using the scattering matrix method, division by zero may happen only in the case of an interface between two media with exactly opposite permittivity and permeability (for a Pendry lens [11], for instance).
Let us finally point out that in this method it is necessary to determine which part of the field is considered upward and which one is downward.When left-handed materials are present in the structure then this choice can be based either on the wave vector or on the Poynting vector since these are opposite in left-handed media.
Analytically, these choices are equivalent inside a layer [40] but of course not outside [41] (or the reflection coefficient would not be correctly defined).Mathematically, this means that the choice of the branch cut for the determination of the propagation constant γ j = (ε j µ j k 2 0 − α 2 ) 1/2 inside a layer has no importance for the determination of the field and the reflection coefficient.That is why this branch cut does not appear when the complex plane representation of a reflection coefficient is computed [20].
But numerically, this choice has an importance : a better stability is obtained if the direction of propagation is such that the wave is decaying in this direction.This can be easily done by chosing a determination of the square root inside the layers so that the imaginary part of the constant propagation is always positive.
The Matlab code which corresponds to this method is provided as supplementary files (size 16 kB, format: m, see field.m,cascade.m,intermediaire.m,field octave.m)with this paper.We tried to make its use as simple as possible for the reader.It is especially adapted to the simulation of structures containing a high number of layers.

Reflection on Bragg mirrors
Two examples are here considered to illustrate how useful numerical methods can be to understand physical situations.Figure 3 shows the structure of the field when an incident beam is reflected by a Bragg mirror.The field inside the Bragg mirror is similar to an evanescent Bloch mode although the field in each layer in still propagative.The structure is illuminated by a beam characterized by a wavelength of λ = 500 nm, an incidence angle of 42 • and a waist of 100λ.FIG. 4 A defect buried in a Bragg mirror is illuminated using a gaussian beam (wavelength λ = 0.965, waist 100λ, incidence angle θ 0 = 24.675• ) in E polarization.The Bragg mirror is a stack of alterning medium 1 ( 1 = 2.2 2 ) layers with a thickness of the thickness of a medium 2 layer is doubled.The total number of patterns is about 40.The structure is rather similar to the one proposed in [22].
The latter remark will allow us to interpret another giant lateral shift [22].Figure 4 shows what happens when a Bragg mirror with a defect is illuminated.Obviously, the field is Bloch-evanescent above and under the defect so that everything happens as if a guided mode was excited in a dielectric slab by evanescent coupling -and no "surface wave" is involved, as suggested previously [22].Moreover, a defect corresponds to a single pole of the reflection coefficient in the bandgap of the mirror so that the classical analysis [4] based on leaky modes is undoubtedly valid.The profile of the reflected beam corresponds exactly to this analysis provided the defect is buried deep enough in the Bragg mirror as shown Figure 4.The exponential decay of the profile in the propagation direction is controlled by the imaginary part of the pole corresponding to the guided mode.

Co-directional coupling
Let us consider the case of two coupled waveguides excited by evanescent coupling.Figure 5 shows the lateral displacement of the reflected beam in three situations (i) for a single waveguide, (ii) for two coupled waveguides with the same incidence angle of the beam and (iii) for the two coupled waveguides and an incidence angle corresponding to the symmetric mode.When a leaky mode is excited (situations i and iii) the lateral shift increases with the beam width, tending towards the prediction of Artmann's formula.It is the case in situation ii, but the lateral shift is surprisingly much higher for narrower beam than the asymptotic shift and even than the lateral shift for the other situations (and the same width of the incident beam).All the distances are given in wavelengths.Blue curve : lateral shift associated to the excitation of the guided mode of a single dielectric slab.For narrow beam the shift is much smaller than for the asymptotic regime, reached for a width of 600λ.Green curve : coupled waveguides and an incidence angle of 33.85 • .The asymptotic regime is very difficult to reach and higher than for one waveguide.Red curve : coupled waveguides and an incidence angle of 33.9  6 The phase of the reflexion coefficient in two different cases of coupling.Red curve : co-directional coupling between two waveguides.The symetric and antisymetric modes are excited for different angles of incidence.The oscillation of the energy between the two waveguides can be obtained for an intermediate incidence angle of 33.9 • .Blue curve : contra-directional coupling between two waveguides.In both cases the modulus of the reflexion coefficient is equal to one over the whole angular range.
Figure 6 shows the phase of the reflection coefficient for the coupled waveguides structure.This is particularly interesting because the ratio between the lateral shift and the beam width is particularly high in this case -and thus particularly easy to measure.Of course, this phenomenon occurs for narrow beams -when the asymptotic regime is not reached.But in this case, Artmann's formula is even misleading -which is not common.The physical explanation of such a phenomenon is given Figure 7.The picture clearly shows that the energy is oscillating between the two waveguides.As long as the energy is "hid-FIG.7 Top: The incoming beam (waist of 100λ, 33.9 • incidence angle) propagating in the upper medium (ε = 5, µ = 1) excites two coupled waveguides (h = 0.285391λ, surrounded by air and separated by λ.The distance between the prism and the first waveguide is of 0.65λ.Bottom: The same, but with one waveguide only.
In both cases the simulation domain is 6000λ large and about 6λ high.den" in the deeply buried waveguide, it does not leak out.
That is what explains the unexpected large lateral shift for small widths.
The profile of the reflected beam is so peculiar (see Figure 8) in this case that it cannot be accurately described by Tamir's approach.The latter would be relevant for a single waveguide, but the exponential decay of the profile is faster in this caseso that Tamir's approach cannot account for the whole lateral shift.

Contra-directional coupling: the light wheel
A new concept has recently been developed in the context of multi-layered structures: when a conventional dielectric waveguide is coupled with a backward waveguide (a layer of LHM with carefully chosen characteristics [42]) then the light forms a localized mode called a "light wheel" [13].This particular mode can be excited using evanescent coupling (i.e. a prism).Due to the evanescent coupling, the light wheel leaks out of the structure, leading to an unusual beam profile.
Here the structure constists in two coupled waveguides which can be excited using an incident beam propagating in a high index medium.The whole structure is presented in Figure 9.
The upper guide is a mere dielectric slab while the lower waveguide is made of a LHM.The thicknesses are chosen so that the two guides present a guided mode for the same propagation constant.
Figure 10 shows the modulus of the field in the structure when a gaussian beam with the correct incidence angle illuminates the structure.
The reflected beam is actually composed of the beam reflected by the first interface and of the leakage of the light wheel as shown Figure 11.This explains why the beam is enlarged (the leakage being as large as the light wheel) but and only slightly displaced (since both the primary reflected beam and the light wheel are centered).
Here, we are not in the asymptotic regime so that the results given by Eqs. ( 1) and ( 3) are not expected to hold.The excitation of the light wheel produces a negative lateral shift of 0.64λ.According to Artmann's formula, this shift should be negative (see the phase of the reflection coefficient Figure 6) but of 6.49λ.And according to Eq. ( 3), there should be no enlargement of the reflected beam.Using an asymptotic analysis in this context could thus be misleading.
Since the reflection coefficient presents two poles with identical real parts and opposite imaginary parts [13] the assumption that it can be approximated with a single zero and a single pole is not valid.This explains why the reflected beam's profile cannot be described using Tamir's approach [4,7].
Therefore, when a light wheel is excited, all the classical analytical analyses fail to describe what happens : the reflected beam is hardly displaced but it is enlarged by the light wheel.This is a case for which reliable numerical tools must be employed -or new analytical tools should be developed, different from what has been proposed until now.Finally, this shows that it would probably be difficult to extend the analysis of Tamir to the case of co-and contra-directional coupling: although in both cases two poles are responsible for what happens, the two phenomena don't have much in common.

CONCLUSION
This work is a contribution to the study of beam reflection on multi-layered structures.We have shown that some peculiar phenomena could lead to exotic reflected beam's profiles and unexpectedly high lateral displacements (for narrow beams).This underlines the limits of the classical approaches of non-specular phenomena through the asymptotic analysis or through analytical models of the reflected beam's profile.
Artmann's formula gives clues that something is happening in a structure, and useful indications on the sign of the lateral shift.But its predictions are not fully reliable (the asymptotic limit is sometimes hardly reached) and it cannot explain lateral shifts.In the cases we have explored, it may even prove misleading.Maybe the community should use asymptotic analysis more carefully and, as far as possible, compare to numeri-cal simulations.Extremely narrow beams are nowadays easy to produce, so that phenomena occuring before the asymptotic regime should not be neglected.
Tamir's work, although of considerable importance, is not definitive : it is for instance unable to describe the profile of a reflected beam when a light wheel is excited.An analytical approach of these exotic phenomena is still lacking.
That is the reason why, with this article, we have given our simulation tools for multi-layered structures and paid much attention to the stability of the methods we have employed.
We have made these programmes as simple to use as possible.Our hope is that our work will be useful to the community -for we are convinced that the physics of multi-layered structures is far from being exhausted.

FIG. 2
FIG. 2 Transmission coefficient of a Bragg mirror computed using transfer matrices (red curve) and scattering matrices (blue curve).The unstability of the transfer matrix method causes a non nul transmission in the bandgap.The Bragg mirror, illuminated in normal incidence, consists of a succession of ε = 2 and ε = 6 media of thickness 88 nm and 51 nm respectively.This curve has been obtain for around 90 patterns of the Bragg mirror.

FIG. 3
FIG. 3 On this image the Bragg mirror (lower part) has been represented 10 times bigger than it really is, in order to see the incoming beam as well as what happens in the photonic crystal.The Bragg mirror is a succession of ε = 2.25 and ε = 4 media of thickness 100 nm and 75 nm respectively.The number of periods is of 30.

FIG. 5
FIG.5The lateral shift versus the width of the incident beam (four times the waist).

FIG. 8
FIG.8The figure shows the profile of the reflected beam in the case of a single waveguide (blue curve) and of two coupled waveguides (red curve).The parameters are given in the caption of Figure7.
FIG.11Profile of the reflected beam when a light wheel is excited (red curve).The profile of the incoming beam has been plotted for comparison (blue curve).The reflected beam is obviously enlarged but not particularly shifted.