## Abstract

An optical system for the generation of partially coherent beams with genuine cross-spectral density functions from spatially modulated globally incoherent sources is presented. The spatial intensity modulation of the incoherent source is achieved by quasi-planar metasurfaces based on spatial-frequency modulation of binary Bragg surface-relief diffraction gratings. Two types of beams are demonstrated experimentally: (i) azimuthally periodic, radially quasi-periodic beams and (ii) rotationally symmetric Bessel-correlated beams with annular far-zone radiation patterns.

Published by Chinese Laser Press under the terms of the Creative Commons Attribution 4.0 License. Further distribution of this work must maintain attribution to the author(s) and the published article’s title, journal citation, and DOI.

## 1. INTRODUCTION

We explore implementations of an experimental scheme, proposed by Gori and Santarsiero [1,2], for the generation of a wide range of spatially partially coherent optical beams. The scheme is based on a globally incoherent field with a specific spatially varying intensity distribution $p(\mathit{v})$ at the input plane of an optical system that performs a linear transformation described by an arbitrary spatial kernel function $H(\mathit{\rho},v)$. Partially coherent fields generated in this way are guaranteed to have genuine (non-negative definite) cross-spectral density functions (CSDs) $W({\mathit{\rho}}_{1},{\mathit{\rho}}_{2})$.

In our implementation the required spatial variation of the intensity distribution is created by first reducing the spatial coherence of a laser beam with a random rotating diffuser [3–5] and then imaging the scattered field onto an intensity-modulating metasurface. The kernel function of the optical system is taken to be of Fourier form, which is sufficient here since the fields to be demonstrated follow the Schell-model [6]. The generation scheme and its practical limitations are discussed in Section 2.

Among numerous possible options, we concentrate on a recently introduced class of azimuthally periodic and radially quasi-periodic partially coherent beams [7]. In a special case these beams have annular far-zone distributions, and in a certain limit they reduce to fundamental Laguerre–Gaussian correlated beams [8,9]. To make the present paper self-contained, we cover the main properties of these fields in Section 3.

The metasurfaces are realized as binary diffractive structures based on Bragg diffraction by surface-relief gratings as information carriers [10]. The spatial intensity information for coherence control is included by spatial-frequency modulation with a continuous variation of the local fill factor of the carrier grating [11,12]. The operating principle, design, and fabrication of the metasurfaces are the subject of Section 4.

The experimental results are presented in Section 5, along with simulations based on the theory in Section 3. Results of beam propagation over several distances $\mathrm{\Delta}z$ in the half-space $z>{z}_{0}$ are presented, and the complex degrees of coherence at the plane $z=0$, which give rise to the propagation properties, are measured with a mirror-based wavefront folding interferometer (WFI) [13].

While only particular examples of partially coherent beams are demonstrated here, the generation scheme is applicable to any beam covered by the model in Refs. [1,2]. It can also be extended to vector beams. These issues are discussed in Section 6 before conclusions.

## 2. DEVISING GENUINE PARTIALLY COHERENT BEAMS

Gori and Santarsiero [1] introduced a class of genuine spatially partially coherent fields, for which the CSD has a non-negative definite form

The schematic setup for generating spatially partially coherent beams with genuine CSDs is illustrated in Fig. 1(a). One starts from a globally incoherent field at the input plane $\mathrm{A}$ of the linear system $\mathrm{T}$, where the intensity distribution is proportional to the kernel function $H(\mathit{\rho},\mathit{v})$. The weight function $p(\mathit{v})$ is embedded in the metasurface at $\mathrm{A}$ to modulate the intensity profile; then the CSD given by Eq. (1) is produced at the output plane $\mathrm{B}$. If desired, the output field can be modulated by adding a complex-amplitude transmission filter with transmission function $t(\mathit{\rho})$ at plane $\mathrm{B}$, which transforms the CSD as

The system we use to generate the incoherent field at plane $\mathrm{A}$ is illustrated schematically in Fig. 1(b). In our implementation a linearly polarized HeNe laser beam (wavelength ${\lambda}_{0}=632.8\text{\hspace{0.17em}}\mathrm{nm}$) is focused to a width ${w}_{\mathrm{D}}$ on a standard rotating diffuser to reduce its spatial coherence area. The scattered field is collimated by lens ${\mathrm{L}}_{1}$ and passes through a Gaussian apodization filter $\mathrm{G}$ to produce a Gaussian Schell-model (GSM) beam [5]. The transmitted beam is then demagnified by an afocal system consisting of lenses ${\mathrm{L}}_{2}$ and ${\mathrm{L}}_{3}$ to illuminate the metasurface with a GSM beam of radius ${w}_{\mathrm{A}}$.

Admittedly, the use of a rotating diffuser is a bulky and somewhat cumbersome approach to reduce spatial coherence. For future needs we therefore put forward an alternative technical solution: moving the diffuser randomly in the $xy$ plane with a piezoelectric device. Random movement over a sufficient area is possible since piezoelectric devices allow fast 2D scanning over distances larger than the required coherence area. Strictly speaking, the field on $\mathrm{A}$ is not completely incoherent, but it is quasi-homogeneous. Therefore the conditions for the validity of Eq. (1) do not hold exactly, but non-negative definiteness is nevertheless ensured since the field is physically realizable.

## 3. PARTIALLY COHERENT MODEL BEAMS

We demonstrate the scheme in Fig. 1 by applying it to a recently introduced class of radially quasi-periodic and azimuthally periodic (RAP) beams [7]. These are generally Bessel-correlated beams with a Gaussian intensity profile at the source plane $z={z}_{0}$, which, in a special case, can have an annular distribution of radiant intensity. The main characteristics of RAP beams are briefly covered to make the paper self-contained and to summarize the formulas that enable numerical simulations of experimental results.

#### A. Radially Quasi-Periodic, Azimuthally Periodic Beams

We introduce RAP beams by assuming that the function $H(\mathit{\rho})$ in Eq. (2) has a Gaussian form $H(\mathit{\rho})=\mathrm{exp}(-{\mathit{\rho}}^{2}/{w}_{0}^{2})$. The weight function $p(\mathit{v})$ is written in polar coordinates $\mathit{v}=(v,\varphi )$ as

The complex degree of coherence (DOC) for solutions $m\ge 1$ is

A closed-form expression for the propagated CSD in the paraxial domain can be found with the standard Fresnel formula [7]. In particular, the spectral density, which we will measure in Section 5, evolves as

#### B. Beams with Annular Far-Zone Radiation Pattern

Let us consider RAP beams with $m=0$ and ${C}_{0}=\pm 1$ in Eq. (5), which gives a weight function

The spectral density at the source plane is still given by Eq. (7), but the complex DOC takes the form

The spectral density at $z>0$,

In the limiting case ${\rho}_{\mathrm{s}}\to 0$, and when ${C}_{0}=-1$, the complex DOC given in Eq. (15) reduces to

If ${C}_{0}=1$ and we again let ${\rho}_{\mathrm{s}}\to 0$, Eq. (15) becomes the DOC of a GSM beam. Figure 2 illustrates some properties of the $m=0$ solutions not discussed in Ref. [7]. Cross sections of the DOC plotted in Fig. 2(a) show central peaks and side lobes of expanding widths when the value of ${\rho}_{\mathrm{s}}$ in the ${I}_{0}$ term in Eq. (15) increases. The far-field intensity in Fig. 2(b) is also affected: the central minimum is close to zero at small values of ${\rho}_{\mathrm{s}}$, which also provides the widest central dark region.

## 4. METASURFACE DESIGN AND REALIZATION

In Ref. [7] we presented an on-axis encoding technique of the metasurface at plane $\mathrm{A}$, which can be realized in the axially symmetric geometry of Fig. 1. However, our experiments showed that this scheme has axial noise, which proved difficult to avoid. We therefore opted for an alternative, off-axis coding scheme. Our method is effectively similar to that used in Ref. [25], with the difference that we modulate the field amplitude with a more compact metasurface, instead of filtering the diffracted light from a spatial light modulator. In addition, we can obtain extremely high diffraction efficiency with the off-axis encoding. The drawback of our technique is the somewhat trickier alignment of the optical components compared to the on-axis method.

#### A. Operating Principle and Carrier Design

The metasurface is realized using a high-frequency carrier grating as introduced in Ref. [10]. This is a binary surface-relief grating with depth $h$ and fill factor $f=c/d$ as design parameters, as defined in Fig. 3(a). The grating operates at large angles of incidence $\theta \sim 30\xb0$. Hence the geometry of Fig. 1 is folded at $\mathrm{A}$ by twice this angle. The design angle of incidence is the Bragg angle $\theta ={\theta}_{\mathrm{B}}$ given by

Assuming that the incident field is linearly polarized in the direction of the grating grooves, we design the carrier grating for maximum efficiency of order $l=-1$ in TE polarization for ${\lambda}_{0}=633\text{\hspace{0.17em}}\mathrm{nm}$, $n=1.49$, and $d={\lambda}_{0}$.

The efficiency of the carrier grating as a function of fill factor $f=c/d$ was calculated using the Fourier modal method (FMM) [26]. The results are shown in Fig. 3(b). The efficiency has two maxima reaching almost unity, at $f\sim 0.16$ and $f\sim 0.62$. We employ the second since the efficiency drops smoothly to zero when $f$ approaches unity. Since the incident field is not fully collimated, we plot in Fig. 3(c) the efficiency as a function of the angle of incidence $\theta $ (Bragg selectivity of the grating) for several fill factors over the range $0\xb0<\theta <40\xb0$. It turns out that the angular selectivity is not a serious limiting factor for operation of the grating.

The intensity coding is performed as local variations of the fill factor $f$ as in Refs. [10,11]. Fill factors between 0.62 and 0.9, marked with dashed vertical lines in Fig. 3(b), are used for efficiency encoding since they are practically feasible. The target grating depth was chosen to be $h=1350\text{\hspace{0.17em}}\mathrm{nm}$, giving ${\eta}_{-1}\approx 0.95$ at the second maximum $f=0.62$.

#### B. Fabrication

The sample was fabricated by a standard electron beam lithography procedure. A 1-inch-diameter and 0.5-mm-thick ${\mathrm{SiO}}_{2}$ substrate was cleaned and spin-coated with a 1.4-μm-thick positive-tone resist layer (PMMA), which was coated with a conductive 30-nm-thick copper layer via thermal evaporation. After creating the grating patterns with electron beam lithography, the copper layer was removed in ${\mathrm{HNO}}_{3}$, followed by resist development in MIBK:IPA (1:3 solution) for 75 s and 30 s rinsing with IPA.

Scanning electron microscope images of one of the samples are shown in the Fig. 4. In Fig. 4(a) we show a cross section around a location where the fill factor changes abruptly from 0.8 to 0.62, and in Fig. 4(b) a top view of the region is shown where $f$ is fairly constant. The fill factor is fairly well reproduced, but in view of Fig. 4(a) the groove depth is somewhat reduced especially at large $f$. This could be improved by increasing the electron beam dose. However, its main effect is to reduce ${\eta}_{-1}$, which should be low at large $f$ according to Fig. 3. In fact the modulation of $h$ helps us to achieve the low values of ${\eta}_{-1}$ required beyond the second vertical line in Fig. 3(b).

## 5. EXPERIMENTAL RESULTS

The focal lengths of the lenses in Fig. 1(b) before plane A are flexible, but we chose ${\mathit{L}}_{1}=150\text{\hspace{0.17em}}\mathrm{mm}$, ${\mathit{L}}_{2}=200\text{\hspace{0.17em}}\mathrm{mm}$, and ${\mathit{L}}_{3}=40\text{\hspace{0.17em}}\mathrm{mm}$. We measured ${w}_{\mathrm{D}}=31\text{\hspace{0.17em}}\mathrm{\mu m}$ and ${w}_{\mathrm{A}}\approx 880\text{\hspace{0.17em}}\mathrm{\mu m}$. The GSM was used as in Ref. [5] to estimate the coherence radius at the filter $\mathrm{G}$ (Thorlabs NDYR20A) to be ${\sigma}_{\mathrm{G}}\approx 210\text{\hspace{0.17em}}\mathrm{\mu m}$, which means that the field is quasi-homogenous with ${\sigma}_{1}/{w}_{\mathrm{G}}\approx 1/21$. Measurement with a WFI gives the results in Fig. 5, where we show cross sections of the spectral density and the absolute value of the spatial DOC at plane $\mathrm{A}$, from which ${w}_{\mathrm{A}}\approx 880\text{\hspace{0.17em}}\mathrm{\mu m}$. While $|\mu (\mathrm{\Delta}x)|$ is not completely Gaussian but has some side lobes, we can estimate ${\sigma}_{\mathrm{A}}\approx 42\text{\hspace{0.17em}}\mathrm{\mu m}$. The field at plane $\mathrm{A}$ is definitely quasi-homogeneous and, since the coherence area is larger than $d$ by a factor of 60, the metasurface works well as a Bragg deflector.

The Gaussian spot at plane $\mathrm{A}$, with ${w}_{\mathrm{A}}\approx 880\text{\hspace{0.17em}}\mathrm{\mu m}$, could be used directly to approximate the Gaussian factor in $p(v)$. Hence, to encode $p(v)$, only the residual RAP factor needed to be realized with the metasurface. Explicitly, this factor is

for beams with $m>0$ and for $m=0$. Here, ${M}_{0}$ is chosen such that the maximum does not exceed unity. Sampled $100\times 100$ versions of the metasurfaces with spatially varying transmission efficiencies are illustrated in Fig. 6. In particular, we show the $m=3$, ${\rho}_{0}=3\text{\hspace{0.17em}}\mathrm{mm}$ solution in Eq. (20) and the $m=0$, ${\rho}_{0}=2\text{\hspace{0.17em}}\mathrm{mm}$ solution in Eq. (21), which are demonstrated experimentally. The modulation contrast of $p(v)$ is much weaker in the former case than in the latter. The linear system $\mathrm{T}$, which does the Fourier transform from plane $\mathrm{A}$ to plane $\mathrm{B}$, is effectively a positive lens of focal length ${f}_{\mathrm{H}}=50\text{\hspace{0.17em}}\mathrm{mm}$. Therefore, the spatial frequency $v$ at plane $\mathrm{A}$ is connected to the real-space coordinate $\rho $ via $v=\rho /{\lambda}_{0}{f}_{\mathrm{H}}$.#### A. Source-Plane Complex Degree of Coherence

The measured amplitude and phase of the DOC at plane $z={z}_{0}$ are shown in Fig. 7 for the two beams considered. The distributions of $\mu (\mathrm{\Delta}x,\mathrm{\Delta}y)$ are evaluated using a mirror-based WFI as in Ref. [27]. We show the measured phases only across a circular area of radius $\mathrm{\Delta}\rho <75\text{\hspace{0.17em}}\mathrm{\mu m}$ for $m=3$ and $\mathrm{\Delta}\rho <42\text{\hspace{0.17em}}\mathrm{\mu m}$ for $m=0$, where the absolute value of the DOC is significant and reliable phase measurements could be made.

The DOC of the $m=3$ solution has an azimuthal periodicity of 60° in the side lobe pattern, whereas the phase in the central region $\mathrm{\Delta}\rho <35\text{\hspace{0.17em}}\mathrm{\mu m}$ is constant. These features, as well as the sharp variations of $|\mu (\mathrm{\Delta}x,\mathrm{\Delta}y)|$, are evident in the measurements even though there are some distortions due to experimental imperfections particularly in the phase measurements within the central region of Fig. 7(e). The abrupt $2\pi $ phase transitions in the central region of Fig. 7(f) are due to wrapping into the $[-\pi ,\pi )$ interval.

#### B. Propagation of Beam Profiles

The propagation of beams that the DOCs in Fig. 7 generate is illustrated in Figs. 8 and 9, where the experimental and simulated results are shown in top and bottom rows, respectively. The experimental results show some speckle-like granularity that could be contributed by the finite coherence area at plane $\mathrm{A}$, where the field is not strictly speaking globally incoherent as assumed in the model. In particular, the discrepancies caused by the granularity and additional side lobes in the coherence function of the $m=3$ RAP beam in the Fig. 7 upper panel are visible in the measured intensity profiles in Fig. 8 at $z=90\text{\hspace{0.17em}}\mathrm{mm}$ and $z=140\text{\hspace{0.17em}}\mathrm{mm}$.

The evolution of the RAP beam with $m=3$ is shown in Fig. 8 at the source plane and three different propagation distances. We assume parameters $R=1/4.5$ and $q=1/40$, which gives ${z}_{\mathrm{G}}=567\text{\hspace{0.17em}}\mathrm{mm}$ and ${z}_{\mathrm{R}}=14.2\text{\hspace{0.17em}}\mathrm{mm}$. At $z=0$ the beam shape is Gaussian with a beam waist ${w}_{0}=338\text{\hspace{0.17em}}\mathrm{\mu m}$. As the propagation distance grows the beam acquires its RAP features. The annular $m=0$ beam in Fig. 9, with ${w}_{0}=210\text{\hspace{0.17em}}\mathrm{\mu m}$, $R=1/7.5$, $q=1/30$, ${z}_{\mathrm{G}}=218\text{\hspace{0.17em}}\mathrm{mm}$, and ${z}_{\mathrm{R}}=7.3\text{\hspace{0.17em}}\mathrm{mm}$, begins to show its dark center at distances $z\sim 10\text{\hspace{0.17em}}\mathrm{mm}$. This property is preserved throughout further propagation.

## 6. DISCUSSION

The RAP beams demonstrated here are not the only possible classes of partially coherent fields that we can generate with the experimental technique discussed in this paper. While Fourier-type kernels allow the demonstration of any Schell-model beam, the use of more general kernels opens the way to produce more general fields. For instance, metasurfaces acting as vortex phase filters with transmittance $t(\mathit{\rho})=t(\rho ,\varphi )=\mathrm{exp}(il\varphi )$, where $l$ is an integer, can be placed at plane $\mathrm{B}$ to produce orbital angular momentum. Likewise, polarizing optical components can contribute to the linear transformations, allowing the realization of vector beams. Space-variant imaging systems with beam shaping or polarization-control metasurfaces (such as polarization gratings [28,29]) in the intermediate Fourier plane of a $4F$ imaging system are options, as are transformations involving interferometers. For example, a WFI transforms any Schell-model beam into a corresponding specular beam [30,31]. We have recently demonstrated such specular transformations using a mirror-based WFI [27].

## 7. CONCLUSIONS

We have presented an implementation of an experimental system for the generation of spatially partially coherent beams with genuine correlation functions, using quasi-planar metasurfaces. The method was applied to particular radially quasi-periodic, azimuthally periodic beams. However, it is applicable to the generation of any physically realizable beams, whether scalar or vectorial.

## Funding

China Scholarship Council; Natural Science Foundation of Hunan Province (2023JJ30383); National Natural Science Foundation of China (61805080); Itä-Suomen Yliopisto; Academy of Finland (320165, 320166, 333938).

## Disclosures

The authors declare no conflicts of interest.

## Data Availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

## REFERENCES

**1. **F. Gori and M. Santarsiero, “Devicing genuine spatial correlation functions,” Opt. Lett. **32**, 3531–3533 (2007). [CrossRef]

**2. **F. Gori and M. Santarsiero, “Devicing genuine twisted cross-spectral densities,” Opt. Lett. **43**, 595–598 (2018). [CrossRef]

**3. **W. Martiensen and E. Spiller, “Coherence and fluctuations in light beams,” Am. J. Phys. **32**, 919–926 (1964). [CrossRef]

**4. **J. D. Farina, L. M. Narducci, and E. Collett, “Generation of highly directional beams from a globally incoherent source,” Opt. Commun. **32**, 203–208 (1980). [CrossRef]

**5. **Q. He, J. Turunen, and A. T. Friberg, “Propagation and imaging experiments with Gaussian Schell-model beams,” Opt. Commun. **67**, 245–250 (1988). [CrossRef]

**6. **A. C. Schell, “A technique for the determination of the radiation pattern from a partially coherent aperture,” IEEE Trans. Ant. Propag. **15**, 187–188 (1967). [CrossRef]

**7. **M. Luo, M. Koivurova, and J. Turunen, “Azimuthally periodic and radially quasi-periodic Bessel-correlated fields,” Opt. Express **30**, 11754–11766 (2022). [CrossRef]

**8. **Z. Mei and O. Korotkova, “Random sources generating ring-shaped beams,” Opt. Lett. **38**, 91–93 (2013). [CrossRef]

**9. **M. Luo and D. Zhao, “Elliptical Laguerre Gaussian Schell-model beams with a twist in random media,” Opt. Express **27**, 30044–30054 (2019). [CrossRef]

**10. **J. Turunen, P. Blair, J. M. Miller, M. R. Taghizadeh, and E. Noponen, “Bragg holograms with binary synthetic surface-relief profile,” Opt. Lett. **18**, 1022–1024 (1993). [CrossRef]

**11. **E. Noponen and J. Turunen, “Binary high-frequency-carrier diffractive optical elements: electromagnetic theory,” J. Opt. Soc. Am. A **11**, 1097–1109 (1994). [CrossRef]

**12. **E. Tervonen, J. Turunen, and J. Pekola, “Pulse-frequency-modulated high-frequency-carrier diffractive elements for pattern projection,” Opt. Eng. **33**, 2579–2587 (1994). [CrossRef]

**13. **A. Halder, H. Partanen, A. Leinonen, M. Koivurova, T. K. Hakala, T. Setälä, J. Turunen, and A. T. Friberg, “Mirror-based scanning wavefront-folding interferometer for coherence measurements,” Opt. Lett. **45**, 4260–4263 (2020). [CrossRef]

**14. **A. T. Friberg and R. Sudol, “Propagation parameters of Gaussian Schell model beams,” Opt. Commun. **41**, 383–387 (1982). [CrossRef]

**15. **J. Durnin, “Exact solutions for nondiffracting beams: I. The scalar theory,” J. Opt. Soc. Am. A **4**, 651–654 (1987). [CrossRef]

**16. **J. Durnin, J. J. Miceli Jr., and J. H. Eberly, “Diffraction-free beams,” Phys. Rev. Lett. **58**, 1499–1501 (1987). [CrossRef]

**17. **R. H. Jordan and D. H. Hall, “Free-space azimuthal paraxial wave equation: the azimuthal Bessel–Gauss beam solution,” Opt. Lett. **19**, 427–429 (1994). [CrossRef]

**18. **D. H. Hall, “Vector–beam solutions of Maxwell’s wave equation,” Opt. Lett. **21**, 9–11 (1996). [CrossRef]

**19. **V. Bagini, F. Frezza, M. Santarsiero, G. Schettini, and G. S. Spagnolo, “Generalized Bessel–Gauss beams,” J. Mod. Opt. **43**, 1155–1166 (1996). [CrossRef]

**20. **J. Turunen, A. Vasara, and A. T. Friberg, “Propagation-invariance and self-imaging in variable-coherence optics,” J. Opt. Soc. Am. A **8**, 282–289 (1991). [CrossRef]

**21. **S. A. Ponomarenko, W. Huang, and M. Cada, “Dark and antidark diffraction-free beams,” Opt. Lett. **32**, 2508–2510 (2007). [CrossRef]

**22. **A. Halder, M. Koivurova, H. Partanen, and J. Turunen, “Paraxial propagation of a class of Bessel-correlated fields,” Opt. Express **26**, 11055–11067 (2018). [CrossRef]

**23. **M. Woerdemann, C. Alpmann, M. Esseling, and C. Denz, “Advanced optical trapping by complex beam shaping,” Laser Photon. Rev. **7**, 839–854 (2013). [CrossRef]

**24. **S. W. Hell and J. Wichmann, “Breaking the diffraction resolution limit by stimulated emission: stimulated-emission-depletion fluorescence microscopy,” Opt. Lett. **19**, 780–782 (1994). [CrossRef]

**25. **Y. Cai, Y. Chen, and F. Wang, “Generation and propagation of partially coherent beams with nonconventional correlation functions: a review [Invited],” J. Opt. Soc. Am. A **31**, 2083–2096 (2014). [CrossRef]

**26. **H. Kim, J. Park, and B. Lee, *Fourier Modal Method and Its Applications in Computational Nanophotonics* (CRC Press, 2017).

**27. **D. Das, A. Halder, H. Partanen, M. Koivurova, and J. Turunen, “Propagation of Bessel-correlated specular and antispecular beams,” Opt. Express **30**, 5709–5721 (2022). [CrossRef]

**28. **F. Gori, “Measuring Stokes parameters by means of a polarization grating,” Opt. Lett. **24**, 584–586 (1999). [CrossRef]

**29. **J. Tervo and J. Turunen, “Paraxial-domain diffractive elements with 100% efficiency based on polarization gratings,” Opt. Lett. **25**, 785–786 (2000). [CrossRef]

**30. **F. Gori, G. Guattari, C. Palma, and C. Padovani, “Specular cross-spectral density functions,” Opt. Commun. **68**, 239–243 (1988). [CrossRef]

**31. **H. Partanen, N. Sharmin, J. Tervo, and J. Turunen, “Specular and antispecular light beams,” Opt. Express **23**, 28718–28727 (2015). [CrossRef]