Heterodyne 3D ghost imaging

Optics Communications 368 (2016) 1–6

Contents lists available at ScienceDirect

Optics Communications journal homepage: www.elsevier.com/locate/optcom

Invited Paper

Heterodyne 3D ghost imaging Xu Yang, Yong Zhang, Chenghua Yang, Lu Xu, Qiang Wang, Yuan Zhao n Department of Physics, Harbin Institute of Technology, Harbin 150001, China

art ic l e i nf o

a b s t r a c t

Article history: Received 1 September 2015 Received in revised form 14 January 2016 Accepted 18 January 2016

Conventional three dimensional (3D) ghost imaging measures range of target based on pulse fight time measurement method. Due to the limit of data acquisition system sampling rate, range resolution of the conventional 3D ghost imaging is usually low. In order to take off the effect of sampling rate to range resolution of 3D ghost imaging, a heterodyne 3D ghost imaging (HGI) system is presented in this study. The source of HGI is a continuous wave laser instead of pulse laser. Temporal correlation and spatial correlation of light are both utilized to obtain the range image of target. Through theory analysis and numerical simulations, it is demonstrated that HGI can obtain high range resolution image with low sampling rate. & 2016 Elsevier B.V. All rights reserved.

Keywords: Heterodyne 3D ghost imaging Range resolution Low sampling rate

1. Introduction Conventional 2D ghost imaging can obtain intensity images through intensity correlation between light in the test path and reference path. The imaging method uses spatial correlation of light to present intensity images of target [1–9]. Due to ghost imaging with high lateral resolution [10] and resistant to atmosphere turbulence [11], researchers apply ghost imaging to many different fields, such as fluorescent imaging [12], polarization imaging [13], optical security [14–16], and computational imaging [17]. Han's research group proposes a 3D ghost imaging lidar system to measure range information of targets [18,19]. In their system, a single pixel detector with time-resolution is used to record intensity of reflected light from targets in different time. Correlating the intensity signal of different echo time and speckle distribution, range image of target can be obtained through 3D ghost imaging system. However, a high-speed data acquisition system should be used to ensure a high range resolution result in conventional 3D ghost imaging. Sampling rate of data acquisition system is usually not high enough to satisfy the needs of range measurement. In this study, a heterodyne 3D ghost imaging (HGI) is proposed to measure range image of target. This system utilizes a continuous wave laser as source. Light intensity is modified by an electro-optical crystal. After intensity modification, a spatial light modulator (SLM) is used to modify speckle distribution. A bucket detector receives echo light from targets. Target range information can be n

Corresponding author. E-mail address: [email protected] (Y. Zhao).

http://dx.doi.org/10.1016/j.optcom.2016.01.049 0030-4018/& 2016 Elsevier B.V. All rights reserved.

obtained by temporal correlation calculation of recorded intensity and modified signals on the electro-optical crystal. Through spatial correlation of speckle distribution and echo light intensity from targets, target intensity information can be calculated. Utilizing this method, a high range resolution result can be obtained under low sampling rate condition. This study provides a new 3D ghost imaging method, which promotes development of ghost imaging in laser radar field.

2. System description and theory analysis HGI system is shown in Fig. 1. As Fig. 1 shows, an electro-optical crystal modulates the frequency of laser. Following, a SLM is used to modulate the spatial distribution of light. An optical transmitting system sends the modulated light out, and an optical receiving system is used to receive light reflected from target. In the system, planar objects are chosen as targets. Reflected light is detected by a bucket detector. Detector signal and modulation signal are processed and sampled by a post-processing circuit. Circuit output is processed by a Fast Fourier Transform (FFT) and reconstructed target range image in the computer. An electro-optical crystal is firstly utilized to modify the intensity of laser, and modulating signal A from the electro-optical crystal is:

⎡ ⎛ 1 ⎞⎤ A = sin ⎢ 2π ⎜ ω 0 + kt ⎟ ⎥ ⎣ ⎝ 2 ⎠⎦

(1)

where ω0 is original frequency of the modulating signal, k is frequency change ratio, and t is modulation time. Supposing I0 is

2

X. Yang et al. / Optics Communications 368 (2016) 1–6

Fig. 1. System schematic of HGI system.

laser intensity, light intensity, after modulated by electro-optical crystal, can be expressed as:

⎡ ⎛ 1 ⎞ ⎤ I = I0 sin ⎢ 2π ⎜ ω 0 + kt ⎟ t ⎥ ⎣ ⎝ 2 ⎠ ⎦

D away from our system. Hence, time delay t′ caused by the target 2D , c

where c is the speed of light. Reflectivity distribution of target is expressed as T (⇀ r ). The echo light intensity received by bucket detector is: T

Ir(i) (t ) =

⎛⎜ ⇀⎞⎟

⎛⎜ ⇀⎞⎟

⎡ ⎢

⎛

∑ T ⎝ r ⎠ I0 I s(i) ⎝ r ⎠ sin ⎢ 2π⋅ ⎜ ω0 + ⇀ r

⎢⎣

⎛ ⎛ ⎞⎞⎤ ⎛ 2D ⎞ ⎟ t + ϕ2 ⎟ ⎟ ⎥ − cos ⎜ 2π ⎜ kt 2 + ⎜ 2ω 0 − t0 − ⎝ ⎠ ⎠ ⎥⎦ c ⎠ ⎝ ⎝

(2)

Next, the light passes through a SLM, which modulates speckle distribution of light. Modulating mode is selected as binary distribution. In ith measurement, the distribution of speckle is Is(i) (⇀ r ), where ⇀ r is horizontal coordinate. Target is located at the distance distance is t′ =

⎛ ⎛ ⎛ 2D ⎞⎞ ⎞ χ⎡ (i ) I˜ (t ) = ⎢ cos ⎜ 2π ⎜ k ⎜ − t0 ⎟ t + ϕ1⎟ ⎟ ⎠ ⎠⎠ 2 ⎢⎣ ⎝ ⎝ ⎝ c

⎝

⎤ 1 ⎛ 2D ⎞ ⎞ ⎛ 2D ⎞ ⎥ ⎟ ⎥. ⎟⎟⎜ t − k⎜t − 2 ⎝ c ⎠⎠⎝ c ⎠⎥ ⎦

ϕ1 =

2 1 2 1 ⎛ 2D ⎞ 2D ⎟ + ω0 kt0 − ω 0 t02 − k ⎜ . 2 2 ⎝ c ⎠ c

ϕ2 =

2 1 2 1 ⎛ 2D ⎞ 2D kt0 − ω 0 t0 2 + k ⎜ ⎟ − ω0 − 2 2 ⎝ c ⎠ c

(

(i ) I I˜ (t )

⎡ ⎛ (i ) 1⎛ 2D ⎞ ⎞ ⎛ 2D ⎞ ⎤ I˜r (t ) = αI0 sin ⎢ 2π ⎜ ω0 + ⎜ t − ⎟ ⎟⋅ ⎜ t − ⎟⎥ c ⎠⎠ ⎝ c ⎠⎦ 2⎝ ⎣ ⎝

T

r ) I s(i) ( ⇀ r ). ∑ T (⇀ ⇀ r

(4)

As Eq. (4) shown, for ith measurement, current intensity varies with time t , and the current signal contains distance information of target. Target distance D can be calculated through a heterodyne (i ) process between current signal I˜ (t ) and modulation signal A . r

After the heterodyne process, the signal can be expressed as:

⎡ ⎛ ⎤ ⎞ 1 (i ) I˜ (t ) = χ sin ⎢ 2π ⎜ ω 0 + k ( t − t0 ) ⎟ ( t − t0 ) ⎥ ⎠ ⎣ ⎝ ⎦ 2 , ⎡ ⎛ 1 ⎛ 2D ⎞ ⎞ ⎛ 2D ⎞ ⎤ ⎟⎥ ⎟⎟⎜ t − ⋅ sin ⎢ 2π ⎜ ω 0 + k ⎜ t − c ⎠⎠⎝ c ⎠⎦ 2 ⎝ ⎣ ⎝

(5)

T where χ = αI0 ∑⇀r T (⇀ r ) Is(i) (⇀ r ) is considered as a constant independent of time, and t0 is time delay of modulating signal. Eq. (5) can be converted into:

0

(6)

(7)

−

2D c

− t0

2

). (8)

k

In Eq. (6), the second item of right side is a high frequency item, which can be filtered out by a low pass filter. The filtered signal can be sampled by data acquisition system and then dealt with Fourier transform. The signal after Fourier transform can be shown as:

(3)

Supposing photoelectric conversion factor of the bucket detector is α , current signal is shown as:

(ω

,

)

⎞⎤ ⎛ 2D χ ⎡ δ⎢ω − k⎜ − t0 ⎟ ⎥. ⎝ c ⎠⎦ 4 ⎣

=

(9)

Eq. (9) is spectrum information of echo signal reflected from target. Through Eq. (9), if the distance of target is D , there is a peak 2D in position ω = k ( c − t0 ) of the signal spectrum. Therefore, target distance is expressed as:

D=

(

ω k

)

+ t0 c 2

.

(10)

At this point, target distance is obtained through heterodyne process. Next, intensity distribution is reconstructed from the value of peak. In the spectrum, peak value corresponds to total intensity of echo light reflected from target at distance D . Through ghost imaging method, intensity distribution of target at distance D can be calculated. Suppose there are M measurements to reconstruct intensity distribution of target. According to ghost imaging theory, the intensity distribution of target is [20]:

1 G(i) (⇀ r)= M where

M

r ) Ib(i) − ∑ Is (i) (⇀ i=1

1 M

M

r )⋅ ∑ Is (i) (⇀ i=1

1 M

M

∑ Ib(i) , i=1

Ib(i) is current intensity of bucket detector in

(11) ith

X. Yang et al. / Optics Communications 368 (2016) 1–6

3

measurement, Is(i) (⇀ r ) is the distribution of speckle, and M is total number of measurements. Ib(i) equals to the peak value in signal spectrum. Therefore, Ib(i) can be expressed as:

Ib(i) =

χ 4

(12)

Submitting Eqs. (12) into (11), target intensity distribution can be calculated. Through Eqs. (10) and (11), target intensity distribution and range information can be calculated by mean of ghost imaging and heterodyne process respectively. Combining these two information, range images of target can be obtained by the HGI system. In HGI system, range information is obtained by heterodyne process. Range resolution is related to resolution of spectrum which is influenced by frequency change ratio and modulation time of modulating signal A . Therefore, sampling rate of data acquisition system does not affect range resolution of HGI system. HGI system proposed in this study can obtained high quality range images with low sampling rate.

Fig. 2. Distribution of target. Table 1 The values of noise parameters.

3. Numerical simulations of HGI

Parameters

To verify the performance of the HGI scheme, we present some numerical simulations in this section. In numerical simulations, targets are selected as three characters “HIT”, which are all planar targets. Fig. 2 shows targets distribution of characters. The distances between three characters and detection system are different. Character “H” is nearest and located at 2900 m away from detection system. Character “I” is farther than character “H” and vertical spacing between two characters is 2 m. Character “T” is farthest of all vertical spacing between “I” and “T” is 2 m as well. In Fig. 2, we use different colors to represent different range values. In the numerical simulations, the modulating signal of electro-op1 tical crystal in Eq. (1) is A = sin [2π (ω0 + 2 kt ) t ], where initial frequency value ω0 = 0Hz, frequency change ratio k¼2.5  104 GHz/s, modulation band width B = 500MHz , modulation time τ is 20μs . In the HGI system, modulation signal should be delayed some time, which is selected as t0 = 1.92 × 10−5s. Echo signal detected by a bucket detector is sampled by data acquisition system, which data sampling rate is only Fs = 40MHz. In practical experiment, noise has a certain degree of influence on HGI system. It is because that HGI system obtains range information through peak location in spectrum. Noise will cause peak jitter and actual peak location cannot be obtained precisely. Besides, intensity image will be also influenced by detector shot noise. In order to ensure that HGI system can work in practical experimental, it is supposed that there is noise in bucket detector which is chosen as an APD detector, and noise can be expressed as: 2 In ¯(t ) = 2eΔf ⎡⎣ Is ¯(t ) + Ib ¯(t ) + Id ¯(t ) F + IL ⎤⎦ + I¯T2 ,

(

)

Values

Ib (t )

4 × 10−9A

Id (t ) F IL

2 × 10−8A 5

Δf

k T R e

10−8A 1GHz 1.38 × 10−23J/K 293K 1kΩ 1.6 × 10−19C

signal after FFT, and vertical axis represents intensity of Fourier spectrum. As Fig. 3 shown, there are three peaks in Fourier spectrum. Horizontal coordinates of the three peaks are corresponding to target ranges and values of peaks represent echo light intensity from the three targets. Due to noise influence, peak location may change in different measurement. Therefore, peak locations of different measurements are averaged. According to average values of three peak horizontal coordinates and Eq. (10), the range information of targets can be obtain and shown in Table 2. In Section 2, we have demonstrated that targets distributions can be obtained through peaks values in Fourier spectrum graph. In Eq. (11), Ib(i) equals to peak values. Simulated targets distributions are shown in Fig. 4.

(13)

where Is (t ), Ib (t ), and Id (t ) are the echo-signal, background, and dark currents of APD detector respectively. e is the charge of electron, Δf is the receiver transmission band, F is the coefficient of noise magnification, and IL is the leakage current. First term of Eq. (1) on the right side 2eΔf [(Is (¯t ) + Ib ¯(t ) + Id ¯(t )) F + IL ] is mean 4kT Δf square value of APD shot noise. Besides, IT¯2 = R is thermal noise current of APD detector, where k is the Boltzmann's constant, T is the effective noise temperature, and R is the resistance equivalent to output resistance. In the simulation, the values of parameter relating to detector noise are shown in Table 1. Effective modulation pixels of SLM are 128 × 128, and there are 10,000 measurements to reconstruct targets range image. Fig. 3 is the frequency information graph of bucket detector in a certain measurement. In Fig. 3, horizontal axis represents frequency of

Fig. 3. Frequency information of bucket detector in a certain measurement.

4

X. Yang et al. / Optics Communications 368 (2016) 1–6

Table 2 Ranges and range errors of different targets. Target Average values of horizontal coordinates Calculated range values of of every peaks ω targets D “H”

3.437 × 106Hz

2900.622m

“I”

3.561 × 106Hz

2901.366m

“T”

6

3.878 × 10 Hz

2903.268m

Combining range information in Table 1 and intensity information in Fig. 4, range image of three targets is obtained. Fig. 5 shows numerical simulation range images with HGI method and conventional method respectively. Fig. 5(a) and (b) are original range image result and post-processed result with HGI method. Fig. 5(c) and (d) are original range image result and post-processed result with conventional ghost imaging method. Different colors are responding to different ranges. In ghost imaging method, backgrounds of all intensity distributions do not equal to 0. Therefore, in Fig. 5(a) and (c), there is amount of image noise existing in the background. In order to get rid of noise in the range image, intensity distribution images are carried on median filter processing. When any pixel values in the filtered images is below a chosen threshold value, the pixel values are set as 0. At last, range image is reconstructed using the processed intensity distribution images as Fig. 5(b) and (d) shown. In Fig. 5(a) and (b), the numerical simulation results verify the feasibility of HGI system and the validity theoretical derivation. Through numerical simulations, three character targets, which are located at 2.9 km away and spaced apart from each other at 2 m distance, can be distinguished with 40 MHz sampling rate of data acquire system and 500 MHz modulation bandwidth. However, when range image of targets is obtained with conventional 3D ghost imaging method as Fig. 5(c) and (d) shown, targets at intervals 2m cannot be distinguished effectively with same sampling rate. Character “I” and “H” share the same color deep yellow, which means that both located at about 2902.5 m away. The results do not conform to reality. Meanwhile, the colors of character “H” located at 2900 m away in Fig. 5(c) and (d) are corresponding to a certain value less than 2899 m. Range resolution of conventional 3D ghost imaging is relatively low. Numerical simulation results verify that HGI has higher range resolution than conventional ghost imaging.

4. Analysis of HGI and conventional method range image quality The HGI system utilizes modulated continuous wave laser as

source instead of pulse laser in conventional 3D ghost imaging. Continuous wave laser in the HGI system possesses high energy and can be used in long distance detection. Besides, the most advantage of HGI is that high quality of range image can be obtained with low sampling rate data acquired system. In this section, comparison of HGI and conventional 3D ghost imaging is performed through theory analysis. Range root-mean-square-error (RMSE) is utilized to evaluate the quality of range image with each method. RMSE is expressed as: 2

n

RMSE =

∑i, j ( Gi, j − Oi, j ) i×j

(14)

where Gi, j is the range image with through each method, and Oi, j is the real range value of the target. (i, j ) represents the coordinate of range images. According to Eq. (14) the RMSE of Fig. 5(a) through HGI system is RMSEH ¼0.7674 m; and the RMSE of Fig. 5 (c) through conventional 3D ghost imaging system is RMSEC ¼1.1094 m. By comparing the RMSE of range image with different methods, it can be concluded that the range error of range image through HGI system is smaller than that through conventional 3D ghost imaging system. Next, range resolution will be also compared between two systems. According to Eq. (10), we have:

ω 2D = − t0. k c

(15) B

In Eq. (15), frequency change ratio can be expressed as k = τ , where B is bandwidth of modulation signal and τ is modulation time. Through differential operation on both sides of Eq. (14), we have:

δD =

τcδω , 2B

(16)

where δD is range resolution, δω is resolution of spectrum. In HGI, sampling rate of data acquired system is Fs and numbers of samF pling is N . Resolution of spectrum can be expressed as δω = Ns . Modulation time, which equals to sampling time, is expressed as N τ = F . Therefore, Eq. (16) is simplified to: s

δD =

c . 2B

(17)

As Eq. (17) shown, range resolution of HGI is independent of data acquired system sampling rate Fs . Modulation bandwidth determines range resolution value. When increasing modulation bandwidth, range resolution is improved. Although the value of sampling rate Fs dose not influence range resolution shown in Eq.

Fig. 4. Results from left to right represent reconstructed intensity distributions of characters “H”, “I” and “T” respectively.

X. Yang et al. / Optics Communications 368 (2016) 1–6

5

Fig. 5. Numerical simulation results: (a) original range image with HGI method. (b) Post-processed range image with HGI. (c) Original range image with conventional ghost imaging. (d) Post-processed range image with conventional ghost imaging.

(17), it is not unlimited for HGI. According to Nyquist's theory, sampling rate Fs should be greater or equal to double maximum value of signal spectrum ωmax . As Fig. 3 shown, the maximum value of signal spectrum is ωmax = 20MHz. Therefore, sampling rate Fs needs to be meet the requirement Fs ≥ 2ωmax = 40MHz . In the numerical simulations, sampling rate is 40 MHz, which satisfies condition Fs = 2ωmax . Modulation bandwidth B is 500 MHz. Submitting the value of modulation bandwidth into Eq. (17), range resolution of HGI can be calculated as 0.3 m. When utilizing conventional 3D ghost imaging to measure range information of target, sampling rate of data Fs acquired system is the major factor affecting range resolution δD′, which can be expressed as:

δ D′ =

c . 2Fs

(18)

As Eq. (18) shown, range resolution and sampling rate have an inverse relationship. When sampling rate is selected as Fs = 40MHz, which is same to the sampling rate used in HGI numerical simulations, range resolution of conventional 3D ghost imaging δD′ is 3.75 m. The calculated result through Eq. (17) is theoretical limit of range resolution by HGI system. When signal-noise-ratio (SNR) is high enough, practical range resolution can approach the limits. If SNR is low, noise is strong. Noise can cause jitter of peak location in spectrum. There is a deviation between measured peak location and actual peak location. Peak location cannot be obtained precisely. Therefore, when SNR is higher, range resolution of HGI system becomes better. Same to HGI system, the calculated result through Eq. (18) is theoretical limit of range resolution by conventional 3D ghost imaging system. The range resolution is influenced by SNR in conventional 3D ghost imaging as well. In order to compare range resolution with the two systems under different SNR, Fig. 6 is drawn. Fig. 6 clearly shows range resolution of the two systems changing along with SNR. The graph’s horizontal axis shows the

SNR range, and the vertical axis shows range resolution. The blue line represents the range resolution changing with SNR of HGI system, and the red line represents that of conventional 3D ghost imaging. As Fig. 6 shown, range resolutions of two systems both decrease when SNR change larger. And under any noise condition, range resolution of HGI system is always smaller than that of conventional 3D ghost imaging. Therefore, it can be concluded that the quality of range image through HGI is better than that through conventional 3D ghost imaging under same noise condition. In this study, only noise effect is analyzed. However, in practical experiment, there may are some potential influences which can reduce range image quality. For instance, the actual location of peak cannot be obtained in spectrum owe to spectrum broadening and jitter of peak location. The reason for this potential influence is that the surface of practical target might be roughness. Atmosphere and noise can cause jitter of peak location. In order to reduce the influence and improve the range image quality, the centroid method can be used to obtain the accurate range

Fig. 6. Range resolution of each method changes along with SNR.

6

X. Yang et al. / Optics Communications 368 (2016) 1–6

information. In addition, influence of atmosphere turbulence can affect the quality of range image in practical experiment, because range image is obtained through intensity image in HGI system. Atmosphere turbulence can lower the quality of intensity image, although at a lower level. Adaptive ghost imaging method [21] and ghost imaging with shaped sources method [22] can be used in HGI system to reduce the impact of turbulence. More other influences of practical experimental setup will be studied in the future.

5. Conclusion In the study, we propose a HGI system, which utilizes electrooptical crystal and SLM to modulate intensity and distribution of light respectively. A single-pixel detector is used to collect echo light from targets. After multiplying Signals from the single-pixel detector and modulation signal of electro-optical crystal by the circuit, output is sampled by data acquired system and dealt with FFT. Range information and total echo light intensity can be obtained through Fourier spectrum. Reflected intensity and distribution of speckle are used to reconstruct target intensity image. Range image can be obtained through intensity image and range information. Numerical simulations are performed to verify feasibility of HGI. Three targets are located at distance of 2900 m and at intervals 2 m. HGI can successfully reconstruct range image of the targets. Under same sampling rate condition, range image of targets are obtained with conventional ghost imaging as a comparison. Conventional method fails to distinguish the targets at intervals 2 m. In addition, RMSEs of two systems are calculated and

range resolutions of two systems are analyzed in analysis section. The calculated results and simulations confirm that HGI can obtain high range resolution image with low sampling rate under noise contamination. HGI overcomes the limitation of sampling rate and has contributed to development of 3D ghost imaging.

References [1] A. Gatti, E. Brambilla, L.A. Lugiato, Phys. Rev. Lett. 93 (2004) 093602. [2] A. Gatti, M. Bache, D. Magatti, E. Brambilla, F. Ferri, A. Lugiato, J. Mod. Opt. 53 (2006) 739–760. [3] A. Valencia, G. Scarcelli, M. D'Angelo, Y. Shih, Phys. Rev. Lett. 94 (063601) (2005) 1–4. [4] A. Gatti, D. Magatti, M. Bache, E. Brambilla, F. Ferri, L.A. Lugiato, Phys. Rev. Lett. 94 (183602) (2005) 1–4. [5] B.I. Erkmen, J.H. Shaprio, Phys. Rev. A 79 (023833) (2009) 1–11. [6] K.W.C. Chan, M.N. O’Sullivan, R.W. Boyd, Phys. Rev. A 79 (2009) 033808. [7] M. Bache, E. Brambilla, A. Gatti, L. Lugiato, Opt. Express 12 (2004) 6067–6081. [8] J. Cheng, S. Han, Y. Yan, Chin. Phys. Soc. 15 (2006) 2002. [9] W. Chen, X. Chen, Eur. Phys. Lett. 110 (2015) 44002. [10] L. Basano, P. Ottonello, Appl. Phys. Lett. 89 (091109) (2006) 1–3. [11] J. Cheng, Opt. Express 17 (2009) 7916–7921. [12] N. Tian, Q. Guo, A. Wang, D. Xu, L. Fu, Opt. Lett. 36 (2011) 3302–3304. [13] D. Shi, S. Hu, Y. Wang, Opt. Lett. 39 (2014) 1231–1234. [14] W. Chen, X. Chen, Appl. Phys. Lett. 103 (2013) 221106. [15] W. Chen, X. Chen, Eur. Phys. Lett. 109 (2015) 14001. [16] S. Yuan, J. Yao, X. Liu, X. Zhou, Z. Li, Opt. Commun. 365 (2016) 180–185. [17] Y. Bromberg, O. Katz, Y. Silberberg, Phys. Rev. A 79 (053840) (2009) 1–4. [18] W. Gong, C. Zhao, J. Jiao, E. Li, M. Chen, H. Wang, W. Xu, S. Han, arXiv:1301.5767 [quant-ph], 2013. [19] H. Yu, E. Li, W. Gong, S. Han, Opt. Express 23 (2015) 14541–14551. [20] B.I. Erkmen, J.H. Shapiro, Phys. Rev. A 77 (043809) (2008) 1–3. [21] D. Shi, C. Fan, P. Zhang, J. Zhang, H. Shen, C. Qiao, Y. Wang, Opt. Express 20 (2012) 27992–27998. [22] C. Luo, J. Cheng, Opt. Lett. 28 (2013) 5381–5384.