Comprehensive numerical design approach for refractive laser beam shapers to generate annular irradiance profiles

Abstract. A refractive laser beam shaper typically consists of either two plano-aspheric lenses or one thick lens with two aspherical surfaces. Ray mapping is a general optical design technique for irradiance reshaping based on geometric optics. Although ray mapping, in principle, allows generating any rotational-symmetric irradiance profile, in the literature this technique is mainly used to transform a Gaussian irradiance profile to a uniform rotational-symmetric profile. For more complex profiles, especially with low intensity in the inner region (such as annular profiles), a high sampling rate is required to ensure an accurate calculation of the surfaces. In practice, the high sampling rate increases the numerical effort to calculate the aspherical surfaces and the simulation time to verify the design considerably. In this work, we evaluate different sampling approaches and surface construction methods. This allows us to propose and demonstrate a comprehensive numerical approach to efficiently design refractive laser beam shapers to generate rotational-symmetric collimated beams with annular irradiance profiles. Ray tracing analysis for several annular irradiance profiles demonstrates the excellent performance of the designed lenses and the versatility of our design procedure.


Introduction
Laser beam shaping is widely used in many industrial and medical applications.In most cases, a laser beam with a Gaussian irradiance profile is reshaped to have a uniform output profile, usually called a flat-top profile.In some applications, such as optical trapping, 1 atom guiding, 2 and laser additive manufacturing, an annular irradiance profile is required in order to ensure the optimal system performance.One example of such an annular irradiance profile is a dark hollow Gaussian (DHG) irradiance profile.There are various proposed techniques for laser beam shaping, especially to generate the rotational symmetric flat-top irradiance profiles.The simplest technique is either to truncate an expanded Gaussian beam 3 or to use an apodized absorption filter. 4uch a technique allows one to generate uniform intensity patterns but is rather energy-inefficient.Typical energy-efficient approaches include refractive 5 or reflective systems, 6 diffractive elements, 7 holograms, 8 and microlens arrays. 9mong these techniques, refractive laser beam shapers are commonly used due to their high optical efficiency and simple structure. 10our general configurations of refractive beam shapers are given in Fig. 1.They are formed either by two plano-aspheric lenses or by one single thick lens with two aspherical surfaces.The general optical design approach of such refractive beam shapers is based on ray mapping within the framework of geometric optics.The problem was initially proposed by Frieden in 1965. 11Then, Kreuzer patented a design procedure in 1969 12 which is still widely used today.The ray mapping technique was further investigated and developed by Rhodes 13 and Hoffnagle. 5,14,15Additional research led to the simplification of the practical design steps 16 and optimization of the design algorithm for better performance. 17ith the ray mapping technique, it is possible to generate any irradiance profile.However, published works mainly focus on flat-top profile generation and do not describe the practical steps of how to generate more complex profiles.For generating annular irradiance profiles with low intensity in the inner region, it is difficult to analytically derive a solution; the available numerical design approach 16 to generate the elementary profiles does not work properly.Based on the ray mapping technique, this work presents an efficient numerical approach to design the refractive laser beam shapers to get annular irradiance profiles.With our design approach, both the numerical effort to calculate the aspherical surfaces and the simulation time to evaluate the designed lens system have been dramatically reduced.Also, the designed aspherical surfaces can be more accurately constructed by comparing different surface construction methods.
This article is structured as follows.In Sec. 2, the design procedure of how to design a refractive beam shaper is described.In the first design step for calculating the mapping function between the radial positions on the two aspherical surfaces, we have introduced both the conventional sampling method and our proposed alternative method.In the last step for constructing the designed surfaces, two methods are investigated and compared.One method is polynomial fitting and the other one is spline fitting which describes the surface profile by piecewise low degree polynomials.In Sec. 3, ray tracing results for different annular irradiance profiles are presented and analyzed.Taking a DHG profile as an example, the efficiency of two sampling methods is compared, highlighting the advantage of the alternative sampling method.Furthermore, the accuracy of two different surface construction methods is discussed and compared.In the end, additional simulation results for two other annular irradiance profiles are presented, underlining the versatility of our design procedure.Finally, in Sec. 4, conclusions are drawn and an outlook is given.

Design Procedure of Refractive Beam Shapers
The design of a refractive beam shaper based on geometric ray mapping is typically divided into four steps, which are listed in Table 1.Following these four design steps, this section describes both the conventional approach to transform a Gaussian irradiance profile to a flat-top profile 16 and our alternative approach to generate annular profiles.The core differences between the two approaches for each step are also given in Table 1.Detailed explanations of different sampling methods and surface construction methods can be found later in their corresponding design steps.The two-lens Galilean-type configuration (see Fig. 2) is used for all designs in this article.It should be noticed that the presented approach can easily be adapted to design other configurations in Fig. 1 by only modifying the slope function Eq. ( 11) in Step 2.

2.1
Step 1: Mapping r to R Based on the law of conservation of energy, the power at the input plane is equal to the power at the output plane assuming there is no energy loss through the lens system.In order to map every radial coordinate r to R, the power inside each circle of the input radial coordinate r should be the same as the power at the inside of each circle of the output radial coordinate R.
Both the input and output irradiance profiles are rotationally symmetric, and can be described by functions F in ðrÞ and F out ðRÞ.The encircled power within r at the input plane is a function of r, given by the integral function G in ðrÞ: Similarly, at the output plane, the function G out ðRÞ calculates the encircled power within R as There are two possible ways to calculate the mapping relationship between r and R from Eq. (1) and Eq.(2).Conventionally, the total power P is directly sampled into equal power bins and, for each power bin, the corresponding radial positions on the first aspherical surface r and on the second aspherical surface R are calculated accordingly.We call this method "sampling P." For annular irradiance profiles, we have found that the numerical calculation effort can be considerably reduced by sampling the output radial coordinate R with alternative equidistant points.More precisely, the complexity of the input and output irradiance profiles has   to be properly compared and the radial coordinate at the plane with the more complex profile should be sampled.
As the desired annular output profile (e.g., a DHG profile) has a null intensity in the inner region, it is considered to be more complex than the input Gaussian distribution.The output radial coordinate R should be sampled.Therefore, we name this alternative sampling method as "sampling R."

Sampling P
The sampling P method is commonly used in the literature to generate a flat-top irradiance profile.Here, both the total input and output power are normalized to 1 or 100%.First, the maximum radii at the input and output planes r max and R max should be properly chosen, so that sufficient power P max (e.g., 99.9%) can be transferred from the input to the output.Then, the total power P max is sampled by equal bins.Supposing the sampling number is N, the sampling index i becomes f1;2; : : : ; Ng.There are a total of (N − 1) bins.For each bin, the encircled power is then given as For each power bin, the corresponding radii at the input plane and output plane r i and R i can be derived from the inverse functions of G in ðrÞ and G out ðRÞ: In this way, r is mapped to R and the generated data pairs fr i ; R i g are used for further calculation.

Sampling R
In order to generate annular irradiance patterns, we propose an equidistant sampling of the output radial coordinate R as an efficient alternative.In this case, we work as follows: 1.A maximum value for R has to be chosen large enough to transmit sufficient power P (e.g., P max ¼ 99.9%).
The value R max should be verified to ensure a feasible aperture size for the designed lens.R max is calculated by 2. R max is equidistantly sampled by a sufficient number of points N. Similarly to previous calculations, the sampling index i is f1;2; : : : ; Ng.Then 3. The encircled power within R i is calculated by 4. The corresponding r i at the input plane is calculated by The calculated data pairs fr i ; R i g provide the mapping relationship between r and R.

Step 2: Calculating the Surface Slopes
In Fig. 2, the slopes v i and V i are of two aspheric lens profiles which are defined as Due to the fact that both input and output beams are collimated, it follows from symmetry that the angles θ 1 and θ 2 are equal and as a result, By applying Snell's law at each aspherical surface and keeping the optical path length between the input plane and the output plane constant, the slope values v i or V i can be described as a function of r i , R i , the separation between two lenses D, and the refractive index of the lenses n 2 . 12The local slope value at every radial position is calculated by

Step 3: Calculating the Surface Sags
After the calculation of the local slopes at each surface, the surface sags z i and Z i can be numerically calculated by integrating the slopes.The specific numerical integration method should be carefully selected in order to minimize the numerical calculation error.In this article, the trapezoidal numerical integration is used, available within MATLAB 18 via the command "TRAPZ".

Step 4: Constructing the Aspherical Surfaces
In order to evaluate the optical performance of the beam shaper, it is necessary to construct smooth surfaces from the calculated discrete data sets of sag and slope values.We evaluate two different surface construction methods: polynomial fitting and spline fitting.

Polynomial fitting
With polynomial fitting an analytical polynomial function is fitted to each set of sag points that describe each surface.To generate a flat-top profile, it is sufficient to use a polynomial function with only even exponents: 16 For more complex irradiance profiles such as annular profiles, we have found that a polynomial function with only even exponents is not adequate no matter how many coefficients are used.Instead, a polynomial function with both even and odd exponents should be used: Here, it should be noticed that r is always positive and the pattern is rotationally symmetric.
The coefficients fA 1 ; A 2 ; : : : ; A n g can be determined by either fitting the polynomial function zðrÞ to the calculated sags or by fitting its first derivative z 0 ðrÞ to the calculated slopes. 16In this article, the polynomial function is fitted to the calculated sags in MATLAB 18 using the command "POLYFITN".This algorithm calculates the coefficients of the polynomial function that fits the data set best by minimizing the least-squares.

Spline fitting
As an alternative, spline fitting uses piecewise low order polynomial functions to match both calculated sag and slope values in order to describe the surface profile. 19We have used a continuous piecewise cubic spline to construct the designed aspherical surfaces in the optical ray tracing software ASAP 20 with the "SPLINE" command.

Ray Tracing Results and Discussions
Our designed refractive beam shapers are evaluated by ray tracing in ASAP and the simulation results are presented and discussed in this section.
The input and output irradiance profiles of our example to transform a Gaussian profile to a DHG profile are shown in Fig. 3.The input beam has a circular Gaussian irradiance profile.The irradiance of a Gaussian profile I G is expressed by where r is the radial coordinate, I a is the maximum (maximum) irradiance at r ¼ 0, and w G is the beam waist.The irradiance at the beam waist is equal to I a × e −2 .In our example, I a is calculated by normalizing the total power to 1, and w G is 2.366 mm.The output beam has a DHG irradiance profile; therefore, the beam can be described by hollow Gaussian beams (HGBs). 21Then, the irradiance I DHG is calculated by n ¼ 0;1; 2; : : : ; where I b is a constant, n is the order of the HGB, and r is the radial coordinate.In the case where n is equal to 0, the DHG is identical by the fundamental Gaussian with beam waist w DHG .As energy is conserved, I b is calculated by also normalizing the total power to 1. n is chosen as 1 and w DHG is 2 mm.

Ray Tracing Results for Designs Obtained with Different Sampling Methods
Two sampling methods have been introduced in Step 1 of Sec.2: sampling P and sampling R. In order to evaluate these two sampling methods, we have used both methods with different sampling numbers to design beam shapers that generate the same DHG profile.After constructing the calculated aspherical surfaces by spline fitting, both designs are simulated in ASAP and ray tracing results are compared.
With the sampling P method, several beam shapers are designed to generate the DHG output profile by using different numbers of sampling points.For each design, a grid of 2001 × 2001 rays has been traced through the lens system with ASAP to directly generate the output intensity profile after the second lens.A comparison of the detected output profiles for different designs is given in Fig. 4. In order to match the required DHG profile, the number of sampling points N has to be as large as 10,000.With fewer sampling points, the output profile departs from the expected DHG profile, especially for the center part.The reason for this is that the original sampling of equidistant power bins in combination with a null/very low intensity at the optical axis results in an equally low sampling rate of the second aspherical lens surface in that region.Therefore, a much higher overall sampling rate (thousands of bins) is necessary to achieve a sufficiently high sampling of the second aspherical lens surface close to the optical axis.This problem is solved by our proposed sampling R method, as it samples the second aspherical surface with equidistant sampling points.
The ray tracing result for a design with sampling R by only 100 sampling points is shown in Fig. 5.It clearly shows that 100 sampling points are already sufficient, which is much less than when we sample the encircled power P.

10,000
Fig. 4 Ray tracing results for designs obtained with sampling P method for different number N of sampling points: the shown simulation results (solid lines) in comparison with the required DHG profile (dashed line) reveal the high sampling rate that is necessary to achieve good agreement of the intensities especially in the central region.With fewer sampling points, the numerical efforts to calculate the lens profiles and the simulation time to evaluate the designed beam shaper can be drastically reduced, as we will explain in more detail in the next section.This finding highlights the benefits of the sampling R method to efficiently generate annular irradiance profiles.

Ray Tracing Results for Designs Obtained with Different Surface Construction Methods
We consider two possible ways to construct the surfaces from the calculated surface sages and/or slopes: polynomial fitting and spline fitting.Designs obtained with these two methods are simulated in ASAP.For each simulation, a grid of 2001 × 2001 rays is traced and evaluated not only at the plane directly after the second lens, but also at planes with certain distances from the second lens along the z axis.
In this way, it is possible to verify how well the generated output profile is maintained after propagation, which is a good indication of the collimation quality of the output beam.
For the method of polynomial fitting, a polynomial function in Eq. ( 13) with 21 coefficients is used.Ray tracing results for the design with a polynomial fitting are presented in Fig. 6.The dashed line is the required DHG profile, whereas the different solid lines correspond to the simulation results for the output profile after certain propagation distances along the z axis.For z ¼ 0.1 mm, the output profile corresponds well to the expected DHG profile.The irradiance profile starts to deviate from the DHG profile when it propagates over distances beyond 100 mm, and after 1000 mm, the difference is very clear.For the design of transforming a Gaussian irradiance profile to a DHG profile, a recently published article has proposed that the first aspherical surface should be described by a polynomial function with fractional instead of integer exponents. 22This proposed polynomial function with fractional exponents has provided a very low mapping error, yet the generated DHG profile seems to not be maintained well enough after propagating by 1000 mm.
Figure 7 shows the simulation results for the design obtained with spline fitting.In this case, the irradiance profile agrees well with the required DHG profile over the complete propagation trajectory, which proves perfect collimation of the output beam.It should be noticed that the simulation time for surfaces constructed by spline fitting is strongly dependent on the number of piecewise polynomial functions, which is defined by the number sampling points.Fewer sampling points significantly reduce the simulation time.

Further Examples of Generated Annular Irradiance Profiles
The presented design approach allows the design of refractive laser beam shapers for various complex irradiance profiles, which are rotational symmetric.As a further example for annular profiles, a Gaussian beam is transformed into an output beam with a cosine irradiance profile, which has both a center peak and an outside ring.The output cosine irradiance profile is a cosine function where a constant is added so that the minimum irradiance is 0. Similarly, we used 100 points to sample the output radial coordinate R and spline fitting for the surface construction.The ray tracing result with a grid of 2001 × 2001 rays is given in Fig. 8.The simulation result is in good agreement with the expected cosine profile.
Fig. 5 Ray tracing results for designs obtained with sampling R method for a moderate number of only 100 sampling points: the simulated DHG irradiance profile (solid line) matches the required DHG output profile (dashed line) very well.
Fig. 6 Ray tracing results for designs obtained with polynomial fitting: z is the propagation distance from the beam shaper.After propagating only 0.1 mm, the simulated output irradiance profile still fits well with the required DHG profile (dashed line); however, the irradiance profile begins to deviate from the required DHG profile with increased propagation distances.
Fig. 7 Ray tracing results for designs obtained with spline fitting: z is the propagation distance from the beam shaper.The simulated irradiance profiles (solid lines) are in excellent agreement with the required DHG profile (dashed line) over the complete propagation trajectory.
As another example, the laser beam shaper is designed to generate a sine irradiance output profile which has two concentric rings.The output sine irradiance profile is a sine function shifted upward to have zero minimum irradiance and moved rightward to have a rotational symmetry axis at 0. Applying the same design approach and simulation procedure, the ray tracing result is shown in Fig. 9.It also agrees well with the expected sine profile.
These two examples conform the versatility of our design procedure and highlight its potential use for generating different annular irradiance profiles with single or multiple rings.

Conclusion
Ray mapping-based optical design is a powerful technique to design the refractive laser beam shapers for irradiance redistribution.If the desired profile is more complex than a flat-top profile, especially with low intensity in the inner region such as annular profiles, the available design procedures have to be modified.In the first step to calculate the mapping function, we have identified an alternative efficient sampling method.Our method by sampling the output radial coordinate R requires much fewer sampling points than the conventional sampling method by sampling the encircled power P. Normally, only 100 sampling points are sufficient.Depending on the particular beam shaping problem, the optical designer should carefully identify and select the appropriate sampling method that allows an accurate and fast design and reconstruction of the lens surfaces.In the last step, to construct the calculated aspherical surfaces, both polynomial and spline fittings can provide satisfactory results.For a polynomial fitting, the polynomial series should have both even and odd exponents rather than only even exponents.Compared with the polynomial fitting, spline fitting has better precision because the generated annular profile can propagate unchanged though a much longer distance along the optical axis.
The described design approach in this article can be easily followed to design a refractive lens system that transforms a Gaussian irradiance pattern to an annular pattern with rotational symmetry.It should be noticed that both input and output beams are collimated, and are thus described by plane wave fronts.A generalization to cope with arbitrary wave fronts will be addressed in future work.
Fig. 8 The output cosine profile is a cosine function shifted upward, so that the minimum irradiance is 0. The simulation result (solid line) of the designed beam shaper to transform a Gaussian profile to a cosine profile shows good agreement with the expected output profile (dashed line).
Fig. 9 The output sine irradiance profile is a sine function shifted upward to have zero minimum irradiance and moved rightward to have a rotational symmetry axis at 0. The simulation result (solid line) of the designed beam shaper to transform a Gaussian profile to a sine profile shows good agreement with the expected output profile (dashed line).

Fig. 1
Fig.1Refractive laser beam shapers can be divided into four general configurations.

Fig. 2
Fig.2Illustration of the two-lens Galilean-type configuration that transforms a collimated input beam to a collimated output beam: an arbitrary ray (from A to B) intersects the input plane at radius r and arrives at the output plane at radius R; z and Z describe the sag values of the first and second aspherical surfaces; t 1 and t 2 define the thickness of the first and second aspherical surfaces; D is the separation between the two lenses.

Fig. 3
Fig.3Irradiance profiles: input Gaussian profile and output DHG profile.r is the radial coordinate and I is the irradiance.

Table 1
Design steps of refractive beam shapers.