Sculpting light freely into user-defined shapes is a central key aspect in many areas of photonics. The list of fields with a strong interest in beam shaping is extensive and ranges from laser material processing [1–4], mode conversion [5], near-eye displays [6], lithography, telecommunications, and practical lighting systems to fundamental research in biology and physics, such as investigations towards revealing the functioning of the brain via patterned photostimulation [7], or the trapping of microparticles [8] and atoms [9]. Complex wavefronts can be shaped by static components such as diffractive optical elements (DOEs), freeform elements [10], and metasurfaces [11], or dynamic devices such as spatial light modulators (SLM) based on liquid crystal, microelectromechanical systems (MEMS) technology, or even soundwaves [12].

Despite the prowess of today’s beam shaping technology, there remain important challenges. One of them is tailoring fields with a lack of coherence, i.e., partially coherent or incoherent light, which is of central importance in many lighting and high power laser applications. In particular configurations, where light sources satisfy the geometric optical regime, some solutions based on ray-optical approaches have proven to be efficient [13,14]. The main advantage of such approaches is that they enforce the design of surfaces with curl-free gradients that are perfectly adapted to freeform elements. However, these methods inherently neglect diffraction and are therefore ill-suited for applications involving miniaturized integrated optics [15] or multicolor beam shaping [6,16,17], as well as for applications requiring more complex partially coherent sources.

Usual representations of partially coherent light introduce four-dimensional second-order correlation functions [18]. Such functions admit, owing to Mercer’s theorem, a mode expansion consisting of spatially coherent and mutually uncorrelated fields. Although computing such mode expansions is in general a difficult task, as it involves solving the eigenvalues of a Fredholm integral equation, some ubiquitous light sources such as multimode lasers or fibers already exhibit a basis of distinct eigenmodes. Moreover, it is often possible to find some phenomenological approximate source models that provide convenient analytical modal descriptions, for instance, a set of spherical waves, which we use here to model a light-emitting diode (LED) source.

In this letter, we present a novel approach for designing multiple phase DOEs to shape user-defined intensity distributions of partially coherent light, which is described by a set of monochromatic, spatially coherent but mutually uncorrelated modes. Exemplary light sources of this category are narrow-band LEDs or sources employing multimode fibers for light delivery.

Our method is enabled by the fact that each mode has unique intensity distributions on each DOE, which results in a unique light transformation for each mode. Our algorithm takes these differences into account and finds optimal phase patterns so that the integrated intensity of all modes at the output plane appears in the desired form. A set of DOEs can also be considered as a simple realization of a volume hologram, which is known to perform different tasks depending on the input mode. Conceptually similar approaches have recently been presented, but using different algorithms and approximately two orders of magnitude fewer input modes [17,19].

Our approach employs numerical wave propagation and is not limited by any geometric constraints such as rotational symmetry. It is, however, based on the thin element approximation (TEA) which must hold for the experimental realization. For the sake of generality, we provide an algorithm that accounts for an arbitrary number of DOEs, although we show that two phase patterns are sufficient in all the situations we considered, to obtain nearly perfect intensity shaping. Moreover, we experimentally demonstrate the concept by redistributing a part of the light from an LED into various shapes using two diffractive patterns displayed on a liquid crystal SLM. Given the widespread use of partially coherent sources in science and technology, we anticipate that our work will be helpful for the advancement of research fields concerned with tailored light fields.

The algorithm we developed consists of an optimization procedure relying on gradient descent to minimize the error between the time-averaged transverse intensity profile $I_S$ of the partially coherent source after crossing the optical system and a user-defined target intensity profile $I_T$. Because we assume that the modes defining the partially coherent source are mutually uncorrelated, their time-averaged intensity $I_S$ is simply the sum of the contributions of all the individual modes constituting the source. The optimization procedure is iterative and consists of three stages. First, a forward pass propagates $N$ input modes $U_0[1:N]$ from the partially coherent source plane to the detection plane where the desired intensity profile $I_T$ is to be detected, through the optical system consisting of $P$ axially separated phase patterns $\Phi [1:P]$. A pseudo-code description of the forward model is presented in algorithm stage 1. Second, the time-averaged input intensity $I_S$ is computed in the detection plane and serves to define the error $err$ based on the squared $l_2$ distance between $I_S$ and $I_T$. Because the goal of the main algorithm is to minimize this error with respect to the phase parameters $\Phi$, the gradient $\nabla _\Phi$ of this error has to be computed. To this end, we develop a backward pass, as is commonly done in computer science to numerically compute the gradients of complicated neural network architectures with gradient backpropagation algorithms. Here we used a few rules from the set defined in [20], leading to the error and gradient computation algorithm stage 2. Third, we need to update the phase patterns $\Phi$. Because we managed to compute a gradient, any first-order optimization method can be employed, the simplest being adding a small gradient step of opposite direction, $\Phi \gets \Phi - \alpha \nabla _\Phi$ with $\alpha > 0$. These three stages need to be repeated iteratively until a convergence criterion regarding the error is achieved, $err < \epsilon$ with $\epsilon > 0$.

In algorithm stage 1, the inputs of the propagate_forward procedure are $U_0$, a list of $N$ transverse modes mutually uncorrelated constituting the partially coherent source, and $\Phi$, a list of $P$ phase patterns initially set to flat phases. The outputs are $U$, a list of modes propagated to the detection plane, and $S$, a list of modes stored just after crossing each phase pattern which will be essential to compute the backpropagated gradient later on. More precisely, $S$ is a list of a list, the outer list being of length $P$ and storing inner lists of modes of length $N$. The data structures present in the algorithms are defined to simplify the notations and provide conciseness, but do not prejudge how they should be implemented in practice to get the best performance. To avoid verbose loops on the mode index, we present the code in a vectorized form by appending a dot ($.$) after the functions that we want to broadcast over the different elements of a list type parameter. Similarly, the vectorized point-wise multiplication of a collection of modes with another collection of modes of the same length or a single phase pattern is defined by a dot operator ($\cdot$). To describe mode propagation from the current plane to the $p$th phase pattern, we introduce a generic prop function, taking both a tranverse mode and an index $p$ as parameters, which can be defined according to the user’s choice. A simple model that we also use in the current study is the angular spectrum (AS) method for free-space propagation and fast Fourier transform (FFT) for Fourier collimation or focusing. In the forward model description 1 as well as in our following experimental demonstration, we choose the detection plane to be Fourier-conjugated to the last phase pattern of the optical system, but this is again arbitrary, up to the user’s choice.

The algorithm stage 2 describes the error and gradients computation which are returned by the backpropagate_gradient procedure. The list of propagated fields $U$ returned by the forward propagation procedure is used to compute the error $err$ and is backpropagated in a mirrored fashion through the same optical system. To achieve gradient backpropagation, all the propagation operators need to be replaced by their adjoint (Hermitian conjugate) operators denoted by an asterisk ($^{*}$).

We experimentally demonstrate the performance and accuracy of our beam shaping algorithm by implementing the setup shown in Fig. 1. A single SLM (Hamamatsu X10468, $16$ mm $\times$ $12$ mm) and a mirror (square with side length of 5 mm) form a tandem configuration of two diffractive patterns. Similar layouts have been used in the past for shaping coherent light and approaching properties of volume holograms, such as wavelength and angular selectivity [17,19,21–25]. Our approach is conceptually related, however, our algorithm is different and the number of modes it manages to handle to accurately describe partially coherent light is substantially higher compared with the previous state of the art (up to 1600 individual coherent modes shown in this letter), which is an essential prerequisite for realizing practical application cases.

 

Fig. 1. Folded setup for holographic beam shaping with two reflections on a single SLM.

Download Full Size | PDF

In this proof of concept study, we employ a green LED (Thorlabs M530L4) as a primary source, with a peak emission wavelength around $522$ nm and bandwidth of $35$ nm. The emitter is square-shaped with a side length of 1 mm with a plastic dome affixed on the top that modifies the spatial emission characteristics. Following a pragmatic approach, we incorporate the optical effects of the plastic dome by modeling the LED according to its image through a 1:1 Keplerian telescope (NA = 0.25), where the emitter appears to be roughly circular with a radius of $r_0 = 900$ µm, exhibiting several dark spots arranged in a square pattern. In our algorithm, we use this image as a replacement for the true LED plus the dome lens. The LED light is filtered by a pinhole of radius $r_{pin} = 75$ µm, which is positioned at a distance $z_{pin}$ from the LED such that the numerical aperture (NA) of the light passing through it is approximately $\textrm {NA} \simeq 0.1$. After collimation by a first lens (focal length $f_1$ = 15.3 mm), the beam radius is measured to be $r_{col} = 1625$ µm. After collimation, the light traverses two diffractive patterns, $\Phi _1$ and $\Phi _2$, at distances of 165 mm and 325 mm from the collimating lens, respectively. Finally, it is directed onto a camera by a second lens ($f_2$ = 15 cm) at a distance $f_2$ from the sensor. Therefore, the camera captures a magnified image of the aperture ($\approx \times 10$) if no diffractive patterns are shown. A linear polarizer is placed in front of the camera, to discard light polarized in a direction insensitive to the SLM.

In the simulation, we treat the LED as a monochromatic, spatially incoherent light source [26]. We discretize the LED emitter model into a square grid of size $N_S\times N_S$, which determines the number of uncorrelated modes representing our source. The pixel size $dx_S = dy_S$ in each dimension is chosen such that the total grid side length is $2r_0$, and the source profile is idealized to a disc of radius $r_0$ with constant intensity. Then, each pixel acts as an independent spherical wave source with an amplitude proportional to the square root of the local intensity, $\sqrt {I}e^{i k r}/r$, that propagates to the pinhole plane where it is clipped by the 150 µm diameter aperture. Thus, the analytical expression of each transverse mode is known in the pinhole plane and serves to define $N_S^2$ input distributions stored in grids of size $512\times 512$ with a pixel size $dx_{pin} = dy_{pin} = 1$ µm. The collimation of the input modes by the lens following the pinhole is modeled via a scaled FFT (Chirp Z-transform), and the resolution after collimation is set to $dx = dy = 20$ µm, which directly matches the SLM pixel size. Obviously, given our disc shape model of the LED, some of the generated modes have zero intensity and could be removed for optimizing memory and computation time (the number of effectively used modes is approximately $\pi N_S^2/4$).

The results of the beam shaping algorithm with two phase patterns are displayed in Fig. 2. We chose $N_S=40$ corresponding to 1600 modes. To demonstrate the flexibility of the algorithm, the LED emission is sculpted into three different shapes: a square with side length of 1.5 mm, a Gaussian ring with radius of $675$ µm and waistof $150$ µm, and a smiley containing fine details with an outer diameter of 2.3 mm. The simulations show an excellent agreement with the theoretical target intensity profiles, which are obtained after a few hundreds of iterations associated with an error bound $\epsilon = 1.10^{-8}$, for targets of total intensity normalized to unity. The typical computation time for such precision is approximately 10 minutes on a mid-range GPU (Geforce RTX 2070 Super). The experimental results obtained with the setup shown in Fig. 1 are likewise excellently matching the theoretical targets, even if we sometimes observe some weak modulations reminiscent of the LED substructure that we did not fully take into account. Moreover, we experimentally observe that our multimode holographic transformations are much less sensitive to misalignment than similar shaping of coherent light, where noticeable speckle usually appears even for small misalignments.

 

Fig. 2. Beam shaping results obtained with two phase patterns for three different target intensity profiles (square, ring, smiley). The first column represents the target intensity profiles set as inputs of the simulation. The second column shows the simulated intensity profiles, and the third column shows the intensity profiles recorded on a camera when the phase patterns computed by the simulation are displayed on the SLM.

Download Full Size | PDF

We would like to emphasize that the high reconstruction quality shown is only possible owing to the mode-dependent readout of the DOE tandem. Each mode is shaped differently to contribute optimally to the target pattern. For more information to better understand this concept, please refer to Supplement 1 and  Visualization 1.

Our simulations and experiments demonstrate that two phase patterns are sufficient to obtain nearly perfect intensity shaping of an incoherent source such as an LED. We performed further simulations considering different mode bases suitable for describing several other partially incoherent light sources, such as Hermite–Gaussian (HG) or linearly polarized (LP) modes, obtaining similar quality levels for all cases (data not shown).

Next we show that these complicated beam transformations cannot be achieved by a single phase pattern. For this purpose, we fix $\Phi _1$ to a flat phase and run an optimization procedure relying only on $\Phi _2$ to achieve the same transformations as above. The results displayed in Fig. 3 confirm that a one phase pattern is not sufficient to ensure a good convergence of the algorithm to the defined targets. Notably, we still observe a good agreement between simulations and experimental results.

 

Fig. 3. Results of the beam shaping algorithm with one phase pattern for the same targets as in Fig. 2.

Download Full Size | PDF

Finally, we discuss the impact of the number of modes used for source modeling on the reconstruction quality. Regardless of the number of modes, the simulation always converges almost perfectly to the given target intensity, but severely undercutting the number of modes required to accurately describe the light source significantly degrades the result. Figure 4 visualizes the influence of the mode number on the experimental reconstructions. An upper bound for the number of modes required to represent our light source can be estimated using Shannon’s sampling theorem. The maximally allowed spatial sampling distance for any object imaged at a given NA is $dx_{min} = \lambda /(2\textrm {NA})$. We can determine the maximum number of modes to be $n_{max} = \pi r_{pin}^2/dx_{min}^2$, which is close to $n_{max} = 2600$ in our case. For our three profiles, we managed to obtain a high fidelity with only approximately 1257 modes, but finer details or sharper edges in the target profile may require an even higher number.

 

Fig. 4. Experimental intensity profiles for varying numbers of modes used for modeling the input beam.

Download Full Size | PDF

In conclusion, we presented a general algorithm designed to shape the intensity of partially coherent light, given its expansion in terms of coherent and mutually uncorrelated modes. We demonstrated the efficiency of this algorithm and evaluated it experimentally with a simple optical system consisting of two phase-only diffractive patterns, which appears sufficient to achieve a very good fidelity. Given the performances of the algorithm, we can envision extending this work to more complex situations such as the combination of several partially coherent beams, eventually with multiple colors, by adapting this algorithm to the design of multilayer or volumic DOEs.