## 3 Numerical Simulations and Results
### Methods Used for Numerical Simulations  
The spectrum of plane-waves method is well known, easy to implement, and efficient for beams with nearly plane wavefront in the near-field. It involves Fourier transforming the input signal and adding a phase, in accordance with the dispersion relation of the waves and distance between input and output planes, to determine the output signal. Mathematically, it is described in Equation (8). \[\textbf{V}(\textbf{r},z)\ =\ \iint\textbf{V}_{0}(\textbf{r}_{0})\ \exp\left[ \text{i}\ \sqrt{k_{0}^{2}-k^{2}}\,z+\text{i}\,\textbf{k}\cdot(\textbf{r}-\textbf{r}_{0}) \right]\text{d}\textbf{k}\text{d}\textbf{r}_{0}, \tag{8}\]where \(k_{0}=\omega/c\) and \(\mathbf{k}\) is the transversal wave-vector of a plane-wave. We use this formula in the case of variation 1 to simulate the propagation of the beam over the distance \(z=1\) m, between SPP and MC, starting with a Gaussian beam with diameter 16 mm, after we multiplied the field with the complex transmission function of SPP (topological charge \(\updownarrow=-1\)). Ray tracing is used everywhere in optics. Therefore, here, we discuss only the way in which we use it to compute the wavefront aberration produced by the combination MC + lens and how we reconstruct the field from ray information traced through an optical component. Firstly, let us define the ray vector (\(\mathbf{t}\), \(t_{z}\)) which points in the direction of the ray and has the property that \(k_{0}\)\(\mathbf{t}=\mathbf{k}\). Let us imagine a divergent bundle of rays that originates from a single point. One ray from the bundle intersects a plane surface placed at distance \(z\) from the starting point, at coordinate \(\mathbf{r}_{s}(\mathbf{t})=z\,\boldsymbol{\mathsf{y}}t_{z}\), while the optical path length along the ray is \(s(\mathbf{t})=z/t_{z}\). We can observe now that the complex amplitude produced by the plane-wave indexed by \(\mathbf{k}\) in the plane \(z\), which usually can be written as \(\exp[\mathrm{i}\,\,k_{z}\,(\mathbf{k})\,z\,+\,\mathrm{i}\,\,\mathbf{k}\cdot \mathbf{r}]\), where \(k_{z}\,(\mathbf{k})=(k^{2}_{0}-k^{2})^{1/2}\), could also be written as \(\exp[\mathrm{i}\,\,k_{0}\,s(\mathbf{k})+\mathrm{i}\,\,\mathbf{k}\cdot \mathbf{r}-\mathrm{i}\,\,\mathbf{k}\cdot\mathbf{r}_{s}(\mathbf{k})]\). Therefore, we can conclude that information about the ray intersection with the surface and optical path length along the ray can be used to model the propagation of a plane-wave. Now, we have to determine the additional phase that each plane-wave acquires in the focal plan due to the aberrations of the system. \[\Phi(\mathbf{k})\,=\,\exp[\mathrm{i}\,k_{0}\,s_{1}(\mathbf{k})+\mathrm{i}\,k_ {0}\,s_{2}(\mathbf{k})-\mathrm{i}\,\mathbf{k}\cdot\mathbf{r}_{s}(\mathbf{k})] \approx\exp[\mathrm{i}\,k_{0}\,s_{1}(\mathbf{k})+\mathrm{i}\,\,\sqrt{k_{0}^{2} -k^{2}}\,f]. \tag{9}\]  
Equation (9) evaluates the aberrations in the focal plane, usually evaluated in the plane tangent to the rear face of the last component or on the surface of a sphere centered in the focal point. \(s_{1}(\mathbf{k})\) is the optical path length along the ray, between the input plane and the plane tangent to the last surface, \(s_{2}(\mathbf{k})\) is the optical path length along the ray indexed by \(\mathbf{k}\) between this plane and the focal plane, and \(\mathbf{r}_{s}\) is the point of intersection between the ray and the focal plane. The approximate value from Equation (9) is equivalent to the wavefront aberration as is usually employed for the Debye-Wolf integral since it is given up to a constant by the optical path length between the input plane and the sphere of radius \(f\) centered in the focal point. We used the correct value from Equation (9), which also takes into consideration that the plane-waves are not directed exactly to the focal point.