Integrable solution for light shaping based on a Fourier-pair mapping

Abstract

In far-field light shaping, one of the design methods is based on a one-to-one map between the irradiance of the source and target. However, an integrability issue may occur in this kind of algorithms, either in the ray mapping method for designing a freeform surface or in those geometric-optics-based methods for achieving a required output phase. We introduce another mapping-type algorithm to tackle the integrability problem, which instead of establishing a mapping between both the source and target irradiance in the space domain, the mapping is assumed on electric fields of a Fourier pair between the space domain and the spatial-frequency domain. By solving the mapping from the Fourier pair, the gradient of the output phase is achieved, that the gradient is equivalent to the obtained mapping function. Moreover, the existence and the characterization of the mapping guarantees the integrability of the gradient so that a smooth output phase can be directly integrated. Based on the obtained smooth output phase, a freeform surface can then be designed for the light-shaping task. Numerical examples are demonstrated for the comparison of the approaches with different mapping assumptions.

Introduction

The design of an optical element for spatial energy redistribution is a fundamental problem of light shaping. In order to redistribute the irradiance, one straightforward thinking is to find a mapping between the irradiance distribution of the source and the target. Therefore, the local energy conservation is usually assumed in the design algorithms, either for the design of diffractive optics [1–6] or refractive optics [7–13].

The optimal mass transport (OMT) model [14–16] is widely used for solving a 2D mapping from the irradiance relation, which is usually formulated as a Monge–Ampère equation. Based on the obtained mapping, different strategies are taken for searching a proper optics to realize the mapping.

One strategy is to directly design the structure of the optical element, which is usually found in those algorithms for the freeform surface design. By using Snell’s law, the gradient of the surface is derived from the mapping and used to construct the freeform surface. This method is often referred to as the “ray mapping method” [8–10, 12, 13]. However, if the mapping is solved from the irradiance relation itself which is decoupled from the surface information, the obtained surface gradient, in general, does not satisfy the integrability condition [11, 17, 18]. Bruneton et al. [17] and Bösel et al. [18] have given mathematical proof that only in paraxial situations, the ray mapping method can provide an integrable gradient for the freeform surface. Therefore, enforcing a direct integration procedure with the gradient data usually causes errors in the designed surface, which may result in a deformed irradiance distribution from the target one. Further optimization processes should be used to improve the performance of the design [8, 10].

An alternative strategy first tries to determine a suitable output phase by using the obtained mapping, where the output phase dominates the light propagation in free space and realizes the mapping. Next, the optical element is then designed based on the obtained output phase. This strategy can usually be seen in the design of a computer-generated hologram (CGH), where the element function of the CGH is simply calculated by the subtraction between the output and input phase [3–5]. Besides, the thin element approximation (TEA) method [19] and inverse local plane interface approximation (LPIA) method [20, 21] also used to design a freeform surface with the retrieved output phase.

The derivation of the output phase from the mapping is under the stationary phase approximation [22], that the gradient of the output phase is connected with the mapping. For example, Feng et al. [19, 23] coupled the phase gradient and the mapping in a complex non-linear partial differential equation to solve the required output phase. However, again, if the mapping is concluded only from the irradiance relation without considering the output phase information, like in [3, 5, 6], the gradient data of the phase derived from the mapping is not necessarily integrable. Similar to the ray mapping method for designing the freeform surface, the integrability issue also happens here for the output phase retrieval. A mathematical proof is provided in this article that the obtained gradient data of the phase is integrable only in paraxial cases.

In order to obtain a proper output phase without solving a complex differential equation, a mapping-type Fourier pair synthesis method is introduced in the previous work of the authors [24], where instead of finding a mapping between the irradiance in the space domain, the mapping is found between the electric fields in the space domain and its corresponding one in the spatial-frequency domain.

In this work, we attempt to prove that the mapping between the Fourier pair gives an integrable solution for the required output phase, which can further be applied for the design of a CGH or a freeform surface. We present and analyze both mapping assumptions from a physical-optics point of view. The integrability condition in both algorithms is discussed mathematically. Numerical examples are shown for the demonstration and comparison of both approaches.

Problem analysis

Field tracing for the far-field light-shaping system

We analyze the far-field light-shaping system by physical optics, and illustrate the notation for the later discussion. The considered light-shaping system is shown in Fig. 1(a). The optical element is to shape the input field from the source to its far-field zone in order to achieve a desired irradiance distribution. The field tracing diagram in Fig. 1(b) illustrates the modeling algorithm for the light-shaping system (Fig. 1(a)). In the field tracing diagram, the electric fields, E=(E_x,E_y,E_z), are indicated through the light path, either in the spatial domain (ρ domain) or in the spatial-frequency domain (κ domain), with ρ=(x,y) is the transversal coordinate on a certain plane and κ=(k_x,k_y) the x- and y-component of the wave vector. $\boldsymbol {\mathcal {B}}$ is an operator that indicates the functionality of the optical element, acts on the field in the ρ domain. The free space propagation step behind the optical element contains two Fourier transform $\boldsymbol {\mathcal {F}}$ and $\boldsymbol {\mathcal {F}}^{-1}$ and a propagation operator $\tilde {\boldsymbol {\mathcal {P}}}$ in the κ domain. Therefore, the field tracing diagram builds up an algorithm that demonstrates how the target field is calculated from the input,

$$ \begin{aligned} \boldsymbol{E}^{\text{tar}}(\boldsymbol{\rho}^{\prime}) &= \boldsymbol{\mathcal{F}}^{-1}\left\lbrace \tilde{\boldsymbol{\mathcal{P}}}\boldsymbol{\mathcal{F}}\left\lbrace \boldsymbol{E}^{\text{out}}(\boldsymbol{\rho}) \right\rbrace \right\rbrace\\ &= \boldsymbol{\mathcal{F}}^{-1}\left\lbrace \tilde{\boldsymbol{\mathcal{P}}}\boldsymbol{\mathcal{F}}\left\lbrace \boldsymbol{\mathcal{B}}\boldsymbol{E}^{\text{in}}(\boldsymbol{\rho}) \right\rbrace \right\rbrace {,} \end{aligned} $$

(1)

Fig. 1

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

where Eⁱⁿ,E^out and E^tar are the electric fields defined on the input plane in front of the optical element, the output plane behind the optical element and the target plane on which the signal for the task is given.

Here, we first assume the functionality of the optical element $\boldsymbol {\mathcal {B}}$ as a phase-only response function to achieve the required output phase. The extra physical effect from the structure of the optical element will then be considered and compensated in the structure design step. Therefore, for the output phase retrieval step, it is assumed that,

$$ \left| {E}_{\ell}^{\text{out}}(\boldsymbol{\rho}) \right| = \left| {E}_{\ell}^{\text{in}}(\boldsymbol{\rho}) \right| \text{,} $$

(2)

where E_ℓ represents any component of the field E, with ℓ=x,y,z. $\left | {E}_{\ell }^{\text {in}}(\boldsymbol {\rho }) \right |$ and $\left | {E}_{\ell }^{\text {out}}(\boldsymbol {\rho }) \right |$ are the amplitudes of the input and output fields, respectively.

It is noted in [25], that if all the operators in Eq. (1) are pointwise operators, the field tracing algorithm establishes a one-to-one map or so-called homeomorphism between the input field Eⁱⁿ(ρ) and the target one E^tar(ρ^′). This is indeed the typical mapping assumption included in all the geometric-optics-based design algorithms.

In fact, due to the homeomorphism is through the whole system, the mapping can be assumed between any two fields in the field tracing diagram (Fig. 1(b)). In the mapping-type Fourier pair synthesis method, the mapping between E^out(ρ) and E^out(κ) is chosen, which is not the same as typical mapping assumption in the literature.

In the following, we compared the phase retrieval methods with these two different mapping assumptions, that one is on ρ^′(ρ) and the other one is on κ(ρ). And the integrability issue of each method is mathematically discussed.

Output phase retrieved from the mapping in the spatial domain

For the output phase retrieval based on the mapping ρ^′(ρ), the algorithm is presented as followed.

Considering a lossless system, the flux is conserved at different positions though the system. Therefore,

$$ \iint {E_{e}}^{\text{in}}(\boldsymbol{\rho}) {\,\mathrm{d}}\boldsymbol{\rho}= \iint {E_{e}}^{\text{out}}(\boldsymbol{\rho}) {\,\mathrm{d}}\boldsymbol{\rho}= \iint {E_{e}}^{\text{tar}}(\boldsymbol{\rho}^{\prime}) {\,\mathrm{d}}\boldsymbol{\rho}^{\prime} \text{,} $$

(3)

where E_eⁱⁿ(ρ),E_e^out(ρ) are the input and output irradiance distribution. Here, since the optical element is considered in its functional embodiment as a phase-only response function, E_eⁱⁿ(ρ) and E_e^out(ρ) is defined on the same reference plane of the optical element and E_e^out(ρ)=E_eⁱⁿ(ρ), where E_eⁱⁿ(ρ) is obtained by propagating the source field to the reference plane. E_e^tar(ρ^′) is the irradiance distribution on the target plane.

Typically in the literature, the homeomorphism is assumed between E_e^out(ρ) and E_e^tar(ρ^′). Therefore, Eq. (3) is derived into its differential form which is also named as the local energy conservation law [11],

$$ \text{det}[J(\boldsymbol{\rho}^{\prime}(\boldsymbol{\rho}))] = \frac{{E_{e}}^{\text{out}}(\boldsymbol{\rho})}{{E_{e}}^{\text{tar}}(\boldsymbol{\rho}^{\prime}(\boldsymbol{\rho}))} {,} $$

(4)

where det[J(ρ^′(ρ))] is the determinant of the Jacobian matrix J(ρ^′(ρ)).

In general, solving the 2D mapping ρ^′(ρ) from Eq. (4) is specified by a mathematical model, the L² Monge-Kantorovich problem, or the so-called “Optimal Mass Transport” (OMT) problem. Several numerical algorithms had been proposed to solve the OMT problem [15, 16].

After the mapping ρ^′(ρ) is solved, by the method of stationary phase according to Bryngdahl [22], the gradient of the output phase is connected with the mapping in the space domain, where their relation is written as:

$$ \nabla \psi^{\text{out}}(\boldsymbol{\rho})=k_{0}n \frac{\boldsymbol{\rho}^{\prime}(\boldsymbol{\rho})-\boldsymbol{\rho}}{\sqrt{\left\| \boldsymbol{\rho}^{\prime}(\boldsymbol{\rho})-\boldsymbol{\rho} \right\|^{2} + L^{2}}} {.} $$

(5)

Here, k₀ is the wave number, n is the refractive index in the free space and L is the propagation distance between the optical element and the target plane.

The existence and the curl-free characterization of the solution for the L² Monge-Kantorovich problem is addressed in the theorem by Brenier [26]. Therefore, a curl-free mapping function ρ^′(ρ) can be solved from Eq. (4). However, due to the nonlinear relation in Eq. (5), ∇ψ^out(ρ) is not necessary a conservative vector field, even though ρ^′(ρ) is. The integrability of ∇ψ^out(ρ) only can be preserved in the paraxial approximation that ∥ρ^′(ρ)−ρ∥≪L. In general, with the gradient data obtained from the mapping in spatial domain, the required output phase ψ^out(ρ) cannot be reconstructed by direct integration.

Output phase retrieved from the mapping in the Fourier pair

For the algorithm of the mapping-type Fourier pair synthesis [24], the homeomorphism is assumed between the field of E^out(ρ) and E^out(κ), which are connected in the Parsevel’s equation.

$$ \iint \left\| \boldsymbol{E}^{\text{out}}(\boldsymbol{\rho}) \right\|^{2} {\,\mathrm{d}}\boldsymbol{\rho}= \iint \left\| \boldsymbol{E}^{\text{out}}(\boldsymbol{\kappa}) \right\|^{2} {\,\mathrm{d}}\boldsymbol{\kappa} {.} $$

(6)

If we assume κ(ρ) is a one-to-one map, Eq. (6) can also be derived to a differential equation.

$$ \text{det}[J(\boldsymbol{\kappa}(\boldsymbol{\rho}))]= \frac{\left\| \boldsymbol{E}^{\text{out}}(\boldsymbol{\rho}) \right\|^{2}}{ \left\| \boldsymbol{E}^{\text{out}}(\boldsymbol{\kappa}(\boldsymbol{\rho})) \right\|^{2}} \text{,} $$

(7)

where det[J(κ(ρ))] is the determinant of the Jacobian matrix J(κ(ρ)).

∥E^out(ρ)∥ is concluded by using Eq. (2). E^out(κ(ρ)) can be calculated by an inverse procedure of the field tracing from E^tar(ρ^′), where E^tar(ρ^′) is usually derived from the given target irradiance E_e^tar(ρ^′), with the far-field assumption.

We denote the target field explicitly as,

$$ \boldsymbol{E}^{\text{tar}}(\boldsymbol{\rho}^{\prime}) = \left[\begin{array}{c} E_{x}^{\text{tar}} (\boldsymbol{\rho}^{\prime}) \\ E_{y}^{\text{tar}} (\boldsymbol{\rho}^{\prime}) \\ E_{z}^{\text{tar}} (\boldsymbol{\rho}^{\prime}) \end{array}\right] \exp[\mathrm{i} \psi^{\text{tar}}(\boldsymbol{\rho}^{\prime})] {.} $$

(8)

Since the target plane is at the far-field zone, the phase of the target field is an approximated spherical phase which is common for all its three components per definition.

$$ \psi^{\text{tar}}(\boldsymbol{\rho}^{\prime}) = k_{0} n \sqrt{\left| \boldsymbol{\rho}^{\prime} \right|^{2}+L^{2}} {,} $$

(9)

where L is the propagating distance from the output plane to the target one.

The amplitude of the target field can be determined by using the given irradiance distribution. Under the far-field assumption, the relation of the irradiance and the electric field is formulated as followed:

$$ \begin{aligned} E_{e}^{\text{tar}}(\boldsymbol{\rho}^{\prime}) &= \frac{n}{2 \mu_{0} c} \hat{\vec{s}}(\boldsymbol{\rho}^{\prime}) \hat{N}(\boldsymbol{\rho}^{\prime}) \left\| \boldsymbol{E}^{\text{tar}}(\boldsymbol{\rho}^{\prime}) \right\|^{2} {,} \end{aligned} $$

(10)

where n is the refractive index, μ₀ the vacuum permeability, and c the speed of light in vacuum. $\hat {\vec {s}}$ is the unit vector of the Poynting vector, which can be calculated by $\hat {\vec {s}} = (x^{\prime }/r, y^{\prime }/r, L/r)$, with (x^′,y^′) is the coordinate on the target plane, and $r=\sqrt {x^{\prime 2}+y^{\prime 2}+L^{2}}$. $\hat {N}$ is the normal vector of the target plane. $\hat {N}=(0, 0, 1)$ if the target plane is perpendicular to the optical axis. ∥E∥ denotes the L² norm of the field, with

$$ \begin{aligned} \left\| \boldsymbol{E} \right\|^{2} &= E_{x}^{2} + E_{y}^{2} + E_{z}^{2} \\ &= E_{x}^{2} + E_{y}^{2} + \left(\frac{s_{x} E_{x} + s_{y} E_{y}}{s_{z}}\right)^{2} {.} \end{aligned} $$

(11)

Therefore, by selecting a polarization state, the target field can be calculated from the given target irradiance.

E^out(κ) is then concluded from E^tar(ρ^′) via an inverse field tracing, that

$$ \begin{aligned} \boldsymbol{E}^{\text{out}}(\boldsymbol{\kappa}) &= \tilde{\boldsymbol{\mathcal{P}}}^{-1} \boldsymbol{\mathcal{F}} \left\lbrace \boldsymbol{E}^{\text{tar}}(\boldsymbol{\rho}^{\prime}) \right\rbrace \\ &= \exp[-\mathrm{i} k_{z}(\boldsymbol{\kappa}) L] \boldsymbol{\mathcal{F}} \left\lbrace \boldsymbol{E}^{\text{tar}}(\boldsymbol{\rho}^{\prime}) \right\rbrace {,} \end{aligned} $$

(12)

with k_z(κ) the z-component of the spatial frequency vector.

After E^out(ρ) and E^out(κ) are prepared, by using the same mathematical model, the L² Monge-Kantorovich problem, the curl-free mapping κ(ρ) can also be solved from Eq. (7).

The stationary phase method is also used to calculate the gradient of the output phase behind the element,

$$ \nabla \psi^{\text{out}}(\boldsymbol{\rho}) = \boldsymbol{\kappa}(\boldsymbol{\rho}) {.} $$

(13)

However, since this time the right side of Eq. (13) is a conservative vector field, ∇ψ^out(ρ) is integrable in any case. Therefore, once κ(ρ) is obtained, ψ^out(ρ) can be calculated by direct integration, regardless of paraxial or non-paraxial, on-axis or off-axis situations.

So far, we present both the schemes with different homeomorphic assumptions in order to reconstruct a required smooth output phase. The mathematical derivation shows that the gradient information of the output phase calculated with the mapping from the irradiance in the spatial domain cannot satisfy the integrability condition. Instead of introducing extra constraints to the method in the spatial domain, the proposed mapping-type Fourier pair synthesis directly results in integrable data for the output phase reconstruction.

Numerical examples of the output phase retrieval

A simple off-axis example with a high divergence angle, where the classic approach that assumes the mapping of the irradiance in spatial domain fails, is demonstrated for the comparison of both approaches, as shown in Fig. 2. We first show the design by the method with the mapping ρ^′(ρ), and illustrate how the non-integrable solution influences the result. We then show the algorithm with mapping κ(ρ) that tackles the problem with an integrable solution.

Fig. 2

(a) An off-axis light-shaping task, shaping a plane wave into a uniform off-axis rectangle pattern. (b)(c) The calculated output phase gradient data, with x- and y-component respectively (Unit: 10⁶rad/m). (d) The calculated output phase by direct numerical integration method with the gradient data (Unit: 10³rad). (e) Irradiance (normalized): simulation result by field tracing, with the designed phase response function used

In the example, the input is a plane wave propagating along the optical axis, with a square aperture of 1×1 mm size. The required optical element is to shape the plane wave into a rectangle homogeneous pattern with a size of 1.5×0.7 m, located on a vertical plane with 1 m away from the shaper, and the target pattern is with 0.5 m offset distance in the y-direction.

Both the irradiance of the input and target are constant values, therefore it is straightforward to conclude that the mapping of the irradiance is between two regular uniform grids, with the shape of a square and a rectangle respectively. We sampled both grids with 401×201 points in x- and y-dimension, and by using Eq. (5), the output phase gradient data is calculated, as shown in Fig. 2 (b) and (c).

The integrability of the resulting gradient data is determined by ψxyout=ψyxout[27, p. 458], where $\psi _{xy} := \frac {\partial ^{2}\psi }{\partial x \partial y}$. Therefore, we use the relative root-mean-squared deviation (RMSD) between $\psi ^{\text {out}}_{xy}$ and $\psi ^{\text {out}}_{yx}$ to determine the integrability of the data, which is evaluated by,

$$ \delta = \sqrt{\frac{\sum_{i}\sum_{j} [\psi^{\text{out}}_{xy}(\boldsymbol{\rho}_{ij}) - \psi^{\text{out}}_{yx}(\boldsymbol{\rho}_{ij})]^{2}} {\sum_{i}\sum_{j} [\psi^{\text{out}}_{xy}(\boldsymbol{\rho}_{ij})]^{2}}} \cdot 100\% \text{.} $$

(14)

The data from Fig. 2 (b) and (c) gives δ=40%.

In order to reconstruct the output phase function based on the data in Fig. 2 (b) and (c), the B-spline technique [28] is applied here for the numerical integration. It is done by introducing B-spline functions with unknown control points to represent the phase function. The control points of the B-spline function is then determined by finding the best fit between the derivative formulas of the B-spline functions and the phase gradient data. By doing so, an optimal approximate smooth phase is obtained and represented by B-spline functions, as shown in Fig. 2 (d).

For testing if the resulting output phase is the required one for the light shaping, we consider the element function as a phase response function Δψ(ρ), which is simply calculated by Δψ(ρ)=ψ^out(ρ)−ψⁱⁿ(ρ). ψⁱⁿ(ρ) is the phase of the input field which is a constant function in this case. A functional component that stores the phase response function of Δψ(ρ) is defined in the software VirtualLab Fusion [29], therefore the obtained output phase from the design can be recalled in the simulation. Simulation is performed with the techniques shown in the field tracing diagram in Fig. 1. A detector set on the target plane gives the irradiance distribution shown in Fig. 2 (e). The shape and the inner irradiance distribution indicate that it deviates quite strongly from the targeted uniform rectangle pattern. Therefore, the approach starts with the mapping of irradiance in the spatial domain cannot provide an accurate output phase.

For the other algorithm that the mapping is established between the Fourier pair, the amplitude of the fields E^out(ρ) and E^out(κ), in both domains, is first prepared from the given source and target irradiance of the task. Fig. 3 (a) and (b) shows both the obtained amplitude of the Fourier pair. We implement the algorithm proposed by Prins [16] to calculate the mapping κ(ρ) in Eq. (7). Their homeomorphism ρ⇔κ is illustrated by two meshes shown in Fig. 3 (c) and (d). The same number of sampling points 401×201 as in the previous approach is used for the design, however, for better visualization, we just show 41×21 mesh nodes in Fig. 3 (c) and (d).

Fig. 3

The squared norm of the amplitude for the Fourier pair: (a) ∥E^out(ρ)∥² in ρ domain, and (b) ∥E^out(κ)∥² in κ domain. (c)(d) The solved one-to-one map between the squared norms of the amplitude in both domains, with the mesh (c) for (a) and its mapped one (d) for (b)

Since the output phase gradient is equivalent to the mapping κ(ρ) according to Eq. (13), it is obtained by applying interpolation technique on the data of κ(ρ). Fig. 4 (a) and (b) show the gradient of ψ^out(ρ). The relative RMSD between $\psi ^{\text {out}}_{xy}$ and $\psi ^{\text {out}}_{yx}$ is δ=0.4%, which is quite smaller than the one from the previous method. The small value indicates the obtained gradient data is integrable.

Fig. 4

The gradient data of the output phase calculated with the mapping type Fourier pair synthesis method, with (a) the x component of the gradient and (b) the y component of it (Unit: 10⁶rad/m). (c) The output phase integrated directly from its gradient data (Unit: 10³rad). (d) The simulated irradiance (normalized) by using the designed phase response function

By direct integration, ψ^out(ρ) is shown in Fig. 4 (c). Similarly, we calculate the phase response function Δψ(ρ) for further investigation. Simulation with the resulting Δψ(ρ) gives the pattern on the irradiance detector as shown in Fig. 4 (d). Here, the detected irradiance clearly indicates that a uniform rectangle pattern which proves the required output phase is obtained accurately.

Demonstration of structure design

The smooth output phase reconstructed by the mapping-type Fourier pair synthesis is also a promising result that can be used for the further structure design of the optical elements, e.g. CGH or freeform surface. Here, we show a demonstration of a freeform lens that is designed based on the obtained output phase.

The freeform lens is set with a predefined planar surface as its front surface and the back one is a freeform surface to be designed. The design of the freeform surface starts with an initial surface. Both the input and output fields are decomposed into local plane waves, with their wave vectors concluded from the input phase and the obtained output phase by the stationary phase function. Both the input and output local plane waves are propagated onto the initial surface and the normal of the surface is corrected by using Snell’s law. The freeform surface then can be concluded from the obtained normals. However, the construction of the freeform lens will introduce amplitude modulation to the input field and Eq. (2) is no longer true. The required output phase ψ^out(ρ) has to be updated for the modulated amplitude behind the freeform lens so that E^out(ρ) and E^out(κ) satisfy the Fourier relation. The mapping-type Fourier pair synthesis illustrated in “Output phase retrieved from the mapping in the Fourier pair” section can again be applied for the output phase retrieval. However, the amplitude of E^out(ρ) should be recalculated with the current freeform lens. The whole algorithm is an iterative approach alternatively performing the output phase retrieval and the freeform surface construction until both the amplitude and the phase of E^out(ρ) meets the Fourier relation with E^out(κ). Since the detailed design algorithm of the freeform lens is not the focus of this article, we just show the results for demonstration. The structure of the freeform lens is shown in Fig. (5) (a) in 3D view and the height profile of the freeform surface is shown in Fig. (5) (b). With the freeform lens, the simulation result in Fig. (5) (c) displays the irradiance distribution detected on the target plane.

Fig. 5

(a) The designed freeform lens. (b) the height profile of the freeform surface (Unit: mm). (c) the simulated irradiance (normalized) with the freeform lens

The above numerical example is a simple one for illustrating the integrability issue in the output phase retrieval with the mapping from the irradiance relation, and in contrast, the phase gradient obtained from the Fourier pair relation does not have an integrability problem. In fact, the latter approach for light shaping is general with no restriction from the phase and shapes of the input field. Another example in the following demonstrates a light shaping problem for a circular shape spherical wave input.

As shown in Fig. 6 (a), the task is to shape a spherical wave to a rectangle pattern in the far-field zone. The full divergence angle of the spherical wave is 30^∘, with a wavelength of 532nm and a circular aperture. The target pattern is located on a vertical plane with 1 m away from the optical element, with a size of 1.2×0.8 m. By the mapping-type Fourier synthesis method, the calculated gradient of the output phase is shown in Fig. 6 (b) and (c). The relative RMSD between the obtained $\psi ^{\text {out}}_{xy}$ and $\psi ^{\text {out}}_{yx}$ reaches δ=0.4%, which indicates the gradient data of the output phase for this example is again integrable. The integrated smooth output phase is shown in Fig. 6 (d).

Fig. 6

(a) An non-paraxial light-shaping task, shaping a spherical wave with circular aperture into a uniform rectangle pattern. (b)(c) The calculated output phase gradient data, with x- and y-component respectively (Unit: 10⁶rad/m and 10⁵rad/m respectively). (d) The calculated output phase by direct numerical integration method with the gradient data (Unit: 10⁴rad)

With the obtained output phase, a freeform lens is again designed for light shaping. We also show the 3D view of the designed freeform lens in Fig. 7 (a). The freeform surface of the lens is shown in Fig. 7 (b). A field tracing simulation is performed with the spherical wave and the designed freeform lens. The irradiance distribution illustrated in Fig. 7 (c) indicates that the lens designed based on the output phase solves the shaping problem in a good way.

Fig. 7

(a) The designed freeform lens. (b) the height profile of the freeform surface (Unit: mm). (c) the simulated irradiance (normalized) with the freeform lens

Conclusion

Homeomorphism through the light path in light-shaping system can be assumed between different fields if all the operators can be considered pointwise. The field tracing diagram provides a comprehensive view for the design algorithms with different mapping assumptions. The retrieval of the required output phase for the design can be proceeded by using the mapping solved between different field pairs. By the mapping-type Fourier pair synthesis, due to the special relation between the output phase gradient and the mapping of the Fourier pair, the integrable gradient data is obtained, which guarantees the output phase can be accurately reconstructed by a direct integration procedure. A smooth freeform surface can be then designed for the shaping based on the retrieved output phase.