## 2 Theoretical Analysis  
Based on Fresnel-Huygens theory and considering a mask with \(Q\) apertures arranged on a ring that is illuminated with collimated laser beams, the optical field near the focal region \(z=f(1+\Delta)\) with \(\Delta\ll 1\) can be expressed as [30]:  
\[\Phi_{Q}(\rho,\phi;\theta_{n})=e^{ik(f+\rho^{2}/2\Delta)}\sqrt{\frac{2}{\pi \omega^{2}}}e^{-\rho^{2}/\omega^{2}}\Psi_{Q}(\rho,\phi;\theta_{n}) \tag{1}\]  
with  
\[\Psi_{Q}(\rho,\phi;\theta_{n})=\sqrt{\frac{1}{Q}}\sum_{s=0}^{Q-1}e^{i\frac{ik \phi}{2f}\rho\cos\left(\phi-\frac{2\pi s}{Q}\right)}e^{is\theta_{n}} \tag{2}\]  
where (\(\rho\),\(\phi\)) is the polar coordinate of (\(x\),\(y\)), \(f\) is the focal length of lens, \(a\) is the radius of the ring, \(k=2\pi/\lambda\), \(\lambda\) is the wavelength of laser, \(\omega=\lambda f/(\pi\omega_{0})\), \(\omega_{0}\) is the radius of aperture, and \(\theta_{n}=2\pi n/Q\) is the fixed phase difference between adjacent incident beams. Here, \(Q\) and \(n\) are both integer numbers. For a phase difference between adjacent beams greater than \(\pi\), the structure of the optical field will be equal. Therefore, the phase parameter \(n\) is restricted to a range of \(0\leq n\leq Q/2\) for a given \(Q\). The wave function \(\Psi_{Q}(\rho,\phi;\theta_{n})\) represents various two-dimensional crystalline or quasicrystalline structures with different topologic charges. The wave patterns \(\left|\Phi_{Q}(\rho,\phi;\theta_{n})\right|^{2}\) near the focal region can be confirmed to be nearly propagation invariant over a specific distance range. At the limit of \(Q\rightarrow\infty\), a finite sum can be a certain kind of approximation of an integral with a new variable \(\alpha=2\pi s/Q\) and \(d\alpha=2\pi/Q\). Consequently, the wave function \(\Psi_{N}(\rho,\phi;\theta_{n})\) can be derived as:  
\[\Psi_{Q}(\rho,\phi;\theta_{n})=J_{n}\left(\frac{ka}{2f}\rho\right)=\sqrt{\frac {1}{2\pi}}\int_{0}^{2\pi}e^{i\frac{ka}{2f}\rho\cos\left(\phi-\alpha\right)}e^ {in\alpha}\;d\alpha \tag{3}\]  
The result indicates that high-order Bessel Gaussian beams can be achieved by interference of multiple Gaussian beams. From Equation (1), the wave patterns of \(\Phi_{N}(\rho,\phi;\theta_{n})\) can be numerically calculated with \(\omega_{0}=0.13\) mm, a = 3 mm, and \(f=1000\) mm. Figure 1 shows the calculated intensity distributions of \(\left|\Phi_{Q}(\rho,\phi;\theta_{n})\right|^{2}\) for \(Q=6\), \(8\), \(12\), \(21\), \(30\), and \(50\), with different relative phase \(\theta_{n}\) from \(n=0\) to \(n=8\). For \(n=0\), the calculated patterns with equal phase (\(n=0\)) are shown in the first column. The patterns exhibit crystalline, quasicrystalline, and kaleidoscopic structures that are consistent with the previous results [28]. For the case of \(Q=50\), the central part of the intensity distribution has some resemblance to a Bessel beam of zeroth order. This feature comes from the limiting behavior of \(Q\rightarrow\infty\). On the other hand, for the case of \(n\
eq 0\), the intensity distributions show that a variety of rich quasicrystalline wave patterns can be obtained by manipulating the relative phase \(\theta_{n}\). It can be seen that the wave patterns are repeated as \(n>Q/2\), which is consistent with the previous prediction. The minimal intensity can be discovered in the center of the wave patterns and the region increases with the increasing value of \(n\) for all cases of the \(Q\)-multiple beams. Moreover, the central area of the wave patterns resembles a high-order Bessel beam at larger values of \(Q\). While, increasing the value of \(n\), a similar behavior can be found in terms of a Bessel Gaussian beam with the order \(n\). The numerical results indicate that the relative phase \(\theta_{n}\) can benefit the generation of quasicrystalline beams. In the following, we experimentally constructed an optical system to demonstrate the theoretical predictions.