Freeform imaging systems: Fermat’s principle unlocks “first time right” design

For more than 150 years, scientists have advanced aberration theory to describe, analyze and eliminate imperfections that disturb the imaging quality of optical components and systems. Simultaneously, they have developed optical design methods for and manufacturing techniques of imaging systems with ever-increasing complexity and performance up to the point where they are now including optical elements that are unrestricted in their surface shape. These so-called optical freeform elements offer degrees of freedom that can greatly extend the functionalities and further boost the specifications of state-of-the-art imaging systems. However, the drastically increased number of surface coefficients of these freeform surfaces poses severe challenges for the optical design process, such that the deployment of freeform optics remained limited until today. In this paper, we present a deterministic direct optical design method for freeform imaging systems based on differential equations derived from Fermat’s principle and solved using power series. The method allows calculating the optical surface coefficients that ensure minimal image blurring for each individual order of aberrations. We demonstrate the systematic, deterministic, scalable, and holistic character of our method with catoptric and catadioptric design examples. As such we introduce a disruptive methodology to design optical imaging systems from scratch, we largely reduce the “trial-and-error” approach in present-day optical design, and we pave the way to a fast-track uptake of freeform elements to create the next-generation high-end optics. We include a user application that allows users to experience this unique design method hands-on.


Introduction
Optical imaging systems have been playing an essential role in scientific discovery and societal progress for several centuries 1,2 . For more than 150 years scientists and engineers have used aberration theory to describe and quantify the deviation of light rays from ideal focusing in an imaging system, and to develop methods to design diffractionlimited imaging systems 3 . Until recently most of these imaging systems included spherical and aspherical refractive lenses (dioptric systems) or reflective mirrors (catoptric systems) or a combination of both (catadioptric systems). The last decennia, with the introduction of new ultra-precision manufacturing methods, such as singlepoint diamond turning and multi-axis polishing 4 , twophoton polymerization 5,6 or other additive manufacturing technologies 7 , it has become possible to fabricate lenses and mirrors that have at least one optical surface that lacks the common translational or rotational symmetry about a plane or an axis. Such optical components are called freeform optical elements 8 and they can be used to greatly extend the functionalities, improve performance, and reduce volume and weight of optical imaging systems that are principal parts of spectrometers 9 , telescopes 10,11 , medical imaging systems 12,13 , augmented and virtual reality systems 14 , or lithography platforms 15,16 . As such imaging systems including freeform optical elements will be key to tremendously advance science and engineering in a wide range of application domains such as astronomy, material research, chip fabrication, visualization, metrology and inspection, energy production, safety and security, biotechnology and medical imaging and diagnosis 17 .
Today the design of optical systems largely relies on efficient ray-tracing and optimization algorithms for which a variety of commercial software 18 and optimization algorithms 19 are available. During an optical design cycle, different parameters of the optical system such as the optical material parameters, radii, coefficients, and positions of the optical surfaces are varied to optimize a defined merit function that indicates the image quality for a given field of view 20 . These merit functions are typically "wild" with many local minima, and there is no guarantee that local or global optimization algorithms will lead to a satisfactory design solution 21 . A successful and widely used optimization-based optical design strategy therefore consists of choosing a well-known optical system as a starting point (e.g. from literature) and steadily achieving incremental improvements. Such a "step-and-repeat" approach to optical design, however, requires considerable experience, intuition and guesswork, which is why it is sometimes referred to as "art and science" 22 . This applies especially to freeform optical systems. With the potential of freeform optics widely recognized and its advantages essential to the creation of next-generation ultra-performing optical systems, it has become a top priority to overcome the laborious and hardly reproducible trial-anderror approaches currently used for their design.

Freeform optical design strategies
So far four main design strategies for freeform optical systems have been proposed and developed. All of them target to guide researchers and optical designers to a sufficiently good freeform design that, if needed, can then be used as starting point for subsequent design optimization.
A first strategy involves the so-called direct design methods, such as the Simultaneous Multiple Surface (SMS) design method 23 . They rely on solving geometrical or differential equations describing the freeform optical system under study to achieve a well-performing initial design [23][24][25][26][27][28][29][30][31] . Although these methods show clear merits, so far, they lack a straightforward path to increase the number of optical surfaces that can be calculated 25,29,31 . Because of this limitation in scalability, their applicability remained limited until today.
A second strategy focuses on the automated design and advanced optimization methods to generate highperformance freeform systems with reduced efforts by the optical designer [32][33][34][35][36][37] . Although this strategy aims to provide practical design tools to achieve systems with better optical performance, they do not offer valuable insights in the optical design process. Therefore, the cause for not reaching a satisfactory solution often remains unclear, with the only option to restart the process repeatedly.
A third strategy consists in calculating an initial design that is free of certain aberrations and then rely on optimization techniques to balance all aberrations and yield the best overall imaging performance. One can for example start with an initial rotationally symmetric on-axis design, that has been corrected for several aberrations [38][39][40][41] , then introduce freeform surfaces and iterate from there, while unobscuring the light path in the system by introducing tilts for the optical surfaces. Alternatively, first-order unobstructed, plane-symmetric systems of three or four spherical surfaces can be calculated as starting points 42,43 . So far, with this design strategy, only a limited number of distinct low-order aberration terms have been canceled or controlled.
A fourth strategy is based on nodal aberration theory 44-46 that has been extended to freeform surfaces [47][48][49] . This approach allows us to predict and visualize the contribution of Zernike terms of an optical surface to the aberration fields and provides information to the designer which Zernike terms are required to correct an aberration of the system 50,51 . So far, nodal aberration theory is one of the more systematic freeform design strategies. It can guide the experienced designer towards a successful design provided a substantial and in-depth understanding of the underlying theory.
In this paper, we present a novel hybrid direct design method (patent no. WO/2019/129872 52 ) for optical imaging systems that also allows us to systematically expand the range of aberration terms that need to be controlled, suppressed, or canceled. The method is holistic because it can be used to design imaging systems that include catoptric and/or dioptric spherical, aspherical, and/or freeform surfaces and components. It allows us to match user-defined conditions such as minimal blurring for each individual order of aberrations by calculating the corresponding coefficients of the optical surfaces, only neglecting order-crossing aberration correlations. As a result, our proposed method allows us to calculate initial imaging systems "first time right" for a given geometry. In that sense, it provides a local and deterministic solution rather than a global solution. Moreover, the method is scalable, as it allows us to add additional optical surfaces in an unrestricted manner. The method also allows the calculation of aberrations of sufficiently high orders to immediately and accurately estimate the overall imaging performance of the system under design. We highlight the substantial added value for the optical design community by demonstrating the systematic, scalable, deterministic, flexible, and holistic character of the design approach with distinct high-end examples, both for catoptric and catadioptric systems. Furthermore, we provide the opportunity for hands-on experience with a design application for imaging systems based on three freeform mirrors.

Materials and methods
Freeform optical systems, aberrations, and their mathematical description There are various possibilities to classify aberrations in freeform optical systems 53,54 . For the design method presented in this paper we distinguish three different categories, according to the overall system symmetry present: (1) systems with two orthogonal planes of symmetry, (2) systems with one plane of symmetry, and (3) systems without any symmetry. Here, we will only consider systems with one plane of symmetry as it is the most common category. The two other symmetry categories can be treated similarly.
To describe the transverse ray aberrations (see Fig. 1) in the case of an optical system with the x-z plane of symmetry we follow the systematic description of aberrations of non-rotationally symmetric systems by Barakat and Houston 53 . Following an ansatz similar to that used in most aberration theories, we assume an arbitrary ray that is emitted from an object at infinity with a field direction vector H = (sin(H x ), sin(H y ), (1sin(H x ) 2 -sin(H y ) 2 ) 1/2 ) for field angles H x and H y that passes through the pupil of the system (first surface in Fig. 1) at (x p , y p , f1(x p ,y p )), and that finally intersects the image plane at h = (h x , h y , D). This intersection can be written in vector form as a series expansion h = h 0 + ϵ (1) denotes the ideal image location as a function of H x and H y . Note that H x and H y can describe an infinitely distant as well as a finite object where rays are emitted from an object with coordinates H = (H x , H y , z 0 ). The deviation from the ideal image location is then described by the ray aberration polynomials ϵ (i) = ðϵ x ; ϵ y ; 0Þ as the sum over the transverse ray aberration coefficients ϵ x;j;k;l;m and ϵ y;j;k;l;m ϵ x ðx p ; y p ; H x ; H y Þ ¼ in the x-and y-directions, respectively, where the sum of the indices j þ k þ l þ m ð Þdenotes the aberration order. It is important to remark that the number of independent and non-vanishing aberration coefficients per order depends on the considered symmetry of the freeform optical system. The transverse aberration series expansions are very similar to "standard" aberrations and thus related to Nodal or vectorial aberration theory so that a conversion of aberration coefficients is possible. For example, all ϵ j;k;l;m coefficients with l ¼ m ¼ 0 can be identified as spherical aberrations, coefficients with l ¼ 1, m ¼ 0 or l ¼ 0, m ¼ 1 as comatic aberrations, and so on.
We now consider a sequence of N refractive and/or reflective optical surfaces f i (i =1…N) (in Fig.1 these are three mirrors) aligned along a principal ray path. To describe their surface functions, we use a power series representation with coefficients f i;st : In the case of one plane of symmetry (here the x-z-plane), all coefficients in y with odd t vanish. Next, we introduce the ray mapping functions (u i (x p ; y p ; H x ; H y ), v i (x p ; y p ; H x ; H y )) in the x-and y-directions that describe where an arbitrary ray, described by variables (x p ; y p ; H x ; H y ), will intersect each optical surface f i . The ray mapping functions (u i , v i ) for the ith surface are power series in (x p ; y p ; H x ; H y ) that we write in a similar way as the ray aberration expansions in Eqs. (1) and (2): The chosen principal ray path then determines all vertices of the optical surfaces in a 3D geometry through the series coefficients (u i;0000 , v i;0000 , f i;00 ). As the surfaces may be An arbitrary incident ray with field direction vector H and field angles H x and H y Pupil coordinates of the system Real intersection with image plane Fig. 1 Schematic layout of a typical plane-symmetric freeform optical system, here, consisting of three mirrors. The plane of symmetry is indicated by the white dotted lines on each surface. The ray path of an arbitrary ray is shown including the introduced mathematical description of the initial ray direction and intersections with each mirror and the local transverse ray aberrations ϵ x and ϵ y in the image plane tilted and thus not share one common optical axis, the principal ray path direction can change from surface to surface. We can define the respective orientation (tilts) of all elements, namely the object, the individual surfaces of the optical elements, and the image plane, through appropriate rotation matrices M H , M i and M h , applied to H, pi = (u i , v i , f i ) and h, respectively. Consequently, the linear surface coefficients f i;10 and f i;01 are zero for all surfaces. The ray path from the object to the image plane then consists of (N+1) segments d 1 …d N+1 that describe the intermediate optical path length distances for an arbitrary ray. Vector geometry and the Pythagorean theorem enable us to express these distances weighted by the refractive indices n i,i+1 of the surrounding materials. For example, the distance d i is calculated as i denotes the refractive index of the material between two consecutive surfaces and ||…|| is the Euclidean norm. Let us consider a fixed arbitrary point on the object (finite or infinite object distance), and a fixed but arbitrary point on the second surface p 2 =M 2 • (u 2 , v 2 , f 2 (u 2 , v 2 )). A ray emerging from the object and passing through p 1 towards p 2 must be such that the total optical path length d 1 + d 2 is an extremum according to the modern formulation of Fermat's principle. With points at the boundaries kept fixed, the only remaining variables that can be changed to reach an extremum are u 1 and v 1 at the point Following similar arguments, we can derive two sets of differential equations for all defined distances d i pairwise from object to image space: An optical system consisting of N surfaces is thus fully described by N differential equations D i,x and N differential equations D i,y for i = 1…N, for a given but arbitrary pupil plane. Suppose that all functions h x , h y , u i , v i , and f i are analytic and smooth solutions of the differential equations D i,x and D i,y , then Taylor's theorem implies that the functions must be infinitely differentiable and have power series representations as defined in Eqs. (1)-(5). We can employ a power series method 55 to find solutions to the derived differential equations D i,x and D i,y . This method substitutes the power series of Eqs. (1)-(5) into the differential equations in Eqs. (6) and (7) in order to determine the values of the series coefficients. To calculate certain coefficients, we differentiate Eqs. (6) and (7) with respect to x p ; y p ; H x ; H y evaluated at It is important to remark that Eqs. (8) and (9) correspond to the x-components and the y-components of the mapping functions and ray aberrations in the local image plane.

Deterministic freeform optical design and system evaluation
A preparatory step in the design of a freeform optical system is to specify the layout, number and types of surfaces to be designed and the location of the stop. If we define for example the ith surface as the stop, then (9). The layout of the optical system is defined by the path of the chief ray for the central field. The tilts of optical surfaces are entered using rotation matrices. We can now derive the differential Eqs. (6) and (7) and evaluate Eqs. (8) and (9) for the zero-order case for j = k = l = m = 0. This condition is mathematically equivalent to the defined chief ray path and is as such automatically fulfilled. At this stage we can include optical specifications such as the entrance pupil diameter (ENPD), focal length (FL), field of view (FOV) and design wavelength as required by the targeted application.
With the differential equations established and the overall system specifications introduced, two design steps need to be taken: (1) solve the nonlinear first-order case using a standard nonlinear solver or by making use of equivalent first-order optics tools that will provide structurally similar nonlinear equations 42,43 ; (2) solve the linear systems of equations in ascending order by setting unwanted aberrations to zero or by minimizing a combination thereof as required by the targeted specifications of the imaging freeform system. These two steps are identical for all freeform optical designs and are implemented as follows: 1. We evaluate Eqs. (8) and (9) for all (i, j, k, l, m) for the first-order case with j +k +l +m = 1. This results in a nonlinear system of equations for the second order surface coefficients f i;st with s+t = 2, the mapping coefficients u i;1;0;0;0 , u i;0;0;1;0 ,v i;0;1;0;0 , v i;0;0;0;1 and the aberration coefficients ϵ x;1;0;0;0 , ϵ x;0;0;1;0 ,ϵ y;0;1;0;0 , ϵ y;0;0;0;1 . The latter are four firstorder ray aberration coefficients that can be set to zero or that can be minimized to calculate the unknown surface and the mapping coefficients. Extra conditions for the second order surface coefficients can be imposed as additional equations if desired. In cases where the number of equations will exceed the number of unknowns, the first-order optical powers of some surfaces can serve as additional initial degrees of freedom. The equations can now be solved using a standard nonlinear solver. 2. For each of the higher orders with Ω = j + k +l +m = 2, 3…, we can determine the exact relationship between the surface, mapping and aberration coefficients by evaluating Eqs. (8) and (9) for all corresponding (i, j, k, l, m). This becomes a linear process once the previous order has been solved. Each linear system of order o then relates the surface coefficients f i;st with s+t = Ω+1 to the mapping coefficients u i;jklm , v i;jklm and the aberration coefficients ϵ x;jklm; , ϵ y;jklm with j +k +l +m = Ω. For each order, we apply the elimination method for solving linear systems to eliminate the unknown mapping coefficients and to obtain a reduced linear system that expresses the aberration coefficients as linear equations of the unknown surface coefficients of that order. The resulting linear system can be either squared, overdetermined or underdetermined. In all cases, we use the Moore-Penrose pseudoinverse (MATLAB's pinv function) to solve the linear system for the respective surface coefficients to get the minimum-norm least-squares solution. If we would calculate the least-squares solution for an overdetermined linear system, the aberrations of the considered order would be weighted equally. This, however, would not take the system specifications into account. We therefore define a basic set of weighting factors WF j;k;l;m ¼ ENPD=2 ð Þ jþk FOV= p 2 lþm that is multiplied with the reduced linear equations of same index (j, k, l, m). It is important to stress that an exact but weighted linear least-squares solution 56 is thus very different from commonly used weighting factors in optimization. The weighted least-squares solution for the reduced linear system then takes both the maximum entrance pupil diameter and largest diagonal field diameter into account to simultaneously minimize all properly weighted aberrations for each order, and to calculate the corresponding surface and aberration coefficients. The calculated coefficients are now substituted into the original linear system to calculate the remaining unknown mapping coefficients of that order. The set of surface coefficients f i;st that have as such been obtained in a deterministic way, now fully describe the N reflective and/or refractive optical surfaces. As such they form a satisfactory "first time right" solution of the freeform system-under-design, while the aberrations and mapping coefficients can be used to evaluate the imaging quality and the geometry of the optical system. Alternative solutions can be further explored, for example by introducing different weighting factors to rebalance the corrected aberrations, e.g. in favor of a wider field of view.
In addition, the proposed power series method enables a function-based analytic ray-tracing evaluation by accurately calculating higher order mapping and aberration coefficients in ascending order and for all combinations of j þ k þ l þ m ¼ 1; 2; 3 up to a desired order that is typically higher than that of the preceding surface coefficient calculations. If not stated otherwise, we always calculate the surface coefficients up to the sixth order, and the mapping and aberration coefficients up to the eighth order, and this throughout the paper. Thus, the full physical behavior of the optical system can be immediately interpreted by the optical designer. The obtained "first time right" solution thus provides an initial imaging system that can be further fine-tuned, for example by using a classical ray-tracing software and applying the embedded advanced optimization algorithms. Figure 2 summarizes the overall design process in a flowchart, which highlights the difference between input parameters for the optical designer (upper row in blue) and the calculated series coefficients (lower row in orange).

Two-dimensional design example for illustration purpose
First, we apply the outlined general process of the previous section to design a Cassegrain-type aplanatic onaxis (all rotation angles are zero) two-mirror system (see Fig. 3a). The two distances between mirror 1 and 2, and mirror 2 and the image plane are 250 and 300 mm, respectively. The focal length (f L ) is 2000 mm and the ENPD is 150 mm. Due to the rotational symmetry of the problem, we can calculate the solution in the x-z-plane. As such the given Eqs. (1)-(5) simplify with partially vanishing sums and indices k=m=t=0. The chosen principal ray path determines the constant terms of the polynomial series in Eqs. (3) and (4) to bef 1;00 ¼ u 1;0000 ¼ u 2;0000 ¼ 0 and f 2;00 ¼ À250. The optical path length segments are ; H x Þ and D = 50. Next, we can derive two differential equations D 1;x and D 2;x according to Eq. (6) for the defined path length segments and afterwards replace u 1 by x p to position the stop at the first mirror. The two subsequent calculation steps are executed as follows: 1. We evaluate ∂x j p ∂H l x D i;x (cf. Eq. (8)) for i=1,2 and for j ¼ 1; l ¼ 0 and j ¼ 0; l ¼ 1 at x p ¼ H x ¼ 0, respectively, which results in four nonlinear algebraic equations. By setting the first-order aberrations ϵ x;1;0;0;0 and ϵ x;0;0;1;0 equal to zero, there are an equal number of four unknown coefficients f 1;20 , f 2;20 , u 2;1000 , and u 2;0010 to solve for by using MATLAB's fsolve function. 2. For each higher order Ω ¼ j þ l ¼ 2; 3; 4; 5 we continue to evaluate ∂x j p ∂H l x D i;x for i = 1,2 and for both j ¼ Ω; l ¼ 0 (spherical aberrations) and j ¼ Ω À 1; l ¼ 1 (comatic aberrations) at x p ¼ H x ¼ 0, respectively. Thus, we obtain a constant number of four linear equations for every order. By setting the aberration coefficients ϵ x;Ω;0;0;0 and ϵ x;ðΩÀ1Þ;0;1;0 equal to zero for each order, there is an equal number (squared system) of four unknown coefficients of the surfaces f 1;ðΩþ1Þ0 , f 2;ðΩþ1Þ0 and the mapping function u 2;Ω000 , u 2;ðΩÀ1Þ010 to immediately solve for. As expected from the imposed on-axis symmetry, all polynomial series coefficients for Ω ¼ 2; 4 vanish. The aplanatic solution that we calculated is shown in Fig. 3a and the respective spot diagrams after creating a rotationally symmetric system out of the calculated mirror profiles are shown in Fig. 3b.
As expected from our imposed design conditions, both third-order spherical aberration and third-order coma vanish (calculated from Seidel coefficients using Zemax, the file is made available).

Advanced catoptric and catadioptric freeform design examples
To further illustrate the potential of our deterministic direct optical design method, we apply the "first time right" method to two highly advanced state-of-the-art optical freeform systems that were recently developed by To allow readers to follow and evaluate the design method hands-on and in realtime we have programmed and compiled a graphical user interface (GUI) enhanced user application in C++ for deployment to MATLAB Web Application Server. In the shown screenshot in Fig. 4, partial obscuration is caused by mirror 2, which is directly reflected in a penalty term that is added to the figure of merit function. The Supplementary Information provides all necessary additional information on this user application.
The first reference design example is a catoptric imaging system by Bauer et al. 50 for the visible spectrum with x-zplane symmetry and the object at infinity. It consists of three freeform mirrors with the stop at the first mirror. The targeted system volume is 60 L. The layout of the optical components can be described by three distances (between mirror 1 and 2, mirror 2 and 3, and mirror 3 and the image) and four rotation angles (one for each of the three mirrors, and one for the image), which provides seven degrees of freedom for the optical designer. These distances and angles can now be adjusted to create an unobscured starting geometry that meets certain geometrical constraints such as the given target system volume. This can be done manually and in real time with the user application. Details on geometry selection procedures are given in the Supplementary Information. With the entrance pupil diameter, focal length, field of view and design wavelength specified, we can initiate the approach as described in the previous section: 1. Evaluating Eqs. (8) and (9) for all (j +k +l +m) = 1 and i = 1…3 results in a nonlinear system of 12 nonvanishing equations with 18 unknowns. Setting the four first-order ray aberration coefficients ϵ x;1;0;0;0 , ϵ x;0;0;1;0 , ϵ y;0;1;0;0 , ϵ y;0;0;0;1 to zero leaves eight mapping coefficients u i;1;0;0;0 , u i;0;0;1;0 , v i;0;1;0;0 , v i;0;0;0;1 (i=2,3) and six second order surface coefficients f i;st (s +t = 2) as unknowns. We can define two of the three mirrors to have a base curvature, for example f 1;20 ¼ f 1;02 ¼ c 1 and f 3;20 ¼ f 3;02 ¼ c 3 , reducing the number of unknown surface coefficients from six to four. The nonlinear system can now be solved. The calculated design that we obtained by optimizing the seven degrees of freedom of the geometry (three distances and four angles) is shown in Fig. 5a, where the system layout cross-section is combined with the full 3D peak-to-valley freeform departures (PV) from the best-fit base sphere for each mirror, respectively. Figure 5b shows the spot diagrams for six selected fields based on our aberration calculations. With an average RMS spot radius of about 5 μm, our directly calculated system provides an already well-corrected "first time right" solution that can be readily further optimized. Next, all forty previously calculated surface coefficients and the initial seven degrees of freedom are used as variables for further optimization, e.g. using MATLAB's lsqnonlin solver. After ten iterations, which only take a few minutes a Cross-section of the directly calculated initial system combined with peakto-valley freeform departures (PV) from the base sphere for the primary, secondary and tertiary mirror. b Corresponding spot diagrams for six selected fields based on aberration calculations (blue triangles) and ray tracing (red crosses) in comparison. c Cross-section of the subsequently optimized system combined with peak-to-valley freeform departures (PV) from the base sphere for each mirror. d Corresponding spot diagrams for the same six fields based on aberration calculations (blue triangles) and ray tracing (red crosses) in comparison of calculation, the system already reaches diffractionlimited performance close to the reported performance 50 for almost the full field of view. The optimized design shows slightly increased and moderate freeform departures distributed among the three mirrors. The results are shown in Figs. 5c, d accordingly. The aberration-based performance estimation of our method was found to be in excellent agreement with spot diagram data from classical ray tracing (calculated using Zemax, overlaid red cross symbols in Figs. 5b, d).
It is important to emphasize that these results mainly highlight the very effective nature of this new initial freeform design and evaluation tool. Any further final optimization can be performed by using a suitable multivariable optimizer or a classical ray-tracing software with embedded advanced optimization algorithms. Thus, the directly calculated system solution provides an excellent starting point that can be used for subsequent optimization. By using a non-uniform field weighting, and/or a wave front-error-based optimization merit function for the "first time right" design, a diffraction-limited performance can be readily reached.
As stated in the introduction, increasing the number of calculated surfaces is a major issue for most current direct design methods. As a second example from literature 57 , we therefore have selected a monolithic freeform objective for a very compact infrared camera with four optical surfaces, two of which are refractive and two reflective. We chose this example to illustrate the scalability of our design method, well beyond the capabilities of most present-day direct freeform design approaches. We furthermore highlight its holistic character, since in this case we are dealing with a catadioptric system with freeform and aspherical surfaces. Here we design an objective consisting of three aspheres and one freeform at the second surface, with the stop placed at the first surface. It was originally designed by the company Jenoptik AG Jena and discussed by Kiontke et al 57,58 . An all-freeform redesign of the monolithic objective has been proposed and discussed by Reshidko et al 59 . We follow these references as closely as possible with the following system requirements: an F/1.4 design covering a 37×25-degree field-of-view, an 8.4 mm entrance aperture, made of optical-grade germanium and operating in the long-wave infrared region (LWIR) from 8-12 μm. With the object at infinity, the layout of the system can be described by four distances and five angles, defining the principal ray path from object to image and the respective positions and orientation of all elements and corresponding rotation matrices. It is important to emphasize that the following calculation steps are almost identical to the previous three-mirror example, whereas only the nonlinear first design step is slightly altered: A comparison with the catoptric three surfaces design example makes clear that the only differences with the catadioptric system designed here are the altered set of underlying Eqs. (6)- (9) and the adapted solution of the nonlinear first-order equations. Everything else remains the same no matter what types or how many surfaces are used. Important to remark here is that all orders higher than the first-order result in linear systems of equations that can be fast and very efficiently transformed and solved. This example clearly demonstrates again the highly effective nature of our proposed design and evaluation method. In addition, it highlights two further important features: (1) our method allows us to simply combine refractive and/or reflective surfaces of spherical, aspheric or freeform shapes; (2) the method also straightforwardly enables to increase the number of calculated surfaces of an optical system without considerably increasing the computational complexity of the problem. This direct

Discussion
Equipping lens-and/or mirror-based optical systems with freeform optical surfaces makes it possible to deliver highly original imaging functionalities with superior performance compared to their more traditional (a)spherical counterparts, such as enhanced field-of-view, increased light-collection efficiencies, larger spectral band, and higher compactness. Until now mathematical models and design strategies for freeform optics remained limited and failed to provide deterministic solutions. In particular, the identification of a suitable initial design has often proven to be a painstaking and time-consuming trial-and-error process.
In this paper, we reported on the first deterministic direct design method for freeform optical systems that is not restricted by the aberration terms that can be controlled. The method relies on Fermat's principle and allows a highly systematic generation and evaluation of directly calculated freeform design solutions that can be readily used as starting point for further and final optimization. As such, this new method allows the straightforward generation of "first time right" initial designs that enable a rigorous, extensive and real-time evaluation in solution space when combined with available local or global optimization algorithms.
We can summarize the main features of our new method as follows: 1. Holistic: The method can be used to design imaging systems that include catoptric and/or dioptric spherical, aspherical, and/or freeform surfaces and components. It works for all present aberrations in freeform optical systems up to a desired order. 2. Deterministic: As many current design strategies strongly rely on multi-step optimization routines, final design results often differ considerably for slightly different starting points. The here presented method allows a fast and deterministic solution as well as a very systematic evaluation of the solution space while providing more detailed insights into the fundamental physics of the optical system than most conventional design approaches. 3. Scalable: So far, several direct design methods have been reported for freeform optical systems. However, these methods have in common, that they are developed for a specific system layout, e.g. two or three freeform mirrors. Once a designer would like to add one or more mirrors to the design, the methods do not typically scale with the number of calculated optical surfaces. In contrast, the here presented method follows a clear set of rules that allows us to add additional refractive and/or reflective surfaces to a system. 4. All-round: the method can be beneficial for "junior" and "senior" optical designers alike. The in-depth insights into aberrations with this approach are very valuable for skilled senior optical designers. Inexperienced optical designers with a less advanced understanding of aberration theory on the other hand can equally benefit from it, as this in-depth insight is not a must-have to make excellent use of the method. Finally, the drastic reduction in parameter space, i.e. the large ratio of calculated coefficients over input coefficients, make this method very attractive for designers with a strong background in optimization. The deterministic, holistic, and scalable nature of this all-round method therefore has the full potential to create a true paradigm shift in how freeform imaging systems will be designed and developed from now on to suit a wide range of applications.
Author contributions F.D. conceived the idea and developed the 'first time right' method and algorithms. F.D. and H.T. wrote the paper. Both authors contributed to the discussion and preparation of the manuscript.