On anomalous optical beam shifts at near-normal incidence

We develop the theory of optical beam shifts (both Goos-Hanchen and Imbert-Fedorov) for the case of near-normal incidence, when the incident angle becomes comparable with the angular beam divergence. Such a situation naturally leads to strong enhancement of the shifts reported recently [ACS Photonics 6, 2530 (2019)]. Experimental results find complete and rigorous explanation in our generalized theory. In addition, the developed theory uncovers the unified origin of the anomalous beam shifts enhancement via the Berry phase singularity. We also propose a simple experimental scheme involving quarter-wave plate that allows to observe the giant transverse and longitudinal, spatial and angular beam shifts simultaneously. Our results can find applications in spin-orbit photonics, polarization optics, sensing applications, and quantum weak measurements.

The reflection and refraction of a plane electromagnetic wave at a dielectric interface are the basic physical processes inherent to all-optical systems and devices. They are rigorously described by the Snell's law and Fresnel equations. Nevertheless, in practice, we usually deal with optical beams of a finite width for which the plane wave approximation is oversimplified. This results in the deviation from geometrical optics that manifests itself as spatial shifts of the beam known as the Goos-Hänchen (GH) and Imbert-Fedorov (IF) shifts. These shifts can be explained in terms of weak material-mediated interaction of the beam spectrum (orbit state) and its polarization (spin state) [1][2][3][4].
Today, the GH and IF shifts, including the photonic spin Hall effect (PSHE) as a particular example of IF shift, are well-studied in different material systems including atomic optics, optical sensors, graphene, metasurfaces, polarizers, uniaxial crystals, etc [5][6][7][8]. Usually, the spatial (angular) beam shifts are very small -typically of the order of the light wavelength (the beam angular spectrum variance) -that limits their application [1][2][3][4]. The shifts can be enhanced under several specific conditions including the near-Brewster incidence [9][10][11][12], material resonances [13,14], exceptional points [15], and output beam polarization post-selection [16]. In all these cases, the enhancement occurs at large angles of incidence (typically, > 10 • ). However, the anomalous enhancement of the PSHE in hyperbolic metamaterial slabs has been experimentally observed recently at a near-normal angle of incidence [17,18]. This is very counter-intuitive as at the strictly normal incidence the PSHE completely disappears. Therefore, at first sight, there are no rigorous arguments for the enhancement of PSHE at small angles of incidence. Moreover, the standard theory of optical beam shifts [4] fails to explain the observed effect.
In this Letter, we develop a generalized theory of beam shifts addressing small incident angles. We show that anomalously large spatial and angular shifts of both types (GH and IF) could be simultaneously observed for the beam near-normally transmitted or reflected at an uni-axial slab. Usually, the longitudinal (GH) and transverse (IF) shifts are caused by different origins and considered independently. Here, we show that anomalous GH and IF shifts have the same nature caused by the Berry phase singularity when the incident angle becomes comparable to the angular divergence of the beam. While the partial cases of near-normal IF shift enhancement have been observed [17][18][19][20], we predict another intriguing feature, namely the anomalous GH shift which is commonly known to exist under total internal reflection conditions [4]. In the frame of the developed generalized theory, we analyze necessary and optimum conditions for the giant beam shifts of all types (GH and IF, linear and angular) and reveal that a conventional lowbirefringent quarter-wave plate (QWP) can be used in such a remarkably simple experiment. Finally, we analyze the relevant beam parameters, polarization structure and intensity profiles of the shifted transmitted beam.
Standard theory of optical beam shifts.-We consider first the oblique incidence of a monochromatic optical beam on a flat vacuum-medium interface, see Fig. 1. The values of the shifts could be calculated within the Jones matrix formalism in paraxial approximation. The Jones matrixT a relates the incident and reflected/transmitted plane wave amplitudes in the beam coordinate frame |E a =T a (ϑ, µ, ν)|E , where |E ∝ |e · f (µ, ν) is the incident constituent plane-wave Jones vector, index a = r, t denotes reflected and transmitted waves, respectively, |e = (e X , e Y ) T is the Jones vector of the central plane wave in the incident beam, and f (µ, ν) is the incident beam Fourier spectrum expressed through in-plane (µ) and out-of-plane (ν) deflections of the non-central wave vectors [4]. The linear ( X a , Y a ) and angular ( P a X , P a Y ) GH ( X a , P a X ) and IF ( Y a , P a Y ) shifts are then obtained using their original definition meaning the quantum-mechanical expectation values of the coordinate and momentum operators [4]: 1. Schematic geometry of the problem. The long black arrow marks the non-central constituent plane wave with incident wavevector k i , angular in-plane (µ) and out-of-plane (ν) deflection components. Its plane of incidence forms a polar deflection angle φ with the one of a central plane wave. The yellow arrows are collinear with the wavevectors of the corresponding central plane waves. The incident angle of the central plane wave is marked as ϑ. The inset schematically shows the incident wavevectors of the non-central plane waves of a beam forming the large polar deflection angles.
where the integration is taken over the transverse momentum components k a X = γ a k 0 µ and k a Y = k 0 ν in the beam coordinate frame, k 0 = 2π/λ is the light wavevector in vacuum, λ is the light wavelength, and factor {γ t,r } = {1, −1} accounts for the inversion of the xcomponents of the wavevectors with reflection. Here, N a = E a |E a is the normalization factor proportional to the beam intensity. At large incidence angles, Jones ma-trixT a is slightly inhomogeneous across the beam spectrum (∂T a /∂k a X,Y ∼ λ), which leads to small shifts of the order of the light wavelength and angular spectrum variance for spatial and angular shifts, respectively (see Supplemental Material [21]).
Geometric Berry phase singularity.-The obvious absence of an IF shift under strictly normal incidence on one hand, and its singularity in approximate expressions for large incidence angles (∝ cot ϑ where ϑ is the incidence angle [4]) on another hand, assume the sharp IF shift peak in the vicinity of normal incidence. This motivates us to resolve explicitly this geometric singularity by analyzing the distribution of polar deflection angle φ, see Fig. 1. Under large incident angles the polar deflection angle φ of any constituent plane wave is small Another situation takes place for the near-normal incidence, where different plane-wave components constituting the beam have planes of incidence oriented arbitrarily with respect to the central plane wave and cover a large range of polar deflection angles (at strictly normal inci-dence, φ ∈ [−π/2, π/2], see the inset of Fig. 1). In an anisotropic medium, the wide spreading of the incidence planes results in the significant difference of Fresnel's amplitudes for different plane-wave components. Thus, under some specific conditions (anisotropy degree, optical axis orientation, retardation phase of anisotropic slab), the optical beam incident at sufficiently small angles can experience anomalously large beam shifts. To show this explicitly, we note that at near-normal incidence, Jones matrixT a (ϑ, µ, ν) depends on incident angle ϑ, in-plane µ and out-of-plane ν deflections only implicitly, through the polar deflection angle φ [T a (ϑ, µ, ν) ≡T a (φ)], where the latter could be approximated as follows [21] φ(ϑ, µ, ν) ≈ tan −1 ν ϑ + µ .
When ϑ is of the order of beam divergence, one can notice from Eq. (3) that the polar deflection angle φ may be large enough for a finite portion of the beam spectrum, and even close to ±90 • at ϑ −µ. Thus, Jones matrixT a (φ) will be substantially inhomogeneous across the beam spectrum (∂T a /∂k a X,Y ∼ w 0 , where w 0 is the beam waist), which, in turn, will cause anomalous beam shifts (see Fig. S2 in Supplemental Material [21]).
In fact, the polar deflection angle at near-normal incidence (ϑ → 0) is equivalent to the geometric Berry phase, Φ B = −φ cos ϑ −φ [4]. Thus, the considered geometric singularity of the polar deflection angle simultaneously means the Berry phase singularity. The large IF shift arises from large transverse momentum derivatives of the Berry phase across the beam spectrum according to its definition.
Another surprising consequence of Eq. (3) is that it leads to the anomalous GH shift. In conventional situation (large angles of incidence), the GH shifts are caused by the spatial dispersion of the scattering coefficients [4]. In contrast, at near-normal incidence GH shift arises from large longitudinal momentum gradient of the Berry phase across the lateral parts of the beam spectrum. Note that such longitudinal gradients caused by strong spin-orbit coupling are negligible at large incident angles.
Optical beam shifts at near-normal incidence.-To exhibit the near-normal Berry phase singularity we consider transmission through an anisotropic slab -a waveplate. Following the Jones matrix formalism we calculate analytically the GH and IF shifts of the Gaussian beam transmitted through the uniaxial waveplate when the beam incidence plane is parallel to the waveplate optical axis, which lies at the interface (Fig. 1). The Gaussian beam angular spectrum is defined by the Fourier spectrum f (µ, ν) = exp[−κ 2 (µ 2 +ν 2 )/2], where κ = k 0 w 0 / √ 2 is the inverse beam divergence. Results for the four transmitted beam shifts read (the detailed derivation is done in Supplemental Material [21]) where S 1 = |e x | 2 − |e y | 2 , S 2 = 2 [e * x e y ], and S 3 = 2 [e * x e y ] are the corresponding Stokes parameters of the incident beam, N t = (τ + + τ − S 1 Λ Y (κϑ))/2 is the squared norm of the transmitted beam, t − = t e − t o , τ +,− = |t e | 2 ± |t o | 2 , and τ × = 2 (t e t * o ) are real coefficients depending on the ordinary t o and extraordinary t e transmission amplitudes for the plane wave under normal incidence [21]. Hereinafter, we use the dimensionless shifts values defined as P X,Y = P X,Y · κ 2 /k and X, Y = X, Y · k. The nonlinear geometric resonant factors Λ Y and Λ X are the following Note that the squared norm of transmitted beam state N t (ϑ) plays a role of a slight renormalization, in sharp contrast to conventional incident angles cases, where it causes singularity in the shifts by approaching zero [9][10][11][12][13][14][15][16]23]. For the reflected beam, we arrive at the same results as in Eqs. (4)-(7) with substitutions {t e , t o } → {r e , r o }, γ t → γ r , and N t → N r . Equations (4)- (7) are the main result of this work applicable for any cases except κϑ 1. One can see that each shift in Eqs. (4)- (7) is the product of three factorsincident beam polarization (Stokes parameters), material properties of the uniaxial slab (τ +,−,× ) and geometric resonance [f (κϑ)]. The second term in the expression (7) for the spatial IF shift is a counterpart of the anomalous PSHE connected with the circular polarization of the incident beam. Moreover, due to the presence of the first term in (7) the IF shift can be obtained with the arbitrary polarization state except strictly TM-or TEpolarization.
We also note the remarkable connection of the four shifts to the geometric resonant terms Λ X,Y /(ϑN t ). First, there is a proportionality between spatial and angular shifts, which is only broken by the spin-Hall term in Eq. (7). The simultaneous appearance of all possible optical beam shifts (namely, spatial and angular, GH and IF shifts) at the same near-normal angle of incidence is unique. At large incident angles (κϑ 1), spatial and angular GH shifts can be simultaneously observed Nonlinear geometric resonant factors (8) appearing in shifts (4)- (7): ΛY (κϑ)/ϑ (dark blue) and ΛX (κϑ)/ϑ (light blue), along with their asymptotes at large incidence angles 1/ϑ (dashed black) and 2/(κ 2 ϑ 3 ) (dotted black), respectively. Here, we used the inverse beam divergence κ = 700, namely λ = 630 nm, w0 100 µm. Inset shows the dependence of function ΛY (ϑ)/ϑ on the inverse beam divergence κ. The arrow indicates the sequence of increasing beam waist, κ = 300, 700, 2000, respectively.
only in special cases such as lossy media [24] or vortex beams [25]. Second, there is a further similarity of geometric resonant terms along Y and X directions, whereas Λ X (κϑ) = e −κ 2 ϑ 2 Λ Y (iκϑ). This reflects the fact that all anomalous near-normal incident shifts have the same origin, namely the geometric Berry phase singularity. In addition, the near-normal anomalous shifts form a connection to superoscillation-related physics [26].
Impact of beam width.-The optical beam shifts under near-normal incidence depend substantially on the beam waist w 0 via the inverse beam divergence κ, as it is expressed by nonlinear functions Λ X,Y (κϑ) [Eq. (8)]. The IF shifts are anomalous and have asymptotes ∝ 1/ϑ at large incidence angles (κϑ 1) which is hinted by the standard theory [4]. The GH shifts have higherorder asymptotes at large incident angles ∝ 1/(ϑ 3 κ 2 ) [ Fig. 2]. This explains why the anomaly in GH shifts was not anticipated by the standard theory for uniaxial slab, in which the spatial GH shift is absent except the case of total internal reflection. Moreover, one can notice that the amplitude of the optical beam shifts peak is proportional to the beam waist, while the corresponding critical angle is inversely proportional to it (see the inset in Fig. 2). In general, the critical angle corresponding to the maxima of GH and IF shifts is about a few tenths of a degree for a beam width of several dozen microns at λ = 630 nm. This result absolutely agrees with the recent experimental observations of the anomalous PSHE under near-normal incidence at hyperbolic metamaterials [17,18].
Dependence on slab material parameters.-From a rigorous analysis of near-normal shifts in isotropic media we found that the geometric Berry phase singularity at nearnormal incidence is indeed manifested in the shifts only when the difference of Fresnel coefficients for TM-and TE-polarizations is finite at normal incidence [21]. Then, we focus on the conventional waveplate with optical axis lying at the interface within a beam incidence plane. In the general case, the proposed generalized theory could be applied to arbitrary optical axis orientation and extended to biaxial and bianisotropic media. A particular case of an uniaxial slab with the optical axis perpendicular to the interface has been recently considered in Refs. [19,20], where the enhanced PSHE combined with the spin-controlled vortex generation was observed and explained by the divergence of the Berry phase. All material dependences in Eqs. (4)- (7) are factored out in four dimensionless coefficients, τ + , τ − , τ × , and t − . Both angular shifts are proportional to τ − , while the spatial shifts are proportional to τ × except for the PSHE term in the spatial IF shift which is proportional to |t − | 2 . Therefore, we search numerically for the maximum values of these parameters as functions of the waveplate retardation Γ = (n e − n o )k 0 δ z , where δ z is the slab thickness, and n = (n o + n e )/2 is the average refractive index of the uniaxial slab, see Fig. 3. First, we find that both |τ + | and |τ × | achieve maximum for the quarter-wave plate condition, i.e. for Γ = ±π/2 (mod 2π), while |t − | 2 /2 reaches maximum value for a half-wave plate (HWP). This may also be shown analytically for small birefringence |n e − n o | 1, see the Supplemental Material [21]. The dependence of all parameters on Γ is oscillatory due to the Fabry-Pérot resonances in the slab [21,27]. Secondly, we consider a QWP and show that max(|τ + |) grows as a function of average refractive index n up to the saturation limit max(|τ + |) ≈ 1 at n ≈ 10, while max(|τ × /2|) reaches the same maximum, but in a distinct limit n → 1. Hence, for realistic QWP with n ≈ 3, |τ + | and |τ × | can be around 0.6 and 1.2, respectively. Thus the low-birefringence QWP is a perfect candidate for observing all four types of anomalous shifts simultaneously (although PSHE has maximum for a HWP). Fig. 4 shows the optical beam shifts dependencies on the small incidence angle for the beam transmitted through the QWP with optical axis lying in the beam plane of incidence. It is important to note that a linear IF shift can change the shift direction to the opposite one at some angle under elliptically-polarized beam illumination (Fig. 4c). This is caused by the interplay between differently polarized terms of the linear IF shift [Eq. (7)]. The switching of the IF direction could be also observed for strongly anisotropic systems due to the polarization mixing, especially when the optical axis is arbitrarily oriented [28].
Intensity and helicity profiles of the anomalously shifted beam.-The intensity and S 3 (helicity) patterns simulated with Jones matrixT t are shown in Fig. 5 for the corresponding input eigenpolarizations. The helicity dynamics for an ordinary TM-polarized beam (S 1 = 1) provides additional insight into the nature of the anomalous shifts at near-normal incidence. When the incidence angle ϑ exceeds the beam divergence 1/κ, PSHE naturally occurs owing to the so-called circular birefringence Panels a-d: normalized helicity (Stokes parameter S3) patterns in k-space (a1-d1) and in real space (a2-d2) for the ordinary polarized (S1 = 1) beam near-normally transmitted through QWP. The ellipses in panels (a1-d1) indicate the corresponding polarization profile in k-space. Panels e-h: intensity profiles at maximum linear and angular GH and IF shifts (shown by white arrows) of transmitted beam for their respective eigenpolarizations. Parameters are the same as in Fig. 4. The plots in the panels are the squares with the side 10w0 for all real-space panels, and 10/κ for all k-space panels.
( Fig. 5a) [7]. At near-normal incidence ϑ 1/κ, spin-Hall doublets from lateral and opposite parts of beam Fourier spectrum become prominent, and an octupole helicity profile is built up in k-space with a phase singularity at a distance ∼ k 0 ϑ from the geometric-optics beam axis (Figs. 5b-5c). In real space, at ϑ = 0, the helicity is carried by a mode resembling a Bessel-Gaussian mode with m = 4 [29] and intensity null on the beam axis (Fig. 5d). This mode comes from the extraordinary component induced by the waveplate and deforms along with the k-space profile, accompanying the anomalous spatial shifts. Transverse helicity imbalance accompanies the angular IF shift for incident beam polarization different from S 1 . The presence of such a mode and the polarization singularity provide another close connection to the superoscillations [26]. Finally, we demonstrate explicitly the beam shifts of all four kinds in Figs. 5e-5h. The anomalous GH angular shifts are accompanied by longitudinal beam spectrum deformations around this polarization singularity (Fig. 5e).
To conclude, we have developed a theory of spatial and angular shifts of optical beams undergoing near-normal incidence. The theory is able to explain experimental peculiarities of strongly enhanced IF shifts obtained when the angle of incidence was less than one degree, and predicts the additional anomalous GH shifts. The principal point of the theory is that anomalously large shifts of all types are shown to have the same fundamental originthe singularity in the Berry phase appearing in the beam Fourier spectrum. The results have been illustrated for an analytically treated example of a uniaxial waveplate. We propose a very simple configuration with a quarterwaveplate, which supports anomalously large simultaneous near-normal incidence shifts of all kinds in manifold times exceeding shifts observed separately in conventional cases. Thus, the anomalous near-normal shifts are completely feasible for direct experimental verification. The results obtained could arise for other types of waves (acoustic, electron beams, etc.) material anisotropy such as biaxial and bianisotropic media.
O. Y. and A. B. acknowledge support from the Foundation for the Advancement of Theoretical Physics and Mathematics "BASIS". The authors thank Konstantin Bliokh from RIKEN for fruitful discussions. In order to demonstrate the keynote point distinguishing our work from all other studies devoted to GH and IF shifts, we derive the strict expression for the polar deflection angle φ [Eq. (3)], responsible for the deflection of noncentral constituent plane wave from the central one. The corresponding relation between the polar deflection angle φ and the incident angle ϑ, small in-plane µ and out-of-plane ν deflections can be derived from Fig. S1:

Supplemental Material: On anomalous optical beam shifts at near-normal incidence
where θ = ϑ + µ is the total incident angle of the non-central wave in the plane of incidence. Importantly, the explicit form of sin φ and cos φ could be expressed as follows: FIG. S1. Schematic geometry for the calculation of polar deflection angle φ for a component plane wave with wave vector depicted by a gray arrow. The gray vertical plane marks the incidence plane of central plane-wave in the beam.
Within the previous studies, e.g. [S1], it was assumed that ν, µ ϑ, so the polar deflection angle was determined as φ = ν/ sin ϑ being always negligibly small in comparison to the angle of incidence φ ϑ. This statement is absolutely relevant for any not near-zero angle [S1].

arXiv:2107.09738v1 [physics.optics] 20 Jul 2021
However, in the case of near-normal incidence, when ϑ becomes comparable with ν and µ, one can notice from Eq. (S1) that the polar deflection angle may be extremely large. Under condition ϑ + µ = 0 it could be even close to 90 • . Another appealing feature of Eq. (S1) is the fact that φ becomes µ-dependent under near-normal incidence. This arbitrary beam-width-dependent polar angle deflection under the special case of near-normal incidence is the keynote fundamental insight of this work leading to the unconventional behaviour of GH and IF shifts described in the main text. In particular, beam waist can strongly affect not only the amplitude but also the sign of φ. Therefore, the planes of incidence for different plane waves constituting the beam can be substantially different.
Under near-normal incident angles, the geometric Berry phase of constituent non-central plane wave is induced by the azimuthal rotation of the plane of incidence and is completely determined by the polar deflection angle Φ B = −φ cos ϑ −φ [S1]. So, the intriguing properties of the optical beam shifts under near-normal incidence are defined by the singularity of the geometric Berry phase.

SUPPLEMENTAL NOTE 2: JONES MATRIX FORMALISM FOR OPTICAL BEAM INCIDENT ON UNIAXIAL SLAB ALONG OPTICAL AXIS AT NEAR-NORMAL INCIDENCE
We follow the formalism of the Jones matrix describing the beam interaction with the interface including spatial dispersion effects in the momentum space representation [S1] T a =Û a † (ϑ a , k a )F a (ϑ, k a )Û a (ϑ a , k a ). (S3) Here, index a denotes reflected (r) and transmitted (t) beams, respectively;Û a (ϑ a , k a ) is the rotational transformation of the fields from the beam coordinate frame to the spherical basis of the TE and TM modes;F a (ϑ, k) = (f a pp , f a ps ; f a sp , f a ss ) is the Fresnel matrix for the fields in the spherical basis of the TE and TM modes; ϑ is the angle of incidence, ϑ a is the reflected/refracted angle; k a k a c + (k a X u a X + k a Y u Y ) is the wave vector of constituent plane wave of the beam, k a c = n a k 0 (sin ϑ a , 0, cos ϑ a ) is the wave vector of the beam central plane wave, n a is the refractive index of the corresponding medium, k 0 = 2π/λ is the free-space wave vector, λ is the wavelength of incident light, u a X and u Y are the unit vectors in the geometric-optics beam coordinate system (beam propagates along z-axis), k a X = k 0 γ a µ and k a Y = k 0 ν are the transverse wave vector components of the constituent plane waves, γ a = cos ϑ/ cos ϑ a , µ and ν are the in-plane and out-of-plane deflections of non-central wave vectors as it is shown in Fig. S1. Hereinafter, we consider only square 2 × 2 matrices corresponding to the transverse components of the fields which are essential for a paraxial beam. Next, we derive the rotation matrixÛ a (ϑ, k). As already mentioned, this matrix corresponds to the rotational transformation of the electric field of the particular constituent plane wave in a beam from the beam coordinate frame to the (local) spherical basis of the TE and TM modes, where the Fresnel matrix would be diagonal, which involves three consequent rotationsÛ a (ϑ a , k a ) =R y (θ a )R z (φ a )R y (−ϑ a ) [S1]. Here, the first (in acting order) is the rotation from beam coordinate frame to laboratory frameR y (−ϑ a ), followed by rotation from the laboratory frame to the local spherical coordinate frame (which is aligned with the constituent plane-wave component plane of incidence and has axis z normal to the interface), obtained with rotation by a polar deflection angleR z (φ a ), and the rotationR y (θ a ) back to the beam coordinate frame. The part of matrixÛ a (ϑ, k) for transverse fields has the following form:Û a (ϑ, k) = cos ϑ cos θ cos φ + sin ϑ sin θ cos θ sin φ − cos ϑ sin φ cos φ . (S4) Under near-normal incidence, we may set in Eq. (S4) that cos ϑ ≈ 1, cos θ ≈ 1 in the leading order, and neglect sin ϑ sin θ ≈ ϑ (ν) 2 + (ϑ + µ) 2 cos φ, leaving µ-dependence only in sin φ and cos φ (which leads to neglecting terms ∝ µ 2 and higher-order terms in ϑ). Therefore, the rotation matrix in the leading-term approximation is just the matrix of rotation by polar angle φ in transverse beam coordinates: The ordinary (o) and extraordinary (e) transmission and reflection coefficients for an uniaxial plate could be deduced using 4 × 4 transfer matrix method [S2] where n o and n e are the refractive indices for ordinary and extraordinary wave in an uniaxial crystal, respectively, η o,e = (1 − n o,e )/(1 + n o,e ), andδ z = k 0 δ z is the dimensionless normalized slab thickness, δ z is the slab thickness.
Here, we have neglected O(ϑ 2 ) terms in all transmission coefficients since ϑ → 0.
In the general form, Fresnel matrices for transmitted (F t ) and reflected (F r ) plane waves can be expressed as follows:F t = t pp t ps t sp t ss ,F r = r pp r ps r sp r ss , where t ij and r ij are the complex transmission and reflection coefficients. The corresponding transmittance and reflectance are defined as T ij = |t ij | 2 and R ij = |r ij | 2 , respectively. Indices p and s correspond to TM and TE polarization, respectively. The components of the transmissive Fresnel matrix (S8) are t pp = t e cos 2 ϕ + t o sin 2 ϕ, t ss = t o cos 2 ϕ + t e sin 2 ϕ, where ϕ is the angle between the incidence plane for constituent beam plane wave and the optical axis of an uniaxial crystal. The expression (S9) shows that at nonzero angles ϕ off-diagonal components t sp = t ps of the Fresnel matrix mix s and p polarizations of a plane wave, with maximum mixing at ϕ = π/4. Hereinafter, we consider the special case, when the optical axis of an uniaxial crystal lies at the intersection of the planes of the interface and central beam plane wave, so that φ ≡ ϕ. Therefore, in this approximation ϑ enterŝ F t,r (φ) only implicitly through the deflection angle φ. This is the least polarization-mixing and counter-intuitive case for polarization dependence of beam shifts, while one can generally consider the arbitrary values of ϕ leading to the increase of the polarization mixing effects. Moreover, this case also describes the optical axis lying within the interface of uniaxial crystal, but oriented perpendicularly to the central plane wave of a beam, with the corresponding substitution t e,o → t o,e .
The Jones matrices for transmitted light in an uniaxial waveplate are readily derived from their definition in the main text, and at any plate thickness we find the following compact form for them: where complex coefficients τ +,− and ρ +,− are defined as follows The product of Jones matrices appearing in Eq. (1) could be expressed as followŝ where In the case of reflection, the product of Jones matrices (S12) remains the same with proper substitution t o,e → r o,e .

SUPPLEMENTAL NOTE 3: OPTICAL BEAM SHIFTS FOR THE QUARTER-WAVE PLATE UNDER NEAR-NORMAL INCIDENCE
In this section, we derive the explicit form of Jones matrices for a uniaxial waveplate and analyze the impact of materials parameters (namely, anisotropic refractive index tensor and thickness) of uniaxial slab on the magnitude of optical beam shifts.
According to Eqs. (4)-(7) of the main manuscript, we find that the angular and linear optical beam shifts are proportional to τ − and τ × , respectively. Here, we perform the analytical analysis for low-birefringent (|n o − n e | 1) uniaxial plate in order to determine the waveplate retardation corresponding to maximum of τ − and τ × . In the case of transmitted beam, one can write where Γ = (n o −n e )δ z is the waveplate retardation, x = (1−η 2 e ), y = η 2 e ·e 2iδzne . Fabri-Pérot condition for constructive interference in a waveplate implies that the phase 2δ z n e ≡ 0 mod 2π, which ensures that y should be a real number. In addition, since |η e,o | < 1 we can state that 0 < (1 − y) < 1. Then, the maximum is achieved for the most different denominators in Eq. (S14), i.e. for e 2iΓ = −1, which gives finally Γ max ≡ ±π/2 mod 2π corresponding indeed to the quarter-wave plate retardation.
Similarly, we can prove that maximum of τ × = 2 (t e t * o ) also appears for the quarter-wave retardation. First, we can write explicitly . (S15) Then, after noting that y should be real-valued (see discussion for τ − above) it takes the form .
(S16) Therefore, maximum τ × is achieved for Γ max ≡ ±π/2 mod 2π, which is again the quarter-wave plate retardation. Analytic formulas for maximum values of |τ − | and |τ × | at near-quarter-wave plate retardation thus read where η = (1 − n)(1 + n), and n = (n o + n e )/2 is the average refractive index. Analytical results (S17) are plotted in the lower panel of Fig. 3 in the main text (solid lines) along with the numerical results (points), where we find a good agreement.
Moreover, taking into account the form of Jones matrix (S10), we may see why anomalous GH shifts arise in the case of normal incidence at the QWP. It is well known that the GH shifts are proportional to the µ-derivatives of Jones matrix elementsT t 11 andT t 22 [S1]. In Fig. S2, we plot |T t 11 | for near-normal (left panel) and not-near-normal (right panel) incidence in constituent plane wave deflection coordinates µ and ν. We see that at not-near-normal incidence, the derivative ∂T t 11 /∂µ κ within the beam angular width and does not depend on ν to the first order. However, at near-normal incidence, the derivative ∂T t 11 /∂µ ∼ κ for wavevectors with ν = 0. Thus, the contribution from the lateral plane waves is indeed the source of anomalous GH shifts at near-normal incidence.

SUPPLEMENTAL NOTE 4: NEAR-NORMAL SHIFTS IN ISOTROPIC MEDIA
In this section, we prove that the near-normal shifts in isotropic homogeneous media do not exhibit anomalous shifts of any kind. We find that all shifts are factorized into two parts: (i) a resonant geometric part and (ii) materialand polarization-dependent part, where the latter part destroys the geometric resonance. The only exception is the angular IF shift which has a maximum for strictly circular input polarization, which is however subwavelength and thus not anomalous. These results further suggest that a finite difference of Fresnel amplitudes for TE and TM waves at normal incidence (as it happens for any anisotropic medium, including low-birefringent uniaxial crystals) would be sufficient condition for the observation of the anomalous shifts. In this case, the geometric resonance term overcomes the material-dependent one and plays a crucial role for the observed optical beam shifts peaks.
Following Eq. (3.13) from Ref. [S1] we write the Fresnel matrix for a homogeneous isotropic medium under small incidence angles ϑ :F . (S18) HereX a p,s = (1/f a p,s )∂ 2 f a p,s /∂ϑ 2 are the normalized second-order derivatives of Fresnel coefficients. Then we derive the Jones matrix relating the incident and secondary fields in the beam coordinate frames according to Eq. (S3): where the "relative Fresnel coefficients' splittings"Ỹ a p,s = 1 − (γ a ) −1 (f a s,p /f a p,s ) are the non-divergent at small angles analogues of Y a p,s =Ỹ a p,s cot ϑ functions in Ref. [S1]. IF shift could be expressed as follows according to Eq. (1): where integration is taken over the transverse momentum components k a X = γ a k 0 µ and k a Y = k 0 ν in the beam coordinate frame. For the light beam incident at the air-material flat interface, we use the Gaussian field distribution with angular spectrum f (µ, ν) = exp[−κ 2 (µ 2 + ν 2 )/2], where κ = k 0 w 0 / √ 2 and w 0 is the beam waist. First, we focus on the matrix integral by ν in the numerator. After matrix and Gaussian differentiation using the product rule, resulting matrix multiplication, and discarding terms linear in ν which bring no contribution to the integral , and the shift is again not peaked, being cubic in small angles (see Fig. S3, panel (b), red curve).
The analogous procedure could be conducted for another types of shifts, namely X , P X and P Y , leading to the similar suppression of the geometric resonance by the material-dependent term proportional to [ϑ 2 + O(ϑ 2 )] at near-normal incidence. Thus, the anomalous optical beam shifts cannot be observed for isotropic medium and require the finite difference between the Fresnel amplitudes at strictly normal incidence.

SUPPLEMENTAL NOTE 5: NEAR-NORMAL SHIFTS IN UNIAXIAL WAVEPLATE
In this section, similarly to the previous one, we briefly go through the calculations of all the shifts for the anisotropic waveplate. The Jones matrix for the waveplate has already been discussed in detail in Supplemental Note 3. As in the main text, we calculate the shifts for the transmitted beam, while the calculations for reflected beam could be done analogously.
Linear IF shift: Y · k 0 . The matrix product under the integral: iT t † (∂T t /∂ν) − iκ 2 ν ·T t †T t . The second term of this expression (taken apart from −iκ 2 ) has already appeared in the calculation of angular IF shift P Y , and thus leads to the proportional term, −iS 2 · τ − · Λ Y (κϑ)/ϑ. The first part leads to the non-vanishing