Two-dimensional Talbot effect of the optical vortices and their spatial evolution

We report on the experimental and theoretical study of the near-field diffraction of optical vortices (OVs) at a two-dimensional diffraction grating. The Talbot effect for the optical vortices in the visible range is experimentally observed and the respective Talbot carpets for the optical vortices are experimentally obtained for the first time. It is shown that the spatial configuration of the light field behind the grating represents a complex three-dimensional lattice of beamlet-like optical vortices. A unit cell of the OV lattice is reconstructed using the experimental data and the spatial evolution of the beamlet intensity and phase singularities of the optical vortices is demonstrated. In addition, the self-healing effect for the optical vortices, which consists in flattening of the central dip in the annular intensity distribution, i.e., restoring the image of the object plane predicted earlier is observed. The calculated results agree well with the experimental ones. The results obtained can be used to create and optimize the 3D OV lattices for a wide range of application areas.

Structured light is a key issue of photonics, whose potential is coming to be implemented 1 . Usually, structured light represents the nontrivial intensity, polarization, and phase distributions of the coherent light, which can result in extraordinary characteristics capable of transforming the light-matter interactions. One of the brightest manifestations of structured light is the beams containing phase singularities also called optical vortices(OVs) 2 . Their feature is the orbital angular momentum (OAM) arising from the helical wavefront imposed by an azimuthal phase dependence of exp (ilϕ) [3][4][5] , where l is the topological charge (TC) of the phase singularity 6 . The orbital angular momentum of a light beam can be transferred to matter 7 and vice versa. Since the date of their discovery 2,3 , the OAM-carrying beams have found wide application in optical manipulations [8][9][10][11] , imaging 12,13 , and quantum key distribution 14 and made it possible to implement the Tbit s −1 data transmission rates in optical communications [15][16][17][18] . The recent progress in this field and future directions were reviewed in 5,19 . Among these directions are complex two-and three-dimensional (2D, 3D) lattices of optical vortices [20][21][22] . The formation of well-ordered spatial arrays of optical vortices is highly attractive, e.g., for micromanipulation in biology, medicine, and materials science 23 , two-dimensional optical trapping (creating arrays of microscale sculpted atom optical states in two and three dimensions) 1 and sorting 24 , photolithography 25,26 , micromachining and structuring, etc. 27 The most obvious approach to creating a spatial array of optical vortices is the use of division of a wavefront of the OAM-carrying beams. In this connection, it is reasonable to examine the near-field diffraction of the OAM-carrying beams at a 2D amplitude diffraction grating. It should be noted that the diffraction of the OAM-carrying beams was studied for the cases of a single slit 28,29 , sector-shaped diaphragm 30 , and 1D gratings 31 . It is well-known that a plane wave incident onto a periodic grating is diffracted such that the Talbot effect occurs in the near-field 32 . The physics behind this phenomenon is interference of diffracted waves, which results in the intensity distribution (image) self-reproduction at distances multiple of Talbot length Z T = 2� 2 / , where is the grating period, is the wavelength of the incident wave. This phenomenon is well-studied, both theoretically and experimentally, for plane waves 33 and theoretically investigated for the beams comprising phase singularities 20,34,35 . In the experiment, the Talbot effect of the OAM-carrying beams was observed in the 1D case 31 , as well as for the two-dimensional configuration at the THz radiation 35,36 and near-infrared single photons 37 . In the 1D case 31 , the experimental intensity patterns were in the form of a periodic set of stripes, while for two-dimensional configuration there is an array of beams with a topological charge. The Talbot effect was explored using a superposition of two optical lattices generated by a superposition of two quasi-OAM states in 38 and exploited for generating optical vortex arrays by multiplexing metasurface design 22 . Till now, no consistent Scientific Reports | (2020) 10:20315 | https://doi.org/10.1038/s41598-020-77418-y www.nature.com/scientificreports/ theoretical and experimental study of near-field diffraction of the optical vortices at a 2D amplitude diffraction grating in the visible range has been considered. Here, we report the results of observation of the two-dimensional Talbot effect of the optical vortices in the visible range and formation of three-dimensional lattices of optical vortices consisting of beamlet arrays and demonstrate their spatial evolution.

Experimental and numerical results
In the paper we study a near-field diffraction of the optical vortices at a 2D amplitude diffraction grating. In the experiments, the azimuthal phase modulation was imposed onto a Gaussian beam from a He-Ne laser using a 2D phase-only spatial light modulator (SLM, PLUTO-NIR-011, Holoeye, a pixel pitch of 8 µm ). The laser source operated in a single-mode regime and produced the linearly polarized radiation at a wavelength of 632.8 nm. The laser radiation was expanded by a 3× beam expander before launching onto the SLM. A fork-shaped hologram loaded in the SLM provided the phase modulation in the form cos(G x x + lϕ) , rather than exp(iϕ) , where G x is the reciprocal lattice vector. This allowed us to obtain the high-quality optical vortex in the first diffraction order, as shown in Fig. 1a. Then, the topological charge of the resulting vortex was measured by the method described in 39 . The beam with a specific TC was reduced and directed along the sample normal to the input facet. The amplitude 2D diffraction grating was located at the opposite side of the sample (the output facet) (See Methods for details concerning with the sample fabrication and characterization). Near-field diffraction patterns were imaged by a 100× objective and measured using a monochrome CCD camera mounted on a motorized translation stage. As can be seen in Fig. 2 (top row), the measured intensity distributions (profiles) for the TC l = +1 consist of well-ordered 2D annular beamlet arrays with the ring-shaped maxima for all the Talbot planes, which differs from 1D case, where linear array of stripes is observed 31 . These maxima suffer from  www.nature.com/scientificreports/ distortions caused by the weak intensity at the center of the incident annular beam. The shape of the maxima is different for each unit cells in the beamlet array. The unit cells are defined so that each one contains a single beamlet at its center and it approximately corresponds to a hole in the grating structure. The beamlets undergo deformations, which strengthen with distance from the center of the image plane and take a semicycle (U-shaped) form, as can be clearly seen for the forth Talbot plane. The beamlet positions coincide with the lattice points of the object structure for the integer Talbot planes, while the images corresponding to the fractional mZ T /2 planes (see, for instance, the 3/2Z T plane) are shifted by half a period relative to those corresponding to the integer planes Z T . Although this property is the same as in the classical Talbot effect, the annular beamlet arrays are observed instead of the hole images of the object structure. All the calculated results were obtained using the numerical simulation based on the Fresnel-Kirchhoff formalism described in Sec. Methods. The calculated intensity ( I(x, y, z) ∼ |E(x, y, z)| 2 ) and phase ( ψ(x, y, z) = arg(E(x, y, z)) ) profiles obtained using Eq. (13) are in good agreement with the experimental ones ( Fig. 2). In the calculation, we used a beam radius of w = 150 µm in Eq. (2) to fit the intensity distributions to the experimental ones. Analysis of the phase distributions showed phase singularities in each beamlet array cell (Fig. 2, bottom row), which are distinguished by the vortex and anti-vortex, depending on the sign acquired. Here, we use the term anti-vortex for the vortex of the opposite sign. The phase increases counter-clockwise in the vicinity of the optical vortices with a positive TC. According to 5 , the TC of the optical vortex can be found using the path integral over the closed curve C subtended by a phase surface (over the field cross section) containing the singularities Here, ψ(r, ϕ) is the phase distribution of a scalar wavefield presented in the polar coordinates. Eq. (1) allows us to calculate the TC of all optical vortices in the field cross section taking into account their values and signs, while their positions are defined by intersections of zero-level contours of both real and imaginary parts of the field. Although the optical singularities within each beamlet array cell are grouped in the vortex-anti-vortex pairs, the total TC over a single cell calculated using Eq. (1) is zero. This is not surprising and follows from the deterministic rule also known as the sign principle 40 or nucleation of optical vortices based on the properties of the wavefields topology. In addition, there exists an uncompensated optical vortex at the center of the image (the so-called global vortex), which is related to the TC of the incident beam. The beamlets become broader and smoother with an increase in the propagation coordinate z or the number of a Talbot plane. They fill the entire cell area in the fourth Talbot plane and then begin overlapping. This trend is observed for the beamlet phase profiles as well.
Similar to the case of l = +1 , the respective experiments were carried out for the incident beams with TC values of l = −4, −3, −2, −1, +2, and +3 . As was observed in the experiments and proved by our calculations, the beamlet arrays of the optical vortices were observed in all the cases. Figure 3 illustrates the results for TCs of l = −1, +1, +2, and +3 at the third Talbot plane. All the distributions exhibit the rotational symmetry and have An increase in the TC results in spreading of the beamlets and their overlap. The latter factor will result in the interaction between the beamlets. This complicates the intensity profiles of the beamlet array. It is supposed that the optical singularities will penetrate to the neighboring beamlet array cells, so they cannot be related to the given cell any more. In spite of such a behavior, the well-ordered diffraction pattern is formed in the far-field, as shown in Fig. 4. Individual diffraction orders represent the optical vortices whose TC is equal to the TC of the incident beam; in particular, in this case, we have l = +3 . All the diffractive orders have fork-shaped phase profiles, except for the zero order, which has a spiral (azimuthal) phase profile. The results described above are related to the special case of diffraction at the integer and fractional Talbot planes. Of special interest is the analysis of the 3D spatial structure of the beamlets, i.e., the configuration of the optical vortices inside the area between two neighboring integer Talbot planes. Therefore, we numerically analyzed the spatial distribution of the intensity and phase singularities of the OVs. Figure5a shows the calculated beamlet isointensity volumes between the third and fourth Talbot planes for the 2 × 2 central unit cell array for l = +1 and the respective top and side views (insets on the left). As can be seen, the beamlets undergo spatial evolution; specifically, each annular intensity beamlet is separated into four channels according to the number of neighboring lattice sites and then these channels merge with the formation of the annular intensity structures at the 3 1 2 Talbot plane placed between the initial beamlet positions in the xy plane corresponding to the third Talbot plane. After that, the annular structure is reconfigured back into the four channels and each of them contributes to the neighboring annular intensity maxima at the four Talbot planes. The spatial intensity profile of the beamlets at the fourth Talbot plane is sufficiently close to the profile corresponding to the third Talbot plane. Such a spatial behavior is typical of a rectangular lattice inside the area enclosed by any neighboring Talbot planes. To build an experimental beamlet structure, all the transverse profiles taken along the z axis were stacked in a single array. Then, the experimental beamlet structures were reconstructed using standard isovolume visualization tools. Figure5b shows the reconstructed experimental beamlets isointensity volume taken at the 1/e intensity level. Good agreement between the theoretical and experimental results proves our analytical model. The numerical analysis allows one to visualize the trajectories of OVs, which are shown in Fig. 5a (insets on the right). All optical vortices can be distinguished by the sign of the phase singularity and by its spatial behavior. The optical vortices with the positive TC nutate in the vicinity of array lattice points along the z axis, while the position of the global OV at the center of the lattice remains constant. The optical vortices with the negative A natural (intrinsic) manifestation of the Talbot effect is a spatial dependence of the transverse near-field intensity distribution previously integrated over the second transverse coordinate 32 , which is also known as a Talbot carpet. It provides valuable information on the structural parameters and the light illuminating the structure. The Talbot carpets can be easily obtained for both 1D structures and plane waves. However, an obstacle arises when we refer to the optical vortices diffracted at a 2D diffraction grating. This obstacle concerns the complex spatial intensity distribution over the unit cell of a specific image, leading to the dependence of the results on the integration limits. Nevertheless, we obtained the Talbot carpets for the configuration under study. To do that, we integrated all the intensity profiles corresponding to the coordinates z over two rows of unit cells at the image center in the y direction. The results of the integration are presented in Fig. 6 for l = +1 and l = +3 . These dependences have the structures well-ordered in both the transverse and forward direction. A number of the Talbot planes can be seen in all the cases. The calculated Talbot length ( Z T = 2� 2 / ), i.e., the distance between the neighboring Talbot planes, was 316 µm , which is slightly different from a measured value of 324 µm . The numerically calculated Talbot carpets are in good agreement with the experimental ones.

Discussion
The Talbot carpets shown in Fig. 6 manifest themselves in the self-healing effect consisting in flattening of the central dip in the intensity distribution related to the optical vortex, i.e., restoring the image of the object plane. This effect was predicted for the optical vortices in 34 . It was attributed to a transverse energy flow from the peripheral high-intensity area of the annular beam to the central null intensity area. As can be seen in Fig. 7, the intensity profiles in the first Talbot planes have the low-intensity areas at the center, but become more uniform at 4Z T . Under the experimental conditions, this effect takes place for a limited number of Talbot planes, specifically, for 3, 4, 5 Talbot planes at l = +1 and 4 − 7 Talbot planes at l = +3 . Our study showed that the self-healing effect is manifested in different ways, so the following situations are distinguished by the TC value of the incident beam: (i) if the TC is odd ( ±1, ±3, ±5, . . . ), the intensity profiles in the self-healing area represent a beamlet array of ring-shaped maxima whose positions coincide with the lattice points of the grating structure or intermediate between the lattice points; (ii) if the TC is multiple of 4 ( ±4, ±8, ±12, . . . ), the patterns represent a beamlet array of round spots of the almost uniform intensity, which resembles a perfect self-image of the object structure; in the next Talbot plane, the number of maxima is doubled and a half of them are located at the lattice points of the grating structure, while the others hold intermediate positions between them; (iii) if the TC is unevenly even ( ±2, ±6, ±10, . . . ), the intensity patterns resemble the patterns for the previous case, but they are shifted by a half-period along any single coordinate. Additionally, the measured intensity profiles in the Talbot planes under illumination by the optical vortices experience the angular rotation. In particular, in the case of illumination by the beams with an TC of l = +3 , the intensity profiles taken at 4Z T are rotated clockwise by 3 deg. The change in the TC sign will result in the change in the rotation direction. This effect is only observed for the optical vortices and predicted by the calculations.
The approach used (see Sec. Methods) can be applied to the gratings of other types, including the phase ones, by specifying the grating transmission function T(x 0 , y 0 ) entering Eq. (3).

Conclusions
In summary, we experimentally and theoretically studied the near-field diffraction of the optical vortices at a twodimensional amplitude grating. The Talbot effect for the optical vortices in the visible range was experimentally observed and the corresponding Talbot carpets for the optical vortices were observed for the first time. It was shown that the complex field distributions in the Talbot planes have the phase singularities (optical vortices), whereas the intensity surrounding these singularities has an annular structure. The total TC of these singularities within the object area in the integer and fractional Talbot planes is retained and equal to the TC of the incident beam. It was numerically shown that the phase distributions within the unit cells of a specific image at the integer Talbot planes causes a complex structure of the bounded phase singularities of the opposite signs, i.e., the vortices and anti-vortices, so the resulting TC is zero. It was shown that the spatial configuration of the light field is a complex three-dimensional lattice of the beamlet-like optical vortices. A unit cell of the lattice of the OVs was reconstructed from experimental data and spatial evolution of the beamlet intensity and phase singularities of optical vortices was demonstrated. In addition, we observed the self-healing effect for the optical vortices, which consists in flattening of the central dip in the annular intensity distribution, i.e., restoring the predicted image of the object plane. The calculated results agree well with the experimental ones. The results obtained can be used to create and optimize the 3D OV lattices for a wide range of application areas. where E 0 is the maximum amplitude, l is the topological charge and w is the effective beam radius. The parameter s = ±1 corresponds to the positive and negative TC, respectively. Further, we assume that l is a positive integer. Then, the field amplitude at the entrance of the grating is where T(x 0 , y 0 ) is the transmission function of the grating. It is assumed that the grating is sufficiently thin, so in Eq. (3) we have x ′ ≈ x 0 and y ′ ≈ y 0 . For a 2D grating, T(x 0 , y 0 ) will be a periodic function of the coordinates x and y; therefore, we have T(x 0 , y 0 ) = T(x 0 + m� x , y 0 + n� y ) , where x and y are the grating periods in the x and y directions and m and n are integers. The diffracted field amplitude in the near-field behind the grating at the distance z is given by the Fresnel-Kirchhoff integral 41 Equation (4) can be presented in the form is a periodic function of the coordinates (x, y) with the periods x and y , then it can be expanded into the Fourier series where G x = 1/� x and G y = 1/� y , t mn are the Fourier coefficients and Using the binomial distribution, the term (x 0 + iy 0 ) l in Eq. (5) can be rewritten as where C l q = l! q!(l−q)! are the binomial coefficients. Accounting for Eqs. (6) and (8) where H m are the Hermite polynomials, and the result is T(x 0 , y 0 ) = m,n t mn exp[i2π(mG x x 0 + nG y y 0 )], T(x 0 , y 0 ) exp −i2π(mG x x 0 + nG y y 0 ) dx 0 dy 0 . www.nature.com/scientificreports/ Here, b xm = 2πmG x − (k/z)x and b yn = 2πnG y − (k/z)y , a = kw 2 /2z . It can be shown that Finally, taking into account Eq. (4) the diffracted field can be represented as The intensity and phase profiles of the diffracted field in an arbitrary plane at the coordinate z can be calculated using Eq. (13).
Fabrication and characterization of the sample. The sample was a transparent quartz plate, which was polished to the optical quality and then coated with a silver film with a thickness of about 200 nm. The grating was fabricated by ion etching using a focused ion beam setup (FB-2100, Hitachi). The structure represents a 2D regular array of round holes, as shown in Fig. 1b. The structure contains 40 periods along each axis and has overall sizes of 400 × 400 µm 2 . The respective periods were 10 µm , while the hole diameter was about ∼ 5 µm , so the duty cycles of the structure were 0.5. The use of a nontransparent silver coating allowed us to achieve the high level of the amplitude modulation of the intensity, as can be seen in the optical microscopy image presented in Fig. 1c. License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.