In recent years, nondiffracting beams have attracted much attention due to their peculiar properties^[1,2]. Such beams are able to keep their transverse intensity distribution over a long distance in free space and quickly reestablish their structures when encountering barriers^[3]. Since the pioneering work of Durnin et al., Bessel beams, which are the exact solution of the Helmholtz wave equation in circular cylindrical coordinates, were presented^[4] and have shown outstanding performance in the fields of optical trapping^[5,6], optical microscopy^[7,8], material processing^[9], and so on. Apart from Bessel beams, plane waves^[10], Mathieu beams^[11], and Weber beams^[12] are the families of propagation invariant beams, which are, respectively, the nonparaxial solutions of the Helmholtz equation under Cartesian coordinates, elliptic cylindrical coordinates, and parabolic cylindrical coordinates. In particular, a family of nondiffracting beams forms a set of complete orthogonal solutions. In this case, any linear superposition of the same type of diffraction-free modes with the same radial wave vector keeps the nondiffracting properties^[13]. Recently, Kovalev et al. presented a linear superposition of Bessel modes, whose complex amplitude can be expressed in terms of Lommel functions of two variables. Hence, such a beam is called Lommel beam^[14]. Compared to Bessel beam whose transverse intensity profile is radial symmetry, the Lommel beam has a reflective symmetry with respect to both Cartesian coordinate axes^[15]. In addition, a Lommel beam has an adjustable symmetrical transverse intensity pattern and continuously variable orbital angular momentum (OAM), which can benefit optical manipulation^[16,17] and optical information encoding^[18,19].

When optical beams propagate in materials with inhomogeneous index distributions, they will suffer distortions or scattering. In order to promote the applications of Lommel beams in different materials, the propagation properties of Lommel beams in some aberration-weak media such as gradient index medium^[20] and weak ocean turbulence^[21,22] have been explored. However, the wavefront of the beam could be completely scrambled when it encounters the highly scattering media^[23–25], such as the strong atmosphere or biological tissues, which restricts its applications behind the highly scattering media. To overcome the scattering, wavefront shaping techniques have opened up the possibility to reconstruct the structured light beam through highly scattering media^[26–31]. Up to now, some kinds of the structured light beams such as Bessel beams^[26], donut beams^[26–28], optical needle beams^[29], optical bottle beams^[31], and perfect vortex beams^[30] have been generated through highly scattering media. However, to the best of our knowledge, the construction of the Lommel beam through highly scattering media has not been demonstrated yet.

In this paper, we utilize a digital micromirror device (DMD) to produce Lommel beams through a highly scattering medium by means of a transmission matrix (TM)-based point spread function (PSF) engineering method. All of the experimental results of the constructed Lommel beams agree well with the theoretical predictions. Moreover, the optical field distributions and OAM can be engineered continuously by adjusting the parameters in the angular spectrum. The intensity profiles of the generated beams were recorded, and their phase distributions were measured by the phase shifting method. Furthermore, the raster scanning of Lommel beams was demonstrated, and a number of Lommel beams were constructed simultaneously. The constructed Lommel beams under high scattering could promote the development of optical communication and optical manipulation behind highly scattering media.

Figure 1 illustrates the principle. As shown in Fig. 1(a), when an optical Lommel beam propagates through a highly scattering medium, the beam occurs with multiple scattering, and the wavefront is scrambled, thus producing a random speckle field distribution behind the highly scattering medium. In this situation, the performance of the Lommel beam is completely suppressed. For the mathematical description, the optical field at the input plane and output plane can be, respectively, set as an $N\times 1$ vector and an $M\times 1$ vector, where $N$ and $M$ are the numbers of the input modes and output modes. Due to the scattering process being linear and deterministic, the relationship between the output field ($E_{\mathrm{out}}$) and input field ($E_{\text{in}}$) can be described as^[32]$$E_{\mathrm{out}}={\mathit{TE}}_{\text{in}},$$(1)where $T$ is the TM of the highly scattering medium, which connects the relationship between the input field and output field. The TM can be derived from $4N$ speckle intensity patterns at the output plane with the phase shifting method, as shown in Fig. 1(b). Then, according to the method of TM-based PSF engineering^[26], in order to engineer the PSF distribution at the focus, a numerical filtering operation needs to be performed on the virtual Fourier plane of the output mode of $T$, which is determined as $$T^{\text{filt}}={\Im }^{-1}\{\Im \{T\}\times A\},$$(2)where $\Im \{\}$ and ${\Im }^{-1}\{\}$ are, respectively, the Fourier transform and inverse Fourier transform. $A$ is the Fourier transform of the desired PSF distribution at the focus and is employed as a mask to filter the calibrated TM.

For the diffraction-free Lommel beam^[14], its complex amplitude in cylindrical coordinates can be written as $$E(r,\phi ,z)=c^{-n}\exp \left(iz\sqrt{k^2-{\alpha }^2}\right)U_n[c\alpha r\exp (i\phi ),\alpha r],$$(3)where $c$ is a dimensionless complex coefficient with a general form of $c=c_0\exp (i{\phi }_0)$, the non-negative modulus $c_0$ takes values smaller than unity to guarantee the convergence of the series, and ${\phi }_0$, which ranges in (0, $2\pi$), is the angle of $c$. $k=2\pi /\lambda$ is the wave number of light with wavelength $\lambda$, $n$ is the topological charge, and $\alpha$ is the transverse wavenumber. $U_n$ is the Lommel function of two variables, which is given by $$U_n(\omega ,\xi )=\sum _{p=0}^{\infty }{(-1)}^p{\left(\frac{\omega }{\xi }\right)}^{n+2p}{J}_{n+2p}(\xi ).$$(4)

The corresponding Fourier transform of the Lommel beam is expressed as $$A(\rho ,\phi )=\frac{1}{\lambda \alpha }\frac{{(-i)}^n\exp (\mathit{in}\phi )}{1-{[c\exp (i\phi )]}^2}\delta (\rho -\frac{\alpha }{k}),$$(5)where $(\rho ,\phi )$ are the polar coordinates in the spectral plane, and $\delta (x)$ is the Dirac delta function.

From Eq. (5), we can see that the angular spectrum of the ideal Lommel beam is an infinitely thin ring due to the existence of $\delta (x)$. According to Ref. [26], when the thin ring is directly used as a mask for numerical filtering in the method of TM-based PSF engineering, it will remove most of the low frequencies of the input field that contributes to the focus, leading to the generated PSF distribution at the focal plane with a low signal-to-backgound ratio (SBR). Thus, in order to generate Lommel beams with a high SBR, we set the annular mask with a thickness and modify Eq. (5) as $$A(\rho ,\phi )=\frac{1}{\lambda \alpha }\frac{{(-i)}^n\exp (\text{in}\phi )}{1-{[c\exp (i\phi )]}^2}[\mathrm{circ}\left(\frac{\rho }{{\rho }_1}\right)-\mathrm{circ}\left(\frac{\rho }{{\rho }_2}\right)],$$(6)where $\text{circ}(x)$ is the circle function, i.e., inside the circle, $x<1$, $\text{circ}(x)=1$, whereas outside the circle, $\text{circ}(x)=0$; ${\rho }_1$ and ${\rho }_2$ are the outer and inner radius of the ring, respectively.

After the complex $A(\rho ,\phi )$ is constructed, we employed it as a mask to filter the calibrated TM as Eq. (2), and the resulting operator $T^{\text{filt}}$ can then be used to compute the complex input field with the operation of optical phase conjugation (OPC). Finally, with the calculated input field impinging on this highly scattering medium, the desired Lommel beam can be produced at the output plane, as shown in Fig. 1(c).

Figure 2 sketches the experimental setup. To achieve a rapid wavefront shaping, we utilized a high-speed DMD (Vialux V-7001) as a spatial light modulator, which can switch at a rate of 22.727 kHz. A laser beam with $\lambda =532\mathrm{nm}$ (Cobolt 04-01 Series) was expanded by a telescope constituted by L1 and L2 and reflected by M1 to fully illuminate the surface of the DMD. Then, with the assistance of a $4f$ configuration and a spatial filter, the superpixel method^[15,33] enables the DMD to shape the complex amplitude of the desired light in its first-diffraction-order beam. Each pixel of the DMD is a micromirror, which can be controlled to turn on or off independently. In the superpixel method, the square regions of nearby micromirrors are grouped into superpixels. The key point of this method is to make the phase prefactors of the target plane response of each micromirror within a superpixel distribute uniformly over a circle in the complex plane with a fixed phase step. With this method, we can modulate different target fields by turning on different combinations of micromirrors in a superpixel. The modulated beam then impinged on a ZnO scattering layer, whose thickness is about 280 µm via an objective lens OBJ1 ($10\times$, $\mathrm{NA}=0.25$). Note that the ZnO scattering layer was employed as highly scattering medium. After that, another objective OBJ2 ($20\times$, $\mathrm{NA}=0.4$) was used to collect the transmitted light field and a CMOS camera (D752, PixeLINK) was employed to image the light field. The distance between the ZnO scattering layer and the focal plane was about 200 µm. For the TM measurement, the phase shifting method^[32] was employed, and a plane wave was introduced as a reference beam. In the measurement of TM, each of the input modes was turned on sequentially with the rest turned off, and its corresponding transmitted speckle field was measured from four interferometric measurements where the signal light was phase-shifted by 0, $\pi /2$, $\pi$, and $3\pi /2$. With this procedure, each column of the TM was calibrated in sequence. In experiment, $N=32\times 32$ segments on the DMD and $M=480\times 480$ pixels on the CMOS camera were, respectively, used as the input and output modes. In order to obtain a high signal-to-noise ratio, the Hadamard basis was adopted.

In order to verify the validity of our method, we first constructed a Lommel beam with $n=2$, $c=0.7$ through ZnO scattering layer. The intensity and phase profiles of the desired beam are shown in Figs. 3(a) and 3(b). Its corresponding intensity and phase profiles in the angular spectrum are shown in Figs. 3(c) and 3(d). Note that the outer and inner radii of the annular mask were ${\rho }_1=240$ and ${\rho }_2=40$ pixels of the output plane in order to increase the SBR at the focal plane. Employing this angular spectrum as a filtering mask in the TM-based PSF engineering method, we were able to generate the desired Lommel beam through highly scattering media. Figures 3(e) and 3(f) are the intensity and phase profiles of the generated Lommel beam in simulation. In comparison, the corresponding profiles of the experimentally generated Lommel beam behind the ZnO scattering layer are, respectively, illustrated in Figs. 3(g) and 3(h). In order to make a further comparison among the theoretical distribution, simulation result, and experimental result, their intensity profiles along the white dashed line and the blue dashed line are plotted in Figs. 3(i) and 3(j), respectively. It can be observed that both of the experimental and simulated results agree well with the theoretical distribution. In this case, the Lommel beam was successfully constructed through highly scattering media with our proposed method.

Further, we demonstrated that the distribution of generated Lommel beams can be tailored flexibly by designing their corresponding angular spectrum with appropriate parameters in the TM-based PSF engineering method. In Figs. 4(a)–4(d), it is clearly observed that the constructed Lommel beams tend to be more stretched along the $y$ axis with the increase of $c_0$ ($c_0=0.1$, 0.4, 0.7, 0.9) when other parameters remain the same ($n=2$, ${\phi }_0=0$). Figures 4(e)–4(h) are the corresponding phase profiles. Moreover, the size of the Lommel beams increases with the increasing topological charges. Figures 4(i)–4(l) are, respectively, the intensity patterns of Lommel beams with different topological charges ($n=1$, 2, 3, 4) when other parameters are the same ($c_0=0.7$, ${\phi }_0=0$). The corresponding phase patterns are shown in Figs. 4(m)–4(p). In addition, we investigated the control of the orientation of the Lommel beams generated through the ZnO scattering layer. Figure 5 illustrates the intensity and phase profiles of the generated Lommel beams with different parameters ${\phi }_0$ (${\phi }_0=0$, $\pi /4$, $\pi /2$, $3\pi /4$) and the same parameters ($c_0=0.6$; $n=2$). All of these experimental results conform to the theoretical distribution.

Apart from shaping a single beam at a fixed position, two-dimensional raster scanning of the beam could benefit optical manipulation. Based on the TM method and the fast switching ability of DMD, we were able to achieve the rapid scanning of the generated Lommel beams through the ZnO scattering layer. The corresponding experimental results are presented in Figs. 6(a)–6(f). Note that the parameters of the Lommel beams were $c_0=0.7$, $n=1$, and ${\phi }_0=0$. What is more, except for constructing a single Lommel beam [Fig. 6(g)], we demonstrated that we can also generate multiple Lommel beams simultaneously by setting multiple focusing locations in applying the OPC of the filtered TM. As shown in Figs. 6(h) and 6(i), two and four Lommel beams are produced through the ZnO scattering layer at the same time. The ability of constructing multiple Lommel beams simultaneously could benefit the parallel trapping and manipulation of a number of microparticles through the highly scattering media.

In summary, we have experimentally constructed various Lommel beams through highly scattering media by applying an angular spectrum with appropriate parameters as a filtering mask in the TM-based PSF engineering method. All of the established Lommel beams match well with the theoretical predictions. As expected, the field distributions and OAM were engineered continuously by adjusting the beam parameters in the angular spectrum. In addition, the rapid raster scanning of Lommel beams through the ZnO scattering layer was demonstrated by employing the fast switching ability of the DMD. Moreover, the simultaneous construction of multiple Lommel beams was realized through the ZnO scattering layer. The method can also be extended to the other complex media, such as biological tissues and multimode fibers. We believe that this work will benefit the applications of Lommel beams behind the highly scattering media.