## 2 Materials and Methods
### Numerical Simulations  
In order to determine the far-field amplitude and phase distributions of vortex beams, a computer program was used to simulate the propagation of an incident beam with an arbitrary amplitude and phase distribution through a focusing lens. The code, written in Wolfram Mathematica, sums the Huygens wavelet contributions using the Fresnel approximation from the near field to points in the far field [22]. In this paper, we model a flattop beam passing through a spiral phase plate and a perfect focusing lens, however, the computer code can be used to calculate the effect of a wide range of beam profiles, phase distributions and optical aberrations. Since we use an SLM to imprint phase patterns onto an incident beam, the phase distribution is written modulo \(2\pi\), allowing for the full use of the dynamic range of the SLM in the zero to \(2\pi\) range, while still allowing for the equivalent of arbitrarily high phases to be expressed. For extremely short pulses, where the number of optical cycles in the pulse is comparable to the value of the charge \(m\), using the full phase pattern would result in a temporal \"walk-off\" of the different phase-shifted portions of the beam within the focal region. Using the modulo \(2\pi\) phase has the added benefit of limiting these ultrashort pulse width effects. In Figure 1, we present an example of the code and output of the simulations done in Wolfram Mathematica. The program from Figure 1a computes the radial intensity plot of a vortex through paraxial approximation using Equation (5), and computing \(\left|\psi_{m}(R,z)\right|^{2}\). Figure 1b,c correspond to the 2D and 3D plots of the radial intensity of a vortex with a topological charge \(m=15\), computed with respect to a horizontal rate \(R/\sqrt{z}\), respectively. The superposition of a pair of vortices is simulated by implementing code that allows the visualization of the intensity profile of various charge vortices limited by a Gaussian beam profile and then computes their interference, based on Equations (6)-(8). The simulated intensity patterns are presented in Figures 2 and 3, for equally and opposite charged coplanar and non-coplanar vortices, respectively. In these simulations, no free parameters are used, in order to compare the experimental results directly with the simulations presented in this section. For two equally charged coplanar vortices, there should only be constructive interference since both vortices are identical, as presented in Figure 2a. When both coplanar vortices have opposite charges, there should be destructive interference producing a number of fringes around the vortex equal to twice the charge, as presented in Figure 2b. A phase tilt can be added to the numerical results to better simulate what occurs in the laboratory, where individual vortices intersect at a small angle. For equally charged vortices, the phase tilt causes interference fringes to appear in the vortex. The number of fringes increases as the phase tilt does. These fringes run along one direction, as the typical fringes seen when two coherent and in-phase laser beams intersect. When two oppositely charged vortices are superimposed with a phase tilt, the result is an interference pattern that has varying spacing between the interference fringes around the vortex. This pattern is known as the firebird pattern, becoming more pronounced as the phase tilt increases. Figure 3 shows the described interference patterns. Figure 1: Vortex simulation in Mathematica. (**a**) Code of paraxial approximation to generate (**b**) the intensity plot and (**c**) the 3D pattern of a focused vortex. It is also of interest to obtain the intensity patterns of the superposition of multiple wavelength vortices. Through the use of a blazed fork pattern to diffract vortices into the first order, a wide range of wavelengths can be used to generate quality vortices without the need for multiple spatial light modulators.