1. Introduction

The far-field light-shaping problem is aimed at redistributing the energy of the light source for certain applications. Regarding the design of the optical element to be used for the shaping, mainly two strategies are available in literature. One is a structure-oriented approach, where the optical element structure is developed based on geometric relations, e.g, the freeform surface design algorithms [1–5]. The other considers the functionality of the required optical element first, which is mainly a phase-only response function, and based on that a structure is designed as a second step to realize its functionality, with certain constraints. The latter approach is usually favored in diffractive optics, especially for paraxial beam shaping [6]. Typically, once the element function is obtained, the thin element approximation (TEA) method is applied to design the element structure [7–9]. Here, we follow the latter strategy, and mainly consider the functional embodiment of the required element in this article. However, we generalize the approach and do not restrict ourselves to the paraxial situation.

The functional embodiment can be concluded by an inverse thinking of the field tracing techniques, which shows that synthesizing a Fourier pair of the electromagnetic fields is the main requirement for far-field light shaping. The Iterative Fourier Transform Algorithm (IFTA) is one of the successful numerical methods [10–13] . With IFTA, a required output phase is obtained in an iterative Fourier transform process, and the phase response function for the element is nothing else than the subtraction of the required output phase from the input phase, which can be calculated from the given source. However, IFTA is not a general approach. A mathematical proof given later in Section 3 demonstrates that in certain shaping systems, when the fields of the Fourier pair have a homeomorphic relationship, IFTA will exhibit a stagnation problem.

To tackle the above-mentioned problem, we introduce a mapping-type Fourier pair synthesis approach. Instead of an iterative process, the approach results in a required output phase in a single step, without need of iteration. Moreover, the output phase obtained by the approach is a smooth function, from which a smooth element function or even a smooth element structure can be designed, thus reduction of speckles [14,15] or higher conversion efficiency [16] may result. With a non-paraxial example, the method is demonstrated in Section 4. In Section 5, a simulation with physical optics is performed to investigate the result. It reveals high accuracy for the non-paraxial system. Finally, in section 6, we show that even though the homeomorphic situation cannot be fulfilled for certain systems, the result from that approach is a reasonable initial guess for a subsequent IFTA.

2. Field tracing in light-shaping systems

We first analyze the light-shaping problem under the framework of field tracing [17]. Together with the modeling techniques, a light-shaping system is illustrated abstractly for the following discussion in Fig. 1. In light-shaping tasks, the known information is usually the source field and the signal (typically the target energy distribution defined by electric-field-related quantities, e.g. irradiance, radiant intensity, etc). The optical element to be designed shapes the source field to achieve the target, as shown in Fig. 1(a). A field tracing diagram (Fig. 1(b)) indicates the modeling methods used for the system, which are sequentially applied in the light path. The field solvers for each part of the system are applied in either the spatial domain ($\boldsymbol {\rho }$ domain) or the spatial-frequency domain ($\boldsymbol {\kappa }$ domain).

 

Fig. 1. A light-shaping system with a field tracing diagram illustrating the modeling techniques.

Download Full Size | PDF

In Fig. 1(b), $\boldsymbol {E}(\boldsymbol {\rho })$ represents the vectorial electric field, i.e. $\boldsymbol {E}(\boldsymbol {\rho })=(E_x(\boldsymbol {\rho }), E_y(\boldsymbol {\rho }), E_z(\boldsymbol {\rho }))$, defined in the $\boldsymbol {\rho }$ domain with $\boldsymbol {\rho }=(x, y)$. $\boldsymbol {E}^{\textrm {in}}(\boldsymbol {\rho })$ and $\boldsymbol {E}^{\textrm {out}}(\boldsymbol {\rho })$ are fields defined on the plane in front of the optical element, and the plane behind it; $\boldsymbol {E}^{\textrm {out}}(\boldsymbol {\rho }')$ is the field on the target plane. ${\tilde {\boldsymbol {E}}}(\boldsymbol {\kappa })$ are the fields defined in the $\boldsymbol {\kappa }$ domain, with $\boldsymbol {\kappa }=(k_x, k_y)$, the $x$ and $y$ component of the wave vector. ${\boldsymbol {\mathcal {B}}}$ is the operator used to model the optical element, which acts on the field in the $\boldsymbol {\rho }$ domain. For the system in this article, we assume it is a phase-only response function. ${\tilde {\boldsymbol {\mathcal {P}}}}$ denotes the free space propagation operator in the $\boldsymbol {\kappa }$ domain. And ${\boldsymbol {\mathcal {F}}}$ and ${\boldsymbol {\mathcal {F}}}^{-1}$ are the direct and inverse Fourier transform operators. Therefore, the field tracing diagram builds up an algorithm and demonstrates how the target field is calculated from the input,

3. Fourier pair synthesis

In general, attempting to obtain $\boldsymbol {E}^{\textrm {out}}(\boldsymbol {\rho })$ from the given signal is not straightforward, due to the lack of full field information about the target or sometimes the constraints of a certain design task. In a functional-embodiment design, the element function is assumed to be a phase-only response function, while the amplitude modulation will be taken into account later in the structure design step. Therefore, the amplitude of $\boldsymbol {E}^{\textrm {out}}(\boldsymbol {\rho })$ can be assumed equal to that of the input field $\boldsymbol {E}^{\textrm {in}}(\boldsymbol {\rho })$, and its phase is to be determined with the constraint from the signal, or the required field ${\tilde {\boldsymbol {E}}}^{\textrm {out}}(\boldsymbol {\kappa })$ in the $\boldsymbol {\kappa }$ domain which is concluded from the signal. Therefore, the Fourier pair synthesis between $\boldsymbol {E}^{\textrm {out}}(\boldsymbol {\rho })$ and ${\tilde {\boldsymbol {E}}}^{\textrm {out}}(\boldsymbol {\kappa })$ is the critical point for the design.

3.1 Iterative Fourier transform algorithm (IFTA)

The IFTA is one of the algorithms available for the Fourier pair synthesis, which numerically searches for the required $\boldsymbol {E}^{\textrm {out}}(\boldsymbol {\rho })$, and satisfies the merit functions for the field ${\tilde {\boldsymbol {E}}}^{\textrm {out}}(\boldsymbol {\kappa })$ [10,11]. In practice, the IFTA works on a single component of the field. Here, we use the notation $E_\ell$ to represent any component of the field $\boldsymbol {E}$, with $\ell =x, y, z$.

The iterative process of IFTA repeats the Fourier transform between each component of the fields ${E}_\ell ^{\textrm {out}}(\boldsymbol {\rho })$ and ${\tilde {E}}_\ell ^{\textrm {out}}(\boldsymbol {\kappa })$, along with the amplitude constraints $C_{\boldsymbol {\rho }}$ in the $\boldsymbol {\rho }$ domain and $C_{\boldsymbol {\kappa }}$ in the $\boldsymbol {\kappa }$ domain respectively. The amplitude constraint $C_{\boldsymbol {\rho }}$ replaces the amplitude of ${E}_\ell ^{\textrm {out}}(\boldsymbol {\rho })$ with that of ${E}_\ell ^{\textrm {in}}(\boldsymbol {\rho })$ in each iteration. While the amplitude constraint $C_{\boldsymbol {\kappa }}$ is to replace the calculated amplitude of the field in the $\boldsymbol {\kappa }$ domain with the amplitude of ${\tilde {E}}_\ell ^{\textrm {out}}(\boldsymbol {\kappa })$, which is concluded from ${E}_\ell ^{\textrm {out}}(\boldsymbol {\rho }')$ via an inverse modeling of the free space propagation. The algorithm retrieves the required output phase for ${E}_\ell ^{\textrm {out}}(\boldsymbol {\rho })$ that fulfills these constraints, or at least does so approximately.

In general, the IFTA uses the Fast Fourier Transform (FFT) method for performing the Fourier transform, that implies each field value of ${\tilde {E}}_\ell ^{\textrm {out}}(\boldsymbol {\kappa })$ is contributed to by all the field values from ${E}_\ell ^{\textrm {out}}(\boldsymbol {\rho })$. However, in certain cases, in which a homeomorphism exists between the fields ${E}_\ell ^{\textrm {out}}(\boldsymbol {\rho })$ and ${\tilde {E}}_\ell ^{\textrm {out}}(\boldsymbol {\kappa })$, or in other words, each field value of ${\tilde {E}}_\ell ^{\textrm {out}}(\boldsymbol {\kappa })$ is contributed by a single corresponding field value of ${E}_\ell ^{\textrm {out}}(\boldsymbol {\rho })$, the Fourier transform can be accurately performed by the homeomorphic Fourier transform (HFT) [18], whereas IFTA will fail to converge.

For example, let us consider the case of a non-paraxial system where the homeomorphic assumption is valid. The phase part of the field $\boldsymbol {E}^{\textrm {out}}(\boldsymbol {\rho })$ is only a strong smooth phase, which is a common phase for each field component. The fields of the Fourier pair are written explicitly as:

According to the theory of the HFT [18], the Fourier transform results in a relationship between the wavefront phase in the two domains,

3.2 Mapping type Fourier pair synthesis

On the other hand, if the Fourier pair exhibits a homeomorphic behavior, the stationary phase formula

Finding an appropriate mapping is a critical point. However, with the amplitudes in both domains known, the mapping can be solved by the well-known Parseval’s equation,

If we assume $\boldsymbol {\kappa }(\boldsymbol {\rho })$ is a one-to-one map, Eq. (7) can be derived to a differential equation:

The mathematical model of the $L^{2}$ Monge-Kantorovich problem, or the so-called “Optimal Mass Transport” (OMT) problem is applied here to solve $\boldsymbol {\kappa }(\boldsymbol {\rho })$ from Eq. (8). The curl-free characterization of the solution of the OMT problem also guarantees the mapping $\boldsymbol {\kappa }(\boldsymbol {\rho })$ as a conservative vector field. Therefore, a smooth function ${\psi }^{\textrm {out}}(\boldsymbol {\rho })$ can be obtained by a direct integration of Eq. (6)

So far, another mathematical model, a mapping-type Fourier pair synthesis is introduced for light shaping in the homeomorphic situation. The approach starts with a mapping relation between the fields in both domains and obtains a smooth common output phase for all three field components. Therefore, the mapping-type approach not only tackles the stagnation problem, but also provides a smooth function without dislocations. And in Eq. (8), one can see that the mapping-type Fourier pair synthesis method takes the full vectorial field into account.

In literature, a one-to-one map relation (or homeomorphism) is also assumed in those geometric-optics-based algorithms for light-shaping design, either for the design of holographic optical element [19–21] or freeform surfaces [2,3,22,23]. Usually, the homeomorphism is introduced for the irradiance between the input and target in the $\boldsymbol {\rho }$ domain so that local energy conservation is ensured. However, the homeomorphic assumption in the $\boldsymbol {\rho }$ domain is only valid when all the operators in the field tracing diagram (Fig. 1(b)) are pointwise operators, including the ${\boldsymbol {\mathcal {F}}}$ [24]. Therefore, we like to point out that the HFT is also included in the design logic of those algorithms.

4. Numerical example

In this section, the algorithm with the homeomorphic assumption for the Fourier pair synthesis is demonstrated with an example. A simple light-shaping task is shown in Fig. 1. We attempt to design a smooth phase response function, modulating the input light and achieving a required energy distribution at the far-field zone. The source is a Gaussian wave, with a wavelength of 532 nm and half divergence angle of $1^{\circ }$. The target pattern is the portrait of Mr. Ernst Abbe, which is given as an irradiance distribution with a size of $1\times 1\;\textrm {m}$ (shown in Fig. 2(a)). The irradiance of the pattern is scaled for better visualization. For a straightforward discussion, the target plane is set perpendicular to the optical axis and located $1\;\textrm {m}$ away from the component.

 

Fig. 2. An example for the discussion of the design algorithm. (a) Target irradiance distribution (normalized). The squared norm of the amplitude for the Fourier pair: (b) $\left \| \boldsymbol {E}^{\textrm {out}}(\boldsymbol {\rho }) \right \|^{2}$ in $\boldsymbol {\rho }$ domain , and (c) $\left \| {\tilde {\boldsymbol {E}}}^{\textrm {out}}(\boldsymbol {\kappa }) \right \|^{2}$ in $\boldsymbol {\kappa }$ domain. (d)(e) The designed homeomorphism between the squared norms of the amplitude in both domains, with the mesh (d) for (b) and its mapping one (e) for (c).

Download Full Size | PDF

For the Fourier pair synthesis, the amplitudes of both the field $\boldsymbol {E}^{\textrm {out}}(\boldsymbol {\rho })$ and ${\tilde {\boldsymbol {E}}}^{\textrm {out}}(\boldsymbol {\kappa })$ are prepared first. With the amplitude constraint in the $\boldsymbol {\rho }$ domain, the amplitude of $\boldsymbol {E}^{\textrm {out}}(\boldsymbol {\rho })$ is equivalent to the amplitude of $\boldsymbol {E}^{\textrm {in}}(\boldsymbol {\rho })$ that is obtained by propagating the source field to the component plane. On the other hand, by assuming a polarization state, for instance linearly polarized along $x$, the target field $\boldsymbol {E}^{\textrm {tar}}(\boldsymbol {\rho }^{\prime })$ can be calculated from the given irradiance distribution. $\boldsymbol {E}^{\textrm {tar}}(\boldsymbol {\rho }^{\prime })$ exhibits a spherical phase function since it is located in the far-field zone. The propagation step is modeled inversely to obtain the field ${\tilde {\boldsymbol {E}}}^{\textrm {out}}(\boldsymbol {\kappa })$. The mapping-type method works with the squared norm of both the field $\boldsymbol {E}^{\textrm {out}}(\boldsymbol {\rho })$ and ${\tilde {\boldsymbol {E}}}^{\textrm {out}}(\boldsymbol {\kappa })$, therefore we show their squared norm values in Fig. 2(b) and (c), respectively.

In order to calculate $\boldsymbol {\kappa }(\boldsymbol {\rho })$ from Eq. (8) with the data in Fig. 2(b) and (c), we modify the numerical algorithm proposed by Prins [25] to solve the equation, that in every iteration of the algorithm, an integration process of the mapping function is introduced. Then we differentiate this obtained potential function of the mapping as the input for the next iteration. In this way, the algorithm converge to an integrable mapping fast. Besides, a dummy square-shape uniform density function $D({\boldsymbol {\xi }})$ is introduced, so that the integration over its defined area is equal to both sides of Eq. (7). Instead of calculating $\boldsymbol {\kappa }(\boldsymbol {\rho })$ between $\left \| \boldsymbol {E}^{\textrm {out}}(\boldsymbol {\rho }) \right \|^{2}$ and $\left \| {\tilde {\boldsymbol {E}}}^{\textrm {out}}(\boldsymbol {\kappa }) \right \|^{2}$ directly, we separately solve two mappings, $\boldsymbol {\kappa }({\boldsymbol {\xi }})$ and $\boldsymbol {\rho }({\boldsymbol {\xi }})$, between the dummy function and each squared norm of the amplitude. By doing so, the mapping $\boldsymbol {\kappa }(\boldsymbol {\rho })$ is the combination of $\boldsymbol {\kappa }({\boldsymbol {\xi }})$ and $\boldsymbol {\rho }({\boldsymbol {\xi }})$, while the homeomorphism $\boldsymbol {\rho } \leftrightarrow \boldsymbol {\kappa }$ is still maintained. The same approach has been proposed in $\boldsymbol {\rho }$ domain light shaping algorithm [26]. The resulting homeomorphism is shown in Fig. 2(d) and (e), where the coordinates $\boldsymbol {\rho }$ and $\boldsymbol {\kappa }$ are sampled on the mesh nodes.

From the two meshes, the mapping $\boldsymbol {\kappa }(\boldsymbol {\rho })$ can be directly concluded, with both data $k_x(\boldsymbol {\rho })$ and $k_y(\boldsymbol {\rho })$ shown in Fig. 3(a) and (b) respectively. The phase $\psi ^{\textrm {out}}(\boldsymbol {\rho })$ is calculated with Eq. (6), by integrating the data of the mapping. The B-spline technique [27] is applied here for the numerical integration, and as a result, the integrated phase $\psi ^{\textrm {out}}(\boldsymbol {\rho })$ is represented by the B-spline functions, as shown in Fig. 3(c).

 

Fig. 3. The calculated mapping $\boldsymbol {\kappa }(\boldsymbol {\rho })$ with (a) $k_x(\boldsymbol {\rho })$ and (b) $k_y(\boldsymbol {\rho })$ (unit: $10^{6} \textrm {rad/m}$). (c) The required wavefront phase $\psi ^{\textrm {out}}(\boldsymbol {\rho })$ resulting from the mapping $\boldsymbol {\kappa }(\boldsymbol {\rho })$ (unit: $10^{3} \textrm {rad}$). (d) The phase response function $\Delta \psi (\boldsymbol {\rho })$ for the optical element (unit: rad). (e) The residual phase after a spherical phase function fitting from (d) (unit: rad).

Download Full Size | PDF

So far, the required phase for $\boldsymbol {E}^{\textrm {out}}(\boldsymbol {\rho })$ is obtained by the Fourier pair synthesis with the homeomorphic assumption. And the phase response function for the component is simply calculated by the subtraction of the phase of the input field from that of the output, $\Delta \psi (\boldsymbol {\rho })=\psi ^{\textrm {out}}(\boldsymbol {\rho })-\psi ^{\textrm {in}}(\boldsymbol {\rho })$. The resulting $\Delta \psi (\boldsymbol {\rho })$ is shown in Fig. 3(d), which is mainly a spherical phase function with additional modulation in Fig. 3(e). Therefore, the functionality of the required component can be interpreted as shaping the phase of the input field, giving a spherical phase to dominate the field propagation for achieving the majority of the target area, and introducing an aberrant phase to redistribute the energy after propagation.

5. Physical-optics simulation

For investigating the designed result, simulation with the calculated phase response function is done in the software VirtualLab Fusion [28]. In the field tracing engine, VirtualLab Fusion simulates the system with the techniques shown in the field tracing diagram in Fig. 1. A mathematical criterion in the background decides what kind of Fourier transform algorithm can be applied for the ${\boldsymbol {\mathcal {F}}}$. In this case, the software determines that ${\boldsymbol {\mathcal {F}}}$ is accurately calculated by the HFT. Therefore, the field behind the functional component in this example indeed exhibits a homeomorphic behavior. An irradiance detector is set in the system and the result is shown in Fig. 4(a), which is coincident with the target pattern in Fig. 2(a).

In addition, for more careful testing, we enforced the ${\boldsymbol {\mathcal {F}}}$ to be performed by the Fast Fourier Transform (FFT) method. And the result is shown in Fig. 4(b). The result in Fig. 4(a) is almost equivalent to the one in Fig. 4(b), with only slight deviations at the edge which can be neglected, again revealing that the HFT is a valid choice in this example. Consequently, the homeomorphism assumption between the field $\boldsymbol {E}^{\textrm {out}}(\boldsymbol {\rho })$ and ${\tilde {\boldsymbol {E}}}^{\textrm {out}}(\boldsymbol {\kappa })$ is checked as valid, and the application of the homeomorphic assumption in the Fourier pair synthesis attains an accurate design.

 

Fig. 4. Detector result: irradiance. Simulation with (a) HFT, (b) FFT.

Download Full Size | PDF

Besides, the smooth phase response function in Fig. 3(d) can be used further for the element structure design. For example, a computer-generated hologram (CGH) can directly be designed, since the element function of the CGH is the phase response function itself. A freeform surface may be designed as well based on the phase response function. In this non-paraxial case, TEA is not appropriate for the surface design, however, the inverse local plane interface approximation (LPIA) method may be applied [29]. The local wave vector from the input field and the gradient of the phase response function are used to calculate the normals of the surface, by using Snell’s law. Then the calculated normals are used for constructing the surface. Usually, an iterative process is necessary for achieving the surface that realizes the required phase response [30,31]. However, the design and comparison of different element structures are out of the scope of this article.

6. Discussion

For further investigation of the design algorithm, another example where the homeomorphism does not exist is also demonstrated. Here, the task is similar to the one shown above, with all the parameters in the task kept the same, except the size of the target pattern, which is reduced to $100\times 100\;\textrm {mm}$, so that the output is in a more paraxial situation. Following the same logic, the homeomorphic assumption is again made for the Fourier pair of the field behind the optical element. Applying the same algorithm in Section 4., the phase response function is designed, as shown in Fig. 5(a). Again, the algorithm results in a smooth element function. Figure 5(b) shows the residual phase in addition to a fitted spherical part. Compared with the previous example, one can see that the designed phase response function is similar except for the difference in their gradient.

 

Fig. 5. (a) The phase response function $\Delta \psi (\boldsymbol {\rho })$ for the optical element (unit: rad). (b) The residual phase after a spherical function fitting of (a) (unit: rad). (c)(d) Detector result: irradiance. Simulation with (c) HFT, (d) FFT.

Download Full Size | PDF

With the designed $\Delta \psi (\boldsymbol {\rho })$, the simulation is again performed via physical-optics modeling methods. With different Fourier transform algorithms, HFT v.s. FFT, the detector this time gives quite different results of irradiance distribution, as shown in Fig. 5(c) and (d). For the simulation with HFT, since the same assumption is made in both the design and the modeling, the result evidently coincides with the target pattern. However, by the analysis with rigorous FFT, Fig. 5(d) shows obvious diffraction effects and lower resolution compared to the target pattern, which indicates that the designed function is indeed not valid for the task. The homeomorphism assumption for the Fourier pair is far from a good approximation in this case.

In summary, we bring the homeomorphic assumption into the content of Fourier pair synthesis. With both the examples above, it is demonstrated that the critical point for the mapping-type algorithm is whether the field behind the optical element is in its homeomorphic zone, or from a modeling point of view, whether the HFT can be accurately applied for the following operator. If so, typically in non-paraxial situations where the output field has a strong wavefront phase, the mapping-type design algorithm can provide a smooth phase response function. Otherwise, the designed result is not a valid function for the light-shaping task. Its validity should be investigated by physical-optics modeling. However, when the homeomorphic behavior is broken, IFTA can be applied properly. Moreover, the IFTA can benefit from the result of the homeomorphic algorithm, which is shown in what follows.

Since in the second task the homeomorphic assumption is not fulfilled, the IFTA is appropriate for the design. We start the IFTA with an initial guess given by the result of the homeomorphic algorithm. With about 200 iterations, the IFTA converges. Figure 6 shows the phase response function, with (a) the initial one and (b) the optimized one after IFTA, with both in a $2\pi$ modulo mode. The figures only show the residual part after subtracting the same spherical phase. Figure 6(c) shows the simulation result with the optimized function, which compared to the irradiance distribution in Fig. 5(d), has improved resolution and no diffraction fringes.

 

Fig. 6. The residual phase part after a spherical phase fitting of the phase response function: (a) the initial one and (b) the optimized one (in $2\pi$ modulo, unit: rad). (c) The detector result: irradiance, simulated with the optimized phase response function.

Download Full Size | PDF

It was shown by Aagedal et al [14] that the speckle problem in the IFTA may be resolved by a well-introduced initial phase function. The speckle-free pattern as shown in Fig. 6(c) indeed indicates that the smooth phase response function resulting from the homeomorphic algorithm is a good option to start with. In such a paraxial situation, TEA can be applied on the final function after IFTA optimization for the structure design.

7. Conclusion

Shaping the far-field of light mainly requires a Fourier pair synthesis algorithm. For those systems, the field behind the shaper exhibits a homeomorphic behavior, typically in non-paraxial situations, IFTA may end up with weak convergence. In such cases, the mapping-type algorithm provides a fast and valid design. The required output phase can be directly concluded from the mapping relation between the Fourier pair. And a smooth phase response function can be obtained. For those systems without homeomorphic behavior, one can also start the design with the mapping-type Fourier pair synthesis method, but the algorithm may lead to an inaccurate function. However, the resulting phase response function is a well-suited initial guess for the IFTA. The mapping type algorithm enriches the Fourier pair synthesis class, provides new insights for light shaping, and introduces more options for the design.