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
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) is one of these hollow beams. In our previous work, we have reported a scheme to generate HGBs with spatial filtering on amplitude . 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  and information processing [15, 16, 17, 18, 19, 20, 23, 22, 23, 24]. Yang-Gu algorithm  is an improved version of G-S algorithm and was used to the phase retrieval algorithm in fractional Fourier transform domain . 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 SLM1 and SLM2, are located at the input plane P1 and the Fourier spectrum plane P2, 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 E3(x′, y′) will also be at its waist, with the amplitude A3(x′, y′) and and a constant phase profile ϕ 3(x′, y′) at the plane P3, which can be expressed as [4, 5]
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
where ϕ 1(r)=ϕ1+Δϕ 1(r) and ϕ2(ρ)=ϕ 2(ρ)+Δϕ2(ρ), and ϕ2(ρ) is the phase distribution at the plane P2 just after the phase modulation of SLM2. ℋ 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 A3(r′) and ϕ 3(r′) in Eq. (3), we can compute A 2(ρ) and ϕ2(ρ) by using HT as following
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 A2(ρ) only. In fact, however, we can obtain the numerical solution by employing a phase retrieval algorithm. As we know, Yang-Gu phase retrieval algorithm  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 A3(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 0εt. The corresponding complex amplitude A 1, k(r)exp[iϕ′1, k(r)] can be obtained. (2) The HT of A1, 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,A2)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, A2)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-n0>Δn, the convergence closed to correct value of ϕ 1(r) is considered as stop and εt is decreased to continue the convergence. If mse(A′2, k, A2)min<εt and εt>ε0, the parameters εt and ϕ 1, 0(r) are renewed with better values mse(A′2, k, A2)min and ϕ1, k(r), respectively, and label parameters n 0 and n are also updated. If εt<ε0, 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 A 2(ρ).
After the acquirement of the corresponding recovered phases ϕr1(r) and ϕr2(ρ), the modulation phase profiles on SLM1 and SLM2 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 A3(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 c0 εt of ϕr, k(r) can increase the convergence speed of this phase retrieval algorithm. In our algorithm the value of ε t 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 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′).
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.
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] [PubMed]
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] [PubMed]
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).
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] [PubMed]
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] [PubMed]
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]
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] [PubMed]
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] [PubMed]
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]