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

Vittorio Aita  Department of Physics and London Centre for Nanotechnology,
King’s College London, Strand, London WC2R 2LS, UK  Mykyta Shevchenko
Department of Electronic & Electrical Engineering, University College London
(UCL), Torrington Place, London WC1E 7JE, UK  Francisco J. Rodríguez-Fortuño
Department of Physics and London Centre for Nanotechnology, King’s College
London, Strand, London WC2R 2LS, UK  Anatoly V. Zayats  Department of Physics
and London Centre for Nanotechnology, King’s College London, Strand, London
WC2R 2LS, UK

(28th August 2024)

###### Abstract

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 and manipulation 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 the case of focused vortex beams
propagating through a uniaxial slab, considering both the cases of scalar and
vectorial vortexes in the common framework of a Laguerre-Gaussian mode basis.
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. As a demonstration of its versatility, we apply the developed
approach to describe propagation of high-order vortex beams in uniaxial
materials with various dispersion characteristics, exploring the 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.

##  I 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 ($\ell$): an
integer describing how many times the phase wraps in $[0,2\pi]$ 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 $\ell$ and the handedness by its sign.
Therefore, OVBs possess an orbital angular momentum (OAM) of $\pm\hbar\ell$
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 non-
uniform 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 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 $\mathrm{LG}_{\ell p}$ and the
latter as a superposition of orthogonally polarised OVBs with opposite
topological charges $\pm\ell$ [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.

##  II Vectorial diffraction theory in isotropic multi-layered media

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

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 distribution 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
$\overrightarrow{\mathbf{r}}$ in space originating from the superposition of
an infinite number of plane waves propagating from the aperture to the point
$\overrightarrow{\mathbf{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]:

$\overrightarrow{\mathbf{E}}(\overrightarrow{\mathbf{r}})=-\frac{ik\,e^{-ikf}}{2\pi
f}\iint\limits_{\Omega}\overrightarrow{\mathbf{E}}_{\Omega}\,e^{i\,\overrightarrow{\mathbf{k}}\cdot\overrightarrow{\mathbf{r}}}\,\differential\Omega\,,$
(1)

where $f$ is the distance between the point $\overrightarrow{\mathbf{r}}$ and
the pupil,
$\overrightarrow{\mathbf{k}}=k_{0}\overrightarrow{\mathbf{k}}_{r}=\omega/c\overrightarrow{\mathbf{k}}_{r}$,
with $\overrightarrow{\mathbf{k}}_{r}$ being the “relative” wavevector in the
medium, $\overrightarrow{\mathbf{E}}_{\Omega}$ is the electric field strength
within the solid angle $\Omega$ under which the pupil is seen from
$\overrightarrow{\mathbf{r}}$. The integration over $\Omega$ corresponds to a
summation of all the plane waves directed to the point
$\overrightarrow{\mathbf{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 $\overrightarrow{\mathbf{k}}$,
weighted by an appropriate $\overrightarrow{\mathbf{A}}$ function [22, 18, 23,
24]:

$\displaystyle\overrightarrow{\mathbf{E}}(x,y,z)$
$\displaystyle=\iint\limits_{-\infty}^{\infty}\overrightarrow{\mathbf{A}}(k_{x},k_{y})\,e^{i(k_{x}x+k_{y}y)}\,e^{\pm
ik_{z}z}\differential k_{x}\,\differential k_{y}\,,$ (2a)
$\displaystyle\overrightarrow{\mathbf{A}}(k_{x},k_{y})$
$\displaystyle=\frac{1}{4\pi^{2}}\iint\limits_{-\infty}^{\infty}\overrightarrow{\mathbf{E}}(x,y,0)\,e^{i(k_{x}x+k_{y}y)}\differential
x\,\differential y\,,$ (2b)

where the function $\overrightarrow{\mathbf{A}}(k_{x},k_{y};z)$, which is the
angular spectrum of
$\overrightarrow{\mathbf{E}}(\overrightarrow{\mathbf{r}})$, is the weight of
each plane wave component propagating at the direction
$\overrightarrow{\mathbf{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:

$\displaystyle\overrightarrow{\mathbf{A}}(k_{x},k_{y};z)=\overrightarrow{\mathbf{A}}\\!\left(k_{x},k_{y},0\right)e^{\pm
ik_{z}z}\,,$ (3) $\displaystyle\quad k_{z}\equiv
k_{0}\sqrt{1-k_{\mathrm{r}x}^{2}-k_{\mathrm{r}y}^{2}}\,.$ (4)

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]):

$\displaystyle r^{s,p}_{\mathrm{I}}$ $\displaystyle=\frac{e^{2ik_{1z}z_{0}}\
\left(r_{12}^{s,p}+r_{23}^{s,p}e^{2idk_{2z}}\right)}{1+r_{12}^{s,p}r_{23}^{s,p}e^{2idk_{2z}}}\,,$
(5a) $\displaystyle t^{s,p}_{\mathrm{I}}$
$\displaystyle=\frac{e^{i(d(k_{2z}-k_{3z})+z_{0}(k_{1z}-k_{3z}))}\
t_{12}^{s,p}t_{23}^{s,p}}{1+r_{12}^{s,p}r_{23}^{s,p}e^{2idk_{2z}}}\,,$ (5b)
$\displaystyle r^{s,p}_{\mathrm{II}}$
$\displaystyle=\frac{e^{i(2dk_{2z}+z_{0}(k_{1z}+k_{2z}))}\
r_{23}^{s,p}t_{12}^{s,p}}{1+r_{12}^{s,p}r_{23}^{s,p}e^{2idk_{2z}}}\,,$ (5c)
$\displaystyle t^{s,p}_{\mathrm{II}}$
$\displaystyle=\frac{e^{i(d(k_{2z}-k_{3z})+z_{0}(k_{1z}-k_{3z}))}t_{12}^{s,p}t_{23}^{s,p}}{1+r_{12}^{s,p}r_{23}^{s,p}e^{2idk_{2z}}}\,,$
(5d)

where the numeric subscript labels the medium (two number pairs $ij$ refer to
the interface between medium $i$ and $j$), so that $r_{ij}$ and $t_{ij}$ are
the Fresnel reflection and transmission coefficients, respectively, for the
interface 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 $\overrightarrow{\mathbf{k}}$. The electric field in the
integration volume is piece-wise defined as

$\overrightarrow{\mathbf{E}}(\overrightarrow{\mathbf{r}})=\begin{dcases}\overrightarrow{\mathbf{E}}_{1}=\overrightarrow{\mathbf{E}}_{\mathrm{i}}+\overrightarrow{\mathbf{E}}_{\mathrm{r}1}&\mathrm{for}\quad
z\leq z_{0}\\\
\overrightarrow{\mathbf{E}}_{2}=\overrightarrow{\mathbf{E}}_{\mathrm{t}1}+\overrightarrow{\mathbf{E}}_{\mathrm{r}2}&\mathrm{for}\quad
z_{0}<z\leq z_{1}\\\
\overrightarrow{\mathbf{E}}_{3}=\overrightarrow{\mathbf{E}}_{\mathrm{t}2}&\mathrm{for}\quad
z>z_{1}\end{dcases}\,,$ (6)

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

$\displaystyle\overrightarrow{\mathbf{E}}_{2}=\overrightarrow{\mathbf{E}}_{\mathrm{t}1}+\overrightarrow{\mathbf{E}}_{\mathrm{r}2}=\overrightarrow{\mathbf{E}}_{\mathrm{t}1}^{s}+\overrightarrow{\mathbf{E}}_{\mathrm{t}1}^{p}+\overrightarrow{\mathbf{E}}_{\mathrm{r}2}^{s}+\overrightarrow{\mathbf{E}}_{\mathrm{r}2}^{p}\,,$
(7a)
$\displaystyle\overrightarrow{\mathbf{A}}_{\mathrm{t}1}^{(s,p)}=t^{s,p}_{\mathrm{I}}\,\overrightarrow{\mathbf{A}}_{\mathrm{I}}^{(s,p)}\,,\qquad\overrightarrow{\mathbf{A}}_{\mathrm{r}2}^{(s,p)}=r^{s,p}_{\mathrm{II}}\,\overrightarrow{\mathbf{A}}_{\mathrm{I}}^{(s,p)}\,.$
(7b)

A useful application of the angular spectrum formalism is the study of focused
fields, which has led to the RW theory [26, 21]. 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\theta$ (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 $\theta$ 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):

$\displaystyle\overrightarrow{\mathbf{E}}(\rho,\varphi,z)=-\frac{ikfe^{-ikf}}{2\pi}$
$\displaystyle\times\int\limits_{0}^{\theta_{\mathrm{max}}}\int\limits_{0}^{2\pi}\overrightarrow{\mathbf{E}}_{\infty}(\theta,\phi)\,\underbrace{e^{ikz\cos\theta}e^{ik\rho\sin\theta\cos(\phi-\varphi)}}_{e^{i\overrightarrow{\mathbf{k}}\cdot\overrightarrow{\mathbf{r}}}}\overbrace{\sin\theta\,\differential\phi\,\differential\theta}^{\differential\Omega}\,,$
(8a) $\displaystyle\overrightarrow{\mathbf{E}}_{\infty}$
$\displaystyle=\Bigl{[}t^{s}\bigl{(}\overrightarrow{\mathbf{E}}_{\mathrm{L}}\dotproduct\hat{\mathbf{n}}_{\varphi}\bigr{)}\hat{\mathbf{n}}_{\varphi}+t^{p}\bigl{(}\overrightarrow{\mathbf{E}}_{\mathrm{L}}\dotproduct\hat{\mathbf{n}}_{\rho}\bigr{)}\hat{\mathbf{n}}_{\theta}\Bigr{]}\sqrt{\frac{n_{1}}{n_{2}}}\sqrt{\cos\theta}\,,$
(8b)

where $\overrightarrow{\mathbf{E}}_{\mathrm{L}}$ is the electric field
incident on the surface of the lens. The integration in the reciprocal space
has been limited to the angle $\theta_{\mathrm{max}}$ corresponding to the
half aperture of the lens field of view
($\mathrm{NA}=n\sin\theta_{\mathrm{max}}$), while $\phi$ assumes all the
possible values in $[0,2\pi]$ interval. The solution propagating along
negative $z$ has been discarded for obvious physical reasons. The choice of
spatial coordinates $(x=\rho\cos\varphi,\,y=\rho\sin\varphi)$ 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]

$\int\limits_{0}^{2\pi}\\!\\!\matrixquantity(\cos(m\phi)\\\
\sin(m\phi))e^{ix\cos(\phi-\varphi)}\,\differential\phi=2\pi
i^{m}\,J_{m}(x)\matrixquantity(\cos(m\varphi)\\\ \sin(m\varphi)),$ (9)

which allows for an analytical solution for at least one of the integrals,
once the amplitude of the $\overrightarrow{\mathbf{E}}_{\infty}$ components
are written as a combination of trigonometric functions. In order to solve Eq.
(8), the $\overrightarrow{\mathbf{E}}_{\infty}$ field distribution needs to be
made explicit. Focusing through isotropic multilayered media can easily be
implemented, following the procedure described in [13, 14].

##  III Vector diffraction theory in uniaxial media

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

In solving the boundary problem for an isotropic multi-layered 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 ($\hat{\mathbf{z}}$), 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

$\boldsymbol{\varepsilon}_{2}=\matrixquantity(\varepsilon_{x}&0&0\\\
0&\varepsilon_{x}&0\\\ 0&0&\varepsilon_{z})\,.$ (10)

Each plane wave component in the angular spectrum has to satisfy the wave
equation
$\overrightarrow{\mathbf{k}}\crossproduct\overrightarrow{\mathbf{k}}\crossproduct\overrightarrow{\mathbf{E}}+\mu\varepsilon\overrightarrow{\mathbf{E}}=0$
which, using Eq. (10) and reducing the system to a two-dimensional problem
($k_{\mathrm{r}y}=0$) can be written as:

$\matrixquantity[k_{\mathrm{r}z}^{2}-\mu\varepsilon_{x}&0&-k_{\mathrm{r}x}k_{\mathrm{r}z}\\\
0&(k_{\mathrm{r}x}^{2}+k_{\mathrm{r}z}^{2})-\mu\varepsilon_{x}&0\\\
-k_{\mathrm{r}x}k_{\mathrm{r}z}&0&k_{\mathrm{r}x}^{2}-\mu\varepsilon_{z}]\overrightarrow{\mathbf{E}}=0\
,$ (11)

where the subscript $r$ denotes the “relative” wavevector
($\overrightarrow{\mathbf{k}}=k_{0}\overrightarrow{\mathbf{k}}_{r}=\left(\omega/c\right)\,\overrightarrow{\mathbf{k}}_{r}$),
as previously mentioned. The solutions for the electric field
$\overrightarrow{\mathbf{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
($\overrightarrow{\mathbf{H}}=(1/\omega\mu_{0}\mu)\overrightarrow{\mathbf{k}}\crossproduct\overrightarrow{\mathbf{E}}$).
These two modes and their fields are presented in Table III.

Table 1: Modes of the electromagnetic field in a uniaxial material. Polarisation | Condition | E-field | H-field
---|---|---|---
s | $k_{\mathrm{r}x}^{2}+k_{\mathrm{r}z}^{2}=\mu\varepsilon_{x}$ | $\overrightarrow{\mathbf{E}}^{s}=\mathrm{E}_{0}^{s}\matrixquantity(0\\\ 1\\\ 0)$ | $\overrightarrow{\mathbf{H}}^{s}=\frac{\mathrm{E}_{0}^{s}}{c\mu_{0}\mu}\matrixquantity(\mp k_{\mathrm{r}z}\\\ 0\\\ k_{\mathrm{r}x})$
p | $\frac{k_{\mathrm{r}x}^{2}}{\varepsilon_{z}}+\frac{k_{\mathrm{r}z}^{2}}{\varepsilon_{x}}=\mu$ | $\overrightarrow{\mathbf{E}}^{p}=\frac{\mathrm{E}_{0}^{p}}{\sqrt{\mu\varepsilon_{x}}}\matrixquantity(\pm k_{\mathrm{r}z}\\\ 0\\\ \mp\frac{\varepsilon_{x}}{\varepsilon_{z}}k_{\mathrm{r}x})$ | $\overrightarrow{\mathbf{H}}^{p}=\frac{\varepsilon_{x}\mathrm{E}_{0}^{p}}{\omega\mu_{0}\sqrt{\mu\varepsilon_{x}}}\matrixquantity(0\\\ 1\\\ 0)$

The solutions in Table III 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:

$\displaystyle\hat{\mathbf{e}}_{s}^{\pm}=\frac{1}{\sqrt{\mu\varepsilon_{x}}}\Bigl{(}\mp
k_{\mathrm{r}y},\pm k_{\mathrm{r}x},0\Bigr{)}\,,$ (12a)
$\displaystyle\hat{\mathbf{e}}_{p}^{\pm}=\frac{1}{\mu}\left(\pm\frac{k_{\mathrm{r}x}k_{\mathrm{r}z}}{\varepsilon_{x}},\pm\frac{k_{\mathrm{r}y}k_{\mathrm{r}z}}{\varepsilon_{x}},\mp\frac{k_{\mathrm{rT}}^{2}}{\varepsilon_{z}}\right)\,.$
(12b)

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
$\overrightarrow{\mathbf{k}}$ scale equally. In contrast, there is mixing of
transverse and longitudinal components of $\overrightarrow{\mathbf{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.

##  IV 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
$\overrightarrow{\mathbf{E}}_{\infty}$ 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]:

$\displaystyle|\mathrm{LG}_{\ell,p}|$
$\displaystyle=\frac{w_{0}}{w(z)}\sqrt{\frac{2p!}{\pi(p+|\ell|)!}}\Biggl{(}\frac{2\rho^{2}}{w(z)^{2}}\Biggr{)}^{|\ell|/2}$
$\displaystyle\times\mathbb{L}_{\ell
p}\Biggl{(}\frac{2\rho^{2}}{w(z)^{2}}\Biggr{)}\,,$ (13a)
$\displaystyle\arg\bigl{(}\mathrm{LG}_{\ell,p}\bigr{)}$
$\displaystyle=i(2p+\ell+1)\arctan\biggl(\frac{z}{z_{r}}\biggr{missing})-\frac{\rho^{2}}{w(z)^{2}}$
$\displaystyle-i\,k\frac{\rho^{2}}{R(z)}-i\ell\phi\,,$ (13b)

where $\mathbb{L}_{\ell p}$ represents the generalised Laguerre polynomials
[31], with $\ell$ 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}=\pi w_{0}^{2}/\lambda$ 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(\rho^{2}/R)\rightarrow 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
$f_{w}(\theta)=\exp[-\frac{1}{f_{0}^{2}}\frac{\sin[2](\theta)}{\sin[2](\theta_{\mathrm{max}})}]$
should be introduced. Here, $f_{0}=\frac{w_{0}}{f\sin\theta_{\mathrm{max}}}$
is a geometrical factor representing the ratio between the beam and the lens
lateral sizes. With this choice, the dependence on $\rho$ is also modified so
that
$\rho^{2}=\frac{w_{0}\sin[2](\theta)}{f_{0}\sin[2](\theta_{\mathrm{max}})}$.
With these approximations, the mode amplitude only depends on the variable
$\theta$, while its phase only on $\phi$:

$\mathrm{LG}_{\ell p}=\mathrm{C}_{\ell
p}f_{w}(\theta)\,\Theta(\theta)\,\Phi(\phi)\,,$ (14)

where $\Theta$ and $\Phi$ collect the terms depending on $\theta$ and $\phi$,
respectively, and $\mathrm{C}_{\ell p}$ all the constant factors:

$\displaystyle\mathrm{C}_{\ell p}=\sqrt{\frac{2p!}{\pi(p+|\ell|)!}}$ (15a)
$\displaystyle\Theta(\theta)=\Biggl{(}\frac{2\sin[2](\theta)}{w_{0}f_{0}\sin[2](\theta_{\mathrm{max}})}\Biggr{)}^{|\ell|/2}\mathbb{L}_{\ell
p}\Biggl{(}\frac{2\sin[2](\theta)}{w_{0}f_{0}\sin[2](\theta_{\mathrm{max}})}\Biggr{)}$
(15b) $\displaystyle\Phi(\phi)=\matrixquantity(\Phi_{x}(\phi)\\\
\Phi_{y}(\phi))\,,$ (15c)

where the two-dimensional vector ($\Phi_{x},\,\Phi_{y}$) is needed for the
description of vectorial vortices. The electric field incident on the lens
($\overrightarrow{\mathbf{E}}_{\mathrm{L}}$) can then be written in a general
way as:

$\overrightarrow{\mathbf{E}}_{\mathrm{L}}=\mathrm{C}f_{w}(\theta)\,\Theta(\theta)\,\matrixquantity(\Phi_{x}(\phi)\\\
\Phi_{y}(\phi))\,,$ (16)

where the exact definition of ($\Phi_{x},\,\Phi_{y}$) depends on the type of
vortex.

###  IV.1 Scalar vortex

In the case of a scalar vortex, $\Phi_{x}=\mathrm{E}_{0x}\exp(-i\ell\phi)$,
$\Phi_{y}=\mathrm{E}_{0y}\exp(-i\ell\phi)$, and the dependence on $\phi$ is
the same for both components resulting in uniform SOP. The beam polarisation
can then be completely described by the coefficients
$(\mathrm{E}_{0x},\mathrm{E}_{0y})$, which represent the projections of the
SOP on the polarisation basis $\\{\ket{\mathrm{H}},\,\ket{\mathrm{V}}\\}$.
Denoting these projections as $\bra{\rm H}\ket{\psi_{s}}$ and $\bra{\rm
V}\ket{\psi_{s}}$, the scalar SOP $\ket{\psi_{s}}$ can be obtained as

$\ket{\psi_{s}}=\ket{\rm H}\bra{\rm H}\ket{\psi_{s}}+\ket{\rm V}\bra{\rm
V}\ket{\psi_{s}}\,.$ (17)

###  IV.2 Vector vortex

Vectorial vortexes can be modelled as a superposition of LG beams with
opposite values of $\ell$, the same value of $p$, and orthogonal circular
polarisations. As a consequence of this, the dependence on $\theta$ of both
components of $\overrightarrow{\mathbf{E}}_{\mathrm{L}}$ is the same, while
the phase terms have different dependence on $\phi$. After introducing the
polarisation basis [16]

$\displaystyle\ket{\rm R,\ell}$ $\displaystyle=\frac{\ket{\rm H}-i\ket{\rm
V}}{\sqrt{2}}e^{i\ell\phi}\,,$ (18a) $\displaystyle\ket{\rm L,\ell}$
$\displaystyle=\frac{\ket{\rm H}+i\ket{\rm V}}{\sqrt{2}}e^{-i\ell\phi}\,,$
(18b)

the vectorial SOP $\ket{\psi_{v}}$ can be written as the sum of orthogonal
circular scalar vortices:

$\ket{\psi_{v}}=\ket{\rm R,\ell}\bra{\rm R,\ell}\ket{\psi_{v}}+\ket{\rm
L,\ell}\bra{\rm L,\ell}\ket{\psi_{v}}\,,$ (19)

where the projections on the helicity basis (Eq. 18) are
$\mathrm{E}_{0R}=\bra{\rm R,\ell}\ket{\psi_{v}}$ and $\mathrm{E}_{0L}=\bra{\rm
L,\ell}\ket{\psi_{v}}$. To consistently label SOPs in Eq. 19 and Eq. 17, the
coefficients $\mathrm{E}_{0R}$, and $\mathrm{E}_{0L}$ in $\ket{\psi_{v}}$ can
be expressed in terms of $\mathrm{E}_{0x}$ and $\mathrm{E}_{0y}$, so that

$\ket{\psi_{v}}=\frac{\mathrm{E}_{0x}+i\mathrm{E}_{0y}}{\sqrt{2}}\ket{\mathrm{R},\ell}+\frac{\mathrm{E}_{0x}-i\mathrm{E}_{0y}}{\sqrt{2}}\ket{\mathrm{L},\ell}\,.$
(20)

Within this notation framework, assuming $\ell=0$, both $\ket{\psi_{s}}$ and
$\ket{\psi_{v}}$ return uniform SOPs, expressed with the two-dimensional Jones
vector $(\mathrm{E}_{0x},\mathrm{E}_{0y})$. For example,
$(\mathrm{E}_{0x},\mathrm{E}_{0y})=(1,0)$ gives the SOP $\ket{\rm H}$ and
$(1,-i)/\sqrt{2}$ corresponds to $\ket{\rm R}$. Values of $\ell\neq 0$ will
otherwise produce scalar SOPs from $\ket{\psi_{s}}$ and vectorial ones from
$\ket{\psi_{v}}$.

The examples of various SOPs are shown in Fig. 3. The SOPs corresponding to
$\ell=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 $\ell\neq 0$ is chosen. For
example, the point associated with $\ket{H,0}$, which in a scalar case
corresponds to a horizontal SOP, represents instead a radial beam $\ket{H,1}$
for $\ell$ = 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
$\ell=\pm$1, $\pm$2 and the same values of $(\mathrm{E}_{0x},\mathrm{E}_{0y})$
are also shown.

Consistent with the above polarisation description, the functions $\Phi_{x}$
and $\Phi_{y}$ for vectorial SOPs can also be re-written in terms of
$\mathrm{E}_{0x}$ and $\mathrm{E}_{0y}$:

$\displaystyle\Phi_{x}(\phi)=\frac{1}{\sqrt{2}}\quantity\Big[\mathrm{E}_{0x}\,\cos(\ell\phi)-i\mathrm{E}_{0y}\sin(\ell\phi)]\,,$
(21a)
$\displaystyle\Phi_{y}(\phi)=\frac{1}{\sqrt{2}}\quantity\Big[\mathrm{E}_{0x}\sin(\ell\phi)+i\mathrm{E}_{0y}\cos(\ell\phi)]\,.$
(21b)

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 $\pm$1, are shown in
correspondence of the scalar SOPs having the same values of
($\mathrm{E}_{0x},\mathrm{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 $\mathrm{E}_{0x},\mathrm{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 $\pi/2$
rotation of the pattern. North and south poles represent uniform circular
polarisation regardless of the topological charge.

###  IV.3 Electric field evaluation

Adopting the same piecewise definition of the electric field as in Eq. (6),
the electric field distribution can be calculated as

$\displaystyle\overrightarrow{\mathbf{E}}_{\mathrm{i}}(\overrightarrow{\mathbf{r}})$
$\displaystyle=-\frac{ikfe^{-ikf}}{2\pi}$
$\displaystyle\times\sum\limits_{\alpha=s,p}\int\limits_{0}^{\theta_{\mathrm{max}}}\int\limits_{0}^{2\pi}\overrightarrow{\mathbf{E}}_{\infty}^{(\alpha)}(\theta,\phi)\,e^{i\overrightarrow{\mathbf{k}}^{+}_{\alpha}\dotproduct\overrightarrow{\mathbf{r}}}\sin\theta\,\differential\phi\differential\theta\,,$
(22a)
$\displaystyle\overrightarrow{\mathbf{E}}_{\mathrm{r}j}(\overrightarrow{\mathbf{r}})$
$\displaystyle=-\frac{ikfe^{-ikf}}{2\pi}$
$\displaystyle\times\sum\limits_{\alpha=s,p}\int\limits_{0}^{\theta_{\mathrm{max}}}\int\limits_{0}^{2\pi}r_{\mathrm{j}}^{(\alpha)}(\theta)\,\overrightarrow{\mathbf{E}}_{\infty}^{(\alpha)}(\theta,\phi)\,e^{i\overrightarrow{\mathbf{k}}^{-}_{\alpha}\dotproduct\overrightarrow{\mathbf{r}}}\sin\theta\,\differential\phi\differential\theta\,,$
(22b)
$\displaystyle\overrightarrow{\mathbf{E}}_{\mathrm{t}j}(\overrightarrow{\mathbf{r}})$
$\displaystyle=-\frac{ikfe^{-ikf}}{2\pi}$
$\displaystyle\times\sum\limits_{\alpha=s,p}\int\limits_{0}^{\theta_{\mathrm{max}}}\int\limits_{0}^{2\pi}t_{\mathrm{j}}^{(\alpha)}(\theta)\,\overrightarrow{\mathbf{E}}_{\infty}^{(\alpha)}(\theta,\phi)\,e^{i\overrightarrow{\mathbf{k}}^{+}_{\alpha}\dotproduct\overrightarrow{\mathbf{r}}}\sin\theta\,\differential\phi\differential\theta\,,$
(22c)

where the subscript $j$ labels the region of space and $\alpha$ labels the $s$
and $p$ components of each vector. The term
$\overrightarrow{\mathbf{E}}_{\infty}$ is obtained substituting Eq. (16) into
Eq. (8)b, once the functions $\Theta$ and $\Phi$ in Eq. (16) have been
determined for the chosen beam type. Upon substitution of
$\overrightarrow{\mathbf{k}}^{\pm}$ 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 $\phi$, 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
$\overrightarrow{\mathbf{E}}_{\mathrm{t}1}$, for an incident field given by a
generic VVB, whose SOP is given by $(\mathrm{E}_{0x},\,\mathrm{E}_{0y})$
(assuming it to be normalised to one). In the following, the electric field
taken as an example will be simply labelled as $\overrightarrow{\mathbf{E}}$
for simplicity of notations.

Table 2: Functions $\zeta^{(s,p)}_{i,\pm}(\theta)$ and
$\eta^{(s,p)}_{i}(\theta)$ in the integrals Eq. 26 calculated in the medium
$i$, for either $s$ or $p$ polarisation and propagation along
$\pm\hat{\mathbf{z}}$ direction. Labels of the Fresnel coefficients have the
same meaning as previously introduced.
$\mathbf{\zeta^{(s,p)}_{i,\pm}(\theta)}$
---
Medium 1 | $\zeta^{(s,p)}_{1,\pm}(\theta)=\zeta_{1,\pm}(\theta)=\exp\quantity\Big(\pm i\,k_{0}z\cos\theta\sqrt{\varepsilon_{1}\mu_{1}})$
Medium 2 | $\zeta^{s}_{2,\pm}(\theta)=\exp\quantity\Bigg(\pm i\,k_{0}z\sqrt{\frac{\varepsilon_{1}\mu_{1}\quantity\big(\cos 2\theta-1)}{2}+\varepsilon_{2x}\mu_{2x}}),\quad\zeta^{p}_{2,\pm}(\theta)=\exp\quantity\Bigg(\pm i\,k_{0}z\sqrt{\varepsilon_{1}\mu_{1}\quantity\big(\cos 2\theta-1)\frac{\varepsilon_{2x}}{2\varepsilon_{2z}}+\varepsilon_{2x}\mu_{2x}})$
Medium 3 | $\zeta^{(s,p)}_{3,\pm}(\theta)=\zeta_{3,\pm}(\theta)=\exp\quantity\Bigg(\pm i\,k_{0}z\sqrt{\frac{\varepsilon_{1}\mu_{1}\quantity\big(\cos 2\theta-1)}{2}+\varepsilon_{3}\mu_{3}})$
$\mathbf{\boldsymbol{\eta^{(s,p)}_{i}(\theta)}}$
| Reflection | Transmission
Medium 1 | $\quad r_{\mathrm{I}}^{s}\quad$ | $r_{\mathrm{I}}^{p}\,\cos\theta$ |  - | -
Medium 2 | $r_{\mathrm{II}}^{s}$ | $r_{\mathrm{I}}^{p}\,\sqrt{\quantity\Big[\varepsilon_{1}\mu_{1}\quantity\big(\cos 2\theta-1)]\varepsilon_{2x}}$ | $t_{\mathrm{I}}^{s}$ | $t_{\mathrm{I}}^{p}\,\sqrt{\quantity\Big[\varepsilon_{1}\mu_{1}\quantity\big(\cos 2\theta-1)]\varepsilon_{2x}}$
Medium 3 | - | - | $t_{\mathrm{II}}^{s}$ | $t_{\mathrm{II}}^{p}\,\sqrt{\varepsilon_{1}\mu_{1}\quantity\big(\cos 2\theta-1)+2\varepsilon_{3}\mu_{3}}$

In the first step, the field on the lens plane
($\overrightarrow{\mathbf{E}}_{\infty}$) is divided into its $s$ and $p$
components (Eq. (8b)) and the projections on the basis vectors are calculated:

$\displaystyle\mathrm{E}_{\infty}^{s}(\theta,\phi)$
$\displaystyle=\overrightarrow{\mathbf{E}}_{\mathrm{L}}(\theta,\phi)\dotproduct\hat{\mathbf{n}}_{\varphi}(\theta,\phi)=\bigg{(}\mathrm{E}_{0y}\cos[(1-\ell)\phi]-\mathrm{E}_{0x}\sin[(1-\ell)\phi]\bigg{)}\Theta(\theta)\,,$
(23a) $\displaystyle\mathrm{E}_{\infty}^{p}(\theta,\phi)$
$\displaystyle=\overrightarrow{\mathbf{E}}_{\mathrm{L}}(\theta,\phi)\dotproduct\hat{\mathbf{n}}_{\rho}(\theta,\phi)=\bigg{(}\mathrm{E}_{0x}\cos[(1-\ell)\phi]+\mathrm{E}_{0y}\sin[(1-\ell)\phi]\bigg{)}\Theta(\theta)\,.$
(23b)

They are then multiplied by the unit vectors given in Eq. (12), according to
the studied case, in order to obtain the full vector
$\overrightarrow{\mathbf{E}}_{\infty}^{s,p}$. After that, for each component
of the angular spectrum, the factors solely depending on $\theta$ are
collected together, yielding:

$\displaystyle\mathrm{E}_{\infty,x}^{s}(\theta,\phi)=\mathrm{C_{0}}\,\sqrt{\cos\theta}\sin\theta\,\Theta(\theta)\,\mathbf{E}_{\infty}^{s}(\theta,\phi)\sin\phi\,\Xi(\theta,\phi;\varphi)\,,$
(24a)
$\displaystyle\mathrm{E}_{\infty,x}^{p}(\theta,\phi)=\mathrm{C_{0}}\,\sqrt{\cos\theta}\sin\theta\,\Theta(\theta)\,\mathbf{E}_{\infty}^{p}\sqrt{\frac{\mu_{1}\varepsilon_{1}(\cos
2\theta-1)}{2\mu_{2}\varepsilon_{2z}}+1}\cos\phi\,\Xi(\theta,\phi;\varphi)\,,$
(24b)

with $\Xi(\theta,\phi;\varphi)=\exp[i\sqrt{\varepsilon_{1}\mu_{1}}\rho
k_{0}\cos(\phi-\varphi)\sin\theta\bigr{]}$,
$C_{0}=\mathrm{C}\sqrt{\mu_{1}\varepsilon_{1}}$.In the above expressions, the
term $\sqrt{\cos\theta}$ comes from the intensity law, while the term
$\sin\theta$ (which is in addition to the one contained in $\Theta$) comes
from the Jacobian of the integration variables transformation. In the next
step, the integration over $\phi$ is carried out applying Eq. (9) and upon
simplification of the trigonometric functions of $\phi$, when necessary. With
Wolfram Mathematica, 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 trigonometric functions, returning the corresponding
combination of Bessel functions. The list has been built using all the
possible combinations found in the calculations, to ensure the possibility to
write a general solution for any arbitrary beam LGℓp. Upon completion of this
step, the field spectral components become:

$\displaystyle\mathrm{E}^{s}_{\infty,x}(\theta;\varphi)$
$\displaystyle=\pi\mathrm{C}k_{0}i^{-\ell}\sin(\theta)\sqrt{\cos(\theta)}\Theta(\theta)\sqrt{\mu_{1}\varepsilon_{1}}$
$\displaystyle\times\biggl{[}\biggl{(}\mathrm{E}_{0x}\cos[\varphi(\ell-2)\Bigr{]}-\mathrm{E}_{0y}\sin[\varphi(\ell-2)\Bigr{]}\biggr{)}J_{2-\ell}\Bigl{(}\sqrt{\varepsilon_{1}\mu_{1}}\rho
k_{0}\sin\theta\Bigr{)}$
$\displaystyle+i^{2\ell}\biggl{(}\mathrm{E}_{0x}\cos\ell\varphi-\mathrm{E}_{0y}\sin\ell\varphi\biggr{)}J_{\ell}\Bigl{(}\sqrt{\varepsilon_{1}\mu_{1}}\rho
k_{0}\sin\theta\Bigr{)}\biggr{]}\,,$ (25a)
$\displaystyle\mathrm{E}^{p}_{\infty,x}(\theta;\varphi)$
$\displaystyle=\pi\mathrm{C}k_{0}i^{-\ell}\sin(\theta)\sqrt{\cos(\theta)}\Theta(\theta)\sqrt{\mu_{1}\varepsilon_{1}}\sqrt{\frac{(\cos(2\theta)-1)(\mu_{1}\varepsilon_{1})}{2\mu_{2}\varepsilon_{2z}}+1}$
$\displaystyle\times\biggl{[}\biggl{(}\mathrm{E}_{0y}\sin[\varphi(\ell-2)\biggr{]}-\mathrm{E}_{0x}\cos[\varphi(\ell-2)\Bigr{]}\Bigr{)}J_{2-\ell}\Bigl{(}\sqrt{\varepsilon_{1}\mu_{1}}\rho
k_{0}\sin\theta\Bigr{)}$
$\displaystyle+i^{2\ell}\Bigl{(}\mathrm{E}_{0x}\cos\ell\varphi-\mathrm{E}_{0y}\sin\ell\varphi\Bigr{)}J_{\ell}\Bigl{(}\sqrt{\varepsilon_{1}\mu_{1}}\rho
k_{0}\sin\theta\Bigr{)}\biggr{]}\,.$ (25b)

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:

$\displaystyle\mathcal{I}_{1}$
$\displaystyle=\int_{0}^{\theta_{\mathrm{max}}}\,J_{\ell}\Bigl{(}k_{0}\sqrt{\varepsilon_{1}\mu_{1}}\rho\sin\theta\Bigr{)}\,f_{w}(\theta)\sqrt{\cos\theta}\sin\theta\,\quantity\big(\cos\theta+1)\,\Theta(\theta)\,\zeta^{(s,p)}_{i,\pm}(\theta)\,\eta^{(s,p)}_{i}(\theta)\,\differential\theta$
(26a) $\displaystyle\mathcal{I}_{2}$
$\displaystyle=\int_{0}^{\theta_{\mathrm{max}}}\,J_{2-\ell}\Bigl{(}k_{0}\sqrt{\varepsilon_{1}\mu_{1}}\rho\sin\theta\Bigr{)}\,f_{w}(\theta)\sqrt{\cos\theta}\sin\theta\,\quantity\big(\cos\theta-1)\,\Theta(\theta)\,\zeta^{(s,p)}_{i,\pm}(\theta)\,\eta^{(s,p)}_{i}(\theta)\,\differential\theta$
(26b) $\displaystyle\mathcal{I}_{3}$
$\displaystyle=\int_{0}^{\theta_{\mathrm{max}}}\,J_{2+\ell}\Bigl{(}k_{0}\sqrt{\varepsilon_{1}\mu_{1}}\rho\sin\theta\Bigr{)}\,f_{w}(\theta)\sqrt{\cos\theta}\sin\theta\,\quantity\big(\cos\theta-1)\,\Theta(\theta)\,\zeta^{(s,p)}_{i,\pm}(\theta)\,\eta^{(s,p)}_{i}(\theta)\,\differential\theta$
(26c) $\displaystyle\mathcal{I}_{4}$
$\displaystyle=\int_{0}^{\theta_{\mathrm{max}}}\,J_{1-\ell}\Bigl{(}k_{0}\sqrt{\varepsilon_{1}\mu_{1}}\rho\sin\theta\Bigr{)}\,f_{w}(\theta)\sqrt{\cos\theta}\sin^{2}\theta\,\Theta(\theta)\,\zeta^{p}_{i,\pm}(\theta)\,\eta^{p}_{i}(\theta)\,\differential\theta$
(26d) $\displaystyle\mathcal{I}_{5}$
$\displaystyle=\int_{0}^{\theta_{\mathrm{max}}}\,J_{1+\ell}\Bigl{(}k_{0}\sqrt{\varepsilon_{1}\mu_{1}}\rho\sin\theta\Bigr{)}\,f_{w}(\theta)\sqrt{\cos\theta}\sin^{2}\theta\,\Theta(\theta)\,\zeta^{p}_{i,\pm}(\theta)\,\eta^{p}_{i}(\theta)\,\differential\theta\,.$
(26e)

The explicit form of $\Theta$ depends on the beam type. The functions
$\zeta(\theta)$ and $\eta(\theta)$ represent quantities whose explicit
expression depends on the medium (labelled by $i$), polarisation ($s$ or $p$)
and propagation direction ($\pm\hat{\mathbf{z}}$) they are calculated for
(Table 2). The former ($\zeta$) is the propagation factor and the latter
($\eta$) is a general representation of the Fresnel coefficients needed in
each case. It is worth noting that in the case of the integrals
$\mathcal{I}_{4}$ and $\mathcal{I}_{5}$, only the contribution of the $p$
component of the field is important. Finally, the numerical integration is
performed in MATLAB.

Figure 4: Maps of the integrands $\xi_{\ell\pm 1}$ in Eqs. (26) plotted for
$\ell=-2,0,2$. The colour scale is the same for the three maps.

Similar procedure can be applied to obtain the magnetic field, retrieving its
angular spectrum from the electric field one.

Table 3: Expression of the electric field $\overrightarrow{\mathbf{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 $(\mathrm{E}_{0x},\mathrm{E}_{0y})$ and topological
charge $\ell$. 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 $\eta$ represent transmission coefficients,
according to Table 2. Electric field in medium 2, propagating along positive
$\mathbf{z}$
---
|  | $\mathrm{C_{0}}\,\Big{[}\mathrm{E}_{0x}\quantity\big(f_{\varphi}^{10}-if_{\varphi}^{20})\quantity\Big(\sqrt{2}\mathcal{I}_{1}^{1p}+2\sqrt{\varepsilon_{2x}\varepsilon_{2z}\mu_{2}}\mathcal{I}_{1}^{1s})+\frac{i^{-2\ell}}{2}(i\mathrm{E}_{0x}+\mathrm{E}_{0y})\quantity\big(if_{\varphi}^{12}+f_{\varphi}^{22})\quantity\Big(\sqrt{2}\mathcal{I}_{2}^{1p}-2\sqrt{\varepsilon_{2x}\varepsilon_{2z}\mu_{2}}\mathcal{I}_{2}^{1s})+$
| $\,x\,$ | $-\frac{1}{2}(\mathrm{E}_{0x}+i\mathrm{E}_{0y})\quantity\big(f_{\varphi}^{13}-if_{\varphi}^{23})\quantity\Big(\sqrt{2}\mathcal{I}_{3}^{1p}-2\sqrt{\varepsilon_{2x}\varepsilon_{2z}\mu_{2}}\mathcal{I}_{3}^{1s})\Big{]}$
|  | $\mathrm{C_{0}}\,\Big{[}\mathrm{E}_{0y}\quantity\big(f_{\varphi}^{10}-if_{\varphi}^{20})\quantity\Big(\sqrt{2}\mathcal{I}_{1}^{1p}+2\sqrt{\varepsilon_{2x}\varepsilon_{2z}\mu_{2}}\mathcal{I}_{1}^{1s})+\frac{i^{-2\ell}}{2}(i\mathrm{E}_{0x}+\mathrm{E}_{0y})\quantity\big(f_{\varphi}^{12}-if_{\varphi}^{22})\quantity\Big(\sqrt{2}\mathcal{I}_{2}^{1p}-2\sqrt{\varepsilon_{2x}\varepsilon_{2z}\mu_{2}}\mathcal{I}_{2}^{1s})+$
| $\,y\,$ | $-\frac{1}{2}(i\mathrm{E}_{0x}-\mathrm{E}_{0y})\quantity\big(f_{\varphi}^{13}-if_{\varphi}^{23})\quantity\Big(\sqrt{2}\mathcal{I}_{3}^{1p}-2\sqrt{\varepsilon_{2x}\varepsilon_{2z}\mu_{2}}\mathcal{I}_{3}^{1s})\Big{]}$
OVB | $\,z\,$ | $\mathrm{C_{0}}\,\frac{2\,i^{-2\ell}\varepsilon_{2x}\sqrt{\varepsilon_{1}\mu_{1}}}{\sqrt{\varepsilon_{2z}}}\quantity\big(f_{\varphi}^{10}-if_{\varphi}^{20})\quantity\Big[\quantity\big(\mathrm{E}_{0x}-i\mathrm{E}_{0y})\quantity\big(-if_{\varphi}^{01}+f_{\varphi}^{02})\mathcal{I}_{4}^{1p}+i^{2\ell}\quantity\big(-i\mathrm{E}_{0x}+i\mathrm{E}_{0y})\quantity\big(f_{\varphi}^{01}-if_{\varphi}^{02})\mathcal{I}_{5}^{1p}]$
| $\,x\,$ | $\mathrm{C_{0}}\quantity\Big[\,\quantity\big(\mathrm{E}_{0x}f_{\varphi}^{10}-\mathrm{E}_{0y}f_{\varphi}^{20})\quantity\Big(\sqrt{2}\mathcal{I}_{1}^{1p}+2\sqrt{\varepsilon_{2x}\varepsilon_{2z}\mu_{2}}\mathcal{I}_{1}^{1s})-i^{-2\ell}\quantity\big(\mathrm{E}_{0x}f_{\varphi}^{12}-\mathrm{E}_{0y}f_{\varphi}^{22})\quantity\Big(\sqrt{2}\mathcal{I}_{2}^{1p}-2\sqrt{\varepsilon_{2x}\varepsilon_{2z}\mu_{2}}\mathcal{I}_{1}^{1s})]$
| $\,y\,$ | $\mathrm{C_{0}}\,\quantity\Big[\quantity\big(\mathrm{E}_{0x}f_{\varphi}^{10}+\mathrm{E}_{0y}f_{\varphi}^{20})\quantity\Big(\sqrt{2}\mathcal{I}_{1}^{1p}+2\sqrt{\varepsilon_{2x}\varepsilon_{2z}\mu_{2}}\mathcal{I}_{1}^{1s})+i^{-2\ell}\quantity\big(\mathrm{E}_{0x}f_{\varphi}^{12}+\mathrm{E}_{0y}f_{\varphi}^{22})\quantity\Big(\sqrt{2}\mathcal{I}_{2}^{1p}-2\sqrt{\varepsilon_{2x}\varepsilon_{2z}\mu_{2}}\mathcal{I}_{1}^{1s})]$
VVB | $\,z\,$ | $-4\mathrm{C_{0}}\,i^{1-2\ell}\sqrt{\varepsilon_{1}\varepsilon_{2x}\mu_{1}}{\sqrt{\varepsilon_{1}\varepsilon_{2z}\mu_{1}}}\quantity\big(\mathrm{E}_{0x}f_{\varphi}^{11}+\mathrm{E}_{0y}f_{\varphi}^{21})\mathcal{I}_{4}^{1p}$
Factors appearing in the expressions above
$\mathrm{C_{0}}=i^{\ell}\pi\mathrm{C}k_{0}\frac{\sqrt{\varepsilon_{1}\mu_{1}}}{2\sqrt{\varepsilon_{2x}\varepsilon_{2z}\mu_{2}}}$
$f_{\varphi}^{01}=\cos\varphi\qquad f_{\varphi}^{02}=\sin\varphi$
$f_{\varphi}^{10}=\cos\ell\varphi\qquad
f_{\varphi}^{11}=\cos\quantity\Big[\quantity\big(1-\ell)\varphi]\qquad
f_{\varphi}^{12}=\cos\quantity\Big[\quantity\big(2-\ell)\varphi]\qquad
f_{\varphi}^{13}=\cos\quantity\Big[\quantity\big(2+\ell)\varphi]$
$f_{\varphi}^{20}=\sin\ell\varphi\qquad
f_{\varphi}^{21}=\sin\quantity\Big[\quantity\big(1-\ell)\varphi]\qquad
f_{\varphi}^{22}=\sin\quantity\Big[\quantity\big(\ell-2)\varphi]\qquad
f_{\varphi}^{23}=\sin\quantity\Big[\quantity\big(2+\ell)\varphi]$

To provide a practical example, Table 3 presents the expressions obtained for
the electric field $\overrightarrow{\mathbf{E}}_{t1}$ in both the cases of OVs
and VVBs, where the dependence on the coordinates $(\rho,\varphi,z)$ have been
omitted for simplified notations.

It is worth noting that, when calculating the field for a VVB, the three
$\theta$-integrals ($\mathcal{I}_{1}$, $\mathcal{I}_{2}$, and
$\mathcal{I}_{4}$) need to be evaluated, while the five integrals are needed
for a scalar vortex. This is related to the presence of the $\phi$-dependent
phase factor of the LG modes, which is lost in the superposition of orthogonal
circular vortices with opposite topological charge, needed for VVBs.

###  IV.4 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 interesting 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 $\mathcal{I}_{4}$ and $\mathcal{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,
$\mathcal{I}_{4}$ and $\mathcal{I}_{5}$ contain the terms
$\xi_{1\pm\ell}(\rho,\theta)=J_{1\pm\ell}\Bigl{(}k_{0}\sqrt{\mu_{1}\varepsilon_{1}}\rho\sin\theta\Bigr{)}\sqrt{\cos\theta}\sin^{2}\theta$,
which are calculated for each value of $\rho$ over all the range of values of
the integration variable $\theta$, 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 $\mathcal{I}_{4}$ and $\mathcal{I}_{5}$ show a case for
$\ell=\pm 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 $\rho=0$.
Comparing the maps of $\xi_{\ell}(\rho,\theta)$ for $\ell=0,\pm 1$ (Fig. 4),
it is evident how $J_{0}$ is the only case where for small values of $\rho$
the integral is nonzero. Looking at the same maps for higher values of $\rho$,
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 $\ell=-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 $\ell=1$
for a VVB on a higher-order 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 $\varphi$
is made explicit is:

$\displaystyle\mathrm{E}^{(\mathrm{VVB})}_{z}(\rho,\varphi,z)$
$\displaystyle=-4i\varepsilon_{2z}\sqrt{\mu_{1}\varepsilon_{1}}\mathrm{C_{0}}\,J_{1-\ell}\Bigl{(}k_{0}\sqrt{\mu_{1}\varepsilon_{1}}\rho\sin\theta\Bigr{)}\sqrt{\cos\theta}\sin^{2}\theta\,,$
$\displaystyle\times\biggl{(}\mathrm{E}_{0x}\cos[(1-\ell)\varphi\Bigr{]}+\mathrm{E}_{0y}\sin[(1-\ell)\varphi\Bigr{]}\biggr{)}$
(27a) $\displaystyle\biggl{\\{}\ell=1\implies\biggr{\\}}$
$\displaystyle=-4i\varepsilon_{2z}\sqrt{\mu_{1}\varepsilon_{1}}\mathrm{C_{0}}\,J_{0}\Bigl{(}k_{0}\sqrt{\mu_{1}\varepsilon_{1}}\rho\sin\theta\Bigr{)}\sqrt{\cos\theta}\sin^{2}\theta\,\mathrm{E}_{0x}\,,$
(27b)

so that the only contribution to the longitudinal field can come from the
component of the beam on the state $\ket{\mathrm{H},1}$
($(\mathrm{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,\mathrm{E}_{0y})$, which
corresponds to an azimuthal beam, the longitudinal field will be zero in every
point of the real space.

Figure 5: Vortex beam focusing in free space with NA = $0.9$: (a) Gaussian
beam, (b) scalar vortex beam LG15,0, (c) scalar vortex beam LG10,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.

##  V 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)
$\mathrm{E}_{0x}$, $\mathrm{E}_{0y}$, $ell$, $p$, $w_{0}$ and $\lambda$ 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 (LG00); a high-order scalar vortex (LG15,0); a
scalar vortex with a nonzero radial index (LG10,2); and a vectorial vortex
whose constituent vortices are given by LG±10,2 beams, referred to as the
VVB10,2 (Fig. 5). The wavelength of choice is $\lambda$ = 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 $\hat{\mathbf{x}}$ direction), hence only
one of the initial electric field components is nonzero ($\mathrm{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 polarisation: 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 $\mathrm{E}_{0x},\mathrm{E}_{0y}$ and $\ell$ (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 LG15,0 vortex, reaching a $35\%$
for LG10,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.

Figure 6: Propagation of (a) LG10,2 and (b) VVB10,2 beams focused with NA =
0.9 in a three-layered system. The central layer can be a (a, yellow frames)
isotropic ($\varepsilon_{2x}=\varepsilon_{2z}=2.25+i5\times 10^{-2}$), (b,
purple frames) elliptic ($\varepsilon_{2x}=2.25+i5\times 10^{-2}$,
$\varepsilon_{2z}=1.9+i5\times 10^{-2}$), (c, blue frames) hyperbolic
($\varepsilon_{2x}=2.25+i5\times 10^{-2}$, $\varepsilon_{2z}=-1.9+i5\times
10^{-2}$) medium, while the first and third layers are chosen as free space
($\varepsilon_{1}=\varepsilon_{3}$ = 1). The description of the content of the
panels is the same as in Fig. 5. The multiplication factors refer to the panel
they are shown in.

#### 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 LG10,2 and VVB10,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
$\varepsilon_{2x}=\varepsilon_{2z}=2.25+i5\times 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 $\hat{\mathbf{z}}$ direction of the
reference frame, which coincides with the normal to the interfaces. Depending
on the sign of the permittivity components $\real(\varepsilon_{2x})$ and
$\real(\varepsilon_{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 $\varepsilon_{2x}=2.25+i5\times 10^{-2}$,
$\varepsilon_{2z}=1.9+i5\times 10^{-2}$. In this case, the transverse field,
affected by $\varepsilon_{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 $\real(\varepsilon_{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 $\imaginary(\varepsilon_{2x})$ and
$\imaginary(\varepsilon_{2z})$ would produce different damping rates for
transverse and longitudinal components and interplay between them).

#### Uniaxial slab: hyperbolic dispersion.

If eventually $\real(\varepsilon_{2x})\real(\varepsilon_{2z})<0$, the second
layer has a hyperbolic dispersion (Fig. 2, blue frames). The chosen
permittivities for this case are: $\varepsilon_{2x}=2.25+i5\times 10^{-2}$,
$\varepsilon_{2z}=-1.9+i5\times 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 [37, 7].

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

##  VI 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 hyperbolic 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 open-source 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.

## Acknowledgemnts.

This work was supported in part by the ERC iCOMM Project (No. 789340) and the
ERC Starting Grant No. ERC-2016-STG-714151-PSINFONI. All the data supporting
finding of this work are presented in the Results section and are available
from the corresponding author upon reasonable request.

## References

  * Allen  _et al._ [1992] L. Allen, M. W. Beijersbergen, R. Spreeuw, and J. Woerdman, Orbital angular momentum of light and the transformation of Laguerre-Gaussian laser modes, Physical review A 45, 8185 (1992).
  * Shen  _et al._ [2019] Y. Shen, X. Wang, Z. Xie, C. Min, X. Fu, Q. Liu, M. Gong, and X. Yuan, Optical vortices 30 years on: OAM manipulation from topological charge to multiple singularities, Light: Science & Applications 8, 90 (2019).
  * Naidoo  _et al._ [2016] D. Naidoo, F. S. Roux, A. Dudley, I. Litvin, B. Piccirillo, L. Marrucci, and A. Forbes, Controlled generation of higher-order Poincaré sphere beams from a laser, Nature Photonics 10, 327 (2016).
  * Tang  _et al._ [2018] D. Tang, J. Wang, J. Ma, B. Zhou, P. Yuan, G. Xie, H. Zhu, and L. Qian, Theoretical study on second-harmonic generation of focused vortex beams, Laser Physics Letters 15, 035402 (2018).
  * Wang  _et al._ [2023] H. Wang, H. Wang, Q. Ruan, J. Y. E. Chan, W. Zhang, H. Liu, S. D. Rezaei, J. Trisno, C.-W. Qiu, M. Gu, _et al._ , Coloured vortex beams with incoherent white light illumination, Nature Nanotechnology 18, 264 (2023).
  * Piccardo  _et al._ [2023] M. Piccardo, M. de Oliveira, V. R. Policht, M. Russo, B. Ardini, M. Corti, G. Valentini, J. Vieira, C. Manzoni, G. Cerullo, _et al._ , Broadband control of topological–spectral correlations in space–time beams, Nature Photonics , 1 (2023).
  * Afanasev  _et al._ [2023] A. Afanasev, J. J. Kingsley-Smith, F. J. Rodríguez-Fortuño, and A. V. Zayats, Nondiffractive three-dimensional polarization features of optical vortex beams, Advanced Photonics Nexus 2, 026001 (2023).
  * Maurer  _et al._ [2007] C. Maurer, A. Jesacher, S. Fürhapter, S. Bernet, and M. Ritsch-Marte, Tailoring of arbitrary optical vector beams, New Journal of Physics 9, 78 (2007).
  * Zhan [2009] Q. Zhan, Cylindrical vector beams: from mathematical concepts to applications, Advances in Optics and Photonics 1, 1 (2009).
  * Krasavin  _et al._ [2018] A. V. Krasavin, P. Segovia, R. Dubrovka, N. Olivier, G. A. Wurtz, P. Ginzburg, and A. V. Zayats, Generalization of the optical theorem: experimental proof for radially polarized beams, Light: Science & Applications 7, 36 (2018).
  * Shen  _et al._ [2023] Y. Shen, Q. Zhang, P. Shi, L. Du, X. Yuan, and A. V. Zayats, Optical skyrmions and other topological quasiparticles of light, Nature Photonics , 1 (2023).
  * Mansuripur [1986] M. Mansuripur, Distribution of light at and near the focus of high-numerical-aperture objectives, JOSA A 3, 2086 (1986).
  * Quabis  _et al._ [2000] S. Quabis, R. Dorn, M. Eberler, O. Glöckl, and G. Leuchs, Focusing light to a tighter spot, Optics communications 179, 1 (2000).
  * Youngworth and Brown [2000] K. S. Youngworth and T. G. Brown, Focusing of high numerical aperture cylindrical-vector beams, Optics Express 7, 77 (2000).
  * Biss and Brown [2004] D. Biss and T. Brown, Cylindrical vector beam focusing through a dielectric interface: reply to comment, Optics Express 12, 970 (2004).
  * Milione  _et al._ [2011] G. Milione, H. Sztul, D. Nolan, and R. Alfano, Higher-order Poincaré sphere, stokes parameters, and the angular momentum of light, Physical review letters 107, 053601 (2011).
  * Stratton and Chu [1939] J. A. Stratton and L. Chu, Diffraction theory of electromagnetic waves, Physical Review 56, 99 (1939).
  * Wolf [1959] E. Wolf, Electromagnetic diffraction in optical systems-I. An integral representation of the image field, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 253, 349 (1959).
  * Wolf and Li [1981] E. Wolf and Y. Li, Conditions for the validity of the Debye integral representation of focused fields, Optics Communications 39, 205 (1981).
  * Sheppard  _et al._ [1977] C. J. R. Sheppard, A. Choudhury, and J. Gannaway, Electromagnetic field near the focus of wide-angular lens and mirror systems, IEE Journal of Microwaves Optics and Acoustics 1, 129 (1977).
  * Richards  _et al._ [1959] B. Richards, E. Wolf, and D. Gabor, Electromagnetic diffraction in optical systems, II. Structure of the image field in an aplanatic system, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 253, 358 (1959).
  * Whittaker [1902] E. Whittaker,  _A Course of Modern Analysis: An Introduction to the General Theory of Infinite Series and of Analytic Functions, with an Account of the Principal Transcendental Functions_ (University Press, 1902).
  * Novotny and Hecht [2012a] L. Novotny and B. Hecht, Theoretical foundations, in  _Principles of Nano-Optics_ (Cambridge University Press, 2012) p. 12–44, 2nd ed.
  * Mandel and Wolf [1995] L. Mandel and E. Wolf, Some useful mathematical techniques, in  _Optical Coherence and Quantum Optics_ (Cambridge University Press, 1995) p. 92–146.
  * Lekner [1994] J. Lekner, Optical properties of a uniaxial layer, Pure and Applied Optics: Journal of the European Optical Society Part A 3, 821 (1994).
  * Wolf and Gabor [1959] E. Wolf and D. Gabor, Electromagnetic diffraction in optical systems - i. an integral representation of the image field, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 253, 349 (1959).
  * Novotny and Hecht [2012b] L. Novotny and B. Hecht, Propagation and focusing of optical fields, in  _Principles of Nano-Optics_ (Cambridge University Press, 2012) p. 45–85, 2nd ed.
  * Picardi  _et al._ [2017] M. F. Picardi, A. Manjavacas, A. V. Zayats, and F. J. Rodríguez-Fortuño, Unidirectional evanescent-wave coupling from circularly polarized electric and magnetic dipoles: An angular spectrum approach, Physical Review B 95, 245416 (2017).
  * Jackson [2003] J. D. Jackson, Electrodynamics, classical, in  _digital Encyclopedia of Applied Physics_ (John Wiley & Sons, Ltd, 2003).
  * Saleh and Teich [1991] B. E. Saleh and M. C. Teich, Beam optics, in  _Fundamentals of Photonics_ (John Wiley & Sons, Ltd, 1991) Chap. 3, pp. 80–107.
  * Abramowitz and Stegun [1968] M. Abramowitz and I. A. Stegun,  _Handbook of mathematical functions with formulas, graphs, and mathematical tables_ , Vol. 55 (US Government printing office, 1968).
  * Scully and Zubairy [1991] M. O. Scully and M. S. Zubairy, Simple laser accelerator: Optics and particle dynamics, Physical Review A 44, 2656 (1991).
  * Dorn  _et al._ [2003] R. Dorn, S. Quabis, and G. Leuchs, Sharper focus for a radially polarized light beam, Physical review letters 91, 233901 (2003).
  * Veselago [1968] V. G. Veselago, Reviews of topical problems: the electrodynamics of substances with simultaneously negative values of $\epsilon$ and $\mu$, Soviet Physics Uspekhi 10, R04 (1968).
  * Pendry [2000] J. B. Pendry, Negative refraction makes a perfect lens, Physical review letters 85, 3966 (2000).
  * Smith and Schurig [2003] D. Smith and D. Schurig, Electromagnetic wave propagation in media with indefinite permittivity and permeability tensors, Physical Review Letters 90, 077405 (2003).
  * Aita and Zayats [2023] V. Aita and A. V. Zayats, Enhancement of optical spin-orbit coupling in anisotropic enz metamaterials, IEEE Photonics Journal 15, 1 (2023).
  * Aita and Rodríguez-Fortuño [2024] V. Aita and F. J. Rodríguez-Fortuño, Infocus - interaction of focused complex beams with uniaxial slabs 10.18742/25133507.v1 (2024).