Interplay between diffraction and the Pancharatnam-Berry phase in inhomogeneously twisted anisotropic media
Chandroth Pannian Jisha, Alessandro Alberucci, Lorenzo Marrucci, Gaetano Assanto
IInterplay between diffraction and the Pancharatnam-Berry phase in inhomogeneouslytwisted anisotropic media
Chandroth P. Jisha, ∗ Alessandro Alberucci,
2, 3, † Lorenzo Marrucci,
4, 5 and Gaetano Assanto
2, 3, 6 Centro de F´ısica do Porto, Faculdade de Ciˆencias,Universidade do Porto, PT-4169-007 Porto, Portugal Optics Laboratory, Tampere University of Technology, FI-33101 Tampere, Finland N oo EL - Nonlinear Optics and OptoElectronics Lab, University “Roma Tre”, IT-00146 Rome, Italy Dipartimento di Fisica, Universit`a di Napoli Federico II, IT-80100 Naples, Italy CNR-ISASI, Institute of Applied Science and Intelligent Systems, IT-80078 Pozzuoli (NA), Italy CNR-ISC, Institute for Complex Systems, IT-00185 Rome, Italy (Dated: May 11, 2019)We discuss the propagation of an electromagnetic field in an inhomogeneously anisotropic mate-rial where the optic axis is rotated in the transverse plane but is invariant along the propagationdirection. In such a configuration, the evolution of an electromagnetic wavepacket is governed by thePancharatnam-Berry phase (PBP), responsible for the appearance of an effective photonic potential.In a recent paper [A. Alberucci et al., “Electromagnetic confinement via spin-orbit interaction inanisotropic dielectrics”, ACS Photonics , 2249 (2016)] we demonstrated that the effective potentialsupports transverse confinement. Here we find the profile of the quasi-modes and show that thephotonic potential arises from the Kapitza effect of light. The theoretical results are confirmed bynumerical simulations, accounting for the medium birefringence. Finally, we analyze in detail aconfiguration able to support non-leaky guided modes. I. INTRODUCTION
In the last few years, a great deal of attention has beendevoted to investigating light propagation in inhomoge-neously anisotropic materials, featuring a rotation of theoptic axis but no variations in the local refractive indices.In such configuration, a phase distribution proportionalto the local rotation angle of the optic axis can be su-perposed to the beam wavefront, leading to the so-calledplanar photonics [1–4]. Such space-dependent phase de-lay is due to a gradient in the Pancharatnam-Berry phase(PBP) [5, 6], the latter being associated with a changein the beam polarization. For a closed path (same initialand final state) the PBP is provided by the solid anglesubtended by the closed circuit on the Poincar´e sphere[7].Historically, the idea of modulating the phase of an elec-tromagnetic wave via the PBP is due to Bomzon andcoworkers, who used a metallic sub-wavelength grating(nowadays it would be called a two-dimensional metama-terial) [8, 9]. This concept was then applied by Marrucciand collaborators to nematic liquid crystals (NLCs), nat-ural materials where the pointwise orientation of the op-tic axis can be engineered by tailoring the boundary con-ditions via optical alignment techniques [10–12]. Thisfull control over the local NLC director paved the wayto new functionalities, for example the conversion of spinangular momentum (polarization degree of freedom) toorbital angular momentum (spatial degree of freedom)[13, 14] and the realization of polarization gratings in ∗ [email protected] † alessandro.alberucci@tut.fi both amplitude [15, 16] and phase [17–19]. Later on, theinterest on PBP was boosted thanks to the introductionof metasurfaces and the demonstration of polarization-dependent optical devices [3, 20], such as lenses [21, 22],holograms [23], vortex plates in polymerizable NLCs [24],deflectors based on the spin-Hall effect [25, 26]. Recently,the PBP has been used to tailor the beam wavefront uponreflection from a layer of chiral liquid crystals [4, 27–29].In this Paper we investigate the propagation of finite-size electromagnetic wavepackets in inhomogeneouslyanisotropic materials, subject to a pointwise rotation ofthe principal axes (Fig. 1) but invariant along the prop-agation direction z . Previously, using geometric opticsother authors analyzed the dependence of the beam tra-jectory on the photon for odd distribution of the rotationangle, the so called optical Magnus effect or spin-Hall ef-fect [30, 31]. For even rotations, the interplay betweendiffraction and PBP was discussed for wide beams as su-perpositions of plane waves in Refs. [32, 33], whereas theapproach we undertake hereby holds valid until longitu-dinal field components become relevant, i.e., for beams ofsize comparable to or narrower than the wavelength [34].In a recent Letter we showed that light propagates underthe influence of an effective photonic potential, leadingto leaky guided modes for bell-shaped distributions ofthe optic axis rotation [35] (Fig. 1). In this Article weprovide a complete and detailed description of the the-ory -including a quantitative comparison with numericalsimulations- first published in [35]. The paper is struc-tured as follows. In Section II we briefly recall the basicsof optical propagation in twisted anisotropic materials inthe absence of diffraction. In Section III we generalize tothe three-dimensional case the paraxial equations govern-ing light propagation found in Ref. [35] to the case of long(with respect to the Rayleigh length) samples. We also a r X i v : . [ phy s i c s . op ti c s ] M a r FIG. 1. (a) Definition of the rotation angle θ . (b) Gaussiandistribution of the rotation angle θ on the plane xy when max-imum rotation is 90 ◦ : the arrows correspond to the local opticaxis (i.e., the extraordinary axis), superposed to the spatialdistribution of θ represented as a color map. When δ = π ,the distribution plotted in (b) corresponds to a polarization-dependent lens [3, 21]. sketch out the comparison with Pauli equation for quan-tum particles. In Section IV we first discuss the physi-cal reason behind the absence of longitudinally-invariantmodes, and then we find the equations governing thequasi-modes, including the higher-order components ofthe localized waves. In Section V we provide an extendedphysical explanation about the origin of the photonic po-tential as a Kapitza effect, i.e. a periodic modulation ofthe PBP versus propagation providing a z -independenteffective photonic potential [36, 37], by using the plane-wave solution shown in Section II. In Section VI first therelationship between the symmetry of the twisting an-gle and light behavior is recalled, and then the profile ofthe CW (continuous wave) component of the quasi-modeis computed for the first time. In Section VII we checkthe validity of the theoretical results using two differ-ent types of numerical simulations, based on an in-houseBPM (Beam Propagation Method) code written in therotated reference system and on an open-access FDTD(Finite-Difference Time Domain) code. Through numer-ical simulations we check the validity of effective pho-tonic potential and -in the guiding case- the accuracy ofthe associated quasi-modes found in Section VI. In Sec-tion VIII we report the design of a non-leaky (V-shaped)waveguide in a twisted structure and test its propertiesvia FDTD simulations. Finally, in Section IX we summa-rize the contents of the Article, pinpointing the noveltieswith respect to Ref. [35] and illustrating future develop-ments, both on the theoretical and on the experimen-tal/technological sides. II. PANCHARATNAM-BERRY PHASE IN ATHIN SAMPLE
We consider a non-magnetic anisotropic uniaxial crys-tal [38]. In the reference system of the principal axes (subscript D ), the relative permittivity is (cid:15) D = (cid:15) ⊥ (cid:15) (cid:107)
00 0 (cid:15) ⊥ , (1)where (cid:15) ⊥ and (cid:15) (cid:107) are the dielectric constants for fields os-cillating normal and parallel to the optic axis ˆ n , respec-tively. The optic axis ˆ n is normal to the wavevector k andat an angle θ ( x, y ) with respect to the axis y [Fig. 1(a)].Noteworthy, (cid:15) ⊥ and (cid:15) (cid:107) are spatially uniform: this con-dition, together with the previous assumption ˆ n · k = 0,rules out any refractive index changes in the transverseplane [39]. We then introduce the ordinary and extraor-dinary refractive indices as n ⊥ = √ (cid:15) ⊥ and n (cid:107) = √ (cid:15) (cid:107) ,respectively. The local rotation of the optic axis (i.e., ofthe principal axes) around z can be expressed as R ( θ ) = cos θ sin θ − sin θ cos θ
00 0 1 . (2)The matrix (2) operates on a vector in the laboratoryframework xyz and yields its coordinates in the Cartesiansystem aligned with the local principal axes. The relativepermittivity (cid:15) in xyz is then provided by (cid:15) ( x, y ) = R [ − θ ( x, y )] · (cid:15) D · R [ θ ( x, y )] . (3)When diffraction in the anisotropic material can be ne-glected, light propagation can be modeled by using planewaves. In fact, the Jones formalism [40] can be employedaccounting for the transmission dependance on the trans-verse coordinate ( x, y ) through the rotation angle θ . Inthe plane-wave approximation, light propagation throughan anisotropic slab of length z is given by (cid:18) E o ( x, y, z ) E e ( x, y, z ) (cid:19) = (cid:18) e ik n ⊥ z e ik n (cid:107) z (cid:19) · (cid:18) E o ( x, y, E e ( x, y, (cid:19) , (4)where E o and E e are the local ordinary and extraordinarypolarization components, respectively [Fig. 1(a)].Since ordinary and extraordinary directions vary with( x, y ), it is more convenient to refer to the two circularlypolarized beams which are eigensolutions of the rotationoperator (2). We thus introduce the circular polarizationbasis ˆ L = (ˆ x − i ˆ y ) / √ R = (ˆ x + i ˆ y ) / √ P (cid:18) E L E R (cid:19) = P · (cid:18) E o E e (cid:19) = 1 √ (cid:18) − i i (cid:19) · (cid:18) E o E e (cid:19) , (5)where the subscripts L and R correspond to LCP andRCP, respectively. Finally, the overall transmission ofthe anisotropic slab is [39] (cid:18) E L ( x, y, z ) E R ( x, y, z ) (cid:19) = e i ¯ nk z (cid:18) cos (cid:0) δ (cid:1) − i sin (cid:0) δ (cid:1) e iθ − i sin (cid:0) δ (cid:1) e − iθ cos (cid:0) δ (cid:1) (cid:19) · (cid:18) E L ( x, y, E R ( x, y, (cid:19) , (6)where δ ( z ) = k z ∆ n is the retardation between ordi-nary and extraordinary components, ∆ n = n (cid:107) − n ⊥ isthe birefringence and ¯ n = ( n (cid:107) + n ⊥ ) / z = 0 by setting, forexample, E R ( x, y,
0) = 1 and E L ( x, y,
0) = 0. When δ = (2 l + 1) π , with l an integer, the RCP wave turns intoLCP: this change in polarization state is accompaniedby a dynamic phase change ∆ φ dyn = (2 l + 1) π ¯ n/ ∆ n anda further phase shift ∆ φ geo = ± θ purely of geometricorigin, a manifestation of PBP [10]. Planar photonics el-ements based on PBP usually work in the latter regime, i.e., with the length of the anisotropic material designedto achieve a π phase delay between ordinary and extraor-dinary components. III. PANCHARATNAM-BERRY PHASE IN ANEXTENDED SAMPLE
Neglecting anisotropy in the diffraction operator [42],Maxwell’s equations in the paraxial approximation read ∇ (cid:18) E x E y (cid:19) + k (cid:18) (cid:15) xx ( x, y ) (cid:15) xy ( x, y ) (cid:15) yx ( x, y ) (cid:15) yy ( x, y ) (cid:19) (cid:18) E x E y (cid:19) = 0 , (7)where the relative permittivity is given by (3). Theparaxial approximation in Eq. (7) allows one to neglectthe longitudinal electric field whenever the beam size ex-ceeds the wavelength [34].We apply the slowly varying envelope approximationthrough the transformation E o = e ik n ⊥ z ψ o and E e = e ik n (cid:107) z ψ e , i.e., factoring out the dynamic phase respon-sible for polarization rotation versus propagation. Forparaxial wavepackets, the second derivatives along z canbe neglected and Eq. (7) yields [35]2 ik n ⊥ ∂ψ o ∂z = −∇ t ψ o + (cid:34)(cid:18) ∂θ∂x (cid:19) + (cid:18) ∂θ∂y (cid:19) (cid:35) ψ o + (cid:18) ∂ θ∂x + ∂ θ∂y (cid:19) ψ e e ik ∆ nz + 2 (cid:18) ∂θ∂x ∂ψ e ∂x + ∂θ∂y ∂ψ e ∂y (cid:19) e ik ∆ nz , (8)2 ik n (cid:107) ∂ψ e ∂z = −∇ t ψ e + (cid:34)(cid:18) ∂θ∂x (cid:19) + (cid:18) ∂θ∂y (cid:19) (cid:35) ψ e − (cid:18) ∂ θ∂x + ∂ θ∂y (cid:19) ψ o e − ik ∆ nz − (cid:18) ∂θ∂x ∂ψ o ∂x + ∂θ∂y ∂ψ o ∂y (cid:19) e − ik ∆ nz , (9)where ∇ t = ∂ x + ∂ y . Equations (8-9) indicate that thewaves are not subject to refractive index gradients, astransverse modulation is only due to the pointwise rota-tion of the principal axes. For small birefringence, i.e., n ⊥ ≈ n (cid:107) , Eqs. (8-9) resemble Pauli equation for a chargedmassive particle moving in a bidimensional space [43] i (cid:126) ∂ ψ ∂t = − (cid:126) m ∇ xy ψ + U ( x, y ) I · ψ + H LS ( x, y, t ) · ψ , (10)with I the identity operator and H LS the Hermitian ma-trix responsible for spin-orbit coupling [44] [45]; ψ isa two-component spinor with elements ψ o and ψ e , re-spectively. The term containing U ( x ) is a scalar po-tential acting equally on both components; the termwith H LS ( x, y, t ), proportional to the Pauli matrix S ,can be associated to an equivalent time-dependent mag-netic field responsible for spin-rotation (in our case powerexchange between extraordinary and ordinary compo-nents). IV. QUASI-MODES IN (1+1)D
Equations (8-9) do not support z -invariant modes dueto the explicit dependence on the propagation coordinate z . This is due to the unavoidable power exchange be-tween ordinary and extraordinary components when theoptic axis rotates across the transverse plane. Notewor-thy, the coupling between the two components does notvanish for any coordinate transformation. To clarify this,let us take a given point P = ( x P , y P ) in the transverseplane with θ P = θ ( P ) and assume the wavepacket to be,e.g., purely extraordinary. Due to diffraction, waveletsin P will spread towards adjacent points P + dP , where θ ( P + dP ) = θ P + dθ due to the optic axis rotation.In P + dP , light will then be in a superposition of or-dinary and extraordinary polarization states, no matterhow small is the change in θ . This can also be explainedthrough a direct analogy with the quantum mechanicsof spinning particles: according to Eq. (10), photons be-have like spin 1/2 particles [46] subject to a magnetic fieldrotating in the plane xy versus the propagation coordi-nate z . Thus, the commutation rule between orthogonalmagnetic fields forbids the existence of stationary/energyeigenstates [43]. Non-Abelian propagation of light in in-homogeneously anisotropic materials was discussed ear-lier in the context of geometric optics [31].From what stated above, no eigenmodes of Eqs. (8-9)exist, as confirmed by direct numerical investigations.Nonetheless, the system could support quasi-modes, i.e.,periodically-varying solutions of finite lateral extension[47, 48]. Hereafter, for the sake of simplicity we con-sider the one-dimensional case with ∂ y = 0. The pe-riodic quasi-modes of Eqs. (8-9) can be calculated withthe ansatz ψ j = g j ( x ) e i (cid:80) ∞ p = −∞ (cid:82) β ( j ) p ( x,z ) dz ( j = o, e ) , (11) where g j ( x ) are z -independent functions, whereas thecomplex exponentials account for spatial variations -in both phase and amplitude - on longitudinal scalesequal or smaller than the beat-length λ/ ∆ n . The pe-riodicity of the quasi-mode is ensured by β ( j ) p ( x, z ) = β ( j ) p ( x ) exp (cid:16) i πp ∆ nλ z (cid:17) [47].To calculate the quasi-modes of a guiding PBP structure,we start from Eqs. (8-9) with the ansatz (11) and find − k n ⊥ g o (cid:88) p β ( o ) p ( x, z ) = − (cid:40) d g o dx + 2 i dg o dx ∂φ o ∂x + g o (cid:34) i ∂ φ o ∂x − (cid:18) ∂φ o ∂x (cid:19) (cid:35)(cid:41) + (cid:18) dθdx (cid:19) g o + (cid:20) d θdx g e + 2 dθdx (cid:18) dg e dx + ig e ∂φ e ∂x (cid:19)(cid:21) e ik ∆ nz e i ( φ e − φ o ) , (12) − k n (cid:107) g e (cid:88) p β ( e ) p ( x, z ) = − (cid:40) d g e dx + 2 i dg e dx ∂φ e ∂x + g e (cid:34) i ∂ φ e ∂x − (cid:18) ∂φ e ∂x (cid:19) (cid:35)(cid:41) + (cid:18) dθdx (cid:19) g e − (cid:20) d θdx g o + 2 dθdx (cid:18) dg o dx + ig o ∂φ o ∂x (cid:19)(cid:21) e − ik ∆ nz e i ( φ o − φ e ) , (13)where we introduced φ j ( x, z ) = (cid:80) p (cid:82) β ( j ) p ( x, z ) dz . Aquasi-mode defined by (11) conserves its profile in prop-agation provided that dβ ( j )0 dx = 0 is satisfied.Equations (12-13) contain periodic terms on both sides,thus each harmonic can be equalized separately. Thecomputation is much easier if we assume that φ o − φ e is constant versus z : after the derivation we will verify a posteriori that the latter condition is satisfied to firstorder.Focusing on the cw components p = 0, the latter satisfythe eigenvalue problem − β ( o )0 k n ⊥ g o = − d g o dx + (cid:34)(cid:18) dθdx (cid:19) + 1 k (∆ n ) ∞ (cid:88) p =1 p dβ ( o ) p dx dβ ( o ) − p dx (cid:35) g o , (14) − β ( e )0 k n (cid:107) g e = − d g e dx + (cid:34)(cid:18) dθdx (cid:19) + 1 k (∆ n ) ∞ (cid:88) p =1 p dβ ( e ) p dx dβ ( e ) − p dx (cid:35) g e . (15) Thus Eqs. (14-15) show that the terms coupling ex-traordinary and ordinary waves in Eqs. (8-9) act on thecw components g j ( x ) ( j = o, e ) as additional contri-butions (summations over p ) to the photonic potential.Such higher-order contributions, however, from a physi-cal point of view are expected to be relevant only whenthe beat-length L B = λ/ ∆ n exceeds the Rayleigh dis-tance L R = πnw m /λ , with w m the transverse extent ofthe localized mode.When higher-order contributions can be neglected (i.e., L B /L R (cid:28) V ( x ) = 12 n j k (cid:18) dθdx (cid:19) . (16)We need now to evaluate the degree of accuracy in us-ing the potential (16) into Eqs. (14-15). Equations (12-13) also determine the profile of the functions β ( j ) p ( x )for p (cid:54) = 0, yielding the periodic evolution of the quasi-mode along z . The knowledge of β ( j ) p ( x ) allows evaluat-ing the weight of the term ∝ (cid:80) ∞ p =1 1 p dβ ( o ) p dx dβ ( o ) − p dx on thephotonic potential. For the sake of simplicity, assumingnon-negligible β ( j ) p only for | p | ≤
1, for p = 1 we get − k n ⊥ g o β ( o )1 = − λπ ∆ n dg o dx dβ ( o )1 dx − λ π ∆ n g o d β ( o )1 dx + g e d θdx + 2 dg e dx dθdx , (17) − k n (cid:107) g e β ( e )1 = − λπ ∆ n dg e dx dβ ( e )1 dx − λ π ∆ n g e d β ( e )1 dx . (18)Similarly, for p = − − k n ⊥ g o β ( o ) − = λπ ∆ n dg o dx dβ ( o ) − dx + λ π ∆ n g o d β ( o ) − dx , (19) − k n (cid:107) g e β ( e ) − = λπ ∆ n dg e dx dβ ( e ) − dx + λ π ∆ n g e d β ( e ) − dx − g o d θdx − dg o dx dθdx . (20)In the general case the functions β ( j ) p satisfy a second-order ordinary differential equation with coefficientsgiven by the cw component g j ( x ). A drastic approxima-tion is at hand for large ∆ n , when Eqs. (17-20) provide β ( e )1 = β ( o ) − = 0, β ( o )1 = − (cid:16) g e d θdx + 2 dθdx dg e dx (cid:17)(cid:46) (2 k n ⊥ g o )and β ( e ) − = (cid:16) g o d θdx + 2 dθdx dg o dx (cid:17)(cid:46) (cid:0) k n (cid:107) g e (cid:1) . Substitutingback into Eqs. (14-15), we find that the photonic poten-tial is proportional to (cid:0) dθdx (cid:1) + O (∆ n − ); as an immediateconsequence, at this order of approximation [the approx-imation includes to set n ⊥ ≈ n (cid:107) on the LHS of Eqs. (14-15)] g o = g e = g . The modulation of the beam along z isnow given by β ( o )1 = − d θdx + 2 dθdx d ln gdx k n ⊥ (21)and β ( e ) − = d θdx + 2 dθdx d ln gdx k n (cid:107) , (22)i.e., the quasi-mode undergoes periodic modulation ofboth amplitude and phase profiles. Since β ( e ) − ≈− β ( o )1 , φ o − φ e is constant in first approximation, inagreement with the initial hypothesis. According toEq. (21) it is β ( e ) − ∝ (cid:0) w m k (cid:1) − , in turn providing k (∆ n ) (cid:80) ∞ p =1 1 p dβ ( e ) p dx dβ ( e ) − p dx ∝ (cid:16) L B L R (cid:17) ; the latter resultconfirms that the ratio between the beating length andthe Rayleigh length is the smallness parameter in ourapproximated treatment. V. ORIGIN OF THE EFFECTIVE PHOTONICPOTENTIAL
Since we factored out the dynamic phase when intro-ducing the fields ψ o and ψ e , we can speculate that the FIG. 2. Pancharatnam-Berry phase difference ∆ φ geo versusretardation δ between two regions differing only by a relativerotation θ when spatial dispersion is neglected, i.e., in theplane wave limit. The blue solid lines are the exact resultscomputed from Eq. (23); the red dashed lines are a sinusoidalapproximation for ∆ φ geo sharing the same peak 2 θ of theexact solutions. photonic potential stems from the geometric phase. In or-der to prove this, in analogy with standard GRIN guidesand splitting operators in numerical analysis, let us sep-arate Eq. (7) in two portions, one accounting for diffrac-tive spreading (second derivative along x ) and one forthe inhomogeneous dielectric tensor, respectively. Weguess that the transverse phase distribution from thesecond part is able to compensate diffraction after av-eraging along one beat-length λ/ ∆ n . Owing to the ab-sence of diffraction, the solutions of the second part ofEq. (7) correspond to Eqs. (6), i.e., the exact solutionsin the plane wave limit. On the one hand, accordingto Eq. (6), when a purely circular polarization (eitherLCP or RCP) is launched, the dynamic phase accumu-lated in propagation is k nz . This contribution is clearlyconstant in the transverse plane, thus cannot be respon-sible for the appearance of an x -dependent photonic po-tential. On the other hand, the terms between squarebrackets in Eq. (6) provide a phase delay of geomet-ric origin (i.e., due to variations in light polarization),which we will analyze hereafter. We adopt Pancharat-nam’s original approach [50, 51] to calculate the phasedifference ∆ φ geo between two different polarizations as∆ φ geo = arg [ E ( θ , z ) · E ∗ ( θ , z )] (the superscript ∗ indi-cates the complex conjugate). For a purely RCP wave in z = 0, Eqs. (6) provide∆ φ geo = arg (cid:20) cos (cid:18) δ (cid:19) + sin (cid:18) δ (cid:19) e i ( θ − θ ) (cid:21) . (23)Fig. 2 graphs the geometric phase difference (23) versus δ in the plane-wave limit: ∆ φ geo is periodic with λ/ ∆ n andoscillates between the two extrema 0 and 2 ( θ − θ ). Theoscillation shape versus propagation z depends stronglyon its amplitude: it is sinusoidal for small amplitudes (theexponential can be Taylor-expanded retaining only thelinear term), whereas it is flat-topped when θ − θ = 90 ◦ ;for intermediate values of θ − θ there is a gradual tran-sition between these two limits. Fig. 2 (red dashed lines)plots the sinusoidal approximation of Eq. (23) to illus-trate the discrepancy with the exact form for various rel-ative rotations θ − θ .Equation (23) is the key to understanding the evolutionof a finite-size beam: the beam wavefront acquires a con-tinuous geometric phase-delay ∆ φ geo , which is periodicwith propagation z . Conversely, the amplitude of thesephase-oscillations depends on the optic axis rotation be-tween neighboring points in the transverse plane, i.e., onthe θ distribution and its derivatives across x . Throughthe Kapitza effect of light, a periodic index modula-tion of the form ∆ n = n ( x ) − n = f ( z ) W ( x ) with f ( z ) = (cid:80) l f l exp (2 iπlz/ Λ) yields a z -independent effec-tive photonic potential [37] V Kap = − k n ( x ) − n n = k Λ nπ (cid:32) ∞ (cid:88) l = −∞ f l f − l l (cid:33) (cid:18) dWdx (cid:19) . (24)In essence, V Kap accounts for fast scale modulations onthe continuous wave (cw) component: the rapid oscil-lations in phase distribution entail a modulation of thetransverse momentum k x , the latter providing a net cu-mulative phase across the wavefront due to the effectivekinetic energy k x / (2 k n j ) [36, 37].To further support the interpretation above, we recallthat the PBP is non-transitive [7]: at any given z , thephase difference between two points x and x with ro-tation equal θ = θ ( x ) and θ = θ ( x ), respectively,depends on the polarizations state of the beam in thewhole interval [ x x ] or, otherwise stated, on its pathon the Poincar´e sphere. This is in sharp contrast withthe standard dynamic phase, where the phase differencedepends exclusively on initial and final state. Thus, theaccumulated phase difference needs to be computed tak-ing adjacent points in space. Fig. 2 shows that the pro-file (i.e., the Fourier coefficients f l ) of the periodic func-tion ∆ φ geo along z [given by Eq. (23)] markedly dependson the relative rotation θ = θ − θ , providing distinctKapitza potentials V Kap . In the limit | x − x | → θ , Eq. (23) gives f ( z ) = sin (cid:0) π ∆ nzλ (cid:1) and W ( x ) = 2 n ∆ nθ ( x ). Such ap-proach provides an effective z -invariant potential V = nk (cid:0) dθdx (cid:1) , the correct transverse profile but a factor 2smaller than Eq. (16). Discrepancy arises from the factthat the Kapitza model, elaborated in Ref. [37] for thescalar case, does not account for the complete spin-orbitinteraction occurring in this case [52]. VI. FOCUSING AND DEFOCUSINGPOTENTIALS
According to Eq. (16), the shape of the effective pho-tonic potential strongly depends on the symmetry of θ ( x ). When θ is bell-shaped [solid blue line in Fig. 3(a)],the potential is M-like [dashed green line in Fig. 3(a)] FIG. 3. (a) Rotation angle θ (blue solid line) and corre-sponding photonic potential (dashed green line) versus x fora Gaussian-shaped rotation of principal axes. Here θ = π/ w θ = 5 µ m. (b) As in (a) but when θ has a hyperbolictangent distribution. Here L θ = 5 µ m. When δ = π , (a) and(b) behave like polarization-dependent lens and deflector, re-spectively. and supports leaky modes [37]. Such modes can be com-puted with good accuracy by clipping off the edges of thepotential, i.e., approximating the M-well with a V-profile[53]. For a Gaussian distribution of θθ ( x ) = θ e − x w θ , (25)Fig. 3(a) plots the associated potential V = n j k (cid:104) xθ w θ exp (cid:16) − x w θ (cid:17)(cid:105) . Fig. 4(a) graphs the resultingquasi-modes versus the maximum rotation θ . The quasi-modes computed from (16) hold valid when the higherorder contributions -i.e., the sum over p in Eqs. (14-15)-are negligible, that is, for large enough birefringence. InSection VII A BPM simulations will be employed to checkthe behavior versus the birefringence ∆ n . On the otherhand, these results are valid as long as the mode remainsconfined within the central lobe of the potential [see solidand dashed lines in Fig. 4(b)], regardless of the birefrin-gence ∆ n [37].When θ ( x ) is odd symmetric [solid blue line in Fig. 3(b)],the photonic potential is maximum in the center (around x = 0) and minimum at the edges: light gets repelledfrom the region around x = 0 and no lateral confinementoccurs. For example, if θ is given by θ ( x ) = θ tanh (cid:18) xL θ (cid:19) , (26)the effective potential is V = n j k (cid:104) θ L θ (cid:16) − tanh (cid:16) xL θ (cid:17)(cid:17)(cid:105) [see the dashed green linein Fig. 3(b)]. VII. NUMERICAL SIMULATIONS
In this Section we validate the theory developed aboveby means of numerical simulations. First, we useBPM simulations to verify that the photonic potential
FIG. 4. (a) Intensity profile for a Gaussian distribution of θ and w θ = 5 µ m; each curve is labelled with the correspondingmaximum rotation θ . (b) Corresponding width of the quasi-mode versus maximum rotation angle θ for three different w θ , as labelled (solid lines). The dashed lines graph the halfwidth of the angular distribution: above them the model isnot accurate due to a finite overlap of the mode with the edgesof the photonic potential [37]. Here we solved Eqs. (14-15)considering only the term given by Eq. (16) when computingthe quasi-mode; we assumed n ⊥ = 1 . n (cid:107) = 1 . y and x ,respectively. In (a) and (b) the transverse distribution of θ is given by (25) and (26), respectively. Here θ = 360 ◦ , w θ = L θ = 5 µ m, n ⊥ = 1 . n (cid:107) = 1 .
7, the wavelength is1064 nm. The input is a Gaussian beam of width 3 µ m andflat-phasefront. (16) describes accurately light propagation in a twistedanisotropic medium. Then, we use FDTD simulations toassess the paraxial approximation invoked in going fromEq. (7) to Eqs. (8-9). With FDTD simulations we fi-nally address the accuracy of Eq. (7) in lieu of the exactMaxwell’s equations. A. BPM simulations
The presence of the effective potential given byEq. (16) can be tested by simulating Eqs. (8-9), validin the paraxial limit. To simulate light evolution in the rotated reference system, we use a BPM code based onoperator splitting and the Crank-Nicolson algorithm fordiffraction.First, we check that twisted anisotropic materials ex-hibit a quasi-isotropic response for propagation lengthsmuch longer than the Rayleigh distance and large enoughbirefringence. The results are shown in Fig. 5 for even[Eq. (25) and Fig. 3(a)] and odd symmetry potentials[Eq. (26) and Fig. 3(b)], with n ⊥ = 1 . n = 0 . δ = π , the even and odd distributions correspondto a polarization-dependent lens [21] and a polarization-dependent deflector based on the spin-Hall effect [25], re-spectively. As predicted by our theory, in both cases theelectromagnetic propagation is nearly independent of theinput polarization. If the θ distribution is bell-shaped,diffractive spreading is counteracted and light is laterallytrapped, i.e., waveguiding takes place. If θ is odd, lightis repelled from the symmetry axis x = 0. Summarizing,Fig. 5 confirms qualitatively that light is subject to aneffective potential given by (16).Figure 6 displays light propagation for a Gaussian dis-tribution of θ [see Eq. (25) and Fig. 3(a)] versus themaximum rotation θ , with n ⊥ = 1 . n = 0 . θ increases the beamsundergo a stronger transverse confinement due to the M-shaped potential. At the same time, the figure showshow a purely y -polarized input beam couples a fractionof its power into the x -polarization owing to the point-wise rotation of the optic axis. Regardless of θ , strongcoupling to radiation is observed due to the small spatialoverlap between a linearly polarized Gaussian beam andthe leaky modes.Next, we check that our analytical model is in quanti-tative agreement with the actual solutions, as well. Itis also important to study how electromagnetic propaga-tion depends on the birefringence ∆ n . To that extent,at the input we injected the quasi-modes computed fromthe eigenvalue problem, Eqs. (12-13), in the limit of alarge birefringence ∆ n , i.e., when only the action of thephotonic potential provided by Eq. (16) is considered,see Fig. 4. By means of Eqs. (21-22), our theory predictsthat the second term between square brackets in Eqs. (12-13) goes as (∆ n ) − , and thus can be neglected for largeenough birefringence. Figure 7 illustrates the quasi-modeevolution in the medium for n ⊥ = 1 . n > .
1, guiding is clearly observedfor both ψ o and ψ e : this implies that, for ∆ n > . n ⊥ and n (cid:107) is large enough to be appreciable even on the LHS ofEqs. (12-13), yielding a slightly anisotropic response onlong scales (with respect to the beat length λ/ ∆ n ). Forbirefringence < .
1, higher order terms in Eqs. (12-13)are relevant and the differences between the exact pho-tonic potential and its approximation given by Eq. (16)
FIG. 6. Intensity distribution in the plane xz versus maximum rotation angle θ , computed with BPM simulations when θ isGaussian with w θ = 5 µ m, see Eq. (25); here the refractive indices are n ⊥ = 1 . n (cid:107) = 1 . y -polarized Gaussian beamof waist 3 µ m is launched at the input. In the first and second rows | ψ e | and | ψ o | are plotted, respectively. In the third andfourth rows the corresponding fields | E y | and | E x | are graphed in the laboratory coordinate system, respectively. cannot be neglected anymore: as a net result, guidanceis lost and a polarization-dependent behavior arises. Inparticular, the beam remains partially trapped for bire-fringence down to ∆ n ≈ .
03 and ∆ n ≈ .
01 for ordi-nary and extraordinary inputs, respectively. Finally, thepower coupled to guided modes strongly depends on theoverlap between the input and the quasi-mode, as shownby a direct comparison between Figs. 7 and 6. We notethat the overlap does not depend solely on the spatialdistribution of the optical wavepacket, but also on itslocal polarization: in other words, the quasi-modes arestructured light.
B. FDTD simulations
We now check the previous results by FDTD simu-lations, solving directly the Maxwell’s equations in thetime domain and without the paraxial approximation.We used the open-source code MEEP [54], taking acontinuous-wave excitation at the vacuum wavelength λ = 1 µ m, although our findings hold valid regardlessof the frequency. The source was a Gaussian-shaped col-lection of dipoles, of width 3 µ m across x , infinitesimallynarrow along z and centered in x = z = 0. The emittedradiation from the source was uniformly linearly polar-ized. The inhomogeneous uniaxial medium started in z = 2 µ m.We first verify that the bell-shaped and odd θ distribu-tions correspond to a trapping and a repulsive dynamics,respectively. We will use n ⊥ = 1 . n (cid:107) = 1 . n = 0 .
2) until otherwise specified. Fig. 8(a) illustratesthe wavepacket evolution when the rotation angle is givenby (26) [see Fig. 3(b)]: light is expelled from the centralregion. Fig. 8(b) graphs the confining case correspondingto (25) [see Fig. 3(a)]. In both cases, the FDTD resultsare in excellent agreement with the BPM simulations inFig. 5.We then concentrate on the trend of the guiding effect [ θ satisfying Eq. (25)] versus the maximum rotation angle θ , see Fig. 9. The confinement is strongly enhanced asthe maximum rotation increases, reaching a good degreeof confinement for θ = 360 ◦ . A comparison with Fig. 6is in excellent agreement with the BPM results.Quantitative details on the evolution of the inputwavepacket are gathered by computing its width w η =2 (cid:113)(cid:82) L − L x Π η ( x ) dx ( η = I, E j ) versus propagation z , withΠ I = I/ (cid:82) Idx and Π E j = | E j | / (cid:82) | E j | dx ( j = o, e )the probability densities for intensity and electric fields,respectively. The integrals are limited to the domain[ − L L ] in order to get rid of the diffracting compo-nents generated at the input interface and then radiatedout. Such densities allow evaluating how power shares
FIG. 7. Columns 1 to 5: Intensity evolution in the plane xz when the mode plotted in Fig. 4 is launched at the input z = 0,for different values of n (cid:107) as labelled (the corresponding birefringence ∆ n is 0.001, 0.01, 0.05, 0.2 and 0.4, from left to rightrespectively). Top and bottom rows correspond to local ordinary and extraordinary polarizations, respectively. Column 6:modulus of the field z = 500 µ m versus x and the refractive index n (cid:107) . Here θ is Gaussian with θ = 360 ◦ and w θ = 5 µ m, n ⊥ = 1 . θ obeying Eq. (25)] and focusing [panelb, θ distributions as provided by Eq. (26)] effective potential.The input beam is linearly polarized along y , n (cid:107) = 1 . w θ = L θ = 5 µ m and θ = 360 ◦ . The white dashed line indicatesthe input interface, the red solid lines show beam diffractionfor θ = 0 ◦ . between the two orthogonal polarizations. The secondrow in Fig. 10 shows beam size versus z and versus themaximum rotation θ , respectively. The beam width w I initially broadens due to the strong radiation emittedat the input interface and the coupling into the guide.Over longer propagation distances the wavepacket nar-rows down and stabilizes to a width corresponding tothe leaky mode provided by Eq. (16). We stress that, af- ter the radiation excited at the input interface fades out,another kind of coupling to radiation occurs, inherent toleaky modes in M-shaped waveguides [53]. The corre-sponding losses strongly depend on θ , becoming negligi-ble after several Rayleigh distances for modes narrowerthan the guiding core of the M-guide. The two differentcontributions to radiation can be discriminated more eas-ily by BPM simulations, as in that case the input can betailored to match the effective mode. The role of the over-lap in z = 0 can be evaluated by comparing Fig. 6 and 7for n (cid:107) = 1 . n = 0 .
2) and θ = 360 ◦ :when a Gaussian beam is launched, at the input interfacetwo additional beams are emitted, tilted to left and rightof the impinging wave vector by the same angle. Theseside lobes are strongly dampened when the quasi-modeis launched. The radiation loss in the bulk, related withthe leaky nature of the quasi-mode, can be better appre-ciated for θ = 180 ◦ and n (cid:107) = 1 .
7. In fact, once radiationfrom the interface has moved away from the central re-gion, the amplitude of the central lobe slightly decreaseswith z , both in BPM (Fig. 6) and in FDTD simulations(Fig. 9).We finally investigate light propagation versus in-put polarization and birefringence, the results shown inFig. 11. For ∆ n > .
2, the propagation is almost inde-pendent from the input polarization; at lower birefrin-0
FIG. 9. FDTD simulations versus θ . Evolution of | E x | (first row), | E y | (second row) and the overall intensity | E x | + | E y | (third row) in the plane xz . The input beam is linearly polarized along y .FIG. 10. Direct comparison between analytic theory andFDTD simulations when the input is linearly polarized along y , θ has is Gaussian and w θ = 5 µ m. Upper panels: overall in-tensity cross-sections calculated at z = 120 µ m for a Gaussianwavepacket undergoing diffraction when θ = 0 ◦ (solid blacklines) and the dielectric tensor is transversely rotated (solidblue lines). The red dashed lines are the corresponding quasi-modes from Eqs. (14-15). Lower panels: wavepacket widthversus propagation distance z using Π E y (blue solid line) andΠ I (red solid lines) for the probability densities. The blackdashed line corresponds to a diffracting Gaussian beam in thelimit θ = 0 ◦ . Here L = 10 µ m. gence, the propagation is polarization-dependent. Thesefindings match the BPM results, see e.g., Fig. 7. When n (cid:107) = 1 .
501 and n (cid:107) = 1 .
51 (corresponding to ∆ n = FIG. 11. FDTD simulations versus birefringence for a fixed θ profile. The beam is linearly polarized along x (left column)and y (right column). Birefringence ∆ n is 0.001, 0.01, 0.2 and0.4, from top to bottom row respectively. The red solid linesshow the diffracting case, the dashed white lines correspond tothe initial section of the anisotropic medium. Here n ⊥ = 1 . θ = 360 ◦ and w θ = 5 µ m. .
001 and ∆ n = 0 .
01, respectively), light propagationis polarization-dependent both in BPM and FDTD sim-ulations, owing to the action of the higher order con-1
FIG. 12. V-shaped PBP waveguide. (a) Rotation angle θ (blue solid line) and corresponding potential V from Eq. (27)(green dashed line) for x = 10 µ m and a = 4 × m − . (b)Mode width versus a for (top to bottom curves) x = 5 µ m(blue line), 10 µ m (green) and 20 µ m (red). The starsgraph the width w G of the Gaussian fundamental mode ofthe parabolic potential in the limit x → ∞ . Inset: modeintensity profile when a = 2 × m − for x = 5 µ m (bluesolid line) and x = 20 µ m (red solid line). Blue stars repre-sent the Gaussian mode when x → ∞ . tributions, the latter modifying the shape of the pho-tonic potential with respect to the approximated formulaEq. (16). In both cases, light trapping is stronger forinput beams polarized along y than along x . Specifi-cally, for n (cid:107) = 1 .
51 input wavepackets polarized along y are confined while the orthogonal polarization spreadsout, whereas for n (cid:107) = 1 .
501 wavepacket widens irrespec-tively of the input polarization. Noteworthy, when BPM(Fig. 7) and FDTD (Fig. 11) are compared, it is impor-tant to recall that propagation lengths in the former aremuch larger than in the latter.
FIG. 13. V-shaped PBP waveguide. (a-b) FDTD simulationsin the plane xz for a = 0 , × and 3 × m − , re-spectively (slices from bottom to top) when (a) x = 2 µ mand (b) x = 10 µ m. The red dashed lines correspond to thestraight lines x = ± | x | , i.e., the waveguide edges. VIII. NON-LEAKY WAVEGUIDES
Equation (16) allows tailoring the profile θ ( x ) of thecontinuous rotation to yield a bell-shaped effective wellable to support non-leaky guided modes. Since the low-est value of the potential is zero, for a non-leaky electro-magnetic waveguide the derivative dθ/dx must vanish in x = 0, corresponding to a local extremum of θ . There-fore, θ ( x ) needs to be an even function. The condition V ( x ) > | x | → ∞ must be satisfied to get evanes-cent waves beyond the waveguide edges. Moreover, arealistic waveguide must be (transversely) finite; hence,we have to impose a linear trend θ ( x ) = α (cid:0) | x | − x ♦ (cid:1) )for | x | > x , with x linked to the width of the guide.Finally, for | x | < x we need to select a convex even func-tion θ so that the potential (16) has a relative minimumin x = 0. The simplest choice is a parabolic profile of the2form ax , such that the resulting photonic potential V is V ( x ) = (cid:40) a x n j k for | x | ≤ x , a x n j k for | x | > x , (27)where we satisfied the continuity conditions on θ ( x ) and dθ/dx by setting x ♦ = x / α = 2 ax , respectively.An example is graphed in Fig. 12(a).In the limit x → ∞ the photonic potential is purelyparabolic, thus supporting an infinite number of guidedmodes, the fundamental one being Gaussian with fieldproportional to e − x /w G and width w G = 1 / √ a . The fun-damental mode of the PBP waveguide defined by Eq. (27)corresponds to such solution for widths w G much smallerthan x , as confirmed in Fig. 12(b) graphing the ex-act modes of the potential Eq. (27). We analyze thebeam propagation in such PBP waveguide by means ofFDTD simulations: a synopsis of the results is availablein Fig. 13(a-b). As the coefficient a increases startingfrom zero, the wavepacket undergoes transverse confine-ment. In agreement with Fig. 12(b), narrower waveguides(i.e., smaller x ) yield lower field confinement. IX. CONCLUSIONS AND PERSPECTIVES
In this Paper we examined the propagation of electro-magnetic waves in inhomogeneously anisotropic media,locally twisted in the transverse plane but invariant inpropagation. First, we expanded the results reported inRef. [35], providing a deeper insight on the physical prin-ciples and extending the theoretical computations. Asa matter of fact, we showed that light propagates un-der the influence of an effective photonic potential dueto the spin-orbit interaction, which originates from theperiodic modulation of the Pancharatnam-Berry phaseversus propagation via the Kapitza effect. Noteworthy,this photonic potential features a quasi-isotropic behav-ior, despite the locally anisotropic response of the mate-rial. The attracting or repelling character of the poten-tial strictly depends on the symmetry of the optic axisrotation: in particular, guiding structures correspond tobell-shaped distributions. Due to the inherent power ex-change between the ordinary and the extraordinary com-ponents, the system does not support z -invariant modes[30] but quasi-periodic modes, in agreement with the Flo-quet theorem. While in general the guided modes areleaky owing to an M-shaped potential, we designed a dis-tribution yielding a V-shaped potential and truly guidedmodes.With respect to Ref. [35], all of these results have been(qualitatively and quantitatively) validated by means ofBPM and FDTD numerical simulations. In particular,here we used for the first time a BPM written in the ro- tated framework, specifically developed for this problem.We also demonstrated a very good agreement betweenBPM and FDTD: this is a very important result giventhat BPM codes are much less demanding than FDTDfrom the computational point of view. Our new BPMcode will be thus relevant for the numerical analysis ofq-plates and similar devices, nowadays mainly describedby means of simple plane-wave models. Both the numer-ical techniques confirmed that our theory provides a verygood approximation of the confined mode for a vast rangeof parameters. With respect to the latter statement, wealso explored the limits of our analytical model versusthe medium birefringence and maximum rotation angle,both for weakly and highly anisotropic materials.Our findings pave the way to the design and realiza-tion of a new kind of electromagnetic waveguides, withthe presented mechanism valid regardless of frequency.Twisted structures can be realized, for example, usingliquid crystals [4, 10, 55], laser-written anisotropic glasses[26], metamaterials [3, 5, 56]. With respect to Ref. [39],the geometry proposed hereby is much simpler to realizeas it is structured only in the transverse plane. Fur-thermore, the longitudinal homogeneity of the sampleinhibits back reflections, an important detrimental ef-fect when fabricating high-quality waveguides. Our re-sults can be generalized to more complicated and exoticanisotropic responses, including, e.g., magneto-electriccoupling [57–63]. Several future developments can be en-visaged, including the systematic investigation of lighttrapping versus geometrical and material parameters ofthe structure, and the usage of our model to find the pe-riodic components of the quasi-mode. It will be also in-teresting to investigate the interplay between the Rytov-Vladimirskii-Berry and the Pancharatnam-Berry phaseswhich are both of geometric origin [44], and the connec-tion of our system with the recently introduced gaugeoptics [64, 65].From a more fundamental point of view, this works isan important contribution to the analogy between lightpropagation in anisotropic materials and propagation ofspin 1/2 particles in magnetic fields: in fact, confinementbased on geometric phase can be readily transposed tomatter waves, e.g., charged particles in a Paul trap [66]. ACKNOWLEDGMENTS
C.P.J. gratefully acknowledges Funda¸c˜ao para aCiˆencia e a Tecnologia, POPH-QREN and FSE (FCT,Portugal) for the fellowship SFRH/BPD/77524/2011;she also thanks the Optics Laboratory in Tampere for thehospitality. A.A. and G.A. thank the Academy of Fin-land through the Finland Distinguished Professor grantno. 282858. L.M. acknowledges support from the Euro-pean Research Council (ERC), under grant No. 694683,PHOSPhOR.3 [1] Alexander V. Kildishev, Alexandra Boltasseva, andVladimir M. Shalaev, “Planar photonics with metasur-faces,” Science (2013), 10.1126/science.1232009,http://science.sciencemag.org/content/339/6125/1232009.full.pdf.[2] Nanfang Yu and Federico Capasso, “Flat optics with de-signer metasurfaces,” Nat. Mater. , 139–150 (2014).[3] Dianmin Lin, Pengyu Fan, Erez Hasman, and Mark L.Brongersma, “Dielectric gradient metasurface optical el-ements,” Science , 298–302 (2014).[4] Junji Kobashi, Hiroyuki Yoshida, and Masanori Ozaki,“Planar optics with patterned chiral liquid crystals,” Nat.Photon. , 389–392 (2016).[5] Ze’ev Bomzon, Gabriel Biener, Vladimir Kleiner,and Erez Hasman, “Space-variant Pancharatnam–Berryphase optical elements with computer-generated sub-wavelength gratings,” Opt. Lett. , 1141–1143 (2002).[6] L. Marrucci, C. Manzo, and D. Paparo, “Pancharatnam–Berry phase optical elements for wavefront shaping inthe visible domain: switchable helical modes generation,”Appl. Phys. Lett. , 221102 (2006).[7] M.V. Berry, “The adiabatic phase and Pancharatnam’sphase for polarized light,” J. Mod. Opt. , 1401–1407(1987).[8] Zeev Bomzon, Vladimir Kleiner, and Erez Hasman,“Formation of radially and azimuthally polarized lightusing space-variant subwavelength metal stripe grat-ings,” Appl. Phys. Lett. , 1587–1589 (2001).[9] Erez Hasman, Vladimir Kleiner, Gabriel Biener, andAvi Niv, “Polarization dependent focusing lens by use ofquantized pancharatnamberry phase diffractive optics,”Appl. Phys. Lett. , 328–330 (2003).[10] L. Marrucci, C. Manzo, and D. Paparo, “Optical spin-to-orbital angular momentum conversion in inhomogeneousanisotropic media,” Phys. Rev. Lett. , 163905 (2006).[11] S. Slussarenko, A. Murauski, T. Du, V. Chigrinov,L. Marrucci, and E. Santamato, “Tunable liquid crystalq-plates with arbitrary topological charge,” Opt. Express , 4085–4090 (2011).[12] Jihwan Kim, Yanming Li, Matthew N. Miskiewicz, Chul-woo Oh, Michael W. Kudenov, and Michael J. Escuti,“Fabrication of ideal geometric-phase holograms with ar-bitrary wavefronts,” Optica , 958–964 (2015).[13] Thomas Bauer, Peter Banzer, Ebrahim Karimi,Sergej Orlov, Andrea Rubano, Lorenzo Mar-rucci, Enrico Santamato, Robert W. Boyd, andGerd Leuchs, “Observation of optical polariza-tion M¨obius strips,” Science , 327–332(2016).[15] T. Todorov and L. Nikolova, “Spectrophotopolarimeter:fast simultaneous real-time measurement of light param-eters,” Opt. Lett. , 358–359 (1992).[16] Franco Gori, “Measuring stokes parameters by means ofa polarization grating,” Opt. Lett. , 584–586 (1999).[17] C. Provenzano, P. Pagliusi, and G. Cipparrone, “Elec-trically tunable two-dimensional liquid crystals gratingsinduced by polarization holography,” Opt. Express , 5872–5878 (2007).[18] Nelson V. Tabiryan, Sarik R. Nersisyan, Timothy J.White, Timothy J. Bunning, Diane M. Steeves, andBrian R. Kimball, “Transparent thin film polarizing andoptical control systems,” AIP Advances , 022153 (2011),http://dx.doi.org/10.1063/1.3609965.[19] U. Ruiz, P. Pagliusi, C. Provenzano, and G. Cipparrone,“Highly efficient generation of vector beams through po-larization holograms,” Appl. Phys. Lett. , 161104(2013), http://dx.doi.org/10.1063/1.4801317.[20] Natalia M. Litchinitser, “Photonic multitasking enabledwith geometric phase,” Science , 1177–1178 (2016),http://science.sciencemag.org/content/352/6290/1177.full.pdf.[21] Filippus S. Roux, “Geometric phase lens,” J. Opt. Soc.Am. B , 476–482 (2006).[22] Mohammadreza Khorasaninejad, Wei Ting Chen,Robert C. Devlin, Jaewon Oh, Alexander Y. Zhu,and Federico Capasso, “Metalenses at visible wave-lengths: Diffraction-limited focusing and subwavelengthresolution imaging,” Science , 1190–1194 (2016),http://science.sciencemag.org/content/352/6290/1190.full.pdf.[23] Guoxing Zheng, Holger Mhlenbernd, Mitchell Kenney,Guixin Li, Thomas Zentgraf, and Shuang Zhang, “Meta-surface holograms reaching 80% efficiency,” Nat. Nan-otech. , 308–312 (2015).[24] Sarik Nersisyan, Nelson Tabiryan, Diane M. Steeves, andBrian R. Kimball, “Fabrication of liquid crystal polymeraxial waveplates for uv-ir wavelengths,” Opt. Express ,11926–11934 (2009).[25] Guixin Li, Ming Kang, Shumei Chen, ShuangZhang, Edwin Yue-Bun Pun, K. W. Cheah, andJensen Li, “Spin-enabled plasmonic metasurfaces formanipulating orbital angular momentum of light,”Nano Lett. , 4148–4151 (2013), pMID: 23965168,http://dx.doi.org/10.1021/nl401734r.[26] Xiaohui Ling, Xinxing Zhou, Xunong Yi, Weixing Shu,Yachao Liu, Shizhen Chen, Hailu Luo, Shuangchun Wen,and Dianyuan Fan, “Giant photonic spin hall effect inmomentum space in a structured metamaterial with spa-tially varying birefringence,” Light Sci. Appl. , e290(2015).[27] Mushegh Rafayelyan, Georgiy Tkachenko, and EtienneBrasselet, “Reflective spin-orbit geometric phase fromchiral anisotropic optical media,” Phys. Rev. Lett. ,253902 (2016).[28] Junji Kobashi, Hiroyuki Yoshida, and Masanori Ozaki,“Polychromatic optical vortex generation from patternedcholesteric liquid crystals,” Phys. Rev. Lett. , 253903(2016).[29] Raouf Barboza, Umberto Bortolozzo, Stefania Residori,and Marcel G. Clerc, “Berry phase of light bragg-reflected by chiral liquid crystal media,” Phys. Rev. Lett. , 053903 (2016).[30] K. Yu. Bliokh, D. Yu. Frolov, and Yu. A. Kravtsov,“Non-abelian evolution of electromagnetic waves in aweakly anisotropic inhomogeneous medium,” Phys. Rev.A , 053821 (2007).[31] Konstantin Y. Bliokh, Yuri Gorodetski, VladimirKleiner, and Erez Hasman, “Coriolis effect in optics:Unified geometric phase and spin-Hall effect,” Phys. Rev.Lett. , 030404 (2008). [32] Gabriel F. Calvo and Antonio Pic´on, “Spin-induced an-gular momentum switching,” Opt. Lett. , 838–840(2007).[33] Ebrahim Karimi, Bruno Piccirillo, Lorenzo Marrucci,and Enrico Santamato, “Light propagation in a birefrin-gent plate with topological charge,” Opt. Lett. , 1225–1227 (2009).[34] Melvin Lax, William H. Louisell, and William B. McK-night, “From Maxwell to paraxial wave optics,” Phys.Rev. A , 1365–1370 (1975).[35] A. Alberucci, C. P. Jisha, L. Marrucci, and G. Assanto,“Electromagnetic confinement via spin-orbit interactionin anisotropic dielectrics,” ACS Photonics , 2249–2254(2016).[36] P. L. Kapitza, “Dynamic stability of a pendulum whenits point of suspension vibrates,” Sov. Phys. JETP ,588–592 (1951).[37] A. Alberucci, L. Marrucci, and G. Assanto, “Light con-finement via periodic modulation of the refractive index,”New J. Phys. , 083013 (2013).[38] Mutatis mutandis , all the following results remain validin the more general case of a biaxial crystal.[39] S. Slussarenko, A. Alberucci, C. P. Jisha, B. Piccirillo,E. Santamato, G. Assanto, and L. Marrucci, “Guid-ing light via geometric phases,” Nat. Photon. (2016),10.1038/NPHOTON.2016.138.[40] R. Clark Jones, “A new calculus for the treatment of opti-cal systemsi. description and discussion of the calculus,”J. Opt. Soc. Am. , 488–493 (1941).[41] In this Paper we adopt the the source’s point of view.[42] M. Kwasny, U. A. Laudyn, F. A. Sala, A. Alberucci,M. A. Karpierz, and G. Assanto, “Self-guided beams inlow-birefringence nematic liquid crystals,” Phys. Rev. A , 013824 (2012).[43] P. A. M. Dirac, The Principles of Quantum Mechanics (Oxford Science Publications, Oxford, 1999).[44] K. Y. Bliokh, F. J. Rodr´ıguez-Fortuno, F. Nori, andA. V. Zayats, “Spin-orbit interactions of light,” Nat. Pho-ton. , 796–808 (2015).[45] From the standard analogy between 2D quantum me-chanics and paraxial optics in the monochromatic regime,the propagation coordinate z plays the role of time.[46] Kyle E. Ballantine, John F. Donegan, and Paul R.Eastham, “There are many ways to spin a photon:Half-quantization of a total optical angular momentum,”Science Advances (2016), 10.1126/sciadv.1501748,http://advances.sciencemag.org/content/2/4/e1501748.full.pdf.[47] Jon H. Shirley, “Solution of the Schr¨odinger equationwith a Hamiltonian periodic in time,” Phys. Rev. ,B979–B987 (1965).[48] Lowell S. Brown, “Quantum motion in a paul trap,”Phys. Rev. Lett. , 527–529 (1991).[49] Given the standard paraxial equation 2 ik n ∂A∂z + ∂ A∂x + k ∆ n A = 0, the photonic potential is defined as V = − k n n . [50] S. Pancharatnam, “Generalized theory of interference,and its applications,” Proc. Indian Acad. Sci. A , 0370–0089 (1956).[51] M. V. Berry, “Pancharatnam, virtuoso of the Poincar´esphere: an appreciation,” Current Science , 220–223(1994).[52] K. Y. Bliokh, C. T. Samlan, C. Prajapati, G. Puentes,N. K. Viswanathan, and F. Nori, “Spin-hall effect andcircular birefringence of a uniaxial crystal plate,” Optica , 1039–1047 (2016).[53] Jonathan Hu and Curtis R. Menyuk, “Understandingleaky modes: slab waveguide revisited,” Adv. Opt. Pho-ton. , 58–106 (2009).[54] Ardavan F. Oskooi, David Roundy, Mihai Ibanescu, Pe-ter Bermel, J. D. Joannopoulos, and Steven G. Johnson,“MEEP: A flexible free-software package for electromag-netic simulations by the FDTD method,” Comput. Phys.Commun. , 687–702 (2010).[55] A. D’Alessandro, L. Martini, G. Gilardi, R. Beccherelli,and R. Asquini, “Polarization-independent nematic liq-uid crystal waveguides for optofluidic applications,”IEEE Photon. Techn. Lett. , 1709–1712 (2015).[56] Amir Nader Askarpour, Yang Zhao, and Andrea Al`u,“Wave propagation in twisted metamaterials,” Phys.Rev. B , 054305 (2014).[57] J.B. Pendry, A.J. Holden, D.J. Robbins, and W.J. Stew-art, “Magnetism from conductors and enhanced non-linear phenomena,” Microwave Theory and Techniques,IEEE Transactions on , 2075–2084 (1999).[58] D. R. Smith and D. Schurig, “Electromagnetic wavepropagation in media with indefinite permittivity andpermeability tensors,” Phys. Rev. Lett. , 077405(2003).[59] A. Ghosh, N. K. Sheridon, and P. Fischer, “Voltage-controllable magnetic composite based on multifunc-tional polyethylene microparticles,” Small , 1956–1958(2008).[60] Alexander B. Khanikaev, S. Hossein Mousavi, Wang-Kong Tse, Mehdi Kargarian, Allan H. MacDonald, andGennady Shvets, “Photonic topological insulator,” Nat.Mater. , 233–239 (2013).[61] Carl Pfeiffer and Anthony Grbic, “Metamaterial Huy-gens’ surfaces: Tailoring wave fronts with reflectionlesssheets,” Phys. Rev. Lett. , 197401 (2013).[62] Konstantin Y. Bliokh, Yuri S. Kivshar, and FrancoNori, “Magnetoelectric effects in local light-matter inter-actions,” Phys. Rev. Lett. , 033601 (2014).[63] V. S. Asadchy, Y. Ra’di, J. Vehmas, and S. A. Tretyakov,“Functional metamirrors using bianisotropic elements,”Phys. Rev. Lett. , 095503 (2015).[64] Kejie Fang and Shanhui Fan, “Controlling the flow oflight using the inhomogeneous effective gauge field thatemerges from dynamic modulation,” Phys. Rev. Lett. , 203901 (2013).[65] Qian Lin and Shanhui Fan, “Light guiding by effectivegauge field for photons,” Phys. Rev. X , 031031 (2014).[66] Wolfgang Paul, “Electromagnetic traps for charged andneutral particles,” Rev. Mod. Phys.62