Non-imaging metasurface design for collimated beam shaping

Non-imaging optical lenses can shape the light intensity from incoherent sources to a desired target intensity profile, which is important for applications in lighting, solar light concentration, and optical beam shaping. Their surface curvatures are designed to ensure optimal transfer of energy from the light source to the target. The performance of such lenses is directly linked to their asymmetric freeform surface curvature, which is challenging to manufacture. Metasurfaces can mimic any surface curvature without additional fabrication difficulty by imparting a spatially-dependent phase delay using optical antennas. As a result, metasurfaces are uniquely suited to realize non-imaging optics, but non-imaging design principles have not yet been established for metasurfaces. Here, we take an important step in connecting non-imaging optics and metasurface optics, by presenting a phase-design method for beam shaping based on the concept of optimal transport. We establish a theoretical framework that enables a collimated beam to be redistributed by a metasurface to a desired output intensity profile. The optimal transport formulation leads to metasurface phase profiles that transmit all energy from the incident beam to the output beam, resulting in an efficient beam shaping process. Through a variety of examples, we show that our approach accommodates a diverse range of different input and output intensity profiles. Last but not least, a full field simulation of a metasurface has been done to verify our phase-design framework.


Introduction
Shaping light from incoherent sources, such as light-emitting diodes, has important applications in solar light concentration and indoor and outdoor illumination design.This is traditionally achieved by freeform refractive lenses using design principles from non-imaging optics [1].Non-imaging optics focuses on the optimal transfer of energy from a light source onto a target without requiring the formation of an image of the source.The design of non-imaging freeform lenses is rooted in geometrical optics and is formulated as an inverse problem with the aim to redistribute a light source with a known intensity profile into a desired output intensity.Strategies to tackle this inverse problem can be broadly classified into two approaches [2]: iterative schemes [3][4][5][6] and a nonlinear differential equation approach based on an optimal transport formulation [7][8][9].The iterative schemes are based on an initial lens design that is iterated to achieve the desired beam shaping.This requires a priori knowledge of the lens shape, limiting the complexity of the beam shaping.On the other hand, the optimal transport method makes no a priori symmetry assumptions on the lens shape and has been demonstrated to work for a variety of incident light distributions and output illumination patterns [10].This enables complex illumination patterns [2] akin to those obtained using holographic methods for coherent light sources.Noticeably, realizing such beam shaping requires manufacturing of freeform lenses with surface accuracy and smoothness close to the nanoscale to retain their designed functionality [11].Lens manufacturing based on single point diamond turning and injection molding struggles to meet these requirements [12].In contrast, metasurfaces offer the ability to impart any phase profile with no additional fabrication difficulty [13,14] and large-scale fabrication has been demonstrated using laser printing [15] and nanoimprint lithography [16].
Metasurfaces are nanostructured optical surfaces that are capable of manipulating the properties of light in a judicious manner, which has resulted in a wide range of flat opti-cal devices [17][18][19][20][21][22][23][24][25][26].Their functionality is based on imparting a spatially-dependent phase delay using optical nanoantennas [27,28] to shape the transmitted [29][30][31][32] or reflected wavefront [33].The metasurface design process consists of specifying the phase profile needed to achieve a desired beam shaping functionality.The phase profile can be obtained either by analytical expressions [34], numerical calculations using commercial ray-tracing software [35], as well as holographic methods [36].The latter enables almost arbitrary beam shaping [37], but requires an incident light source that is coherent.While progress has been made recently [38,39], shaping of incoherent light using metasurfaces has remained relatively unexplored [40].Applying the concepts from non-imaging optics to the phase design of metasurfaces is a promising approach to shape incoherent light using flat optics.
Here, we demonstrate a metasurface formulation of the optimal transport approach to determine the phase profile for one-dimensional intensity shaping of a collimated beam.The optimal transport concept ensures that all energy from the incident beam is transmitted to the output beam, resulting in a metasurface phase design with efficient beam shaping.Using the generalized law of refraction, we derive a nonlinear differential equation for the metasurface phase profile, which depends on the intensity profiles of the incident and output beams, as well as the separation distance between the metasurface and target plane.We demonstrate three beam-shaping cases with analytical solutions, and two cases which are solved numerically.The results from these five different combinations of incident and output illuminations serve to demonstrate the versatility of the method in the design of non-imaging metasurfaces.
Finally, we verify our theoretical framework with full-field electromagnetic simulation.We simulate a non-imaging metasurface whose phase is the result of our model, and we compare it to its theoretical output.We quantify the differences between the theoretical solution and the full-field simulation, and we find them to be in good agreement.

Theoretical framework
We consider a normally-incident collimated beam of light with a one-dimensional intensity profile I(x), which is redistributed by a refracting metasurface with phase profile ϕ(x) to a desired target illumination with intensity profile E(t x ) (Fig. 1).The coordinate t x at the target is located a distance L from the metasurface and is connected to the x-coordinate at the metasurface through the refraction caused by the metasurface phase gradient.Our aim is to determine the metasurface phase profile that maps the incident intensity profile to the desired target intensity profile.This constitutes an inverse problem, which we formulate by drawing inspiration from the concept of optimal transport in non-imaging freeform optics [8].
We assume that the metasurface is a refractor, which transmits all incident light and is  lossless.Metasurfaces with near unity transmission and 2π phase coverage can be realized with high-refractive-index dielectrics using Huygens'-type [30,41] and nanopost metaatoms [31].
Energy conservation then dictates that the power incident on the metasurface is equal to the power received at the target [42] W where W and T denote the widths of the metasurface and target plane, respectively, and ∂tx ∂x is the Jacobian describing the change of variable from the target coordinate t x to the metasurface coordinate x.Physically, the Jacobian describes the ray expansion or contraction due to the refraction by the metasurface.Equation (1) requires that the power is locally conserved and results in the differential equation where we have lifted the absolute value of the Jacobian.The proper choice of sign in Eq. ( 2) is connected to the boundary conditions, which we discuss in detail later.
Our aim now is to establish a relationship between the target coordinate t x and the metasurface phase profile ϕ(x).From geometrical considerations (Fig. 1), we obtain where O = (O x , O z ) is the unit vector describing the direction of the refracted ray.The direction of the refracted ray is determined by its wave vector k r = (k x , k z ), which is found from the generalized law of refraction.For a normally-incident beam, the generalized law of refraction in one dimension can be expressed as ∂ϕ ∂x = k x , where ϕ(x) is the phase imparted by the metasurface [43][44][45].The direction of the refracted ray is therefore determined by the metasurface phase gradient as Here, we have introduced a scaled wavelength-independent phase Φ(x) = ϕ(x)/(n 0 k 0 ), where n 0 is the refractive index of the medium between the metasurface and the target plane, and k 0 = 2π/λ is the vacuum wavenumber.Combining Eqs.(3)(4), we arrive at an expression for the phase gradient in terms of the target coordinate Equations ( 2) and ( 5) constitute the basis of our theoretical framework.The solution to Eq. ( 2) provides the target coordinate t x (x) for a desired beam shaping functionality, which is specified through the incident intensity I(x) and output intensity E(t x ).Once the target coordinate is obtained, Eq. ( 5) provides the phase gradient, from which we obtain the metasurface phase profile through integration.

Boundary conditions
The proper choice of sign in the differential equation in Eq. ( 2) is determined by the boundary conditions, which set the trajectories of the light rays impinging on the edges of the metasurface.In non-imaging optics, these boundary conditions are given by the edge-ray principle.The edge-ray principle states that a mapping of all rays from the edge of the source distribution should be directed to the edge of the target illumination distribution, in order to ensure that all rays fall within the target [46].In our case, this yields two possible sets of boundary conditions, as illustrated in the insets of Fig. 2(a,c), 1st set: 2nd set: We illustrate the relation between the choice of boundary conditions and the sign in Eq. ( 2) through the beam shaping of a diverging lens (Fig. 2).In our framework, this is described by a target width larger than the metasurface width (T > W ) as well as incident and output beams of constant intensities, i.e., I(x) = I 0 and E(t x ) = E 0 .From global energy conservation [Eq.( 1)] we obtain E 0 = W T I 0 , and using this to solve Eq. (2) yields Here, we observe that the positive (+) sign can only satisfy the first set of boundary conditions [Eq. ( 6)], while the negative (−) only satisfies the second set [Eq. (7)].Consequently, the choice of sign in the differential equation [Eq.( 2)] comes with a unique set of boundary conditions.The two different approaches achieve the same intensity profile at the target plane, however, the metasurface phase profiles differ substantially.By using the relation in Eq. ( 5) and setting Φ(x = 0) = 0, we find the following analytical expressions for the two phase profiles In the case of a positive sign, the phase profile of the metasurface corresponds to a regular diverging lens, which spreads out the incident beam evenly to achieve a constant intensity at the target plane [Fig.2(a,b)].In contrast, the choice of a negative sign leads to a focusing phase profile with a focal length f smaller than the distance to the target plane (f < L) [Fig.2(c,d)].After the initial focusing at the focal point, the rays spread out to result in the same constant intensity at the target plane as for the diverging lens.It is therefore of value to note that even though the phase profiles and ray trajectories differ, both approaches are equally valid within our theoretical framework.Although both approaches result in the same beam shaping, choosing the negative sign may lead to a lens design with a high numerical aperture, which may be difficult to realize.Thus we choose to focus on the approach described by the positive sign in Eq. ( 2) along with the first set of boundary conditions [Eq. ( 6)], i.e., , with BC:

Beam shaping
Having established the theoretical framework, we now determine the metasurface phase profiles for three different beam shaping examples.An advantage of our theoretical framework is that significant analytical progress can be made for certain beam shaping cases.In particular, shaping a Gaussian input intensity profile to a constant output intensity (and vice versa) can be solved analytically.For the Gaussian-to-constant mapping, we use the input intensity profile I(x) = I 0 exp (−2x 2 /w 2 0 ) with w 0 being the waist radius [47], while the target intensity profile is constant E(t x ) = E 0 , see Fig. 3a.In this case we find that the target coordinate is given by  The constant-to-Gaussian mapping (d) displays a convex phase profile (e) consistent with a focusing effect (f).In both cases the metasurface and target plane have the same width (T = W ). The Gaussian-to-square mapping (g) displays overall defocusing, since the target plane is wider than the metasurface (T = 2W ), as well as local focusing to generate the square target profiles (h,i).For all cases, the metasurface width and metasurface-to-target separation are W = 1 mm and L = 100 mm, respectively.
where erf denotes the error function.We visualize the Gaussian-to-constant shaping in Fig. 3(b,c) with the width of the Gaussian set to w 0 = 0.2 mm and the metasurface and target widths being equal (T = W = 1 mm).We find a concave-type phase profile (Fig. 3b), which is reasonable as the beam shaping taking place is similar to a diverging lens (see Fig. 2a).The rays are directed towards the perimeter of the target plane, resulting in a flattening of the Gaussian beam shape (Fig. 3c).It is however worth noting that the phase gradient differs from a purely diverging lens as the edge rays must be mapped to the edge of the target.Since the target and metasurface widths are equal, this forces the phase gradient to go to zero at the edges of the metasurface (see Eq. ( 5)).
For a constant-to-Gaussian mapping, the input intensity profile is constant I(x) = I 0 , while the target intensity has a Gaussian profile E(t x ) = E 0 exp (−2t 2 x /w 2 0 ) with w 0 = 0.4 mm, see Fig. 3d.Here, we find that the target coordinate is where erf −1 is the inverse error function.The phase profile has a convex shape (Fig. 3e), which leads to focusing of the incident beam to the desired Gaussian target (Fig. 3f).
As a final demonstration, we consider a Gaussian incident beam that is mapped to a target profile consisting of two square bumps with a ratio of 5:1 between the maximum and minimum intensity (Fig. 3g).The minimum output intensity is set to a low but finite value, as zero-valued output intensities violates the local energy conservation dictated by Eq. ( 2).
The width of the target plane is set to twice the width of the metasurface (T = 2W = 2 mm).
To determine the phase profile, we solve the theoretical model numerically.To this end, it is worth noting that the first-order nature of the differential equation in Eq. ( 10) requires in principle only a single boundary condition to determine a solution.To satisfy both boundary conditions dictated by the edge-ray principle [Eq.( 10)], we must limit the solution space to beam shaping of symmetric input and output intensity profiles.This is a consequence of the one-dimensional nature of our theoretical framework.This constraint can be released by extending the framework to two dimensions, since boundary rays can then map to arbitrary points on the edge of the target profile and thereby enable asymmetric intensity profiles, as demonstrated with freeform optics [8].Nonetheless, taking into account the symmetry condition of our framework comes with the advantage that Eq. ( 10) can be solved as an initial value problem.We therefore implemented a numerical solver using a built-in initial value solver in MATLAB, which is based on the Runge-Kutta method.We use the left boundary as the initial condition, i.e., t x (−W/2) = −T /2, and observe that the right boundary condition is automatically satisfied (Fig. 3i).The resulting phase profile and phase gradient contain both diverging and focusing features (Fig. 3h).In particular, the phase gradient consists of an overall positive slope, which spreads the input beam to accommodate the wider target plane (as also seen in Fig. 2a), as well as two local negatively sloped features that serve to focus the input beam into the two square target profiles (Fig. 3h).The ray tracing shows that the boundary conditions are satisfied, and that the designed beam shaping is achieved (Fig. 3i).

Metasurface design and verification
We verify the feasibility of our theoretical model by performing a full-field 2D simulation (COMSOL Multiphysics 6.0) of a miniaturized metasurface.The metasurface is designed to split an incident Gaussian beam into two focused Gaussian profiles (Fig. 4a).Our theoretical model and its numerical implementation determines a wavelength-independent scaled phase profile for the desired beam shaping.We retrieve the actual phase profile of the metasurface at an operating wavelength of λ = 660 nm by the operation ϕ(x) = Φ(x)n 0 2π/λ with n 0 = 1 (Fig. 4b).The phase profile is then discretized according to the look-up table based on silicon nanobeams in air (Fig. 4c).The look-up table is generated using normally-incident plane waves and periodic boundary conditions, and applies for silicon nanobeams with a height h = 0.5λ placed in an array with a period p = 0.5λ.The entire metasurface is constructed from 455 nanobeam elements (located at z = 0 in Fig. 4d), resulting in a metasurface width of 150 µm.We remark that our framework is rooted in geometrical optics and does not take diffraction into account, which is appropriate for shaping of incoherent light.However, in the full-field simulation, diffraction is present due to the interaction of the finite-sized metasurface and the infinite extent of the incident Gaussian beam.To mimic incoherent light and minimize diffraction, the waist radius of the Gaussian beam w 0 is carefully chosen to ensure that the intensity of the beam is near zero at the edges of the metasurface.Simultaneously, the width of the beam also needs to be large enough to ensure a depth of focus that generates a beam profile that is collimated throughout the simulation domain (in the absence of the metasurface).From these considerations, we find that w 0 = 32.5 µm is a suitable waist radius.This trade-off is mainly due to computational limitations on the size of the simulation domain.We believe that an experimental demonstration based on collimated light from, e.g., a light-emitting diode, combined with a significantly larger metasurface would lift this constraint.
The propagation of the output beam and a comparison between the target output intensity and simulated intensity are shown in Fig. 4d and Fig. 4e, respectively.The transmission through the silicon nanobeams dips below unity for beams with a width in the range 70-120 nm (Fig. 4c, dashed lines).The assumption in our theoretical model, that the metasurface has unity transmission is thus not met for most of the metasurface (Fig. 4b, dashed lines).
Nonetheless, the shape of the simulated and target intensities are in good agreement.We observe an intensity ratio of 4.25:1 and a peak to peak distance of 94 µm in the simulated profile, and a slight asymmetry in the peaks.The target profile has a peak separation of 98 µm and a 5:1 intensity ratio, and symmetric Gaussian peaks.Despite the previously discussed losses in the nanobeams, we observe that 91% of the incident power is transmitted in the simulation.This result demonstrates the feasibility of our theoretical model for determining the phase profile for beam shaping applications, as well as how a fairly complex phase profile can be implemented with a metasurface.

Conclusion
In summary, we have demonstrated a theoretical optimal transport framework for obtaining the phase profile of a metasurface for one-dimensional beam shaping with collimated incident light.We show how the resulting nonlinear differential equation can be solved analytically in some cases, and numerically using a Runge-Kutta method in more complex cases.Five specific cases are investigated to demonstrate the versatility of the method within the design of non-imaging metasurfaces, and the method is verified with full-field simulations.We add that although we have only shown examples with collimated beams, our framework can be expanded to describe point-like light sources as well.Thus we believe our work is an important step towards manipulating incoherent light with metasurfaces.
Funding S. R. acknowledges support by the Independent Research Funding Denmark (7026-00117B).

Figure 1 :
Figure 1: Schematic representation of the functionality of the non-imaging metasurface.A normally-incident collimated beam I with intensity profile I(x) is incident on a metasurface with phase profile ϕ(x).The metasurface, positioned at z = 0, has a width of W and refracts the beam onto the target plane with width T and positioned a distance L from the metasurface.The refracted ray O hits the target plane at a specific coordinate t x to produce the desired intensity profile E(t x ).The boundary rays of the beam (blue) map from the edges of the metasurface to the edges of the target plane in accordance with the edge-ray principle.

Figure 2 :
Figure 2: (a,c) Phase profiles and phase gradients, and (b,d) resulting ray plots for a diverging non-imaging metasurface.In (c-d) the boundary conditions are reversed relative to (a-b), as illustrated in the insets.In (b,d) the dashed black lines indicate the boundary rays.The metasurface width, target plane width and metasurface-to-target separation are W = 1 mm, T = 15 mm, and L = 10 cm, respectively.

Figure 3 :
Figure3: Three examples of beam shaping: (a-c) Gaussian-to-constant mapping, (df) constant-to-Gaussian mapping, and (g-i) Gaussian-to-square mapping.The Gaussianto-constant mapping (a) exhibits a concave phase profile (b) that spreads the beam (c).The constant-to-Gaussian mapping (d) displays a convex phase profile (e) consistent with a focusing effect (f).In both cases the metasurface and target plane have the same width (T = W ). The Gaussian-to-square mapping (g) displays overall defocusing, since the target plane is wider than the metasurface (T = 2W ), as well as local focusing to generate the square target profiles (h,i).For all cases, the metasurface width and metasurface-to-target separation are W = 1 mm and L = 100 mm, respectively.

Figure 4 :
Figure 4: (a) Input Gaussian intensity and target double focal spots.(b) Phase profile at the wavelength λ = 660 nm along with corresponding beam widths determined from the look-up table presented in (c).(c) Simulated look-up table for 2D periodic silicon nanobeams in air (n Si = 3.84 + i0.02).The period and height are set to p = h = 0.5λ.(d) Full-field simulation result of the intensity of the output beam.The metasurface is located at z = 0, and centered around x = 0. (e) Simulated intensity evaluated at z = L = 500 µm (green) and the target output intensity (dashed).The metasurface and target planes have widths of W = 150 µm and T = 1.1W , respectively.