## Abstract

We propose a novel phase retrieval algorithm in Hankel (or called Fourier-Bessel) transform domains by using Monte-Carlo method. Based on the proposed algorithm, we investigate the generation of Gaussian-like beams, such as hollow Gaussian beam, Bessel-Gaussian beam and Laguerre-Gaussian beam, with double phase filtering operations. The phase distributions of filters are determined by employing the proposed phase retrieval algorithm. The advantage of the method is that the total energy of the beam is conservative. Numerical simulations are shown to demonstrate the validity of the scheme.

© 2008 Optical Society of America

## 1. Introduction

In recent years, the application and the generation of dark hollow beam (DHB), the intensity of which shows the formation of ring or multiple rings, have been became a focused point in the field of optical physics [1, 2, 3, 4, 5, 6], for the reasons that the DHBs can be used in the fields of optical atom trapping, cold atom guidance and particle controlling. Hollow Gaussian beam (HGB)[4] is one of these hollow beams. In our previous work, we have reported a scheme to generate HGBs with spatial filtering on amplitude [5]. Mathematically this method can generate the exact HGBs, however, the operation of amplitude filtering reduces the total energy of output beam. This drawback will eventually limit the practical application of this scheme as optical tweezers. Therefore, it is obviously significant and necessary to investigate the schemes to generate DHBs without energy loss.

Phase-only filtering will be introduced into the design of new optical scheme because this kind of operator can conserve the total energy of beam in propagation process under paraxial approximation. In the procedure of beam transformation, the information of input and output is known, but detail internal relations are unknown, which is an optical inverse problem. As we known, Gerchberg-Saxton (G-S) phase retrieval algorithm had been applied in beam shaping [7, 8, 9], X-ray system [10, 11, 12], the coherent electron imaging of a carbon nanotube [13] and information processing [15, 16, 17, 18, 19, 20, 23, 22, 23, 24]. Yang-Gu algorithm [25] is an improved version of G-S algorithm and was used to the phase retrieval algorithm in fractional Fourier transform domain [26]. However, to our best knowledge, there have been only few investigations on phase retrieval algorithm in Hankel transform (HT) domain. Moreover, the phase retrieval algorithm among multiple planes is rarely reported and required very high precision between two adjacent planes. In our specific issue of generating HGBs, single phase-only filtering operation cannot complete the task, therefore phase retrieval between multiple planes must be involved. In this article we propose a novel phase retrieval algorithm between three planes. Based on this proposed algorithm, we propose a scheme to generate Gaussian-like beams from a fundamental Gaussian beam with double phase-only filtering both in the input plane and the HT domain. The phase filtering can be implemented by phase-type spatial light modulators (SLMs) and their phase distributions can be calculated with the phase retrieval algorithm. A significant advantage of this method is that the energy is held between input and output planes in this optical system. Moreover, an arbitrary beam with central symmetric amplitude and phase distributions can be generated with this system.

The rest of the paper is organized in the following sequence. In Section 2, the proposed phase retrieval algorithm is introduced. In Section 3, some numerical results are given to demonstrate the validity of the algorithm. Concluding remarks are summarized in the final section.

## 2. Phase retrieval algorithm

Figure 1 illustrates an optical set-up to generate HGBs, in which a double phase filtering process is performed in a classical 4 *f* optical system. Two phase-only filters, exp[*i*Δ*ϕ*
_{1}(*x,y*)]

and exp[*i*Δ*ϕ*
_{2}(*u, ν*)], which can be displayed by spatial light modulators SLM_{1} and SLM_{2}, are located at the input plane P_{1} and the Fourier spectrum plane P_{2}, respectively. The incidence Gaussian beam with fundamental mode on the input plane, *E*(*x, y*), can be generally expressed as

where *G*(*z*), *ω*(*z*) and *R*(*z*) are the functions depending on the variable *z. k*=2*π*/λ is wave vector. For simplicity, we only consider an incidence Gaussian beam at its waist with minimum spot size *ω*
_{0} and a constant phase ϕ_{1}. Therefore the amplitude of the incidence beam *A*
_{1}(*x, y*) can be expressed as

where *ω*
_{0} is the waist of the Gaussian beam and *G*
_{10} is a constant. In such cease the desired output beam *E _{3}*(

*x*′,

*y*′) will also be at its waist, with the amplitude

*A*(

_{3}*x*′,

*y*′) and and a constant phase profile

*ϕ*

_{3}(

*x*′,

*y*′) at the plane P

_{3}, which can be expressed as [4, 5]

$${\varphi}_{3}(x\prime ,y\prime )={C}_{0},$$

where *p* is the order of the HGB.

As we have known, the amplitude and phase distributions of Gaussian beam, Hollow-Gaussian beam, Laguerre-Gaussian beam, Bessel-Gaussian beam and etc, are central symmetric about z direction on the transverse planes. Therefore, in the following analysis we express the amplitude and phase distributions of beam under the polar coordinate system. All the functions of amplitude and phase at the three planes are only radial coordinates dependent, denoted by *r*, *ρ*, *r*′, respectively. From this symmetry, we can express the propagation of light wave in the 4 *f* optical system shown as Fig. 1 by use of the Hankel transform (HT) as follows

and

where *ϕ*
_{1}(*r*)=ϕ_{1}+Δ*ϕ*
_{1}(*r*) and ϕ_{2}(*ρ*)=*ϕ*
_{2}(*ρ*)+Δϕ_{2}(*ρ*), and ϕ_{2}(*ρ*) is the phase distribution at the plane P_{2} just after the phase modulation of SLM_{2}. *ℋ* denotes HT, which of a function *f*(*r*) is expressed as

where *J*
_{0} is the zeroth order Bessel function of the first kind. From the known *A _{3}*(r′) and

*ϕ*

_{3}(

*r*′) in Eq. (3), we can compute

*A*

_{2}(

*ρ*) and ϕ

_{2}(

*ρ*) by using HT as following

$${\phi}_{2}\left(\rho \right)=\mathrm{arg}\left\{\frac{H\left\{{A}_{3}\left(r\prime \right)\mathrm{exp}\left[i{\varphi}_{3}\left(r\prime \right)\right]\right\}}{{A}_{2}\left(\rho \right)}\right\},$$

where the function ‘arg’ denotes the operation of taking phase angle of a complex number. In Eq. (7), HT is used for the inverse HT because the inverse transform of HT is itself yet.

In general case, the following equation

is not satisfied directly, the phase profiles *ϕ*
_{2}(*ρ*) and *ϕ*
_{1}(*r*) are in the black box with known external conditions in the optical system. It is a typical optical inverse problem. Getting analytic solution of Eq. (6) is very difficult from known amplitudes *A*
_{1}(*r*) and *A _{2}*(

*ρ*) only. In fact, however, we can obtain the numerical solution by employing a phase retrieval algorithm. As we know, Yang-Gu phase retrieval algorithm [25] is best one in all iterative phase retrieval algorithms and it can be used in non-unitary transform domain. However, we have tested G-S algorithm and Yang-Gu algorithm for solving Eq. (6) and found out that they cannot deal with this problem because both them do not achieve convergence behavior in this case. Moreover, the error of the amplitude

*A*

_{2}(

*ρ*) will affect the precision of

*A*(

_{3}*r*′) in the transform of next step. Thereby, a convergence phase retrieval with high precision is required.

Here we introduce the idea of Monte-Carlo model to design an iterative phase retrieval algorithm, which is shown in Fig. 2. In the initialization of the process, *ϕ*
_{1, 0}(*r*) is chosen randomly from the interval [-*π, π*]. The amplitudes *A*
_{1}(*r*) and *A*
_{2}(*ρ*) are known and computed possibly, respectively. Other parameters will be pre-defined. The particular iterative process will be continued according to following several steps: (1) Multiple values of the phase *ϕ*
_{1}(*r*) and *ϕ*
_{1, k}(*k*=1, 2, …, *M*), which can introduce parallel computation mechanism into this algorithm, will be generated by adding a random disturbing value *ϕ*
* _{r, k}*(

*r*) located in the interval [-0.5, 0.5]

*c*

*. The corresponding complex amplitude*

_{0}ε_{t}*A*

_{1, k}(

*r*)exp[

*iϕ′*(

_{1, k}*r*)] can be obtained. (2) The HT of

*A*(

_{1, k}*r*)exp[

*iϕ′*(

_{1, k}*r*)] is calculated and get complex amplitudes

*A*′

_{2, k}(ρ)exp[

*i*

*ϕ′*(

_{2, k}*ρ*)]. To express the error between

*A*′

_{2, k}(

*ρ*) and

*A*

_{2}(

*ρ*), a function of mean square error (MSE) is defined as

where *L* denotes the total number of sample point at *ρ*-direction. (3) Here a minimum value of these MSE, mse(*A′ _{2,k},A_{2}*)min will be found out to be regarded as an optimized result and record corresponding value of

*ϕ′*(

_{1, k}*r*), which will substitute

*ϕ*

_{1, 0}when returning to the step (1). If the minimum value satisfies the condition mse(

*A′*)

_{2, k}, A_{2}_{min}>

*εt*, the procedure will return to the step (1). The condition

*n*-

*n*

_{0}≤Δ

*n*is used for checking and controlling the magnitude of random vectors

*ϕ*(

_{r, k}*r*). When

*n*-

*n*0>Δ

*n*, the convergence closed to correct value of

*ϕ*

_{1}(

*r*) is considered as stop and

*ε*is decreased to continue the convergence. If mse(

_{t}*A′*)

_{2, k}, A_{2}_{min}<

*εt*and

*εt*>

*ε0*, the parameters

*εt*and

*ϕ*

_{1, 0}(

*r*) are renewed with better values mse(

*A′*)

_{2, k}, A_{2}_{min}and

*ϕ*(

_{1, k}*r*), respectively, and label parameters

*n*

_{0}and

*n*are also updated. If

*ε*<

_{t}*ε*, we close this recurrence process and consider that pre-defined precision has been achieved. According to the characteristic of this algorithm, it can be refereed as an optimized Monte-Carlo phase retrieval algorithm. The conventional phase retrieval is to explore correct phase information. In this paper, however, the purpose of phase retrieval is to obtain more approximative amplitude

_{0}*A*

_{2}(

*ρ*).

After the acquirement of the corresponding recovered phases *ϕ ^{r}_{1}*(r) and

*ϕ*(

^{r}_{2}*ρ*), the modulation phase profiles on SLM

_{1}and SLM

_{2}can be directly calculated as

Optically these phase distributions can be displayed on the phase-type SLMs. Alternatively, according to Eq. 8 we can fabricate two phase masks like Fresnel zone plate by use of femtosecond laser technology to replace two SLMs.

## 3. Numerical simulation

Figure 3 gives an example of numerical simulation with this algorithm to obtain the secondorder hollow Gaussian beam. The values of all parameters are given as following context: *c*
_{1}=0.99, *c*
_{0}=0.2, Δ*n*=4, *ε*
_{0}=1×10^{-4}, the initial value of *ε _{t}* is equal to 1,

*M*=100, the sample point number and the interval of

*r*(

*ρ*and

*r*′) are equal to 256 and [0, 4]×

*ω*

_{0}, respectively. A group of random vectors is adopted as the initial value of phase

*ϕ*

_{1}(

*r*). Under the condition limited with above mentioned parameters, the total iterative number is 68673 times (This number is also marked with underline and bold in Table 1). The analytical values of amplitudes

*A*

_{2}(

*ρ*) and

*A*(

_{3}*r*′) are depicted for the comparison with their recovered values and are marked by dashed and solid lines, respectively. The recovered phases

*ϕ*′

_{1}(

*r*) and

*ϕ*′

_{2}(

*ρ*) are given in the second row sub-figures in Fig. 3. The precision of the recovered phases

*ϕ*′

_{1}(

*r*) and

*ϕ*′

_{2}(

*ρ*) and hence the corresponding amplitudes

*A*

_{2}(

*ρ*) and

*A*

_{3}(

*r*′) can be increased by choosing a smaller value of

*ε*

_{0}. Figure 4 gives the comparison of the transverse profiles of the HGBs calculated with analytical expression and generated by the phase retrieval algorithm. The results demonstrate well the effectiveness of the proposed algorithm.

It is worth to point out that choosing a proper coefficient *c _{0} ε_{t}* of

*ϕ*(

_{r, k}*r*) can increase the convergence speed of this phase retrieval algorithm. In our algorithm the value of

*ε*is given by the minimum value of the MSE and it will decrease with the convergence accordingly. Such a variable coefficient

_{t}*ε*can enable a high precision phase retrieval results with a reasonable convergence speed. The coefficients

_{t}*c*

_{0}and

*c*

_{1}can also affect the convergence speed. We have also numerically calculated the influence of the coefficients on the iterative steps. The results are shown in the Table 1.

With the same initial values used in Fig. 3, three simulations with G-S, Yang-Gu and the proposed algorithms are tested. Figure 5 shows a comparison of convergence speed between G-S, Yang-Gu and the proposed algorithms, in which the values of MSE are calculated from the analytical and recovered values of amplitude *A*
_{2}(*ρ*). From Fig. 5, one can see that the G-S and Yang-Gu algorithm can not converge to an approximative value of *A*
_{2}(*ρ*), while the proposed algorithm can converge quickly to a desired value. That demonstrates the proposed algorithm is only effective in above mentioned three methods in the issue of generating HGBs. We have also demonstrated that the proposed algorithm is effective to generate other Gaussian-like beams, such as hollow Gaussian beam, Bessel-Gaussian beam and Laguerre-Gaussian beam, by changing the distributions of *A*
_{3}(*r*′) and *ϕ*
_{3}(*r*′).

## 4. Conclusion

We have proposed a novel phase retrieval algorithm based on Monte-Carlo structure in Hankel transform domains. This algorithm can be also expanded into other transform domains to deal with general phase retrieval problem from two amplitude (or intensity) images and it is more effective than the G-S and Yang-Gu algorithm in the behaviors of the convergence. By use of this phase retrieval algorithm, we have designed an optical system with phase-only filtering to generate hollow Gaussian beam or other central symmetric beams. Numerical simulations have demonstrated its validity.

## Acknowledgments

This work was supported by the National Natural Science Foundation of China under Grant No. 10674038 and 10604042, National Basic Research Program of China under Grant 2006CB302901, China Postdoctoral Science Foundation (20080430913), and development program for outstanding young teachers in Harbin Institute of Technology, HITQNJS. 2008. 027. Zhengjun Liu wishes to acknowledge Prof. Amilia Torre, who is working in the ENEA Center, Roman, Italy, for her help in obtaining the references on Hankel transform and some useful discussions.

## References and links

**1. **T. Kuga, Y. Torii, N. Shiokawa, T. Hirano, Y. Shimizu, and H. Sasada, “Novel optical trap of atoms with a doughnut beam,” Phys. Rev. Lett. **78**, 4713–4716 (1997). [CrossRef]

**2. **H. Ito, T. Nakata, K. Sakaki, M. Ohtsu, K. I. Lee, and W. Jhe, “Laser spectroscopy of atoms guided by evanescent waves in micron-sized hollow optical fibers,” Phys. Rev. Lett. **76**, 4500 (1996). [CrossRef]

**3. **Y. K. Wu, J. Li, and J. Wu, “Anomalous hollow electron beams in a storage ring,” Phys. Rev. Lett. **94**, 134802 (2005). [CrossRef]

**4. **Y. Cai, X. Lu, and Q. Lin, “Hollow Gaussian beams and their propagation properties,” Opt. Lett. **28**, 1084 (2003). [CrossRef]

**5. **Z. Liu, H. Zhao, J. Liu, J. Lin, M. A. Ahmad, and S. Liu, “Generation of hollow Gaussian beams by spatial filtering,” Opt. Lett. **32**, 2076 (2007). [CrossRef]

**6. **Y. Zhao, J. S. Edgar, G. D. M. Jeffries, D. McGloin, and D. T. Chiu, “Spin-to-orbital angular momentum conversion in a strongly focused optical beam,” Phys. Rev. Lett. **99**, 073901 (2007). [CrossRef]

**7. **R. W. Gerchberg and W. O. Saxton, “A practical algorithm for the determination of phase from image and diffraction plane pictures,” Optik **35**, 237–246 (1972).

**8. **Z. Zalevsky, D. Mendlovic, and R. G. Dorsch, “Gerchberg-Saxton algorithm applied in the fractional Fourieror the Fresnel domain,” Opt. Lett. **21**, 842–844 (1996). [CrossRef]

**9. **D. F. McAlister, M. Beck, L. Clarke, A. Mayer, and M. G. Raymer, “Optical phase retrieval by phase-space tomography and fractional-order Fourier transforms,” Opt. Lett. **20**, 1181–1183 (1995). [CrossRef]

**10. **K. A. Nugent, T. E. Gureyev, D. F. Cookson, D. Paganin, and Z. Barena, “Quantitative Phase Imaging Using Hard X Rays,” Phys. Rev. Lett. **77**, 2961–2964 (1996). [CrossRef]

**11. **J. Miao, T. Ishikawa, B. Johnson, E. H. Anderson, B. Lai, and K. O. Hodgson, “High resolution 3D X-ray diffraction microscopy,” Phys. Rev. Lett. **89**, 088303 (2002). [CrossRef]

**12. **F. Meng, H. Liu, and X. Wu, “An iterative phase retrieval algorithm for in-line x-ray phase imaging,” Opt. Express **15**, 8383–8390 (2007), http://www.opticsinfobase.org/abstract.cfm?URI=oe-15-13-8383. [CrossRef]

**13. **J. M. Zuo, I. Vartanyants, M. Gao, R. Zhang, and L. A. Nagahara, “Atomic resolution imaging of a carbon nanotube from diffraction intensities,” Science **300**, 1419 (2003). [CrossRef]

**14. **R. Borghi, F. Gori, and M. Santarsiero, “Phase and amplitude retrieval in ghost diffraction from field-correlation measurements,” Phys. Rev. Lett. **96**, 183901 (2006). [CrossRef]

**15. **A. Anand, G. Pedrini, W. Osten, and P. Almoro, “Wavefront sensing with random amplitude mask and phase retrieval,” Opt. Lett. **32**, 1584–1586 (2007). [CrossRef]

**16. **P. Almoro, G. Pedrini, and W. Osten, “Wavefront sensing with random amplitude mask and phase retrieval,” Opt. Lett. **32**, 733–735 (2007). [CrossRef]

**17. **A. K. Ellerbee and J. A. Izatt, “Phase retrieval in low-coherence interferometric microscopy,” Opt. Lett. **32**, 388–390 (2007). [CrossRef]

**18. **J. Zhong and J. Weng, “Phase retrieval of optical fringe patterns from the ridge of a wavelet transform,” Opt. Lett. **30**, 2560–2562 (2005). [CrossRef]

**19. **J. S. Wu, U. Weierstall, J. C. H. Spence, and C. T. Koch, “Iterative phase retrieval without support,” Opt. Lett. **29**, 2737 (2004). [CrossRef]

**20. **N. Nakajima, “Phase and amplitude retrieval in ghost diffraction from field-correlation measurements,” Phys. Rev. Lett. **98**, 223901 (2007). [CrossRef]

**21. **F. Zhang, G. Pedrini, and W. Osten, “Phase retrieval of arbitrary complex-valued fields through aperture-plane modulation,” Phys. Rev. A **75**, 043805 (2007). [CrossRef]

**22. **T. Latychevskaia and H.-W. Fink, “Solution to the twin image problem in holography,” Phys. Rev. Lett. **98**, 233901 (2007). [CrossRef]

**23. **Y. Zhang, G. Pedrini, W. Osten, and H. Tiziani, “Whole optical wave field reconstruction from double or multi in-line holograms by phase retrieval algorithm,” Opt. Express **11**, 3234–3241 (2003), http://www.opticsinfobase.org/abstract.cfm?URI=oe-11-24-3234. [CrossRef]

**24. **Z. Liu and S. Liu, “Double image encryption based on iterative fractional Fourier transform,” Opt. Commun. **272**, 324–329 (2007). [CrossRef]

**25. **G.-Z. Yang, B.-Z. Dong, B.-Y. Gu, J.-Y. Zhuang, and O. K. Ersoy, “Gerchberg-Saxton and Yang-Gu algorithms for phase retrieval in a nonunitary transform system: a comparison,” Appl. Opt. **33**, 209 (1994). [CrossRef]

**26. **B. Dong, Y. Zhang, B. Gu, and G. Yang, “Numerical investigation of phase retrieval in a fractional Fourier transform,” J. Opt. Soc. Am. A **14**, 2709–2714 (1997). [CrossRef]