## 2 Design Method of the Off-Axis Pseudo Bessel Beam and Multiple Pseudo Bessel Beams
Based on these settings, the required compensation phase distribution on the metasurface can be calculated according to the condition of the transverse component of the wave vector and geometric optics. This method can satisfy the condition of forming a standing wave on the plane that is perpendicular to the beam propagation direction, thereby ensuring that all energy radiates along the propagation axis. To get the compensation phase distribution of the reflective metasurface, which can generate the energy-focused, off-axis pseudo Bessel beam, the phase of the E-field at each point on the metasurface can be expressed as  
\[\Phi_{1}(x_{ij},y_{ij})=-\frac{2\pi}{\lambda}\hat{u}_{k}\cdot\overset{\text{ \raisebox{-0.86 pt}{$\rightarrow$}}}{O}\hat{A} \tag{1}\]  
Figure 1: Schematic diagram of **(a)** pseudo Bessel beams generated by a reflective metasurface; **(b)** the proposed method to generate the off-axis pseudo Bessel beam. where \(\hat{a}_{k}\) is the unit vector along to the inward conical wave of each point, as shown in Figure 1b, and the vector corresponding to each point on the metasurface can be expressed as  
\[\overset{\rightarrow}{OA}=\sqrt{x_{ij}{}^{2}+y_{ij}{}^{2}}(\cos\varphi_{ij},\, \sin\varphi_{ij},0) \tag{2}\]  
where \(\varphi_{ij}=\arctan(y_{ij}/x_{ij})\). The vector corresponding to the intersection of the straight line in the direction of any vector \(\hat{a}_{k}\) and the propagation axis is  
\[\overset{\rightarrow}{OB}=\frac{\sin(180^{\circ}-\beta-\alpha)}{\sin\alpha} \sqrt{x_{ij}{}^{2}+y_{ij}{}^{2}}\hat{a}_{k} \tag{3}\]  
where \(\beta\) is the included angle between the vectors \(\overset{\rightarrow}{OA}\) and \(\overset{\rightarrow}{OB}\), and it can be obtained by the vector operation  
\[\cos\beta=\sin\theta_{k}\cos\varphi_{k}\cos\varphi_{ij}+\sin\theta_{k}\sin \varphi_{k}\sin\varphi_{ij} \tag{4}\]  
Finally, the unit vector \(\hat{a}_{k}\) can be expressed as  
\[\hat{a}_{k}=\frac{(\overset{\rightarrow}{OB}-\overset{\rightarrow}{OA})}{ \left|\overset{\rightarrow}{OB}-\overset{\rightarrow}{OA}\right|} \tag{5}\]  
Taking Equation (5) into Equation (1), the phase \(\Phi_{1}(x_{ij},y_{ij})\) at each point of the metasurface can be calculated for the off-axis pseudo Bessel beam, which is also the compensation phase of each unit cell to transform the plane wave to the inward conical wave. However, the feed horn we actually use can only radiate the spherical wave. Therefore, we also need to calculate the compensation phase on the metasurface to convert the spherical wave radiated from the feed horn into the plane wave. The required compensation phase is  
\[\Phi_{2}(x_{ij},y_{ij})=\frac{2\pi}{\lambda}(\sqrt{x_{ij}{}^{2}+y_{ij}{}^{2}+r _{f}{}^{2}}-r_{f}) \tag{6}\]  
where \(r_{f}\) is the distance from the feed horn to the metasurface. Then, the final compensation phase of any off-axis pseudo Bessel beam generated by the feed horn and the metasurface is  
\[\Phi_{k}(x_{ij},y_{ij})=\Phi_{1}(x_{ij},y_{ij})+\Phi_{2}(x_{ij},y_{ij}) \tag{7}\]  
In this way, we can get the compensation phase of each unit cell on the metasurface to generate the \(k\)-th off-axis pseudo Bessel beam, and the compensation phase to generate multiple pseudo Bessel beams can be gotten by using the electromagnetic field superposition principle. The phase compensation of each unit cell whose coordinate is \((x_{ij},y_{ij},0)\) can be calculated by the following equation [34]:  
\[\Phi_{total}(x_{ij},y_{ij})=\arg\{\sum_{k=1}^{n}\exp(j\Phi_{k}(x_{ij},y_{ij}))\} \tag{8}\]  
where \(j\) is the imaginary unit, \(\Phi_{k}\) is the compensation phase that needed to be adjusted on the \(k\)-th beam, and \(n\) is the total number of the multiple pseudo Bessel beams.