1.

Introduction

Optical vortices were discovered in the 19th century.¹ However, they have only been extensively studied in recent decades. Owing to the better understanding of lasers and the development of computers, quantum communication, and other similar fields, optical vortices have been studied in a more detailed and systematic manner. Moreover, due to the characteristics of the diffraction field and its significant potential for applications, optical vortices have gradually become the research focus of many institutions and researchers, especially in countries such as the United States, Australia, Canada, and China.^2,3 It can be applied in multiple research fields, such as optical manipulation, optofluidics, cell cultivation and reproduction, and free-space optical communication based on orbital angular momentum (OAM).^4–7

The methods for generating optical vortex beams have been studied for many decades. They include applying computational holography, mode conversion, spatial light modulator, and spiral phase plates (SPPs).^8,9 Among the above-mentioned methods, the generation method involving the use of an SPP has many advantages, such as a relatively simple and small structure that does not require computers or other active devices. Thus, it is widely applicable. SPPs can generate vortex beams with integer topological charges as well as vortex beams with fractional charges.^10,11 The height of the ideal continuous surface SPP is continuous variable with azimuth, which is used for the theoretical analysis. However, it is difficult to fabricate ideal continuous surface SPPs. A multistaircase SPP, which has an increasing staircase as compared to the continuously spiral and increasing curved surface of the ideal SPPs, is commonly used to fabricate ideal continuous surface SPPs. Therefore, the vortex beams generated by multistaircase SPPs are not ideal. The radial intensity profiles of the Fraunhofer diffraction field of multistaircase SPPs differ from those of the ideal continuous surface SPPs. Although researchers have proposed an approximation method for solving the field by using Fourier expansions and obtained the relative intensities of the different terms in Fourier series,^11,12 the field components have not been solved yet. Through theoretical derivation and modeling simulations, this study is aimed at qualitatively and quantitatively exploring the influence of multistaircase SPPs, with different $M\mathrm{s}$ and intrinsic topological charges $l$, on the characteristics of the Fraunhofer diffraction field. As a result, this study provides a more comprehensive analysis of the optical vortices generated by multistaircase SPPs. Furthermore, the proposed model provides a reference for designing multistaircase SPPs with different $M\mathrm{s}$ and different $l\mathrm{s}$ and generating various types of optical vortices.

2.

Theory

2.1.

Fraunhofer Diffraction Analysis for Multistaircase Spiral Phase Plates

Initially, the characteristics of an ideal continuous surface SPP and that of a multistaircase SPP are analyzed, as shown in Figs. 1(a) and 1(b), respectively. Among these characteristics, $\theta$ is the azimuth angle of the SPP, $\lambda$ is the wavelength of the incident beam, $n$ is the refractive index of the SPP, $n_0$ is the refractive index of air, $\mathrm{\Delta }n$ is the difference between the refractive index of the SPP and that of air, $h$ is the height of the SPP, and $M$ is the number of staircases. The corresponding theoretical SPP topological charge is $l=\mathrm{\Delta }nh/\lambda$.

Fig. 1

(a) An ideal continuous surface SPP and (b) a multistaircase SPP with $M=16$.

However, in practical applications, multistaircase SPPs are commonly applied. Therefore, the height exhibits the characteristic of a stepped spiral growth with an increase in the azimuth angle, and the incident beam will be subjected to the corresponding phase modulation.¹² In this scheme, an SPP with $M$ staircases is used to generate a vortex beam, and the transfer function of the SPP is as follows:¹²

Eq. (1)

$$T(r,\theta )=\mathrm{circl}\left(\frac{r}{R}\right)\exp [i2\pi \frac{l}{M}H(\theta )],$$

Eq. (2)

$$H(\theta )=\mathrm{ceil}\left(\frac{M\theta }{2\pi }\right)=m;m=1,2,\dots ,M,\theta \in [\frac{(m-1)·2\pi }{M},\frac{m·2\pi }{M}],$$

where $T(r,\theta )$ is the phase transfer function of the phase plate and $\mathrm{circl}(r/R)$ is the aperture function. The independent variable $r/R$ has a value ranging from zero to infinity. Inside the aperture hole, $r/R$ takes the value 1, whereas outside the aperture hole, it takes the value 0. Here $H(\theta )$ is the step phase function that is obtained by compressing the horizontal axis of the ceiling function $\mathrm{ceil}(M\theta /2\pi )$. It behaves as $\theta$ obtains $M$ discrete function values when taking values in $M$ different intervals. By extending the domain of $\theta$ from $[0,2\pi ]$ to $(-\infty ,+\infty )$, the corresponding $H(\theta )$ can be expanded by the Fourier expansion method^1,2 such that an infinitely expanded phase transfer function can be obtained:

Eq. (3)

$$T(r,\theta )=\sum _{p=-\infty }^{+\infty }{C}_p\exp (ip\theta ),$$

Eq. (4)

$$\begin{gathered} C_p=\frac{1}{2\pi }\mathrm{circl}\left(\frac{r}{R}\right)\sum _{m=1}^M\underset{2\pi (m-1)/M}{\overset{2\pi m/M}{\int }}\exp (i2\pi lm/M)·\exp (-ip\theta )\mathrm{d}\theta \\ =\frac{1-\exp (ip2\pi /M)}{-2\pi ip}\mathrm{circl}\left(\frac{r}{R}\right)\sum _{m=1}^M\exp \left(i2\pi \frac{l-p}{M}m\right), \end{gathered}$$

where $C_p$ is the coefficient of the $p$’th item of the expansion. If ${\varphi }_m=\exp [i2\pi l-(p/M)m]$ in Eq. (4) is not equal to 1, $C_p$ will always be 0. The $C_p$ takes a nonzero value if and only if ${\varphi }_m$ is equal to 1, namely $(l-p)/M\in Z$.

Eq. (5)

$$C_p=\{\begin{array}{ll}\exp \left(\frac{ip\pi }{M}\right)·\sin c\left(\frac{p}{M}\right)\mathrm{circl}\left(\frac{r}{R}\right),& (l-p)/M\in Z,p\in l-MZ\\ 0,& p\in \text{else}\end{array}.$$

Considering the uniform beam and the theory of Fraunhofer diffraction, the Fraunhofer diffraction field on the back focal plane can be expressed as

Eq. (6)

$$E(\rho ,\phi )=\frac{1}{i\lambda f}\sum _{l=-\infty }^{+\infty }\underset{0}{\overset{+\infty }{\int }}\underset{0}{\overset{2\pi }{\int }}{C}_p\exp [-\frac{ik}{f}r\rho \cos (\phi -\theta )]r\mathrm{d}r\mathrm{d}\theta .$$

The nonzero value of the Fraunhofer diffraction field, which satisfies $(l-p)/M\in Z$, is

Eq. (7)

$$E(\rho ,\phi )=\frac{1}{i\lambda f}\sum _{p=-\infty }^{+\infty }\exp \left(\frac{ip\pi }{M}\right)·\sin c\left(\frac{p}{M}\right)\underset{0}{\overset{R}{\int }}\underset{0}{\overset{2\pi }{\int }}\exp (ip\theta )\exp [-\frac{ik}{f}r\rho \cos (\phi -\theta )]r\mathrm{d}r\mathrm{d}\theta .$$

According to Bessel’s equation,¹³

Eq. (8)

$$\exp [-\frac{ik}{f}r\rho \cos (\theta -\phi )]=\sum _{p=-\infty }^{\infty }{(-i)}^p{J}_p\left(\frac{k}{f}r\rho \right)e^{ip(\theta -\phi )}.$$

Therefore, it can be derived as

Eq. (9)

$$E(\rho ,\phi )=\frac{-i}{\lambda f}\sum _{p=-\infty }^{+\infty }\exp \left(\frac{ip\pi }{M}\right)·\sin c\left(\frac{p}{M}\right)\exp (ip\phi )\underset{0}{\overset{R}{\int }}{J}_p\left(\frac{k}{f}r\rho \right)r\mathrm{d}r.$$

2.2.

Solution of the Finite Hypergeometric Series Summation Fraunhofer Diffraction Field

The integral of the Bessel function is represented by the hypergeometric series¹⁴ as

Eq. (10)

$$\underset{0}{\overset{R}{\int }}{J}_p\left(\frac{k}{f}r\rho \right)r\mathrm{d}r=\frac{{R^2\left(\frac{k\rho R}{2f}\right)}^p}{(p+2)p!}{}_{1}F_2[\frac{p+2}{2};\frac{p+4}{2},p+1;-{\left(\frac{kR\rho }{2f}\right)}^2],$$

where $R$ denotes the radius of the SPP, and it is equal to the radius of the diffraction aperture. We consider the Fraunhofer diffraction field of a uniform beam incident on an ideal continuous surface SPP (the ideal condition), and the diffraction is described as¹⁴

Eq. (11)

$$\begin{gathered} E_p^R(\rho ,\theta )=\frac{{(-i)}^{p+1}\exp (ip\theta )}{(p+2)p!}\left(\frac{k{R}^2}{f}\right){\left(\frac{kR\rho }{2f}\right)}^p{{}_{1}F}_2 \\ \times [\frac{p+2}{2},\frac{p+4}{2},p+1;-{\left(\frac{kR\rho }{2f}\right)}^2], \end{gathered}$$

where ${}_{1}F_2$ is a hypergeometric function.

According to the relationship between,¹⁵

Eq. (12)

$$J_{-p}(x)={(-1)}^n{J}_p(x).$$

Considering a multistaircase SPP with $M$ staircases, the Fraunhofer diffraction field of a uniform beam can be described as the summation of each component field, which corresponds to the Fourier expansion

Eq. (13)

$$\begin{gathered} E(\rho ,\phi )=\frac{-i{R}^2}{\lambda f}\sum _{p=\mathrm{l}}^{l-M{Z}_1}\frac{e^{\left[ip\right(\frac{\pi }{M}+\phi \left)\right]}{\left(\frac{k\rho R}{2f}\right)}^p\sin c\left(\frac{p}{M}\right)}{(p+2)p!}{}_{1}F_2[\frac{p+2}{2};\frac{p+4}{2},p+1;-{\left(\frac{kR\rho }{2f}\right)}^2] \\ +\frac{-i{R}^2}{\lambda f}\sum _{p=-\infty }^{l-M{Z}_2}\frac{e^{\left[ip\right(\frac{\pi }{M}+\phi \left)\right]}\sin c\left(\frac{p}{M}\right){(-1)}^{-p}{\left(\frac{k\rho R}{2f}\right)}^{-p}}{(2-p)(-p)!}{}_{1}F_2[\frac{2-p}{2};\frac{4-p}{2},1-p;-{\left(\frac{kR\rho }{2f}\right)}^2], \end{gathered}$$

where $Z_1>0$, and $Z_2=-Z_1$. Here $Z_1$ and $Z_2$ are defined as the number of the nonzero positive and negative terms, respectively. For example, in the case where $Z_1=4$, the sum of the terms with $p=l$, $l\pm M$, $l\pm 2M$, $l\pm 3M$, and $l\pm 4M$ are chosen. Therefore, $Z_1$ and $Z_2$ represent the approximation quality of the Fourier expansion method.

The Fraunhofer diffraction field of the more practical situation of the Laguerre–Gaussian (LG) beam is^16,17

Eq. (14)

$$E(\rho ,\theta ,f)=\frac{i{w}_0}{w_f}{A}_{l,0}{\left(\sqrt{2}\right)}^l\frac{\mathrm{\Gamma }[1+(l+l+p)/2]}{\mathrm{\Gamma }(l+p+1)}\exp [-i(l+p)(\theta -\frac{\pi }{2}\left)\right]\exp (-\frac{{\rho }^2}{w_f^2}){\left(\frac{\rho }{w_f}\right)}^{l+p}\times F[\frac{p}{2},l+p+1;{\left(\frac{\rho }{w_f}\right)}^2],w_f=\lambda f/\pi w_0,$$

where $w_0$ denotes the waist of the beam, and it is equal to the radius of the diffraction aperture. Here $F(·,·;·)$ is the Kummer’s function and $p$ is the radial mode index.¹⁸ Similar to Eq. (11) and considering the case of the LG beam, the Fraunhofer diffraction field can be described using the unilateral Fourier expansion as

Eq. (15)

$$E(\rho ,\theta ,f)=\frac{i{w}_0}{w_f}\sum _{k=0}^{2Z_1}\sin c\left(\frac{l}{M}\right)\exp \left(\frac{i\pi l}{M}\right)A_{l,0}{\left(\sqrt{2}\right)}^l\frac{\mathrm{\Gamma }[1+(l+l+p)/2]}{\mathrm{\Gamma }(l+p+1)}\exp [-i(l+p)(\theta -\frac{\pi }{2}\left)\right]\exp (-\frac{{\rho }^2}{w_f^2}){\left(\frac{\rho }{w_f}\right)}^{l+p}\times F[\frac{p}{2},l+p+1;{\left(\frac{\rho }{w_f}\right)}^2],$$

where $l=1+kM$, $k=0,1,2,\cdots$.

Finally, the intensity profile of the above-mentioned fields can be obtained using

Eq. (16)

$$I(\rho ,f)=I(\rho ,\theta ,f)={|E(\rho ,\theta ,f)|}^2.$$

The intensity profile of the optical vortices generated by multistaircase SPPs is essential for applications such as optical manipulation and free-space optical communication based on OAM. In a previous study, multistaircase SPPs were used to generate vortex phases and provided a reference for analyzing intensity profiles.¹⁹ In this study, the rule of the error due to the multistaircase is valued by the error function

Eq. (17)

$$f_{\text{error}}=1-[\sin c^2(1/x)].$$

Here $x=M/l$. According to Eqs. (3) and (4), the components of the optical vortex generated by the multistaircase SPP can be divided into two parts: the useful component with $p=l$ and the error components with $p\ne \mathrm{l}$. Thus, the error function $f_{\text{error}}$ represents the power ratio of the error components to the whole components (sum of the useful components and the error components). The $f_{\text{error}}$ has been calculated as shown in Fig. 2. As shown, it reduces sharply and monotonically with the increase of $x=M/l$ in the interval $[1,+\infty )$. Practically, $M/l$ is always $>1$. Thus, considering the power on the useful component and the error components, the rule of the error can be quantified.

Fig. 2

The error function for the situation of the multistaircase SPPs with different $M/l$.

3.

Results and Analysis

By setting the wavelength of the incident beam $\lambda$ to 633 nm, the focal length of the focusing lens $f$ to 20 cm, the radius of the diffraction aperture $(R,w_0)$ to 0.5 cm, and the radial mode index $p$ to 0,¹⁸ the intensity profile of the Fraunhofer diffraction fields was simulated according to Eqs. (13)–(16). In the simulations, three key parameters were discussed under the conditions of a uniform beam and an LG beam, as shown in Figs. 3Fig. 4–5. These parameters were the approximation quality $Z_1$ of the Fourier expansion method, the intrinsic topological charge $l$, and the number of the staircases $M$. In addition, the radial intensities $I_0$ of the Fraunhofer diffraction field, generated by the ideal continuous surface SPPs, were simulated for comparison, as shown in Figs. 4 and 5.

Fig. 3

Radial intensity profiles of the Fraunhofer diffraction field for the situations of the uniform beam: (a) with different values of $Z_1=0$, 1, 2, and 3 and the LG beams (b) with different values of $Z_1=0$, 1, 2, 3, and 4.

Fig. 4

Radial intensity profiles of the Fraunhofer diffraction field for the situation of the uniform beam with different $l\mathrm{s}$ and $M\mathrm{s}$: (a) $l=1$, $M=4$, 8, 16, 32, and $I_0$; (b) $l=10$, $M=32$, 64, 128, 256, and $I_0$; (c) $l=20$, $M=64$, 128, 256, and $I_0$; (d) $l=30$, $M=64$, 128, 256, 512, and $I_0$.

Fig. 5

Radial intensity profiles of the Fraunhofer diffraction field for the situation of the LG beam with different $l\mathrm{s}$ and $M\mathrm{s}$: (a) $l=1$, $p=0$, $M=4$, 8, 16, 32, and $I_0$; (b) $l=10$, $p=0$, $M=32$, 64, 128, 256, and $I_0$; (c) $l=20$, $p=0$, $M=32$, 64, 128, 256, and $I_0$; (d) $l=30$, $p=0$, $M=64$, 128, 256, 512, and $I_0$.

Figure 3 presents the radial intensity profile of the Fraunhofer diffraction field, based on Eqs. (13) and (15), under the conditions of the uniform beam with $l=1$, $M=4$, and $Z_1=0$, 1, 2, and 3 and the LG beam with $l=1$, $M=4$, and $Z_1=0$, 1, 2, 3, and 4. In this case, $I\_Z_1=n$ and $n=0$, 1, 2, 3, and 4 indicate the intensity profile of the Fraunhofer diffraction field of the vortex beams generated by the multistaircase SPPs. Here $n$ is the value of $Z_1$. Therefore, $n$ represents the accuracy of the finite hypergeometric series summation in Eqs. (13) and (15).

On comparing the conditions of the uniform beam and the LG beam, it was observed that, for the condition of the uniform beam, the normalized root mean square error (NRMSE) for the intensity profiles at $Z_1=0$, 1, 2, and 3 and that at $Z_1=1$, 2, 3, and 4 is 0.126, 0.029, 0.008, and 8.680e-4, respectively. This is less than that observed for the condition of the LG beam, which is 0.357, 0.107, 0.042, and 0.015. When $Z_1$ was further increased to 5, the intensity profile of the LG beam in the situation of the multistaircase SPP mentioned in Eqs. (13)–(16) was approximately identical to the situation of the ideal continuous surface SPP. In addition, the intensity profile of the uniform beam under the same condition was approximately identical to the situation of the ideal continuous surface SPP when $Z_1=4$. This approximate identity has been proved for the uniform beam and the LG beam conditions. In the former condition, the NRMSE between the intensity profile at $Z_1=4$ and that at $Z_1=3$ is 8.6807e-4, whereas in the latter condition, the NRMSE between the intensity profile at $Z_1=4$ and that at $Z_1=3$ is 0.015. Accordingly, as the value of $Z_1$ reached 4, the approximation of the Fourier expansion was sufficiently precise. However, the increase in the value of $Z_1$ will result in additional complexity for the computations. In the simulation, when $l$ and $M$ were set to greater than 1 and 4, respectively, it was found that sufficient accuracy can be obtained even when $Z_1$ was set to $<4$. Therefore, for a specific value of $l$ and $M$, the value of $Z_1$ should be appropriately set according to specific simulation results.

Figure 4 presents the intensity profile for the situation of the uniform beam. The parameters are set as $l=1$, 10, 20, and 30 and $M=4$, 8, 16, 32, 64, 128, 256, and 512. For example, $l_1{M}_4$ represents the intensity profile with topological charge $l=1$ and number of the staircases $M=4$ in Figs. 4 and 5. For comparison, the intensity $I_0$ corresponding to the condition of the ideal continuous surface SPP was also simulated with the same $l$ as the situation of the multistaircase SPPs.

As shown in Fig. 4, when $l$ increases, the radius of the dark core and the sidelobe of the ideal situation $I_0$ also increase, whereas the intensity of $I_0$ decreases. In the application field of optical manipulation, it indicated a larger manipulation range but a weaker manipulation ability. In the condition of the multistaircase SPPs, there was a more significant sidelobe effect in the corresponding intensity profile, especially in the case of the higher $l$ and lower $M$, e.g., $l=20$ and $M=32$. This severe sidelobe effect results in negative influence on the applications of the optical vortex.

In general, the intensity profile was closer to $I_0$ when the number of staircases $M$ increased. However, the value of $M$ should also be as low as possible to reduce the complexity and difficulty during the manufacturing process. Accordingly, for certain values of $l$, the corresponding $M$ was obtained with little difference as compare to the ideal situation. When $l$ is 1, 10, 20, and 30 and $M$ is 32, 256, 256, and 512, the NRMSEs between the situations $l_1{M}_{32}$, $l_{10}{M}_{256}$, $l_{20}{M}_{256}$, $l_{30}{M}_{512}$, and $I_0$ are 0.0143, 0.0142, 0.0445, and 0.0221, respectively. Based on these values, a fine accordance between the conditions of multistaircase SPPs and that of ideal continuous surface SPPs can be achieved.

Similarly, in Fig. 5, when $l$ is 1,10, 20, and 30 and $M$ is 32, 256, 256, and 512, the NRMSEs between the situations $l_1{M}_{32}$, $l_{10}{M}_{256}$, $l_{20}{M}_{256}$, $l_{30}{M}_{512}$, and $I_0$ are 0.0093, 0.0142, 0.0563, and 0.0318 respectively. Though the NRMSE between the situations $l_{10}{M}_{256}$ and $I_0$ is 0.059, which is low, Fig. 5(b) shows that the sidelobes still exist when $M$ is 128. To achieve a better result, $M$ should be set to 256. Thus, a fine accordance between the conditions of multistaircase SPPs and ideal continuous surface SPP can also be achieved. Even though it differs from that in the uniform beam condition, the sidelobe effect was more apparent in the condition of the LG beam when $M$ was not sufficiently high. Once the sidelobe effect was suppressed, the difference between the conditions of the two types of SPPs can be negligible. Moreover, Figs. 4 and 5 also indicated that, when $l$ was larger than 30, $M$ reaches 512 or even higher. This may cause difficulties during the fabrication.

4.

Conclusion

To summarize the study, a Fraunhofer diffraction analysis based on a finite hypergeometric series summation was proposed for the optical vortices generated by practically applied multistaircase SPPs. In practical applications, the fabrication of multistaircase SPPs should consider accuracy and complexity and ensure a balance between these two parameters. As verified by the values of the NRMSE and the intensity profile, the above-mentioned analysis was capable of determining the appropriate value for the number of staircases $M$ under the condition of different topological charges $l$. The proposed method and its results can provide a reference for the manufacturing of multistaircase SPPs with sufficient accuracy and acceptable complexity.