Propagation of focused scalar and vector vortex beams in anisotropic media: A semi-analytical approach

In the field of structured light, the study of optical vortices and their vectorial extension--vectorial vortex beams--has garnered substantial interest due to their unique phase and polarisation properties, which make them appealing for many potential applications. Combining the advantages of vortex beams and anisotropic materials, new possibilities for electromagnetic field tailoring can be achieved in nonlinear optics, quantum and topological photonics. These applications call for a comprehensive modelling framework that accounts for properties of both anisotropic materials and vector vortex beams. In this paper, we describe a semi-analytical model that extends the vectorial diffraction theory to focused vortex beams propagating through a uniaxial slab, considering the cases of scalar and vectorial vortexes in the common framework of a Laguerre-Gaussian modes. The model aims to provide a comprehensive description of the methodology, enabling the implementation of complex beams transmission through, reflection from and propagation in uniaxial anisotropic materials for specific applications. We apply the developed approach to propagation of high-order vortex beams in uniaxial materials with various dispersion characteristics>elliptic, hyperbolic and epsilon-near-zero regimes. We show how variations of the medium anisotropy modify the beam structure due to the vectorial nature of their interaction, which results from the different permittivities of the medium for transverse and longitudinal field components. The applicability of the approach can be extended to artificially structured media if they can described by effective medium parameters. The developed formalism will be useful for modelling of interaction of complex beams with uniaxial materials, allowing a common framework for a large variety of situations, which can also be extended beyond the electromagnetic waves.


INTRODUCTION
Since their first introduction [1], optical vortex beams (OVBs) have been subject of countless investigations and have led to important advancements, including in-cavity generation of OVBs, all-optical encryption techniques, OAM-based particles manipulation and beam shaping devices [2][3][4][5][6][7].The feature that makes optical vortices so interesting is their phase singularity, which is quantified by their topological charge (ℓ): an integer describing how many times the phase wraps in [0, 2π] interval in a closed loop around the beam centre.A nonzero topological charge makes the wavefront of OVBs helicoidal, with the number of helices per wavelength distance determined by the value of ℓ and the handedness by its sign.Therefore, OVBs possess an orbital angular momentum (OAM) of ±ℏℓ per photon.OVBs are important in optical communications, quantum optics and imaging, as they provide an additional degree of freedom for photons, introducing the possibility to encode more information in the same beam using topological charge.The combination of co-propagating OVBs leads to a nonuniform polarisation patterns in the resulting beam.In such, the so-called vectorial vortex beams (VVBs), the phase singularities of OVBs translate in polarisation singularities [8,9].These singularities are in many cases accompanied by strong longitudinal fields and often result in unusual behaviour, e.g., the violation of the optical theorem for scattering [10] or formation of topological structures of light, such as optical skyrmions [11].
As interesting as they are in a context of propagation in free space or uniform media, when interacting with anisotropic materials the behaviour of OVBs and VVBs becomes even more complex, leading to potential new opportunities and applications with new frontiers of beam shaping and polarisation control.The task of characterising their propagation in anisotropic media is far from trivial due to the inherent complexities arising from an interplay of the vectorial nature and vortex structure of the beams and the anisotropic properties of the material.While purely numerical simulations may provide the required information, the restrictions on the use of periodic boundary conditions due to the final size of the beam, results in the significant demands on computational resources and convergence issues.Therefore, the development of accurate and efficient modelling becomes imperative for studying and controlling the behaviour of VVBs in such scenarios, and advancing their applications.
In this work, we develop a semi-analytical approach which can be applied to the propagation of focused scalar and vectorial vortex beams through an anisotropic slab.The approach provides an extension of already examined arXiv:2402.04349v1[physics.optics]6 Feb 2024 Each paraxial ray at a height h = f sin θ from the optical axis corresponds to a refracted ray propagating along the direction θ.The angle θ is limited by the numerical aperture of the system to θmax.
cases of tightly focused beams propagating in free space and isotropic materials [12][13][14][15].Both OVBs and VVBs are considered, with the former modelled as Laguerre-Gaussian beams of general order LG ℓp and the latter as a superposition of orthogonally polarised OVBs with opposite topological charges ±ℓ [16].Following a brief introduction to the vector diffraction theory in multi-layered isotropic media, a detailed description of the methodology is presented for an anisotropic uniaxial medium.
Examples of applications of the approach are given, exploring the cases of Laguerre-Gauss beams of different orders, propagating through various categories of uniaxial media.The developed approach will be a useful tool for modelling the interaction of complex beams with material systems of chosen optical properties.The cases covered by the developed method include both isotropic and anisotropic (limited to uniaxial) materials and allow the exploration of various dispersion regimes, elliptic, hyperbolic, as well as epsilon-near-zero.

VECTORIAL DIFFRACTION THEORY IN ISOTROPIC MULTI-LAYERED MEDIA
The physical problem behind this approach is set to find the electric field of a focused beam, in a given volume divided into three domains, with the central one showing an anisotropic dielectric permittivity.One choice could be to follow a very general approach originally developed in [17] in the framework of the vectorial diffraction theory.It consists of an integral definition of the electric field at a fixed observation point, depending on the boundary electric and magnetic fields at the surface of an arbitrarily shaped aperture.For the majority of optical systems, the general approach can be simplified considering the asymptotic condition of the far-field diffraction, which can physically be thought of as the field distri-bution located sufficiently far from the focusing element of the system.Under this approximation, the Stratton-Chu integral can be understood as the electric field at a point ⃗ r in space originating from the superposition of an infinite number of plane waves propagating from the aperture to the point ⃗ r.If the aperture is considered to be a spherical pupil, this approach can be reduced to the well-known Debye-Wolf integral [18][19][20]: where f is the distance between the point ⃗ r and the pupil, ⃗ k = k 0 ⃗ k r = ω/c ⃗ k r , with ⃗ k r being the "relative" wavevector in the medium, ⃗ E Ω is the electric field strength within the solid angle Ω under which the pupil is seen from ⃗ r.The integration over Ω corresponds to a summation of all the plane waves directed to the point ⃗ r.Considering a decomposition of the unknown field into a superposition of plane waves with propagation properties dictated by the structure of the system, the angular spectrum formalism can be developed, which is the basis of the Richards-Wolf (RW) theory of vectorial diffraction [21].
A general solution of the wave equation can be written as the superposition of a number of plane waves with varying wavevector ⃗ k, weighted by an appropriate ⃗ A function [18,[22][23][24]: where the function ⃗ A(k x , k y ; z), which is the angular spectrum of ⃗ E(⃗ r), is the weight of each plane wave component propagating at the direction ⃗ k.It corresponds to the two-dimensional Fourier transform of the field, calculated at a reference plane (here z = 0).This plane can in principle be a plane orthogonal to an arbitrary direction, but it is usually conveniently chosen to be the wave propagation direction.For uniform and isotropic medium, considering a time-harmonic field, the propagation in real space can be calculated as a product in reciprocal space: This represents a powerful formalism for the computation of an electromagnetic field, assuming its distribution can be written in a reference plane.It can also be easily combined with the description of multilayered media by including the Fresnel coefficients describing the system in Eq. (2).For the simple case of multilayered media where all the interfaces are planes parallel to each other and orthogonal to the propagation direction, the boundary conditions can be solved to find the Fresnel coefficients.
In the scope of this paper, the discussion is limited to a slab with three consecutive materials of uniform dielectric media, but the formalism is not restricted to this simplification.Within this assumption, the Fresnel coefficients are (see, e.g., [25]): face between the media i-th and j-th, subscripts I and II refer to the slab interfaces located at z = z 0 and z = z 1 , respectively, d = |z 0 − z 1 | is the thickness of the slab, and s and p are the polarisation components of the field with respect to the plane of incidence and the direction of ⃗ k.
The electric field in the integration volume is piece-wise defined as where the labels i, r, t stand for "incident", "reflected" and "transmitted", respectively.Each term is then projected onto its s and p field components, relating, according to Eq. ( 2), its angular spectrum to the incident field via the Fresnel coefficients from Eqs. (5).The field inside the slab is then A useful application of the angular spectrum formalism is the study of focused fields, which has led to the RW theory [21,26].In simple terms, the field of a focused laser beam is determined by the effect of the focusing system-a lens-on the incoming beam, which is described by a boundary problem at the interface corresponding to the lens.The focusing system is considered to be aplanatic and light comes from a source at infinity, so that the wavefront is assumed to be planar in the plane of the lens.An implicit approximation of this method comes from the use of the angular spectrum formalism in the asymptotic regime: the reference plane is located at infinity, where the wavefront is planar.Assuming the beam is paraxial before the lens, the reference plane can be chosen to coincide with the lens surface.Within asymptotic approximation, the fields in the proximity of the lens can then be formulated in the frame of geometrical optics with two conditions to be fulfilled [27]: (i) the sine condition and (ii) the intensity law.The latter is related to conservation of energy upon propagation through the lens: the energy flux of each ray needs to be constant and the power has to be the same on both side of the lens surface.The former describes the lens boundary as a sphere centred in its geometric focus, with the focal length being its radius, and ensures a one-to-one mapping of the rays incident on the lens to corresponding refracted rays.The distance of each ray from the optical axis-chosen parallel to the propagation direction-can be written as h = f sin θ (Fig. 1) so that independently of the choice of coordinates in the half-space before the lens, it can be mapped onto spherical coordinates on the sphere surface.The angle θ is the refraction angle of the ray at the distance h from the optical axis.
The above approximations lead to a modification of Eq. ( 1): ϕ) e ikz cos θ e ikρ sin θ cos(ϕ−φ) where ⃗ E L is the electric field incident on the surface of the lens.The integration in the reciprocal space has been limited to the angle θ max corresponding to the half aperture of the lens field of view (NA = n sin θ max ), while ϕ assumes all the possible values in [0, 2π] interval.The solution propagating along negative z has been discarded for obvious physical reasons.The choice of spatial coordinates (x = ρ cos φ, y = ρ sin φ) is made for simplifying the integrals computation.Because of the last exponential factor, it is in fact possible to use the integral definitions of modified Bessel functions of the first kind given by [21] 2π which allows for an analytical solution for at least one of the integrals, once the amplitude of the ⃗ E ∞ components are written as a combination of trigonometric functions.In order to solve Eq. ( 8), the ⃗ E ∞ field distribution needs to be made explicit.Focusing through isotropic multilayered media can easily be implemented, following the procedure described in [13,14].

VECTOR DIFFRACTION THEORY IN UNIAXIAL MEDIA
In solving the boundary problem for an isotropic multilayered medium, the conservation of the transverse component of the wave vector (k x , k y ) must be ensured, and its z component can be represented through k x and k y , once the dispersion relation of each medium is known.If the medium is anisotropic, this is generally no longer the case.For this reason, we limit our discussion to the case of a uniaxial crystal with the optical axis along the propagation direction (ẑ), so that the electric field in the (x, y)-plane sees an isotropic medium.In this instance, the permittivity tensor which describes the slab is given by Each plane wave component in the angular spectrum has to satisfy the wave equation ⃗ k × ⃗ k × ⃗ E + µε ⃗ E = 0 which, using Eq. ( 10) and reducing the system to a two-dimensional problem (k ry = 0) can be written as: where the subscript r denotes the "relative" wavevector , as previously mentioned.The solutions for the electric field ⃗ E are then given by the null space of the k-matrix, which only exists when its determinant is zero.This condition can be fulfilled in two cases, corresponding to two possible polarisation modes.The associated magnetic field can be obtained from the Faraday's law ( ⃗ These two modes and their fields are presented in Table I.
The solutions in Table I suggest the possibility to define a basis with respect to the k-vector onto which the electric and magnetic fields [28,29] may be decomposed: The core difference between the iso- [28] and anisotropic cases stems from the different solutions found for the k-vector: upon propagation from Region 1 to Region 2 (Fig. 2), an isotropic slab will make all the components of ⃗ k scale equally.In contrast, there is mixing of transverse and longitudinal components of ⃗ k for an anisotropic slab.
The decomposition of the electric field into its s and p components allows for an easy implementation of the angular spectrum formalism for an uniaxial slab.Knowing the expressions for k z from the modes found for a uniaxial dielectric, the Fresnel coefficients in Eqs. ( 5) can be adapted to an anisotropic layer.

MODELLING APPROACH
Equation 12 can be applied to an anisotropic slab following the same procedure as in the isotropic case described above, with the Fresnel coefficients of the system included in the integrals in Eq. (8).Each field appearing in the piece-wise definition in Eq. ( 6) can be calculated once the explicit form of ⃗ E ∞ is fixed.We apply this semi-analytical approach to both the cases of scalar optical vortices and vector vortex beams.
Both these type of beams can be described as Laguerre-Gaussian beams with a complex amplitude LG ℓ,p [30]:

Polarisation Condition
E-field H-field where L ℓp represents the generalised Laguerre polynomials [31], with ℓ and p being their azimuthal and radial orders, respectively, w(z) describes the beam lateral size as a function of the z coordinate with w 0 being its minimum value, R(z) is the wavefront radius of curvature, and z R = πw 2 0 /λ is the beam Rayleigh range.Following the prescriptions of the angular spectrum formalism, the reference plane is chosen as z = 0, so that w(z) = w 0 , exp ρ 2 /R → 1, and arctan(z/z r )=0.Finally, to account for the limitation to the field of view of the focusing element imposed by its finite size, the apodisation function sin 2 (θmax) should be introduced.Here, f 0 = w0 f sin θmax is a geometrical factor representing the ratio between the beam and the lens lateral sizes.With this choice, the dependence on ρ is also modified so that ρ 2 = w0 sin 2 (θ) f0 sin 2 (θmax) .With these approximations, the mode amplitude only depends on the variable θ, while its phase only on ϕ: where Θ and Φ collect the terms depending on θ and ϕ, respectively, and C ℓp all the constant factors: where the two-dimensional vector (Φ x , Φ y ) is needed for the description of vectorial vortices.The electric field incident on the lens ( ⃗ E L ) can then be written in a general way as: where the exact definition of (Φ x , Φ y ) depends on the type of vortex.

Scalar vortex
In the case of a scalar vortex, Φ x = E 0x exp(−iℓϕ), Φ y = E 0y exp(−iℓϕ), and the dependence on ϕ is the same for both components resulting in uniform SOP.The beam polarisation can then be completely described by the coefficients (E 0x , E 0y ), which represent the projections of the SOP on the polarisation basis {|H⟩ , |V⟩}.Denoting these projections as ⟨H|ψ s ⟩ and ⟨V|ψ s ⟩, the scalar SOP |ψ s ⟩ can be obtained as Vector vortex Vectorial vortexes can be modelled as a superposition of LG beams with opposite values of ℓ, the same value of p, and orthogonal circular polarisations.As a consequence of this, the dependence on θ of both components of ⃗ E L is the same, while the phase terms have different dependence on ϕ.After introducing the polarisation basis [16] |R the vectorial SOP |ψ v ⟩ can be written as the sum of orthogonal circular scalar vortices: where the projections on the helicity basis (Eq.18) are E 0R = ⟨R, ℓ|ψ v ⟩ and E 0L = ⟨L, ℓ|ψ v ⟩.To consistently label SOPs in Eq. 19 and Eq.17, the coefficients E 0R , and E 0L in |ψ v ⟩ can be expressed in terms of E 0x and E 0y , so that Within this notation framework, assuming ℓ = 0, both |ψ s ⟩ and |ψ v ⟩ return uniform SOPs, expressed with the two-dimensional Jones vector (E 0x , E 0y ).For example, (E 0x , E 0y ) = (1, 0) gives the SOP |H⟩ and (1, −i)/ √ 2 corresponds to |R⟩.Values of ℓ ̸ = 0 will otherwise produce scalar SOPs from |ψ s ⟩ and vectorial ones from The examples of various SOPs are shown in Fig. 3.The SOPs corresponding to ℓ = 0, which can be equivalently obtained as scalar vortices or "degenerate" vectorial ones, are presented as points on the surface of the sphere.Each point (E0x, E0y) on the Poincaré sphere can also represent a vectorial polarisation state, if a value of ℓ ̸ = 0 is chosen.For example, the point associated with |H, 0⟩, which in a scalar case corresponds to a horizontal SOP, represents instead a radial beam |H, 1⟩ for ℓ = 1.SOPs belonging to the same series show the same local polarisation distribution along their horizontal central line.In Fig. 3 for each point on the sphere, the SOP of the vectorial vortices obtained for ℓ = ±1, ±2 and the same values of (E 0x , E 0y ) are also shown.
Consistent with the above polarisation description, the functions Φ x and Φ y for vectorial SOPs can also be rewritten in terms of E 0x and E 0y :

Electric field evaluation
Adopting the same piecewise definition of the electric field as in Eq. ( 6), the electric field distribution can be calculated as where the subscript j labels the region of space and α labels the s and p components of each vector.The term ⃗ E ∞ is obtained substituting Eq. ( 16) into Eq.( 8)b, once the functions Θ and Φ in Eq. ( 16) have been determined for the chosen beam type.Upon substitution of ⃗ k ± and of the Fresnel coefficients, expressions depending on polarisation (i.e., s or p) and layer (i.e., being it isotropic or uniaxial) are obtained.Calculations, limited to the analytical integration over ϕ, have been written in Wolfram Mathematica, while the final steps have been carried out on MATLAB.Below, the procedure followed for each of the fields appearing in the piece-wise definition of Eq.( 6) is described, using as a mock example the x component of the field ⃗ E t1 , for an incident field given by a generic VVB, whose SOP is given by (E 0x , E 0y ) (assuming it to be normalised to one).In the following, the electric field taken as an example will be simply labelled as ⃗ E for simplicity of notations.
In the first step, the field on the lens plane ( ⃗ E ∞ ) is divided into its s and p components (Eq.( 8b)) and the projections on the basis vectors are calculated: They are then multiplied by the unit vectors given in Eq. ( 12), according to the studied case, in order to obtain the full vector ⃗ E s,p ∞ .After that, for each component of the angular spectrum, the factors solely depending on θ are collected together, yielding: with Ξ(θ, ϕ; φ) In the above expressions, the term √ cos θ comes from the intensity law, while the term sin θ (which is in addition to the one contained in Θ) comes from the Jacobian of the integration variables transformation.In the next step, the integration over ϕ is carried out applying Eq. ( 9) and upon simplification of the trigonometric functions of ϕ, when necessary.With Wolfram Mathem-atica, it is indeed possible to integrate Bessel functions of the first kind automatically and by means of Eq.( 9), but this is no longer possible if the order of the Bessel function is kept arbitrary.This would imply the need of specifying the desired order of the Laguerre polynomial before completing the first integration.For this reason, the code also includes a list of substitutions used to apply Eq.( 9) to a wide set of combinations of trigonomet- (θ) in the integrals Eq. 26 calculated in the medium i, for either s or p polarisation and propagation along ±ẑ direction.Labels of the Fresnel coefficients have the same meaning as previously introduced.
ric functions, returning the corresponding combination of Bessel functions.The list has been built using all the possible combinations found in the calculations, to en-sure the possibility to write a general solution for any arbitrary beam LG ℓp .Upon completion of this step, the field spectral components become: The components are then re-arranged collecting the Bessel functions involved.The resulting expressions are used in five integrands with which all the field components can be reconstructed: The explicit form of Θ depends on the beam type.The functions ζ(θ) and η(θ) represent quantities whose explicit expression depends on the medium (labelled by i), polarisation (s or p) and propagation direction (±ẑ) they are calculated for (Table II).The former (ζ) is the propagation factor and the latter (η) is a general representation of the Fresnel coefficients needed in each case.
It is worth noting that in the case of the integrals I 4 and I 5 , only the contribution of the p component of the field is important.Finally, the numerical integration is performed in MATLAB.Similar procedure can be applied to obtain the magnetic field, retrieving its angular spectrum from the electric field one.To provide a practical example, Table III presents the expressions obtained for the electric field ⃗ E t1 in both the cases of OVs and VVBs, where the dependence on the coordinates (ρ, φ, z) have been omitted for simplified notations.
It is worth noting that, when calculating the field for a VVB, the three θ-integrals (I 1 , I 2 , and I 4 ) need to be evaluated, while the five integrals are needed for a scalar vortex.This is related to the presence of the ϕdependent phase factor of the LG modes, which is lost in the superposition of orthogonal circular vortices with opposite topological charge, needed for VVBs.

Longitudinal field example
This semi-analytical approach introduces an extension of the Richards-Wolf theory of vectorial diffraction to an anisotropic medium.The advantage of a semi-analytical model, compared to fully numerical studies, is the possibility to handle close expressions for the fields and retrieve the related fundamental information from the functions describing their components.For example, it is interest-ing to look at the longitudinal (z-) component generated by the focusing.In both cases of OVs and VVBs, it solely depends on the p polarisation component of the angular spectrum, since it only contains the integrals I 4 and I 5 .While the spatial distribution of the the longitudinal field depends on the results of the integration and on all the factors contributing to it, it is possible to understand the general trends directly from the integrands.In particular, I 4 and I 5 contain the terms , which are calculated for each value of ρ over all the range of values of the integration variable θ, so that the field at any point in space, depends on the superposition of all the plane waves entering the field of view of the objective.
The integrals I 4 and I 5 show a case for ℓ = ±1, which leads to the appearance of the Bessel function J 0 .This is in fact the only function of the set being nonzero-valued when its argument is zero, which implies that only the terms containing this function will have a non-zero integral in the spatial points corresponding to ρ = 0. Comparing the maps of ξ ℓ (ρ, θ) for ℓ = 0, ±1 (Fig. 4), it is evident how J 0 is the only case where for small values of ρ the integral is nonzero.Looking at the same maps for higher values of ρ, it is expected that the value of these integrals will show regions of different signs for the increasing distance from the origin.
Considering the case of OVs, the z component of the focused field contains the Bessel function J 0 for ℓ = −1 and this appears together with the function J 2 , making more detailed predictions quite cumbersome at this stage, given all the complicated dependencies appearing in the integrals and the general dependence on the SOP.On the other hand, when assuming ℓ = 1 for a VVB on a higherorder Poincaré sphere, the z component only contains the Bessel function J 0 and two particular cases can be highlighted.The z component of the VVB field when the dependence on the coordinate φ is made explicit is: so that the only contribution to the longitudinal field can come from the component of the beam on the state |H, 1⟩ ((E 0x , 0), in this case).This state of polarisation represents a radial beam, which is known for its nonzero longitudinal field, whose intensity can be strongly increased by tight focusing [13,32,33].On the other hand, if the initial SOP is given by (0, E 0y ), which corresponds to an azimuthal beam, the longitudinal field will be zero in every point of the real space.

APPLICATION EXAMPLES
The approach described in the previous section provides the possibility to simulate the focusing of different types of optical vortices and their propagation through a three-layered medium.The required inputs are: (i) E 0x , E 0y , ell, p, w 0 and λ in order to completely characterise the beam; (ii) the permittivities of each layer, to describe its optical behaviour; and (iii) the numerical aperture and the location of the geometrical focus of the lens with respect to the interfaces, to characterise the focusing features of the system.We consider several archetypal cases in order to demonstrate the versatility of the developed approach: a standard Gaussian beam, obtained as a Laguerre-Gauss beam with both indices set to 0 (LG 00 ); a high-order scalar vortex (LG 15,0 ); a scalar vortex with a nonzero radial index (LG 10,2 ); and a vectorial vortex whose constituent vortices are given by LG ±10,2 beams, referred to as the VVB 10,2 (Fig. 5).The wavelength of choice is λ = 650 nm and the beams are focused with an objective of numerical aperture of 0.9 in all cases.All the scalar beams have been simulated assuming horizontal polarisation (i.e., the electric field parallel to the x direction), hence only one of the initial electric field components is nonzero (E 0y = 0).The vectorial vortex is calculated with the same prescription, so that its polarisation state can be thought of as a higher-order equivalent of the horizontal scalar polarisa-tion: scalar and vectorial are located on the same point of the Poincaré sphere, but obtained for different values of topological charge.Any other state on the Poincaré sphere is obtained by a careful selection of E 0x , E 0y and ℓ (cf.Fig. 3).
Free space.The simplest case to model is the propagation of the focused beam in free space, where all the permittivities are set equal to one (Fig. 5).The size of the spatial mesh used in the calculations has been expanded to fully appreciate the spatial variations of the transverse and longitudinal field components, over a macroscopic distance.Given the strong focusing regime, even a simple Gaussian beam develops a nonzero longitudinal field, which shows a two-lobed shape, aligned to the polarisation direction.Although nonzero, the longitudinal field intensity is considerably smaller than the transverse one: the maximum value of the former is approximately 0.6% of the latter.If a nonzero topological charge is introduced, an immediate rise in the longitudinal field strength is observed.Its maximum becomes in fact the 26% of the transverse field for the scalar LG 15,0 vortex, reaching a 35% for LG 10,2 .Common to all the scalar vortices, not limited to the three presented here, the symmetry of the longitudinal field intensity distribution is affected by a focusing-induced astigmatism, which is manifested in a prolate shape of the beam.This is different from a vectorial beam, where the interference between the co-propagating vortices of opposite topological charges results in a cylindrically symmetric shape.Together with a higher degree of symmetry, the intensity of the longitudinal field relative to its transverse counterpart is also increased by the vectorial nature of this types of beam: the longitudinal field maximum becomes approximately 50% of the transverse one, making these beams appealing for applications where a strong longitudinal field is needed.
Table III: Expression of the electric field ⃗ E t1 as an example in the description of the methodology.The vector components are given for both cases of a scalar (OVB) and vectorial (VVB) vortex beams of generic polarisation state (E 0x , E 0y ) and topological charge ℓ.Some terms appearing in the field components have been gathered in extra factors, whose definition is given in the bottom part of the table.All the integrals are labelled with a superscript showing the medium they are calculated in and the polarisation component they refer to.In each integral, the Fresnel coefficients in η represent transmission coefficients, according to Table II.
Electric field in medium 2, propagating along positive z Factors appearing in the expressions above Isotropic lossy slab.The simulated system can be made slightly more complicated by the introduction of a single interface, which can be obtained using the same values of permittivities for two consecutive layers.For the sake of brevity, from here on we are only reporting the results for a more general case of a double interface and limited to the beams LG 10,2 and VVB 10,2 .Any of the regions of the system (Fig. 2) can be isotropic and lossy (complex values of permittivity are supported in the developed approach), so that the first case considered is a lossy dielectric (Fig. 6, yellow frames) with permittivity ε 2x = ε 2z = 2.25 + i5 × 10 −2 .The results are similar to the previous case of free-space propagation with the main difference being a reduction of the overall intensity of the beam along the propagation direction, equally damping both the transverse and longitudinal components of the beam.However, the reflection at the slab boundaries results in the field intensity redistribution across the beam cross-section for both field components.
Uniaxial slab: elliptic dispersion.The central layer can also allow anisotropy, as long as there is a single optical axis and it is aligned with the ẑ direction of the reference frame, which coincides with the normal to the interfaces.Depending on the sign of the permittivity components Re(ε 2x ) and Re(ε 2z ) of Eq. ( 11)), this case describes different dispersion regimes.When their product is positive, the dispersion regime is a conventional elliptic, inheriting the name from the shape of the k-surface characteristic of this case (Fig. 6, purple frames).The permittivities have been chosen to be ε 2x = 2.25 + i5 × 10 −2 , ε 2z = 1.9 + i5 × 10 −2 .In this case, the transverse field, affected by ε 2x , shows the same behaviour as in the previous case of a lossy isotropic dielectric.The differences from the isotropic case are visible in the shapes of the intensity distributions of the beam profiles, augmented by a variation of the standing wave pattern obtained upon multiple reflections inside the anisotropic slab.On the other hand, a smaller real part of Re(ε 2 z), which introduces anisotropy, produces a stronger longitudinal field inside the slab, as a consequence of the conservation of the longitudinal component of the electric displacement vector.The imaginary parts have been kept the same as in the case of an isotropic dielectric in order to avoid overlapping different effects (different values for Im(ε 2x ) and Im(ε 2z ) would produce different damping rates for transverse and longitudinal components and interplay between them).
Uniaxial slab: hyperbolic dispersion.If eventually Re(ε 2x ) Re(ε 2z ) < 0, the second layer has a hyperbolic dispersion (Fig. 2, blue frames).The chosen permittivities for this case are: ε 2x = 2.25 + i5 × 10 −2 , ε 2z = −1.9+ i5 × 10 −2 , where again the transverse component matches the value used for an isotropic dielectric and all the imaginary parts are kept the same.The model reliably reproduces the effects of negative refraction, phenomenon known to happen in hyperbolic materials [34][35][36].The beam is in fact being re-focused to a point outside of the uniaxial slab, and its lateral dimensions are strongly modified inside of it.The modifications of the field intensity distributions are also clearly visible similar to the previous cases but with stronger damping of the longitudinal field and stronger divergence of the transverse field inside the slab, which are connected because of the inter-relation between longitudinal and transverse fields in the beam [7,37].
The described approach is obviously not limited to the case of natural materials and can be applied to structured materials if the effective medium considerations can be used to describe optical properties through effective permittivity.For example, depending on the constituent materials and geometric parameters, plasmonic nanorod-based metamaterials or metal-dielectric multilayered metamaterials may exhibit elliptic or hyperbolic dispersion regimes, or an epsilon-near-zero regime, when a permittivity tensor component approaches zero.The developed approach can be used for simulations of reflection, transmission and absorption spectra, and reveals how the polarisation and intensity distributions of a beam are modified by a particular dispersion regime of the metamaterial [37].

CONCLUSIONS
The developed semi-analytical model for the vectorial diffraction theory in anisotropic uniaxial media can be applied to a wide range of situations, including focused optical vortices, both scalar and vectorial, and their propagation through a dielectric slab, whose optical behaviour can be either isotropic or uniaxial, exploring physically interesting dispersion cases, such as hyper-bolic and epsilon near-zero regimes.The presented approach allows for a comprehensive investigation of various parameters, including the objective numerical aperture, the material permittivity, the slab thickness, and the beam state of polarisation.Moreover, by exploiting the Laguerre-Gauss basis for the description of the vortices, this approach can addresses crucial aspects of optical wave propagation in complex media in regards of orbital angular momentum physics.
The possibility to model such a system can assist the exploration of a wide range of applications, such as optical communication, imaging systems, and laser beam shaping, where the propagation of an optical vortex through an anisotropic slab can play a crucial role.The simulations described here are provided as an opensource software package called "InFOCUS" (Interaction of FOcused Complex beams with Uniaxial Slabs), available at Ref. [38].
In addition to the variety of cases the presented model can be applied to, there are still many potential extensions that could be implemented for future improvements.For instance, extending the model to include the possibility to change the angle of incidence of the incoming beam or, on the material side, including cases described by a more complicated dielectric tensor like, for example, chiral media.

Figure 1 :
Figure 1: Schematics of a paraxial beam focused by a lens.The incoming beam, considered in cylindrical coordinates (ρ, ϕ, z, shown in green), is mapped onto a spherical co-ordinate frame (k, ϕ, θ, shown in blue) upon refraction through a lens.Each paraxial ray at a height h = f sin θ from the optical axis corresponds to a refracted ray propagating along the direction θ.The angle θ is limited by the numerical aperture of the system to θmax.

Figure 2 :
Figure 2: Geometry of the three-layered system.Light is incident from medium 1, region 2 is filled with a uniaxial material.

Figure 3 :
Figure 3: Poincaré sphere representation of the polarisation states.The Stokes vectors S 1 , S 2 and S 3 are used as the coordinate axes of the polarisation space.Scalar polarisation states are represented on the surface of the sphere as black arrows, while vectorial SOPs, calculated according to Eq. (19) with values of topological charge ±1, are shown in correspondence of the scalar SOPs having the same values of (E 0x , E 0y ).The number below each circle indicates the topological charge of the SOP above it, while the two-dimensional vectors indicate the value of E 0x , E 0y used for the corresponding series of SOPs.Similar SOPs can be obtained for the southern hemisphere of the Poincaré sphere with the only difference of a π/2 rotation of the pattern.North and south poles represent uniform circular polarisation regardless of the topological charge.

Figure 5 :
Figure 5: Vortex beam focusing in free space with NA = 0.9: (a) Gaussian beam, (b) scalar vortex beam LG 15,0 , (c) scalar vortex beam LG 10,2 , (d) vectorial vortex beam LG ±10,2 .Colour maps represent the intensity of the (T) transverse and (L) longitudinal electric field components, calculated in both (left) transverse (x − y) at z = 0 and (right) longitudinal (x − z) at y = 0 planes.The line plots are the cross-sections of the corresponding maps along the x or y) axes.The SOP of each type of beam is shown in brackets, according to the basis used for scalar and vectorial beams.

Table I :
Modes of the electromagnetic field in a uniaxial material.

Table II :
Functions ζ