## Abstract

We extend our previous work [Yang et al., Opt. Express **29**, 3621 (2021) [CrossRef] ] and propose an iterative algorithm to design a freeform surface for far-field light shaping. The algorithm alternately performs a wavefront phase design step and a freeform surface construction step. The smooth wavefront phase is designed by the mapping-type Fourier pair synthesis method, and the freeform surface is constructed by using the obtained wavefront phase. The algorithm provides a solid approach that ensures the introduction of the required wavefront phase manipulation for light shaping. Moreover, the related physical effects such as the Fresnel effect and polarization effect are included in the algorithm. We demonstrate the flexibility of the algorithm by examples.

© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

The far-field light-shaping problem deals with designing an optical element that generates a prescribed irradiance distribution of light in the far-field zone. It is indicated in the authors’ previous publication that the problem mainly demands a Fourier pair synthesis [1]. From a physical-optics perspective, the Fourier pair synthesis method provides the required output phase for the shaping. Under a homeomorphic assumption, the method can be simply reduced to finding a mapping between the Fourier pair and an integration step. The method is then named as the mapping-type Fourier pair synthesis method [1].

In this article, we apply the mapping-type Fourier pair synthesis method for a functional phase design and extend it to a structural design. We aim to design a freeform surface for the light-shaping problem in non-paraxial situations. The algorithm for the freeform surface design is developed based on the obtained phase. We solve the problem by the phase design step and the surface design step in an iterative way.

Similar works in the literature were proposed by Feng et al. [2,3]. Feng et al. have developed a Monge-Ampère equation for solving the required output phase. The freeform surface is then constructed by the thin element approximation (TEA) method in the paraxial situation [2], while in non-paraxial situations, the obtained output phase is used to derive an integrable ray map for generating the freeform surface [3]. In [3], Feng et al. successfully designed the freeform surface also in an iterative manner, that the required output phase is recalculated in every iteration and used for the surface construction.

In the phase design, instead of developing a complex non-linear equation like in [2,3], we use the mapping-type Fourier pair synthesis method [1], which can be performed more practically. Moreover, in the surface design step, because our algorithm is developed based on physical-optics modeling techniques, it provides a more general way that includes all the related physical effects in the design. For example, the Fresnel effect of the freeform surface and the polarization of the incident light is naturally taken into account, while these effects are generally excluded either in the Monge-Ampère equation method [4–10] or the “ray mapping method” [11–17]. In the literature, both methods are based on a mapping assumption with the energy conservation between the input and target. The structure of the freeform surface is then solved from the derived Monge-Ampère equation or constructed with the obtained mapping. The freeform surface is calculated only by geometric relations, without considering the related physical effects.

In the following, we first formulate the modeling of the demanded setup and analyze it via physical optics. Then, the iterative procedure of the freeform surface design algorithm is developed and presented step by step. Finally, with examples, we demonstrate the flexibility of the algorithm and show that it can be applied even for the case with multi-mode source.

## 2. Physical optics modeling

In this section, the far-field light-shaping problem is analyzed with physical optics. We illustrate the light-shaping system and the notations in Fig. 1 for the discussion, as same in [18]. The optical element with a freeform surface to be designed modulates the incident light to achieve a given irradiance distribution on the target plane, as shown in Fig. 1(a). In Fig. 1(b), a field tracing diagram shows how to model the electric field of the light passing through the system.

In Fig. 1(b), $\boldsymbol {E}(\boldsymbol {\rho })$ represents the vectorial electric field, i.e. $\boldsymbol {E}\left (\boldsymbol {\rho }\right )=\left ( {E}_{x} \left (\boldsymbol {\rho }\right ), {E}_{y} \left (\boldsymbol {\rho }\right ), {E}_{z} \left (\boldsymbol {\rho }\right ) \right )$, which is defined in the spatial domain ($\boldsymbol {\rho }$ domain), with $\boldsymbol {\rho }=\left ( x,y \right )$. $\tilde {\boldsymbol {E}}\left ( \boldsymbol {\kappa } \right )$ with an over-tilde “˜” notation denotes the fields defined in the spatial-frequency domain ($\boldsymbol {\kappa }$ domain), with $\boldsymbol {\kappa }=\left ( \kappa _{x},\kappa _{y} \right )$, the $x$ and $y$ component of the wave vector, respectively.

Each component of the electric field in the $\boldsymbol {\rho }$ domain is formulated explicitly as,

The $\boldsymbol {\mathcal {B}}$ in the field tracing diagram is the operator used to model the functionality of the freeform element. In this study, the local plane interface approximation (LPIA) method is applied for the operator [19]. Therefore, the relation of the input field in front of the freeform element and the output field behind is written as,

Here, the polarization state of the input field and the Fresnel effect of the freeform element are included in the calculation of the output field.

In Fig. 1(b), for the free space propagation behind the freeform element, the field tracing sequence from $\boldsymbol {E}^{\textrm {out}}\left ( \boldsymbol {\rho }^{\textrm {out}} \right )$ to $\boldsymbol {E}^{\textrm {tar}}\left ( \boldsymbol {\rho }^{\textrm {tar}} \right )$ is formulated as:

In general, the rigorous fast Fourier transform (FFT) technique is employed for both $\boldsymbol {\mathcal {F}}$ and $\boldsymbol {\mathcal {F}}^{-1}$ operators. However, since the target plane is at the far-field zone, $\boldsymbol {\mathcal {F}}^{-1}$ can be accurately performed by the homeomorphic Fourier transform (HFT) method [20], which is a pointwise calculation. Therefore, a one-to-one map, or the so-called homeomorphism is established between $\tilde {\boldsymbol {E}}^{\textrm {out}}( \boldsymbol {\kappa }^{\textrm {out}} )$ and $\boldsymbol {E}^{\textrm {tar}}\left ( \boldsymbol {\rho }^{\textrm {tar}} \right )$. Moreover, if the HFT can be also approximately applied for the $\boldsymbol {\mathcal {F}}$ operator, the homeomorphism is extended to the field pair between $\boldsymbol {E}^{\textrm {out}}( \boldsymbol {\rho }^{\textrm {out}} )$ and $\boldsymbol {E}^{\textrm {tar}}\left ( \boldsymbol {\rho }^{\textrm {tar}} \right )$. It is noted in [18] that this is indeed the assumption of the geometric-optics-based algorithms in the literature. We also use the homeomorphic assumption in our design.

## 3. Design algorithm

In the design, the given information from the light-shaping task is the source and the target irradiance distribution in the far field. We develop our design algorithm by an inverse thinking of the modeling from the previous section.

#### 3.1 Design of the required output wavefront phase

Both the fields $\boldsymbol {E}^{\textrm {in}}\left ( \boldsymbol {\rho }^{\textrm {in}} \right )$ and $\boldsymbol {E}^{\textrm {tar}}\left ( \boldsymbol {\rho }^{\textrm {tar}} \right )$ are first calculated from the given information. $\boldsymbol {E}^{\textrm {in}}\left ( \boldsymbol {\rho }^{\textrm {in}} \right )$ is calculated by tracing the field of the given source to the input plane. $\boldsymbol {E}^{\textrm {tar}}\left ( \boldsymbol {\rho }^{\textrm {tar}} \right )$ is derived from the target irradiance distribution ${E}_e^{\textrm {tar}}\left ( \boldsymbol {\rho }^{\textrm {tar}} \right )$, with the far-field assumption. The derivation is presented as follows.

We denote the target field explicitly as,

The wavefront phase $\psi ^{\textrm {tar}}$ is a spherical shape per definition of the far field [20]. $\psi ^{\textrm {tar}}$ is also common for all its three components.

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 [21, Chap. 11]:

Therefore, by selecting a polarization state $[E_x, E_y]$, the target field can be calculated from the given target irradiance with Eqs. (5)–(8).

For the design, the freeform element is required to generate a proper output field $\boldsymbol {E}^{\textrm {out}}\left ( \boldsymbol {\rho }^{\textrm {out}} \right )$ from $\boldsymbol {E}^{\textrm {in}}\left ( \boldsymbol {\rho }^{\textrm {in}} \right )$, so that after the free space propagation, the target field $\boldsymbol {E}^{\textrm {tar}}\left ( \boldsymbol {\rho }^{\textrm {tar}} \right )$ can be achieved.

We perform an inverse field tracing of $\boldsymbol {E}^{\textrm {tar}}(\boldsymbol {\rho }^{\textrm {tar}})$ to $\boldsymbol {\kappa }-$ domain and calculate the field $\tilde {\boldsymbol {E}}^{\textrm {out}}(\boldsymbol {\kappa }^{\textrm {out}})$ by,

Please remark that in Eq. (4), under the far-field assumption, the value of the phase term introduced by the propagation operator is much larger than the phase of $\tilde {\boldsymbol {E}}^{\textrm {out}}\left ( \boldsymbol {\kappa }^{\textrm {out}} \right )$, i.e. $k_z \left ( \boldsymbol {\kappa }^{\textrm {out}}\right ) \Delta z \gg \textrm {arg}\left [ \tilde {\boldsymbol {E}}^{\textrm {out}}\left ( \boldsymbol {\kappa }^{\textrm {out}} \right ) \right ]$, and contributes mainly to the phase of $\boldsymbol {E}^{\textrm {tar}}\left ( \boldsymbol {\rho }^{\textrm {tar}} \right )$. Therefore, it is allowed to apply a phase freedom for $\tilde {\boldsymbol {E}}^{\textrm {out}}\left ( \boldsymbol {\kappa }^{\textrm {out}} \right )$ and leave only the amplitude of $\tilde {\boldsymbol {E}}^{\textrm {out}}\left ( \boldsymbol {\kappa }^{\textrm {out}} \right )$ as the constraint.

We start the algorithm with an initial freeform element and calculate $\boldsymbol {E}^{\textrm {out}}(\boldsymbol {\rho }^{\textrm {out}})$ by Eq. (3). Then, under the amplitude constraint of $\tilde {\boldsymbol {E}}^{\textrm {out}}(\boldsymbol {\kappa }^{\textrm {out}})$, we design a wavefront phase $\psi ^{\textrm {out}}\left ( \boldsymbol {\rho }^{\textrm {out}} \right )$ for $\boldsymbol {E}^{\textrm {out}}(\boldsymbol {\rho }^{\textrm {out}})$ to fulfill the Fourier relation. Here, we assume the HFT is employed for the Fourier transform $\boldsymbol {\mathcal {F}}$. Therefore, with both the amplitude of $\boldsymbol {E}^{\textrm {out}}\left ( \boldsymbol {\rho }^{\textrm {out}} \right )$ and $\tilde {\boldsymbol {E}}^{\textrm {out}}\left ( \boldsymbol {\kappa }^{\textrm {out}} \right )$ obtained from Eqs. (3) and (9), the mapping-type Fourier pair synthesis method [1] is applied for designing the required $\psi ^{\textrm {out}}\left ( \boldsymbol {\rho }^{\textrm {out}} \right )$ for $\boldsymbol {E}^{\textrm {out}}\left ( \boldsymbol {\rho }^{\textrm {out}} \right )$.

#### 3.2 Amplitude modulation from the freeform surface

After the wavefront phase $\psi ^{\textrm {out}}\left ( \boldsymbol {\rho }^{\textrm {out}} \right )$ is obtained, we proposed an algorithm with the LPIA for the construction of the freeform surface. The designed freeform surface is aimed at realizing the required $\psi ^{\textrm {out}}\left ( \boldsymbol {\rho }^{\textrm {out}} \right )$. The algorithm will be introduced in detail in Section 3.3.

However, the designed freeform surface may introduce new amplitude modulations to $\boldsymbol {E}^{\textrm {in}}\left ( \boldsymbol {\rho }^{\textrm {in}} \right )$, that gives a new amplitude for $\boldsymbol {E}^{\textrm {out}}\left ( \boldsymbol {\rho }^{\textrm {out}} \right )$ according to Eq. (3). Thus, via a Fourier transform, the amplitude constraint of $\tilde {\boldsymbol {E}}^{\textrm {out}}(\boldsymbol {\kappa }^{\textrm {out}})$ is not satisfied. To overcome the problem, we perform the above procedure in an iterative way, so that the output wavefront phase design and the freeform surface construction is alternatively carried out. In this way, the the freeform surface is adapted to realize the required output field under the amplitude constraint.

#### 3.3 Entire algorithm in an iterative procedure

The entire algorithm for the freeform surface design is presented explicitly in this section. The flow chart of the design procedure is illustrated in Fig. 2. The two block diagrams indicate the wavefront phase design step and the freeform surface construction step. To start the algorithm, we set the light-shaping element as a double surface element composed of a certain material. The element contains a predefined surface and the freeform surface to be designed. Our algorithm allows the predefined surface to be arbitrary shapes, which also can be set as either the front or the back surface of the element.

The details of the procedure are listed as follows.

- 1. The freeform surface is initialized, which is denoted as $S^{0}=\left \lbrace \boldsymbol {r} \in \mathbb {R}^{3} \mid \boldsymbol {r} = \left ( \boldsymbol {\rho },z^{0}\left ( \boldsymbol {\rho } \right ) \right ) \right \rbrace$.
- 2. $\boldsymbol {E}^{\textrm {out}}\left ( \boldsymbol {\rho }^{\textrm {out}} \right )$ and $\tilde {\boldsymbol {E}}^{\textrm {out}}\left ( \boldsymbol {\kappa }^{\textrm {out}} \right )$ are calculated by Eq. (3) and Eq. (9), respectively. Then the amplitudes of both fields are extracted for the design of $\psi ^{\textrm {out}}\left (\boldsymbol {\rho }^{\textrm {out}}\right )$, where the mapping-type Fourier pair synthesis method is applied [1]. Please remark that in Eq. 3, the polarization state of the input field and the Fresnel effect of the surface is included in the operator. Therefore, these physical effects are naturally embedded in the design algorithm.
- 3. An internal iterative process is introduced to correct the freeform surface for realizing the required $\psi ^{\textrm {out}}\left (\boldsymbol {\rho }^{\textrm {out}}\right )$. Here, we assume that both $\boldsymbol {E}^{\textrm {in}}\left ( \boldsymbol {\rho }^{\textrm {in}} \right )$ and $\boldsymbol {E}^{\textrm {out}}\left ( \boldsymbol {\rho }^{\textrm {out}} \right )$ are in their homeomorphic zone [20], that the fields can be decomposed into local plane waves. The LPIA operator is applied for realizing the fields propagation through the surfaces, including the predefined surface.
- (a) In the $j-$th iteration, the height profile of the freeform surface is indicated by $z^{0,j} \left (\boldsymbol {\rho }\right )$.
- (b) On both the input and output plane of the freeform element, the fields $\boldsymbol {E}^{\textrm {in}}\left ( \boldsymbol {\rho }^{\textrm {in}} \right )$ and $\boldsymbol {E}^{\textrm {out}}\left ( \boldsymbol {\rho }^{\textrm {out}} \right )$ are decomposed as local plane waves, and the transversal components of each wave vector are calculated by$$\boldsymbol{\kappa}^{\textrm{in}} = \nabla \psi^{\textrm{in}}\left(\boldsymbol{\rho}^{\textrm{in}}\right),$$$$\boldsymbol{\kappa}^{\textrm{out}} = \nabla \psi^{\textrm{out}}\left(\boldsymbol{\rho}^{\textrm{out}}\right),$$where $\psi ^{\textrm {in}}\left (\boldsymbol {\rho }^{\textrm {in}}\right )$ and $\psi ^{\textrm {out}}\left (\boldsymbol {\rho }^{\textrm {out}}\right )$ are the wavefront phase of both fields. Then the wave vectors are derived with
Here, $n^{\textrm {in}}$ and $n^{\textrm {out}}$ are the refractive indices in front and behind the element, respectively.

- (c) Both sets of the local plane waves are propagated to the freeform surface. The intersection positions are indicated by $\boldsymbol {r} \left ( \boldsymbol {\rho }^{\textrm {in}}\right ) = \left ( \boldsymbol {\rho } \left (\boldsymbol {\rho }^{\textrm {in}}\right ),z^{0,j}\left (\boldsymbol {\rho } \left (\boldsymbol {\rho }^{\textrm {in}}\right )\right )\right )$ and $\boldsymbol {r} \left ( \boldsymbol {\rho }^{\textrm {out}}\right ) = \left ( \boldsymbol {\rho } \left (\boldsymbol {\rho }^{\textrm {out}}\right ),z^{0,j}\left (\boldsymbol {\rho } \left (\boldsymbol {\rho }^{\textrm {out}}\right )\right )\right )$, respectively. The wave vectors of each local plane at the intersection positions are indicated as $\boldsymbol {k}^{\textrm {in}}\left (\boldsymbol {r} \left (\boldsymbol {\rho }^{\textrm {in}}\right )\right )$ and $\boldsymbol {k}^{\textrm {out}}\left (\boldsymbol {r} \left (\boldsymbol {\rho }^{\textrm {out}}\right )\right )$.
- (d) For the position $\boldsymbol {r}=\left (\boldsymbol {\rho },z^{0,j}\left (\boldsymbol {\rho }\right )\right )$, Snell’s law is applied to recalculate the normal vector $\boldsymbol {N} (\boldsymbol {r})$ of the freeform surface with respect to the local plane waves from the front and the back. If the freeform surface is the front surface of the element, the following equation is applied,
If the freeform surface is the back surface of the element,

$$\boldsymbol{N}\left(\boldsymbol{r}\right)=n^{\textrm{out}}\boldsymbol{k}^{\textrm{out}}\left(\boldsymbol{r}\right)-n^{\textrm{FE}}\boldsymbol{k}^{\textrm{in}}\left(\boldsymbol{r}\right) {.}$$$n^{\textrm {FE}}$ is the refractive index in the freeform element. The samples of local plane waves from front and back may intersect on different positions on the surface. Thus, in order to apply Snell’s law, the wave vectors of both sets of local plane waves should be interpolated at the same positions. The B-spline technique is applied here for the interpolation of $\boldsymbol {k}^{\textrm {in}}\left (\boldsymbol {r}\left (\boldsymbol {\rho }^{\textrm {in}}\right )\right )$ and $\boldsymbol {k}^{\textrm {out}}\left (\boldsymbol {r}\left (\boldsymbol {\rho }^{\textrm {out}}\right )\right )$. - (e) The gradient of the freeform surface is derived from its recalculated normal vector by$$\nabla z^{0,j+1}\left(\boldsymbol{\rho}\right) = \left(\frac{N_{x}\left(\boldsymbol{\rho}\right)}{N_{z}\left(\boldsymbol{\rho}\right)},\frac{N_{y}\left(\boldsymbol{\rho}\right)}{N_{z}\left(\boldsymbol{\rho}\right)}\right),$$where $N_{x}, N_{y}, N_{z}$ are the three components of the vector $\boldsymbol {N}$.
- (f) The B-spline technique is applied for the integration of $\nabla z^{0,j+1}\left (\boldsymbol {\rho }\right )$, and constructing an updated surface profile $z^{0,j+1}\left (\boldsymbol {\rho }\right )$. The numerical implementation of the B-spline technique is illustrated in the last part of the current section.
- (g) The root-mean-square (RMS) of the deviation between the height profiles of $j$ and $j + 1$ iteration is calculated for a merit function,$$\varepsilon = \sqrt{ \frac{\sum_{\boldsymbol{\rho}} \left( z^{0,j+1}(\boldsymbol{\rho})-z^{0,j}(\boldsymbol{\rho}) \right)^{2}}{N_{\boldsymbol{\rho}}} } {,}$$where $N_{\boldsymbol {\rho }}$ is the sampling number of $\boldsymbol {\rho }$.
- (h) Step (c) to (g) are performed iteratively until a threshold value of $\epsilon$ is achieved.
The internal iterative procedure mentioned above will converge to a freeform surface that realizes the required output wavefront phase $\psi ^{\textrm {out}}\left (\boldsymbol {\rho }^{\textrm {out}}\right )$. So far, the surface $S^{0}$ is updated as $S^{1}$.

- 4. To determine whether the amplitude constraint is fulfilled by the designed freeform surface $S^{1}$, the field $\boldsymbol {E}^{\textrm {in}}\left ( \boldsymbol {\rho }^{\textrm {in}} \right )$ is traced through the updated freeform element and the output field in the $\boldsymbol {\kappa }$ domain is recalculated by$$\tilde{\boldsymbol{E}}^{\textrm{out'}}\left( \boldsymbol{\kappa}^{\textrm{out}} \right)= \boldsymbol{\mathcal{F}} \left\lbrace \boldsymbol{\mathcal{B}}^{\textrm{LPIA}} \boldsymbol{E}^{\textrm{in}}\left( \boldsymbol{\rho}^{\textrm{in}} \right) \right\rbrace,$$where $\boldsymbol {\mathcal {B}}^{\textrm {LPIA}}$ is the LPIA modeling operator of the freeform element.
- 5. We use the relative RMS deviation between the amplitude of $\tilde {\boldsymbol {E}}^{\textrm {out}}\left ( \boldsymbol {\kappa }^{\textrm {out}} \right )$ and $\tilde {\boldsymbol {E}}^{\textrm {out'}}\left ( \boldsymbol {\kappa }^{\textrm {out}} \right )$ as the merit function.$$\eta = \sqrt{ \frac{ \sum_{\boldsymbol{\kappa}^{\textrm{out}}} \left( \left\| \tilde{\boldsymbol{E}}^{\textrm{out}^{\prime}}(\boldsymbol{\kappa}^{\textrm{out}}) \right\| - \alpha \left\| \tilde{\boldsymbol{E}}^{\textrm{out}}(\boldsymbol{\kappa}^{\textrm{out}}) \right\| \right)^{2}}{ \sum_{\boldsymbol{\kappa}^{\textrm{out}}} \alpha \left\| \tilde{\boldsymbol{E}}^{\textrm{out}}(\boldsymbol{\kappa}^{\textrm{out}}) \right\|^{2}} }{,}$$where $\alpha$ is a scaling factor given by$$\alpha = \frac{ \displaystyle \sum_{\boldsymbol{\kappa}^{\textrm{out}}} \left\| \tilde{\boldsymbol{E}}^{\textrm{out}^{\prime}}(\boldsymbol{\kappa}^{\textrm{out}}) \right\| \left\| \tilde{\boldsymbol{E}}^{\textrm{out}}(\boldsymbol{\kappa}^{\textrm{out}}) \right\| }{ \displaystyle \sum_{\boldsymbol{\kappa}^{\textrm{out}}} \left\| \tilde{\boldsymbol{E}}^{\textrm{out}}(\boldsymbol{\kappa}^{\textrm{out}}) \right\| ^{2} } {.}$$
If the value of $\eta$ is close to zero, the amplitude constraint is achieved and the freeform surface is obtained. Otherwise, the freeform surface is altered by performing the wavefront phase design of step (2) and the surface construction procedure from step (3) to (5) iteratively till a value of $\eta$ close to zero is achieved.

In the surface construction step, the freeform surface is obtained with the calculated gradient data. Here, the integration is numerically implemented by using the B-spline functions [22]. We first formulate the surface function as

We set up a minimization functional with respect to the gradient of $z(\boldsymbol {\rho })$ to determine the coefficients $C_{kl}$.

A least-squares method is applied to minimize $\Delta _z(z, \nabla z^{j})$ and solve the coefficients $C_{kl}$ of $z(\boldsymbol {\rho })$ in Eq. (20). In this way, $z(\boldsymbol {\rho })$ can be interpreted as the optimum potential function of the gradient data of $\nabla z^{j}$.

## 4. Design examples

In this section, numerical examples are shown to indicate the validity of the algorithm. The implementation of the freeform surface design algorithm and the relative simulations are performed in the software VirtualLab Fusion [23].

#### 4.1 Plane wave to a certain pattern

The first light-shaping example is to generate a specific irradiance pattern by an incident plane wave. The plane wave is in $x-$ polarized, with its beam size of $1\times 1~\textrm {mm}$ and its wavelength of $532~\textrm {nm}$. The target irradiance distribution is shown as in Fig. 3, with size of $1.2~\textrm {m} \times 0.8~\textrm {m}$. It is located at a vertical plane $1~\textrm {m}$ distance away from the freeform element. Therefore, the system is in a non-paraxial situation with the outgoing light having the divergent angle of $62~{\circ } \times 44~{\circ }$. The irradiance values in Fig. 3 are scaled for better visualization. The freeform element is set with a predefined planar surface as its front surface, and the back surface is to be designed. The refractive index of the medium in the freeform element is set as $1.52$.

To start the iterative process of the design algorithm, the freeform surface is initialized simply as a planar surface. Following the algorithm presented in Section 3.3, the required output wavefront phase $\psi ^{\textrm {out}}\left ( \boldsymbol {\rho }^{\textrm {out}} \right )$ is obtained by the mapping-type Fourier pair synthesis method [1]. The obtained $\psi ^{\textrm {out}}\left ( \boldsymbol {\rho }^{\textrm {out}} \right )$ is then used for the design of the freeform surface. The internal iterative procedure in step 3 of the entire algorithm is applied to construct the freeform surface from both the input/output wavefront phases. In the design, $\boldsymbol {E}^{\textrm {in}}\left ( \boldsymbol {\rho }^{\textrm {in}} \right )$ is decomposed into $101 \times 101$ local plane waves, and $\boldsymbol {E}^{\textrm {out}}\left ( \boldsymbol {\rho }^{\textrm {out}} \right )$ is decomposed into $601 \times 401$. The freeform surface is constructed by B-spline functions with $101 \times 101$ control points. The procedure converges fast that, within $6$ iterations, the deviation of the height profile between two iterations $\varepsilon$ converges to nanometer magnitude. The initial planar surface is then updated with the new freeform surface, and the whole process starts over. The whole process between the wavefront phase design and the freeform surface design converges to a surface that meets the design constraint with $\eta = 0.06$. Figure 4(a) shows the 3D view of the designed freeform element. The height profile of the freeform surface is shown in Fig. 4(b).

To analyze the performance of the designed freeform element, a physical-optics simulation is performed. The operators used in the simulation are listed in Table 1, which are also compared with the operators in the design. In the simulation, VirtualLab Fusion applies a strict mathematical criterion to analyze the output field $\boldsymbol {E}^{\textrm {out}}\left ( \boldsymbol {\rho }^{\textrm {out}} \right )$, and concludes that the HFT is accurate enough here for the Fourier transform. Therefore, the software automatically selects HFT in the simulation. The inverse Fourier transform at the target plane is also decided as HFT by VirtualLab Fusion.

We perform the first simulation with the same $x-$ polarized input beam as in the design task, the irradiance distribution on the target plane is detected, as shown in Fig. 5(a). The result is coincident with the target pattern in Fig. 3, with a relative root-mean-square deviation (RRMSD) of $0.05$ in the irradiance. Therefore, the proposed algorithm results in an accurate freeform surface that solves the light-shaping task. Please note, that the identity between the detected irradiance and the target irradiance also indicates that the Fresnel effect of the freeform surface has been compensated in the design, in such a non-paraxial situation.

To demonstrate that the polarization information of the incident beam is embedded in the design, we perform two more simulations by using a $y-$ polarized and a $45~{\circ }$ polarized plane wave, respectively. The detected irradiance distribution is illustrated in Fig. 5(b) and (c). Compared to Fig. 5(a), both results show deviations from the target distribution, especially the modulation of the brightness at the background area of the pattern. The RRMSD from the target irradiance of both results are $0.16$ and $0.14$, respectively. Therefore, we can conclude that the designed freeform surface is of better performance for the $x-$ polarized input beam. The design based on physical optics takes the polarization from the task into account.

#### 4.2 VCSEL to a Top-hat pattern

In the design algorithm, the construction of the freeform surface mainly uses the information of both the wavefront phase of the input and output fields. Therefore, our proposed algorithm can also be used for multi-mode sources where the input fields of all the source modes share a same wavefront phase.

In the following, another example that shaping a VCSEL source into a top-hat pattern is shown to demonstrate the algorithm. The VCSEL source is modelled with a combination of two Laguerre Gaussian (LG) modes. The fields of the VCSEL source is represented by the incoherent sum of the LG modes. We place the freeform lens sufficiently far away from the VCSEL to ensure the freeform surface is in the far field of the modes. Therefore, the two LG modes share the same spherical wavefront phase per definition of the far field.

The irradiance of the VCSEL source in front of the freeform element is shown in Fig. 6(a), where the freeform element is set at a distance of $2~\textrm {mm}$ to the source. We use the incoherent sum of the LG$_{00}$ and LG$_{01}$ mode to model the VCSEL source, where their weights factor is optimized to fit the irradiance of the VCSEL in the far field. The optimized weight of each mode with respect to their intensity is $55.981:53.895$. The irradiance of both modes are shown in Fig. 6(b) and (c), respectively. In this task, the VCSEL is in a $x-$ polarized state, therefore, both the LG modes are $x-$ polarized. Their wavelength is set as $1~\mu \textrm {m}$. The Rayleigh length of both modes is $19.32~\mu \textrm {m}$ and their half divergence angle is $7.31~{\circ }$.

In this task, the target irradiance distribution is a Top-hat pattern which is formulated as a super Gaussian function

The freeform element is defined with a predefined planar surface as its front surface, and the back surface of the element is a freeform one to be designed. The thickness of the element is $0.5~\textrm {mm}$ and its refractive index $1.5$. The freeform surface is initialized as a spherical surface with $1.5~\textrm {mm}$ radius. Then, the design algorithm is performed. In the example, the number of local plane waves for both input and output are the same, that is $501 \times 501$. We use B-spline functions with $101 \times 101$ control points to represent the freeform surface. The algorithm converged rapidly to $\eta = 0.012$ with $4$ iterations. Figure 7(a) shows the 3D view of the designed freeform element. The height profile of the freeform surface is shown in Fig. 7(b).

A physical-optics simulation is performed for the investigation of the design. The operators used in the simulation are same as in the previous example (Table 1). With the $x-$ polarized VCSEL source and the designed freeform element, the detected irradiance distribution on the target plane is shown in Fig. 8(a). The cross section along the $x-$ axis indicates that the top-hat profile is achieved. The RRMSD from the target irradiance distribution is $0.03$, which indicates that the designed freeform surface solves the shaping problem in an accurate manner.

Similarly, an additional simulation with $y-$ polarized incident light is performed for the demonstration of the polarization effect. The detected irradiance distribution is shown in Fig. 8(b). It is indicated in the cross section that the result is deviated from a uniform profile. The RRMSD is about $0.09$, which is almost three times of the value in $x-$ polarized case. The results in Fig. 8 again demonstrate that the proposed algorithm includes polarization of the incident light.

For each example presented above, the physical-optics simulation in VirtualLab Fusion on a PC (2.39 GHz Intel Core Quad CPU) takes a few seconds.

## 5. Conclusion

A freeform surface design algorithm is proposed for the far-field light-shaping problem. The design of the freeform surface is combined with a wavefront phase design step in an iterative way. The entire algorithm provides an insight that how a Fourier pair synthesis can be applied as a part of a more general solution for light-shaping. The algorithm solves light-shaping problems flexibly for the single-mode light source and the multi-mode source with a single wavefront phase. Moreover, as the algorithm is developed based on the physical optics, the design takes all the relative physical effects into account. Therefore, a high degree of possibility for the design is provided. Although the algorithm in this article is developed under the far-field assumption, it can be generalized straightforwardly when the far field is not reached, which will be our future work.

## Funding

Wyrowski Photonics GmbH.

## Acknowledgments

The authors would like to thank our colleagues Mr. Vignesh Balaji for assistance with proofreading the manuscript. We are also grateful to the company Wyrowski Photonics GmbH for procuring a VirtualLab Fusion [23] license free of charge for the simulations included in this work.

## Disclosures

The authors declare that there are no conflicts of interest related to this article.

## Data Availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

## References

**1. **L. Yang, I. Badar, C. Hellmann, and F. Wyrowski, “Light-shaping design by a Fourier pair synthesis: the homeomorphic case,” Opt. Express **29**(3), 3621–3630 (2021). [CrossRef]

**2. **Z. Feng, B. D. Froese, and R. Liang, “Composite method for precise freeform optical beam shaping,” Appl. Opt. **54**(31), 9364–9369 (2015). [CrossRef]

**3. **Z. Feng, D. Cheng, and Y. Wang, “Transferring freeform lens design into phase retrieval through intermediate irradiance transport,” Opt. Lett. **44**(22), 5501–5504 (2019). [CrossRef]

**4. **H. Ries and J. Muschaweck, “Tailored freeform optical surfaces,” J. Opt. Soc. Am. A **19**(3), 590–595 (2002). [CrossRef]

**5. **R. Wu, P. Liu, Y. Zhang, Z. Zheng, H. Li, and X. Liu, “A mathematical model of the single freeform surface design for collimated beam shaping,” Opt. Express **21**(18), 20974–20989 (2013). [CrossRef]

**6. **K. Brix, Y. Hafizogullari, and A. Platen, “Designing illumination lenses and mirrors by the numerical solution of Monge Ampère equations,” J. Opt. Soc. Am. A **32**(11), 2227–2236 (2015). [CrossRef]

**7. **C. Bosel and H. Gross, “Single freeform surface design for prescribed input wavefront and target irradiance,” J. Opt. Soc. Am. A **34**(9), 1490–1499 (2017). [CrossRef]

**8. **L. L. Doskolovich, D. A. Bykov, E. S. Andreev, E. A. Bezus, and V. Oliker, “Designing double freeform surfaces for collimated beam shaping with optimal mass transportation and linear assignment problems,” Opt. Express **26**(19), 24602–24613 (2018). [CrossRef]

**9. **D. Bykov, L. Doskolovich, A. Mingazov, and E. Bezus, “Optics mass transportation problem in the design of freeform optical elements generating far-field irradiance distributions for plane incident beam,” Appl. Opt. **58**(33), 9131–9140 (2019). [CrossRef]

**10. **L. B. Romijn, J. H. M. ten Thije Boonkkamp, and W. L. IJzerman, “Freeform lens design for a point source and far-field target,” J. Opt. Soc. Am. A **36**(11), 1926–1939 (2019). [CrossRef]

**11. **A. Bauerle, A. Bruneton, R. Wester, J. Stollenwerk, and P. Loosen, “Algorithm for irradiance tailoring using multiple freeform optical surfaces,” Opt. Express **20**(13), 14477–14485 (2012). [CrossRef]

**12. **A. Bruneton, A. Bäuerle, R. Wester, J. Stollenwerk, and P. Loosen, “High resolution irradiance tailoring using multiple freeform surfaces,” Opt. Express **21**(9), 10563–10571 (2013). [CrossRef]

**13. **Z. Feng, L. Huang, G. Jin, and M. Gong, “Designing double freeform optical surfaces for controlling both irradiance and wavefront,” Opt. Express **21**(23), 28693–28701 (2013). [CrossRef]

**14. **Y. Schwartzburg, R. Testuz, A. Tagliasacchi, and M. Pauly, “High-contrast computational caustic design,” ACM Trans. Graph. **33**(4), 1–11 (2014). [CrossRef]

**15. **Z. Feng, B. D. Froese, and R. Liang, “Freeform illumination optics construction following an optimal transport map,” Appl. Opt. **55**(16), 4301–4306 (2016). [CrossRef]

**16. **K. Desnijder, P. Hanselaer, and Y. Meuret, “Ray mapping method for off-axis and non-paraxial freeform illumination lens design,” Opt. Lett. **44**(4), 771–774 (2019). [CrossRef]

**17. **S. Wei, D. Ma, Z. Zhengbo, and Z. Fan, “Least-squares ray mapping method for freeform illumination optics design,” Opt. Express **28**(3), 3811–3822 (2020). [CrossRef]

**18. **L. Yang, I. Badar, C. Hellmann, and F. Wyrowski, “Light shaping by freeform surface from a physical-optics point of view,” Opt. Express **28**(11), 16202–16210 (2020). [CrossRef]

**19. **R. Shi, C. Hellmann, and F. Wyrowski, “Physical-optics propagation through curved surfaces,” J. Opt. Soc. Am. A **36**(7), 1252–1260 (2019). [CrossRef]

**20. **Z. Wang, O. Baladron-Zorita, C. Hellmann, and F. Wyrowski, “Generalized far-field integral,” Opt. Express **29**(2), 1774–1787 (2021). [CrossRef]

**21. **D. J. Griffiths, * Introduction to Electrodynamics* (Cambridge University, 2017), 4th ed.

**22. **I. Badar, L. Yang, C. Hellmann, and F. Wyrowski, “Antiderivative of gradient data by spline model integration,” J. Opt. Soc. Am. A **38**(8), 1187–1193 (2021). [CrossRef]

**23. **Physical optics simulation and design software “Wyrowski VirtualLab Fusion”, developed by Wyrowski Photonics GmbH, distributed by http://www.lighttrans.com LightTrans International UG, Jena, Germany.