## Abstract

Arbitrary shaping of the on-axis intensity of Bessel beams requires spatial modulation of both amplitude and phase. We develop a non-iterative direct space beam shaping method to generate Bessel beams with high energy throughput from direct space with a single phase-only spatial light modulator. For this purpose, we generalize the approach of Bolduc et al. to non-uniform input beams. We point out the physical limitations imposed on the on-axis intensity profile for unidirectional beams. Analytical, numerical and experimental results are provided.

© 2016 Optical Society of America

## 1. Introduction

Bessel beams constitute a class of solutions to Helmholtz equation that are propagation- invariant, i.e. “diffraction-free”. These beams are produced by the interference between an infinite number of plane waves whose wavevectors are distributed on the generatrix of a cone. The transverse profile is made of an intense central core surrounded by a large number of circular lobes of lower intensity. The central core of the beam defines an extremely long focal line and these beams have found various applications both in the linear and nonlinear regime, such as imaging [1–3], optical trapping [4], nonlinear optics [5–8], laser micromachining [9,10].

Since an ideal Bessel beam carries an infinite energy, only apodized solutions can be experimentally realized [1]. Because of this spatial limitation, the intensity of Bessel beams varies along the propagation. While the length of the interference field, *the Bessel zone*, is usually orders of magnitude larger than the Rayleigh range, the variation of the intensity in the main central core hampers several applications where the exact value of the intensity is a key parameter: nonlinear optics, laser nanoprocessing, light-matter interaction.

In this context, we need to control the variation of the intensity of the central core of the Bessel beam, *i.e.* its on-axis intensity distribution. We aim at creating several longitudinal profiles of intensity: flat-top, linear ramps or arbitrary variations, with a femtosecond laser source in order to open space for further applications in the nonlinear regime. We note that with 100 fs pulses, the pulse bandwidth can be neglected for the shaping operation [11].

Many techniques to engineer the intensity distribution along the propagation of Bessel beams have already been reported. Flattening the on-axis intensity distribution has been realized with passive optical components, creating approximations of Bessel beams [12–14]. Zamboni-Rached [15,16] demonstrated that the superposition of several ideal Bessel beams with different spatial frequencies enables the arbitrary control of the on-axis intensity. This technique was used with continuous superpositions of Bessel beams by Dartora *et al* [17] and Cizmar *et al* [18]. In these works, a direct relation is established between the desired on-axis intensity profile and the spatial spectrum of the beam. The peak of the spectrum, corresponding to the spatial frequency,${k}_{r0}={k}_{0}\mathrm{sin}\left(\theta \right)$can be interpreted as the carrier frequency. ${k}_{0}$ and $\theta $ are respectively the norm of the wavevector and the cone angle, i.e. the angle made by the wavevector with the optical axis. This approach was further generalized to nondiffracting hollow beams carrying a vortex charge and to propagation into absorbing media [19], as well as with a superposition of finite-energy Bessel Gauss beams [20].

Spatial Light Modulators (SLMs) are ideally suited for this application, but three issues have to be solved: i/ for applications in the nonlinear regime, the shaping approach needs preserving high-energy throughput, *i.e.* shaping is required in the direct space [11] by opposition to Fourier space shaping [18], ii/ spatial shaping in the direct space requires the modulation of both spatial amplitude and phase of an incoming Gaussian laser beam while SLMs can apply only either phase or intensity modulation and iii/ the phasemask design needs an efficient technique, if possible non-iterative, to gain both time and accuracy.

Arbitrary longitudinal profiles — yet limited by apodization—were created using Fourier space shaping, with phase-masks designed using the iterative Gerchberg-Saxton algorithm, to encode the whole information of the desired spatial spectrum into a single phase-only SLM [18,21,22]. However, Fourier space beam shaping techniques yield extremely low energy conversion. References [23–25] reported experimental demonstration of arbitrary profiles with direct-space shaping, In these works, SLM are used for amplitude modulation [26] which is detrimental to the final energy throughput of the system ( + 1 and −1 diffraction orders are not distinguished). Here, we aim phase modulation so as to preserve energy throughput.

Phase and amplitude shaping with a single phase-only hologram in direct space has been reported with several different approaches. One technique consists in using interlaced phase holograms in the SLM [27], where the desired field is retrieved at zeroth diffraction order. A second approach creates the target field in the first order of diffraction: the modulation in amplitude is obtained by modulating the level of phase wrapping, which is directly linked to the diffraction efficiency [28–30]. The benefit of this approach is to filter out the unmodulated part of the zeroth order that arises from imperfections and sub-100% fill-factor of the SLM. Bolduc *et al* demonstrated how to exactly encode the amplitude and phase in single phase-only SLM [31], but needed in practice an approximation for the experimental realization. Other techniques, more complex, use combination of multiple elements, such as two SLMs [32] and a Digital Micromirror Device (DMD) [33]

We note that all these techniques were developed for uniform (plane-wave) illumination of the shaping elements. Here, we demonstrate arbitrary shaping of the on-axis intensity of Bessel beam. For this purpose, we expand the analytical framework for exact arbitrary encoding in the first diffraction order [31] to the case of non-uniform beams. We show how to avoid the approximation used in Ref [31]. for experimental beam shaping. We clarify the physical limitation to the arbitrary shaping in terms of a convolution by a longitudional response function and show that targets can be quasi-exactly reached provided they are appropriately chosen. We compare numerically and experimentally Bessel beam shaping with uniform and uniformly growing on-axis intensity profiles with a total energy throughput of ~10%.

## 2. Bessel beams with arbitrary on-axis intensity distribution

In the paraxial and scalar approximations, a radially symmetric optical field *E*(*r*, *z* = 0), such as a Bessel beam, can be viewed as a superposition of an infinity of zeroth order Bessel functions of the first kind *J*_{0} .

where *z* is the distance of propagation, *r* is the transverse radial coordinate, and *k _{r}* is the transverse spatial frequency corresponding to the variable

*r*.

*S*(

*k*,

_{r}*z*= 0) is the amplitude of the spatial spectrum or, in other words, the Hankel transform of the optical field

*E*(

*r*,

*z*= 0). From Eq. (1), Cizmar et al demonstrated that the spatial spectrum is linked to the desired on-axis intensity profile

*I*(

*z*), defined at

*r*= 0, as follows [18]:

where *k _{0}* and

*k*are respectively the wave-vector of the beam and the longitudinal spatial frequency. The conical structure of the beam is taken into account by the Bessel phase term exp(

_{z}*ik*

_{z}_{0}

*z*), where

*k*=

_{z0}*k*cos(

_{0}*θ*) is the longitudinal Bessel frequency.

Therefore, Eqs. (1) and (2) define the initial amplitude field $E\left(r,z=0\right)$ that will produce a Bessel beam of given cone angle and on-axis intensity profile *I*(*z*).

We note that for a Bessel beam sufficiently long, the intensity distribution of its spatial spectrum is generally a thin annulus. In the case of Fourier space shaping, it is this amplitude distribution that is realized by amplitude masking of a Gaussian beam [18]. The energy throughput is very weak in this case. In contrast, direct space shaping requires only to apply the phase identical to the one of an axicon on the input Gaussian beam. The energy throughput is much higher.

In the following, we show how to shape an input imperfect Gaussian beam in direct space to generate the target Bessel field $E\left(r,z=0\right)$.

## 3. Direct space shaping of the on-axis intensity of Bessel beams with a phase-only SLM

We produce the Bessel field from the image of a Spatial Light Modulator by a 4-f telescope.

*Experimental setup:* Our setup is described in Fig. 1(a). The light source is an amplified Ti:Sa laser emitting 120 fs pulses at a wavelength of 800 nm. The beam is shaped by a Liquid Crystal On Silicon phase-only SLM (Hamamatsu LCOS-SLM, 792x600 pixels, 20 µm pitch). The Gaussian beam of diameter 4.7 mm at 1/e^{2} is obliquely incident onto the SLM with a small angle of 7°. Then, a 4f optical arrangement realizes two subsequents Fourier transforms with a total demagnification factor of 1/55. Spatial filtering is performed in the Fourier space of lens L_{1}. It selects out only the first order of diffraction, which is the spatial spectrum expressed in Eq. (2). The second optical Fourier transformation is realized with a microscope objective (MO_{1,} NA = 0.3). The Bessel beam is then formed from the focal plane of MO_{1}. The longitudinal position z = 0 corresponds to the focal point of MO_{1}.

The imaging system has a magnification of 27.7 and consists of a microscope objective MO_{2} (NA = 0.4), an imaging lens L_{2} and a CMOS camera (pixel pitch 5.2 μm), placed on a motorized translation stage. We sequentially reconstruct the beam propagation by recording images along the propagation distance *z*.

*Phase-mask design:* Now, we derive the expression of the phase mask to apply on the SLM to generate the target spectrum *S*(*k _{r}*) in the first order of diffraction, for an imperfect but reasonably smooth input amplitude

*A*. It is important to realize that the SLM mask and the plane

_{inc}*z = 0*are optically conjugated by the 4-f telescope. The amplitude distribution in the focal plane of MO

_{1}is a filtered image of the amplitude diffracted in first order by the SLM, with a demagnification factor (-f

_{2}/f

_{1}). For clarity, we ignore the magnification factor in the analytical derivation below.

The target Bessel beam initial amplitude $E\left(r,z=0\right)$can be expressed in terms of its spatial amplitude *A* and phase *ϕ* . To encode both informations into a single phase-only SLM, Davis *et al* [28] suggested to modulate the phase wrapping in order to locally tune the amount of light energy diffracted into the first order. Indeed, the general expression of the imprinted phase mask is given by [31]:

where *m* and *n* correspond to the pixels of the SLM. *M* is a normalized bounded positive function of amplitude (0 ≤ *M* ≤ 1), *F* contains the phase information of the target field, *ϕ _{ref}* is a linear phase ramp used to separate the different diffraction orders, and mod is the modulo function, used for the phase wrapping operation. If the incident laser beam amplitude is

*A*, the complex amplitude of the field after reflection on the SLM is:

_{inc}After developing in Taylor expansion the exponent and expanding in Fourier series the periodic phase term *F + ϕ _{ref}*, we select out only the first term in the Fourier series. After inverse Fourier transform, this first term corresponds to the amplitude that will be diffracted in the 1st diffraction order [31]:

Here, we included in the derivation the slowly-varying amplitude *A _{inc}* of the beam incident on the SLM. By identifying

*U*

_{1}with the amplitude of the Bessel beam

*E*(

*r*,

*z*= 0) =

*A*e

^{i}

*(ignoring the telescope magnification in the expression) , the expressions of the amplitude and phase terms of the phase mask become simply:*

^{ϕ}^{−1}stands for the inverse function of the sinc function defined in the domain [-π,0]. In principle, this method requires the computation of the inverse sinc function at each point of the hologram. Bolduc

*et al*considered that this operation is time consuming and suggested to simplify the expression of

*M*to

*M*≈

*A/A*. However, this approximation is not suited for quantitative approaches. In Fig. 2, we choose as a target (solid line) the on-axis intensity of a Bessel beam, of a cone angle

_{inc}*θ*= 13°, shaped as flat-top with preceding and following sections shaped as parabolas. We have numerically implemented the expression provided by the approximation suggested in Ref [31]. There is only a qualitative agreement with the target. We note that this cone angle value, used throughout the manuscript, is still well within the scalar approximation.

To keep the precision of the analytic result (Eq. (6)), we have produced a lookup table for the inverse sinc function and in a second step, used linear interpolation to numerically compute the phase mask, which is very fast. As shown in Fig. 2 with red circles, the numerical beam propagation with this phase mask generates an on-axis intensity distribution that matches the target with an excellent agreement.

## 4. Physically realizable targets

Arbitrary on-axis intensity profiles can be realized with the approach described above. However, we stress here that we must restrict to physically realizable targets in order to obtain a quantitative agreement between the experimental/numerical profiles with the initial target. We avoid the apodization used in Ref [18].

The physically realizable spectrum in *k _{z}* space is bounded by

*k*

_{0}

*= 2π/λ*on the high-values side and by the numerical aperture NA ${k}_{z}^{\mathrm{min}}={k}_{0}\sqrt{1-N{A}^{2}}$ on the low value one. It reads: ${S}^{\mathrm{exp}}({k}_{z})={S}^{\text{target}}({k}_{z}).R({k}_{z})$, where R is the response function of the system. R = 1 in the interval ${k}_{z}^{\mathrm{min}}<{k}_{z}^{}<{k}_{0}$ and 0 outside.

In the space of propagation z, the response function impacts on the maximal variations of *I(z)* that are experimentally realizable. This latter can be written as a convolution $I\left(z\right)={\left|{A}_{}^{\text{target}}\left(z\right)\ast H\left(z\right)\right|}^{2}$, where *A*^{target}(*z*) is the desired complex on-axis envelope and *H*(*z*) is the Fourier transform of the response function *R.* It reads:

To minimize the distorsions due to the convolution product with the sinc function, the target intensity should not vary on scales smaller than the full width of the central lobe of *H*(*z*). This latter is given by: $\Delta H=4\pi /\left({k}_{0}-{k}_{z}^{\mathrm{min}}\right)$.

As an example, we consider the case of a Bessel beam with uniformly increasing on-axis intensity and infinite slope at the beam termination. Its axial propagation can be described by the following equation:

where I_{max} is the desired maximal intensity of the beam and *z*_{1} = 240 µm. This abrupt profile is shown as a solid line in Fig. 3(b). In Fig. 3(a), we plot as blue line, the modulus of the amplitude spectrum |*S*(*k _{z}*)|, defined by Eq. (2), corresponding to this on-axis profile. The spectrum of the abrupt profile in blue is truncated. The numerical propagation of the truncated spectrum exhibits significant oscillations of the on-axis intensity, as shown as a blue line with marker in Fig. (3.b).

In order to avoid such undesired features, Cizmar *et al.* [18] suggested to use a Gaussian envelope around the spatial spectrum. To clarify the situation, we choose to adapt the target so that the characteristic distance of intensity variation exceeds the width *ΔH =* 54 µm where the NA is 0.3 and the wavelength is 800 nm. For this, we replace the abrupt decrease by a parabolic decrease of the intensity profile as shown in Fig. 3(c):

The parabolic profile has the benefit of reaching zero values of the intensity at a finite distance (*z*_{2}). Here, we used *z*_{2} = 110 µm>*ΔH*, and we observe in Fig. 3(c) that a quasi-perfect agreement is obtained within less than 2%.

The spectrum of such profile is analytically derived in the appendix. As shown in Fig. 3(a) (red lines), the addition of the parabolic decrease significantly reduced the range of spatial frequencies at which the corresponding spatial spectrum is defined. As a result, the retrieved on-axis intensity exhibits fewer oscillations with weaker modulation depth (Fig. 3(c)).

## 5. Results and discussion

In this section, we present experimental results of the propagation of Bessel beams with arbitrary, physically realizable, on-axis intensity profiles using the setup described in Fig. 1. Figure 4(a) shows the intensity distribution of the beam incident on the SLM, which was used to adapt our phase masks to beams with finite extent as commented in section 3. Figure 4(b) shows the phase mask for plane-wave illumination and Fig. 4(c) shows the mask for our experimental input beam. For this case, the target intensity distribution is the flat-top profile shown in Fig. 5.

Figure 5 shows a comparison between simulations and experimental data in case of Bessel beams of cone angle 13 degrees, with flat-top and uniformly increasing on-axis intensities, with three different values of the intensity ramp. We use the same parameters for the parabolic as in Fig. 3.

In all cases, the intensity distributions of the reconstructed beams are in excellent agreement with the target and numerical results. The energy throughput slightly varies from one mask to the other due to the variation of diffraction efficiency, but the typical value is 10%. We note that the balance between the beam parameters (Bessel zone length, etc.) and the input beam distribution strongly influences the total energy throughput. With 10-100 µJ input pulse energies, the setup allows for reaching peak powers of 10^{14}W/cm^{2}, well suited to nonlinear optics and laser processing applications.

## Conclusion

We reported a method that enables the generation of high-intensity Bessel beams with arbitrarily shaped on-axis intensity profiles. We expressed the limitations imposed by numerical aperture and wavelength in terms of a convolution of the desired longitudinal intensity profile with a sinc-shaped function. We expanded the approach of Ref [31]. to non-uniform input beam profiles to encode in a phase-only SLM both phase and amplitude shaping, to generate target spatial spectra in the first order of diffraction, after Fourier filtering. We numerically and experimentally realized different profiles with excellent agreement between the target, numerical realization and experiment. The total high-energy throughput provided by the direct-space encoding and by maintaining finite-size input illumination profiles of the SLM instead of plane wave illumination, allows for generating high-intensity femtosecond pulses for future applications in nonlinear optics or laser processing.

## Appendix

The resolution of Eq. (2) in case of uniformly increasing on-axis intensity is as follows:

*a*=

*k*

_{z}_{0}-k

_{z}; S(x) and C(x) stand for the sine and cosine Fresnel integrals respectively. They are defined as follows:

## Acknowledgments

The authors acknowledge funding from Region Franche-Comte. This work has been performed in cooperation with the Labex ACTION program, contract ANR-11-LABX-01-01.

## References and links

**1. **D. McGloin and K. Dholakia, “Bessel beams: diffraction in a new light,” Contemp. Phys. **46**(1), 15–28 (2005). [CrossRef]

**2. **M. Zhao, H. Zhang, Y. Li, A. Ashok, R. Liang, W. Zhou, and L. Peng, “Cellular imaging of deep organ using two-photon Bessel light-sheet nonlinear structured illumination microscopy,” Biomed. Opt. Express **5**(5), 1296–1308 (2014). [CrossRef] [PubMed]

**3. **J. Zheng, Y. Yang, M. Lei, B. Yao, P. Gao, and T. Ye, “Fluorescence volume imaging with an axicon: simulation study based on scalar diffraction method,” Appl. Opt. **51**(30), 7236–7245 (2012). [CrossRef] [PubMed]

**4. **L. Li, N. Eckerskorn, R. A. Kirian, J. Küpper, D. P. DePonte, W. Krolikowski, W. M. Lee, H. N. Chapman, and A. V. Rode, “quasi-Bessel hollow beam as optical guide for micro-particles,” Proc. SPIE **8810**, 88100N (2013). [CrossRef]

**5. **S. P. Tewari, H. Huang, and R. W. Boyd, “Theory of third-harmonic generation using Bessel beams, and self-phase-matching,” Phys. Rev. A **54**(3), 2314–2325 (1996). [CrossRef] [PubMed]

**6. **A. Dubietis, P. Polesana, G. Valiulis, A. Stabinis, P. Di Trapani, and A. Piskarskas, “Axial emission and spectral broadening in self-focusing of femtosecond Bessel beams,” Opt. Express **15**(7), 4168–4175 (2007). [CrossRef] [PubMed]

**7. **M. A. Porras, A. Parola, D. Faccio, A. Dubietis, and P. D. Trapani, “Nonlinear unbalanced bessel beams: stationary conical waves supported by nonlinear losses,” Phys. Rev. Lett. **93**(15), 153902 (2004). [CrossRef] [PubMed]

**8. **P. Polesana, M. Franco, A. Couairon, D. Faccio, and P. Di Trapani, “Filamentation in Kerr media from pulsed Bessel beams,” Phys. Rev. A **77**(4), 043814 (2008). [CrossRef]

**9. **F. Courvoisier, J. Zhang, M. K. Bhuyan, M. Jacquot, and J. M. Dudley, “Applications of femtosecond Bessel beams to laser ablation,” App. Phys. A **112**(1), 29–34 (2013). [CrossRef]

**10. **F. Courvoisier, R. Stoian, and A. Couairon, “Ultrafast laser micro and nano processing with nondiffracting and curved beams,” J. Opt. Laser Technol. **80**, 125–137 (2016). [CrossRef]

**11. **L. Froehly, M. Jacquot, P. A. Lacourt, J. M. Dudley, and F. Courvoisier, “Spatiotemporal structure of femtosecond Bessel beams from spatial light modulators,” J. Opt. Soc. Am. A **31**(4), 790–793 (2014). [CrossRef] [PubMed]

**12. **Z. Jaroszewicz, J. Sochacki, A. Kolodziejczyk, and L. R. Staronski, “Apodized annular-aperture logarithmic axicon: smoothness and uniformity of intensity distributions,” Opt. Lett. **18**(22), 1893–1895 (1993). [CrossRef] [PubMed]

**13. **C. J. R. Sheppard and S. Mehta, “Three-level filter for increased depth of focus and Bessel beam generation,” Opt. Express **20**(25), 27212–27221 (2012). [CrossRef] [PubMed]

**14. **V. Osipov, V. Pavelyev, D. Kachalov, A. Zukauskas, and B. Chichkov, “Realization of binary radial diffractive optical elements by two-photon polymerization technique,” Opt. Express **18**(25), 25808–25814 (2010). [CrossRef] [PubMed]

**15. **M. Zamboni-Rached, “Stationary optical wave fields with arbitrary longitudinal shape by superposing equal frequency Bessel beams: Frozen Waves,” Opt. Express **12**(17), 4001–4006 (2004). [CrossRef] [PubMed]

**16. **M. Zamboni-Rached, E. Recami, and H. E. Hernández-Figueroa, “Theory of “frozen waves”: modeling the shape of stationary wave fields,” J. Opt. Soc. Am. A **22**(11), 2465–2475 (2005). [CrossRef] [PubMed]

**17. **C. A. Dartora, K. Z. Nobrega, A. Dartora, G. A. Viana, and H. T. S. Filho, “A general theory for the frozen waves and their realization through finite apertures,” Opt. Commun. **265**(2), 481–487 (2006). [CrossRef]

**18. **T. Čižmár and K. Dholakia, “Tunable Bessel light modes: engineering the axial propagation,” Opt. Express **17**(18), 15558–15570 (2009). [CrossRef] [PubMed]

**19. **M. Zamboni-Rached, “Diffraction-Attenuation resistant beams in absorbing media,” Opt. Express **14**(5), 1804–1809 (2006). [CrossRef] [PubMed]

**20. **M. Zamboni-Rached and M. Mojahedi, “Shaping finite-energy diffraction- and attenuation-resistant beams through Bessel-Gauss–beam superposition,” Phys. Rev. A **92**(4), 043839 (2015). [CrossRef]

**21. **L. Li, W. M. Lee, X. Xie, W. Krolikowski, A. V. Rode, and J. Zhou, “Shaping self-imaging bottle beams with modified quasi-Bessel beams,” Opt. Lett. **39**(8), 2278–2281 (2014). [CrossRef] [PubMed]

**22. **C. G. Durfee, J. Gemmer, and J. V. Moloney, “Phase-only shaping algorithm for Gaussian-apodized Bessel beams,” Opt. Express **21**(13), 15777–15786 (2013). [CrossRef] [PubMed]

**23. **T. A. Vieira, M. R. R. Gesualdi, and M. Zamboni-Rached, “Frozen waves: experimental generation,” Opt. Lett. **37**(11), 2034–2036 (2012). [CrossRef] [PubMed]

**24. **T. A. Vieira, M. Zamboni-Rached, and M. R. R. Gesualdi, “Modeling the spatial shape of nondiffracting beams: Experimental generation of Frozen Waves via holographic method,” Opt. Commun. **315**, 374–380 (2014). [CrossRef]

**25. **T. A. Vieira, M. R. R. Gesualdi, M. Zamboni-Rached, and E. Recami, “Production of dynamic frozen waves: controlling shape, location (and speed) of diffraction-resistant beams,” Opt. Lett. **40**(24), 5834–5837 (2015). [CrossRef] [PubMed]

**26. **V. Arrizón, G. Méndez, and D. Sánchez-de-La-Llave, “Accurate encoding of arbitrary complex fields with amplitude-only liquid crystal spatial light modulators,” Opt. Express **13**(20), 7913–7927 (2005). [CrossRef] [PubMed]

**27. **O. Mendoza-Yero, G. Mínguez-Vega, and J. Lancis, “Encoding complex fields by using a phase-only optical element,” Opt. Lett. **39**(7), 1740–1743 (2014). [CrossRef] [PubMed]

**28. **J. A. Davis, D. M. Cottrell, J. Campos, M. J. Yzuel, and I. Moreno, “Encoding amplitude information onto phase-only filters,” Appl. Opt. **38**(23), 5004–5013 (1999). [CrossRef] [PubMed]

**29. **J. Leach, M. R. Dennis, J. Courtial, and M. J. Padgett, “Vortex knots in light,” New J. Phys. **7**, 55 (2005). [CrossRef]

**30. **V. Arrizón, U. Ruiz, R. Carrada, and L. A. González, “Pixelated phase computer holograms for the accurate encoding of scalar complex fields,” J. Opt. Soc. Am. A **24**(11), 3500–3507 (2007). [CrossRef] [PubMed]

**31. **E. Bolduc, N. Bent, E. Santamato, E. Karimi, and R. W. Boyd, “Exact solution to simultaneous intensity and phase encryption with a single phase-only hologram,” Opt. Lett. **38**(18), 3546–3549 (2013). [CrossRef] [PubMed]

**32. **L. Zhu and J. Wang, “Arbitrary manipulation of spatial amplitude and phase using phase-only spatial light modulators,” Sci. Rep. **4**, 7441 (2014). [CrossRef] [PubMed]

**33. **B. Rodenburg, M. Mirhosseini, O. S. Magaña-Loaiza, and R. W. Boyd, “Experimental generation of an optical field with arbitrary spatial coherence properties,” J. Opt. Soc. Am. B **31**(6), A51–A55 (2014). [CrossRef]