## 2 Materials and Methods
### Analytic Solutions  
Following the discussion from [PERSON] [23] for the propagation of the vortex, the transverse wavevectors follow the Poynting vector as it spirals around the azimuth. The wave can be written as a superposition of plane waves with the transverse wavevectors for a number of phase steps all the way around the azimuth. The wave is then written as a summation over the transverse wavevectors, expressed as  
\[\psi_{m}(R,z)=\sum_{K}A_{m}e^{i(K\cdot R+z\sqrt{1-K^{2}})}. \tag{1}\]  
As the number of phase steps approaches infinity, or in the limit as \(K\rightarrow\infty\) around the azimuth, the summation can be replaced as an integral over all transverse wavevectors, \(K\)  
\[\psi_{m}(R,z)=\lim_{K\rightarrow\infty}\sum_{K}A_{m}e^{i(K\cdot R+z\sqrt{1-K^{ 2}})}=\iint_{K}A_{m}e^{i(K\cdot R+z\sqrt{1-K^{2}})}dK. \tag{2}\]  
Considering the orthogonality relation between the basis functions, the amplitude \(A_{m}\) can be determined as a plane wave with the added vortex phase from the near field, \(\psi_{m}(R,z=0)=e^{im\theta}\). Multiplying both sides of Equation (2) by \(e^{-i(K\cdot R^{\prime})}\), and integrating over \(R\), at \(z=0\)  
\[ \begin{array}{ll}\iint_{R}\psi_{m}(R,z)e^{-i(K\cdot R^{\prime})}dR=\iint_{R} \iint_{K}A_{m}e^{i(K\cdot R+z\sqrt{1-K^{2}}-K\cdot R^{\prime})}dKdR\\ \iint_{R}e^{im\theta}e^{-i(K\cdot R^{\prime})}dR=\iint_{K}A_{m}e^{iz\sqrt{1- K^{2}}}\int_{R}e^{i(K\cdot R-K\cdot R^{\prime})}\ dKdR\\ =\iint_{K}A_{m}e^{iz\sqrt{1-K^{2}}}\delta(R-R^{\prime})&dK=\iint_{K}A_{m}e^{i \sqrt{1-K^{2}}}dK\\ \iint_{R}e^{im\theta}e^{-i(K\cdot R^{\prime})}dR&=A_{m}\iint_{K}dK=A_{m}4\pi ^{2}\end{array}  \tag{3}\]  
The expression from Equation (3) is then solved for the amplitude, and the remaining integral can be solved in terms of Bessel functions,  
\[A_{m}=\frac{1}{4\pi^{2}}\iint_{R}e^{-i(K\cdot R-m\theta)}dR=J_{|m|}(KR)\frac{ 1}{K}e^{im\theta}=\frac{|m|(-i)^{|m|}}{2\pi K^{2}}e^{im\theta}, \tag{4}\]  
where \(J_{|m|}(KR)\) is a Bessel function of the first kind, implemented as \(J_{|m|}(KR)=(-1)^{m}J_{m}(KR)\), where \(m\) is the vortex charge and \(\theta\) is the azimuth angle. Equation (4) represents the amplitude of the wave,which can then be substituted into the propagating paraxial wave, where \(\sqrt{1-K^{2}}\to 1-\frac{K^{2}}{2}\), using the approximation for a paraxial wave. The plane wave can then be rewritten as  
\[ \begin{array}{c}\psi_{m}(R,z)=\iint_{K}\frac{(-i)^{|m|}}{2\pi R}\frac{|m|}{K}e ^{im\theta}e^{iz(1-\frac{K^{2}}{2})}dK\\ =|m|e^{im\theta}\int_{0}^{\infty}J_{m}(KR)e^{iz(1-\frac{K^{2}}{2})}\frac{dK}{K} \\ =|m|e^{i(m\theta+z)}\int_{0}^{\infty}J_{m}(KR)e^{-i(\frac{K^{2}}{2})}\frac{dK}{K }=e^{i(m\theta+z)}P_{m}(R)\end{array}  \tag{5}\]  
where \(P_{m}(R)=\sqrt{\frac{R}{8}}(-i)^{\frac{|m|}{2}}\frac{R}{\sqrt{2}}e^{\frac{i \theta^{2}}{2K}}\times\left[J_{\frac{|m|-1}{2}}\left(\frac{R^{2}}{4z}\right)- iJ_{\frac{|m|+1}{2}}\left(\frac{R^{2}}{4z}\right)\right]\), is the combination of half-integer Bessel functions. Then, the superposition of a pair of vortices leads to the intensity contributions from each vortex \(I_{1}\) and \(I_{2}\), as well as the interference term \(I_{1,2}\) that governs the interference pattern, with \(I=I_{1}+I_{2}+I_{1,2}\). The latter term can be determined by Equation (5), where the phase difference arises from the difference in the Bessel functions. The phase difference will then have a dependence on the charge of the vortex in addition to the length of propagation and initial phase. The interference of two vortices is solved analytically by writing \(\psi_{m_{1}}(R,z)\) and \(\psi_{m_{2}}(R,z)\), as  
\[ \begin{array}{c}\psi_{m_{1}}(R,z)=e^{i(m_{1}\theta+z)}e^{-i(k_{1}Rsin(\theta ))}\sqrt{\frac{R}{8}}(-i)^{\frac{|m_{1}|}{2}}\frac{R}{\sqrt{2}}e^{\frac{i \theta^{2}}{4z}}\\ \times\left[J_{\frac{|m_{1}|-1}{2}}\left(\frac{R^{2}}{4z}\right)-iJ_{\frac{|m |+1}{2}}\left(\frac{R^{2}}{4z}\right)\right]\end{array}  \tag{6}\]  
and  
\[