Topology optimization of multicomponent optomechanical systems for improved optical performance

The stringent and conflicting requirements imposed on optomechanical instruments and the ever-increasing need for higher resolution and quality imagery demands a tightly integrated system design approach. Our aim is to drive the thermomechanical design of multiple components through the optical performance of the complete system. To this end, we propose a new method combining structural-thermal-optical performance analysis and topology optimization while taking into account both component and system level constraints. A 2D two-mirror example demonstrates that the proposed approach is able to improve the system’s spot size error by 95% compared to uncoupled system optimization while satisfying equivalent constraints.

. They are found in hightech applications in disciplines such as (nano)metrology, astronomy, life sciences, aerospace, lithography and communications. A cutting-edge example is the Sentinel 5 UV1 space-based spectrometer, shown in Fig. 1. The challenge of optomechanics lies in maintaining the position and shape of optical elements such that the image quality is guaranteed under all working environments (e.g. Giesen and Folgering 2003;Meijer et al. 2008;Pijnenburg et al. 2012). Therefore, optomechanical design is a tightly integrated process involving many disciplines such as thermal control, structural mechanics, motion control and optics.
The optical performance generally encompasses the systems' image quality, optical resolution, and image position accuracy. Both image quality and optical resolution are typically quantified by the spot size or image blur diameter (e.g. Welsh 1991;Doyle et al. 2012). The spot size of an aberration-free system is limited by the wavelength of light. This diffraction limit determines the minimum blur diameter that is achievable by an optical system and hence provides a reference for image quality. The spot size mainly depends on the wavefront quality due to geometric aberrations, whereas the beam position accuracy depends on how accurate the optical components are positioned/oriented. These performance metrics can be determined using geometric ray tracing, which traces the (a) Structural design of the Sentinel 5 UV1 spectrometer. The mirrors are bolted to the housing, and the housing is semi-kinematically connected to the satellite frame via flexures in order to limit deformations due to imposed rigid body displacements of the frame.
(b) Optical design of the Sentinel 5 UV1 illustrating the mirror mounts and optical rays. Note the different mirror back layouts and consistent connections to the frame. The spot position and size depends on the misalignments and surface deformations of each mirror and grating.

Fig. 1
The European Space Agency's (ESA) Sentinel 5 shortwave UV1 band (270-300 nm) spectrometer system to be carried on the MetOp-SG satellite propagation of light rays through an optical system (e.g. Spencer and Murty 1962). For flat or spherical single-mirror systems the wavefront error scales linearly with the Surface Form Error (SFE) of the deformed surface and the pointing error is directly related to the tilt of the surface (e.g. Genberg 1984Genberg , 1999. The optical performance metric of interest thus depends on application and system composition. Many factors may contribute to the inability of an optical system to produce a perfect image, including chromatic and geometrical aberrations, fabrication and alignment errors, (lack of) self-weight and environmental effects such as temperature fluctuations. This work focuses on the reduction of optical performance errors of reflective optical systems induced by (quasi)-static thermal loads. Therefore, structural deformations and temperature differences of the frame should neither excessively distort the mount nor the optical surface. This implies that a mechanically disconnected frame and optical surface combination would be optimal. However, the presence of thermal loads requires material to abduct the heat from the optical surfaces to the frame. In addition, the optical components also require a stiff design, as the structure must constrain the components such that they are not damaged or irreversibly moved after exposure to external conditions such as vibrations, thermal shocks and gravity. To limit excitation from external vibrations, the fundamental elastic eigenfrequency must be higher than a critical lower limit. Adding more practical constraints such as maximum mass and material usage, often linked to costs, the optomechanical mirror support design clearly involves multiple conflicting structural requirements. To make well-founded and justified design trade-offs, careful consideration of the thermomechanical and dynamic performance of optical mounts is required. Optimization techniques can aid in this process.

Problem definition
The current typical design approach of optomechanical instruments is characterized by the optical discipline creating a performance error budget that defines deformation limits for each optical component to the structural discipline. In an iterative design process, the thermal discipline provides the temperature fields and gradients, after which the structural discipline aims to realize a design which meets the deformation limits during thermal and other environmental load cases. From this point of view, an optical system functions as long as the components remain within allowed tolerances of their nominal locations, orientations and deformations. Thus, in the existing thermomechanical design process, the optical performance is not considered directly. Instead, each component is designed and optimized separately to meet a priori defined deformation limits.
To expand the design space and enable further improvements of the optical performance we introduce: 1. system-level optical performance metrics that drive the thermomechanical design and optimization process, 2. simultaneous optimization of all components for the combined metrics, and 3. system-level constraints that, where possible and applicable, replace component-level constraints.

Approach: simultaneous multicomponent multidisciplinary topology optimization
Model-based structural optimization techniques can aid in further improving the performance of optical instruments. Topology optimization is a systematic, bottom-up structural optimization approach that provides maximum design freedom without any prior knowledge of the design. The procedure optimizes the material layout within given design domains in order to maximize a performance measure, while subjected to a given set of loads, boundary conditions, and constraints (e.g. Bendsøe and Sigmund 2003). The method, in combination with mathematical programming, has shown to be able to solve complex multidisciplinary problems with multiple active nonlinear response functions and is capable of producing innovative solutions. The mentioned design requirements have previously been investigated in topology optimization. This includes thermal loads and prescribed temperature differences on stiffness problems (e.g. Rodrigues and Fernandes 1995;Li et al. 2001a;Gao and Zhang 2010;Deaton and Grandhi 2013) and the application to the design of thermomechanical compliant mechanisms (e.g. Sigmund 2001;Ansola et al. 2012). Recently, Zhu et al. (2016) presented a shape preserving design method for topology optimization. By suppressing the elastic strain energy in the local domain, the shape of the concerned domain can be effectively maintained to satisfy the requirements. This is extended by Li et al. (2017) to achieve desired deformation behavior within local structural domains by distinguishing and suppressing specific deformation in a certain direction. The fundamental eigenfrequency as a response function has been thoroughly investigated (e.g. Ma et al. 1995;Pedersen 2000;Du and Olhoff 2007;Tsai and Cheng 2013).
Previously, topology optimization has shown its benefits in the design of mirror mounts for optical performance. For example the minimization of deformations of optical surfaces under static Sahu et al. 2017) and thermal loads (Kim et al. 2005) and mass minimization while constraining deformations as well as the fundamental eigenfrequency (Hu et al. 2017a). Furthermore, semi-kinematic flexible mirror mounts are topologically optimized by, Hu et al. (2017b) who minimize surface form errors subjected to both static and thermal loads constrained by a minimum natural eigenfrequency, and Van der Kolk et al. (2017), who achieve optimal damping characteristics at source frequencies. System optical performance metrics have not been included in the topology optimization framework yet.
An integrated Multidisciplinary Design Optimization (MDO) approach is often applied to couple all involved physics and profit from the interactions between different disciplines, resulting in superior designs. We apply an integrated Structural-Thermal-Optical Performance (STOP) analysis and optimization procedure for this purpose to utilize simultaneous optimization of the optical, structural and thermal design aspects (e.g. Johnston et al. 2004;Doyle et al. 2012;Kuisl et al. 2016). Prior work using STOP optimization couples various analysis tools to obtain the optical performance. This shows great improvement, although the optimization often includes only a limited number of variables such as facesheet thickness and strut diameters (e.g. Williams et al. 1999;Michels et al. 2005;Bonin and McMaster 2007). The results indicate that device performance can profit from integration of optical knowledge, and that further improvement is accessible when all design parameters are considered. The STOP optimization procedure includes optical knowledge at the thermomechanical design level, but the design freedom is not yet fully exploited.
Coupled multicomponent topology optimization can aid in exploiting the component interactions. Simultaneous topology optimization of multiple components has mainly focused on layout design and combined topology and joint location optimization (e.g. Chickermane and Gea 1997;Li et al. 2001b;Zhu et al. 2009Zhu et al. , 2015. Topology optimization involving component interactions to improve a system performance has previously been investigated by Jin et al. (2016Jin et al. ( , 2017, with the focus on simultaneous optimization of multiple coupled actuator mechanisms to minimize the coupling interaction. However, this study was restricted to a single physical discipline. In this work, we combine the foregoing and additionally apply multidisciplinary optimization to a component assembly and optimize for a system-level optical performance metric. This approach will be referred to as integrated System Design Optimization (SDO). To determine the system optical performance, this study uses a simplified version of geometric ray tracing, the ray transfer matrix analysis. The method uses the paraxial approximation to construct a linear operator that describes the behavior of an optical system (e.g. Nazarathy et al. 1986;Smith 2007;Fischer et al. 2008).
The innovation point of the SDO method compared to existing methods is threefold: it combines topology optimization and a full STOP analysis, uses system-level optical performance metrics to drive the thermomechanical design of multiple components simultaneously, and uses system-level constraints to replace multiple equivalent component-level constraints.
We hypothesize that an integrated structural-thermaloptical thermomechanical design optimization procedure taking into account all system components improves the system optical performance compared to individual component optimization, while subjected to equivalent design constraints.
The formulation of the coupled thermomechanical discretized equilibrium equations and modal analysis, topology optimization formulation and sensitivity analysis of a generalized response function are described in Section 2. The extension to a full STOP topology optimization approach is given in Section 3, which focuses on the optical performance prediction from finite element analysis results. The method is first validated on a singlecomponent system, after which the hypothesis is tested by numerical optimization of a two-mirror example as discussed in Section 5. The results are followed by a discussion, recommendations and conclusions presented in Sections 6 and 7.

Formulation of coupled thermomechanical analysis framework for topology optimization
A schematic of a typical two-mirror opto-thermomechanical system is shown in Fig. 2, introducing the necessary nomenclature for the mathematical model. Both domains consist of a linear-elastic homogeneous isotropic material with conductivity k i , Young's modulus E i , Poisson's ratio ν i , density ρ i and Coefficient of Thermal Expansion (CTE) α i . A coordinate system is asserted for every optical component to express the displacements of the optical surfaces. Additionally, each optical element (including space propagation in a medium with constant refractive index) has a coordinate system to track the position and angle of rays with respect to the optical axis. The optical design includes the optical path lengths d, incident angles φ and initial ray position and angles. The response of a system to a given load case depends on the optical design and properties, as well as the component's domain geometries, boundary conditions, structural layouts and material properties. Any system response function depends on the response of multiple components to a given load case, where the response of each domain directly depends on the layout of material, denoted by the design variables s, or indirectly through other responses.
The design variables following from density based topology optimization, one belonging to every finite element in the domain, are bounded by a lower and upper bound, i.e. 0 < s ≤ s i ≤ 1 with i = 1, 2, ...,n, wherê n is the number of variables. The underbound s has a very small value (to avoid numerical issues) denoting the absence of material and,s = 1 providing the element with the assigned initial material properties. Each design variable can have any intermediate value within the given bounds. The material properties are interpolated using a penalization function as discussed in Section 2.1. Sections 2.2 to 2.4 is shown in grey, with nondesign domain illustrated in black. Prescribed displacements U b and temperatures T d , respectively, represent the housing rigid body effects and temperatures. The beam imposes thermal loads Q c onto the mirror surfaces. The position and deformations of mirror i are locally described using coordinate systems O m,i . An incoming wavefront is described by multiple rays, of which each rays' position and the angle is expressed in local coordinate system O r,i with respect to the principal ray incoming on component i. Each mirror makes an angle φ i with relation to the incoming principal ray. The initial optical path lengths of principal ray i is given by d i discuss the thermomechanical and modal analysis, followed by the respective sensitivity analyses.

Penalization scheme
Intermediate density values are implicitly penalized to gradually force the design variables to approach their bounds and facilitate interpretation and improve manufacturability. Most material properties are a function of the design variables, which varies depending on the applied penalization scheme. Common interpolation functions for the Young's modulus R E (s) and conductivity R k (s), are the SIMP and RAMP functions (Bendsøe and Sigmund 1999).
Low-frequency eigenmodes localized in the void regions of the structure should be avoided. Therefore, this study applies a continuously differentiable penalization scheme for R E (s) and R k (s), as proposed by Zhu et al. (2009). The function extends the SIMP interpolation with a linear term to avoid large ratios of R ρ R E for design variables with near zero densities, that is where p is a penalization parameter and q is a weight factor controlling the relative influence. This study consistently uses a power p = 3 and weight q = 0.8. The densities are scaled linearly, thus the volume and mass of an element depends linearly on the accompanying design variable. The method can effectively penalize intermediate densities and avoid spurious eigenmodes for low-density elements while keeping the penalization factor relatively low. The stress is given by where m is the strain from applied forces, T the applied temperature difference and E(s)α T the resulting thermal stress. The thermal stress is a nonlinear function of the design variables via the penalization of the Young's modulus, but the CTE itself is not directly penalized.

Thermomechanical equilibrium equations
The one-way coupled thermomechanical equilibrium equations, neglecting heat transfer by radiation and convection and discretized using finite elements, are given by where K(s) is the global stiffness matrix, H(s) the global conductivity matrix and A(s) the global thermal expansion matrix converting the temperatures to equivalent thermal loads (e.g. Cook 1981). These global matrices are a function of the material layout. Depending on the considered loading conditions, the forces F(s) and applied heat flux Q(s) can become design dependent. The material properties are assumed to be independent of strain and temperature, thus both linear geometric and linear material models are applied. Furthermore, the loads are assumed to be independent of time, thus only steady state heat transfer and (quasi)-static forces apply. Prescribed and unconstrained Degrees of Freedom (DOFs) sets are generally different for the mechanical and thermal problem. To express this in a convenient notation, one can separate the DOFs in terms of unconstrained and prescribed displacement DOFs a and b, as well as unconstrained and prescribed temperature DOFs c and d. Each domain can be a function of externally applied forces F a (s) and heat loads Q c (s), prescribed displacements U b (s) and temperatures T d (s), unconstrained displacements U a (s) and temperatures T c (s) as well as the reaction forces F b (s) and reaction heat loads Q d (s). For the sake of clarity we abbreviate the multiplications and leave out the design variable dependency, such that for instance, For a single domain one can now write the partitioned discretized heat transfer and elasticity equilibrium equations as The equilibrium equations, (5), are solved to obtain both temperature and displacement fields.

Modal analysis and mean eigenvalue
The natural dynamic properties (neglecting damping) of a structure are found by solving the eigenvalue problem for each unknown eigenvalue λ i and corresponding mode shape φ i . Here i is the mode of interest and n the number of DOFs. In this eigensystem the characteristic value equals the squared radial eigenfrequency, i.e. λ i = ω 2 i . All modes are mass normalized satisfying The mean eigenvalue (Ma et al. 1995) is used as a response function to maximize or constrain the minimum fundamental frequency. The function combines them lowest eigenfrequencies of the structure, calculated as This form ensures that the lowest resonant frequency is the main contributor and the influence of modes with a higher eigenfrequency quickly decreases. This implants a simple and effective natural eigenfrequency constraint and ensures the absence of mode switching.

Sensitivity analysis
For any response f (s), the gradient with respect to s must be determined to enable gradient-based optimization. The focus is first on the responses that depend on the thermomechanical analysis. Adjoint sensitivity analysis is efficient due to the large number of design variables, and therefore we augment the response function as where λ a and λ c are the adjoint vectors related to the mechanical and thermal equilibrium equations. The total derivative with respect to the design variables, using the adjoint method to solve the gradient of the objective function in relation to the variables, is given by which holds if the Lagrange multipliers satisfy: Here, loads and boundary conditions are assumed to be design independent and the responses are assumed to solely depend on the displacement field. Then, (10) can be simplified to with the Lagrange multipliers determined by The sensitivities of the mean eigenvalue, as described in (8), with respect to the design variables are The sensitivities of the global matrices can be derived element wise using direct differentiation. With this, the sensitivities have been determined up until the term ∂f ∂U a , which will be discussed in Section 3, after introducing the considered optimization problem and response functions.

Optical performance metrics and sensitivities
This section describes various optical performance metrics relevant for the analysis of reflective optical systems. First, the SFE response will be discussed, after which the analysis of average positional accuracy and spot size are explained. Most commonly the SFE is expressed by the Root Mean Square Error (RMSE) of the deformed configuration compared to the undeformed or another predefined configuration (Genberg 1984). For 3D unit disk surfaces the surface errors are often expressed in Zernike polynomials, which are directly related to typical optical aberrations. For diffraction limited flat or spherical singlemirror systems there exists a simple relation between the RMSE and the Strehl ratio. This is the peak aberrated image intensity compared to the maximum attainable intensity using an unaberrated system. The wavefront is proportional to surface front error. Though, for complex mirrors or multimirror systems the WFE is not directly related to the SFEs and the image quality can only be determined by ray tracing techniques.
To analyze a multi-component system without using numerical ray tracing, the deformed surfaces can be approximated by a fit to directly obtain contributions to the optical surface misalignments, i.e. rigid body movements and SFEs. Next, the ray transfer matrix analysis can be used to track the position and angle of a paraxial ray though a multi-component system, leading to a measure for the averaged positional accuracy. In order to track the rays, all system properties, i.e. optical paths lengths, incident angles and specific properties of the optical components, as well as component transfer functions and misalignments should be known. Finally, the spot size is quantified by the mean average deviation of all rays with respect to the averaged spot position.

Surface form error
The RMSE is the standard deviation of the distance between the deformed surface and the ideal surface and is a global measure for the SFE, see Fig. 3. The deformed surface is constructed from the out of plane displacements, stored in a Fig. 3 Schematic representation of a generalized surface in undeformed (blue) and deformed (green) configuration of mirror i. Out of plane node displacements are indicated by U z,j for the deformed and byÛ z,j for the fitted displacements, where j is the surface node number. A least-squares polynomial fit through the deformed surface (purple) is used to extract the surface misalignments: the tangential and out of plane misalignments, δ x and δ z as well as the rotational misalignment θ y and curvature κ. The component specific properties and misalignments δ cause the input ray r i , with incident angle φ i , to change into r i+1 vector U z , which is fitted to a smooth surface represented by the fitted out of plane displacementsÛ z using a linear least squares regression scheme (Lay 2006), such that where G is the fit matrix. The fit matrix is constructed as where V is the Vandermonde matrix and Y relates the out of plane displacements U z directly to the coefficients of the least-squares solution, i.e. b = YU z . Vector b contains the coefficients of the least-squares fit (size depends on the order of the fitted polynomial) which will be used in Section 3.2. The Vandermonde matrix consists of the terms of a geometric progression and evaluates a polynomial at a set of points, the surface nodal x-coordinates in this case. The tangential displacements are assumed negligible compared to the out of plane displacement, i.e. ||U x || ||U z ||, which ensures V is constant during the optimization and the least squares fit is linear.
The fitted displacements are related to the coefficients viaÛ z = Vb. The nodal residuals of the linear least squares fit, denoted by R, are calculated as where I is the identity matrix. The surface RMSE equals the absolute residual between the deformed surface and the ideal surface. It is preferred to use the Mean Square Error (MSE) when using the response in gradients based optimization to avoid division by zero when taking derivatives. The MSE is defined as whereñ is the total number of nodes on the mirror face, U z,i is the observed andÛ z,i the fitted out of plane displacement at nodal point i.
The sensitivities of the MSE in relation to the free displacements U a , as required in (13) where J R (U a ) is the Jacobian matrix of R with respect to U a and L z selects the appropriate DOFs of the surface of interest from the vector of free displacements, i.e. U z = L z U a .

Positional accuracy
Depending on the application, it is often essential that the image is kept within certain bounds (e.g. within the boundaries of a sensor). To determine the location of a light ray on the image plane we use paraxial ray tracing of multiple rays. Considering a situation where all optical components are symmetric around the optical axis, the positional error of ray j (i.e. the distance of ray j with respect to the optical axis on the image plane), here denoted by ε j , depends on the radial distance and angle of the ray with respect to the optical axis when entering the system, which will be denoted by vector r 0 . Furthermore, it depends on all misalignments δ 1 , ..., δ N of all reflective optics (N is the number of components) and system specific constant parameters p (initial optical path lengths d and angles of incidence φ), thus ε j r 0,j , p, δ 1 , ..., δ N . The lower order misalignments of a surface i in 2D, are the change in curvature κ i , the axial displacement δ z,i (despace), the rotational misalignment θ y,i (tip/tilt), and the radial displacement δ x,i (decenter). The decenter is directly calculated from the tangential displacement of the surface vertex. Other misalignments can be derived from the coefficients of the surface fit b i , which is calculated by the linear least square regression in (16) and shown in Fig. 3.
The radius of curvature R i = 1 κ i is assumed to be constant (parabolic) over the surface for small angular misalignments, that is dz dx i 1, and defined as the reciprocal of the curvature κ i , which equals κ i ≈ d 2 z Determination of the despace, tip/tilt and radius of curvature from the surface fit do not take into account the radial displacement distribution nor the average radial displacement (decenter) of the mirror surface.

Ray transfer matrix analysis
Ray transfer matrix analysis is used to determine the light path through a system based on paraxial approximations by transforming the vector representing the ray, with the appropriate component transfer matrices, which depend on the properties and misalignments of the component. A general misaligned paraxial transformation for component i, such as the mirror in Fig. 3, is denoted as where r i = [x i dx dz i ] T , with x i and dx dz i the distance and angle with respect to the optical axis i in the undeformed configuration. The ray vector after component i, r i+1 , is a linear transformation to the incoming ray r i via the component specific transfer matrix M i , and the influence of the misalignments as contained in the misalignment matrix E i .
Using first-order optical canonical operator theory (e.g. Nazarathy et al. 1986) one can determine the influence of the misalignments of an optical component δ i on the Table 1 General misaligned ray transfer matrices for generalized misaligned space propagation and mirror surfaces (Yuan et al. 2011). Here R i is the radius of curvature, δ z,i and δ x,i are despace and decenter misalignments, θ y,i the tip/tilt contribution and d i the optical path lengths of a principal ray i with angle of incidence φ i on component i position and angle of an incoming ray r i (Yuan et al. 2011). The generalized transfer matrices for the elements used in this study (space propagation and reflective optics) are shown in Table 1. The space propagation transformation matrix does not contain aberrations and changes of refraction index between optical components. The mirror transformation matrix models reflective optics as perfectly reflecting surfaces without scattering, transmission or absorption.
The positional error ε j afterÑ sequential elements with accompanying misalignments E i , 1 considering an initial ray vector r 0,j , equals the distance from the optical axis after the last component, that is ε j = xÑ . The output ray vector equals where the system transformation matrix and misalignment vector are Thus, the positional error of each ray is a function of all system and component optical properties and misalignments.

Averaged positional error and sensitivities
The squared positional error of a single reflective optic, for the purpose of an uncoupled optimization, is solely a function of the misalignments of that specific component, that is whereε i is the average positional error due to misalignments of component i and the number of rays in the system 1 The number of optical elements in a system is generally larger than the number of optical components, since space propagation is also an optical element, i.e.Ñ ≥ N . In most optomechanical applications N = 2N + 1.
m. For the same reason as in (18), the root is omitted. Statistically, it is unlikely that all positional error contributions have the same direction and hence superposition of the errors will lead to overdesign of the individual components, i.e. the errors are uncorrelated. Therefore, the positional errors of independent sources are combined via the Root Sum Square (RSS). The squared positional error taking into account all components is defined as The sensitivities of the positional error contribution of component i in relation to the unconstrained displacements are and the sensitivities of the positional error using an integrated approach equals In both (24) and (26) Note that ∂fε ∂U a,i U a,1 , ..., U a,N , thus the sensitivities of the positional error with respect to the displacements depends also on the displacements of other components in the system. In contrary, the positional error due to component i only depends on that components' displacements, that is The sensitivities of the misalignments δ i are defined as where δ * i = [R i θ y,i δ z,i ] T . Corresponding sensitivity matrix ∂δ * i ∂b i is a transformation scaling the fit coefficients to the misalignments, this is defined as wherek is the order of the fitted polynomial. For a third order polynomial (parabola)k = 3 and where κ i is considered constant over the surface. The order of the fit determines what type of aberrations can be accounted for. The term ∂b i ∂U a,i = Y i L z,i can be derived from (16) and

Spot size and sensitivities
In order to measure the spot size, ray tracing is performed for multiple rays with different initial distance and angle from the optical axis, see Fig. 2. The average resulting deviation from the averaged positional error is a measure for the spot size. The squared spot size due to the misalignments of component i is defined as where ε j andε i are calculated according to (21) and (23). For an integrated system optimization the spot size response is calculated by whereε is calculated with (24). The sensitivities of the MSE spot size due to the misalignments of component i, with respect to the unconstrained displacements, as required to solve for the Lagrange multipliers in (11), are where ∂ε i ∂U a,i is defined in (25) and ∂ε j ∂U a,i is calculated using (27). Similarly, the sensitivities of the MSE spot size taking into account all components' misalignments are where ∂ε ∂U a,i is defined in (26).

Numerical implementation
This section describes practical considerations of the implementation. The initial conditions are set such that designs are initialized with a uniform density field that exactly satisfies the volume constraint. The optimization problem is solved using the Method of Moving Asymptotes (MMA) (Svanberg 1987). The optical performance measures are relatively sensitive to design changes. Therefore, the algorithm is set more conservative (the move limit move is set to 0.1) in order to avoid large jumps in the design space. The optimization is subjected to termination criteria to avoid endless optimization and is considered to be terminated when the design variables and objective function change less than a threshold and all constraints are met, that is ⎡ wheren is the number of design variables, k is the iteration number and i = 1, 2, ...,m withm the number of constraints.
In order to limit the design complexity and to avoid mesh dependency and checkerboard patterns, a general mesh and element-type independent linear spatial filter is implemented (Bruns and Tortorelli 2001). The simplest spatial filter, as used in the presented study, is the linear filter with weights according to where r is the radius of the filter, taken equal to the size of a single element. The radius r ij is the distance between the centroids of elements i and j andn is the number of active elements (equal to the number of design variables). The filtered variables are calculated aŝ wheres are the design variables updated by the optimization algorithm. The filter does not take account of the element volumes as only structured meshes are used.
To further stimulate fully black-and-white designs, the filtered design variables are projected in the direction of their lower and upper bounds by a smooth Heaviside projection function (e.g. Guest et al. 2004;Sigmund 2007;Xu et al. 2010). The projected design variables are penalized using a projection parameter β around a threshold η, given by All given values are constant during the optimization and no continuation strategy is used.

Results
This section discusses two case studies applying the foregoing theory and demonstrating the validity of the proposed method. First, the focus is on the optimization of a single mirror mount to minimize SFEs. Next, the proposed method including the STOP analysis is tested on a two-mirror case.

Single-component surface form error minimization
The focus of this study is on the optimization of a flat mirror mount design subjected to a uniform temperature increase. The goal is to minimize the resulting spotsize due do both boundary and loading conditions. For a singlemirror system the spotsize is directly related to the SFE of the mirror surface, hence minimizing the SFE will suffice. The aim is to show that the proposed method works well for single component problems under various conditions. Both the boundary conditions and eigenfrequency constraint play a dominant role in the possibility to improve the optical performance. The study verifies the single-component topology optimization procedure and investigates the influence of the eigenfrequency constraint on the resulting topology and performance.

Problem definition
The design domain , as shown in Fig. 4, is subjected to an overall temperature increase T, causing the domain to expand and deform. As the temperature field is known, there is no need to solve the heat equation. It is assumed the housing experiences the same temperature increase as the considered design domain and that the design domain is mounted to a housing with different material properties. The resulting change in expansion is introduced in the design domain by application of known prescribed displacements on both fixtures. The objective is to minimize the MSE of the optical surface due to the specified thermal environment, while constrained by a minimum mean eigenfrequency and a maximum volume. Therefore, the mean eigenvalue (8) is Fig. 4 Optical mount design domain of width w, height h and thickness t with CTE α subjected to a uniform temperature increase T. All non design space is indicated in black. Two regions with prescribed displacements are modelled as the interface with a rigid housing structure, which has CTE α 2 , and thus known prescribed displacements U b , which will induce SFE on surface due to the boundary conditions and hence degrade the image quality of the incoming wavefront adopted as a constraint such that the mean eigenvalue must be higher or equal than the minimum elastic mean eigenfrequency ω 2 n . A parameter sweep is performed over a range of minimum eigenfrequency constraints to investigate the influence of the eigenfrequency constraint. The range spans from minimum eigenfrequency constraints where the constraint is inactive up to values where the structure is unable to satisfy the constraint. Additionally, a constant volume constraint is added, in order to ensure a fair comparison.
The problem is formally defined as In this studym = 6 for all cases, to prevent the occurrence of mode switching. Figure 5 shows the resulting topologies and accompanying performances obtained after optimization of (39) for a range of eleven different values for f n (note that ω 2 n = 4π 2 f 2 n ). The required minimum eigenfrequency determines to what degree the optimizer is able to minimize the MSE. Cases with a higher eigenfrequency constraint generally result in  (39) showing topologies (thermal deformations scaled by a factor 100) and RMS surface form error, in μmK −1 , for a range of minimum eigenfrequency constraints f n higher SFE, as a compromise must be made. Note that above a critical frequency the resulting designs do not satisfy the eigenfrequency constraint. The resulting RMS SFEs are given in μm K −1 because the RMS SFE scales linearly with respect to both T and the CTE α, thus the imposed temperature difference or CTE is irrelevant for the final topology.

Optimization results
The topology of the optimal design subjected to a mean eigenfrequency constraint of 500 Hz is shown in more detail in Fig. 6. Designs with a relatively low eigenfrequency constraint tend to possess compliant mechanism-like structures in order to counteract surface deformations. In general, structures with a lower required dynamic stiffness have clear rotation points (shown in red) and thicker beams to support the mirror. The large amount of material underneath the mirror both stiffens the structure and forces the flexure-like structures, which are fixed to the frame (green), to bend outwards effectively flattening the surface (blue).

Two-mirror system spotsize minimization
This example studies the topology optimization of a twomirror system subjected to thermal loads from a light source, as used in for example high-power laser or EUV applications. The study compares the following design approaches: 1. Uncoupled System (US) optimization, and 2. Coupled System (CS) optimization, where the integrated SDO approach is applied.
Both optimization procedures make use of a full STOP analysis, though only the integrated approach considers the component interactions and applies both system and component level constraints. Note that the uncoupled optimization problem is not a typical design approach used in practice, since it does consider system optical performances. The separate components are however artificially decoupled with Fig. 6 Detailed view of the resulting topology of optimization problem (39) with an eigenfrequency constraint of 500 Hz both in undeformed and deformed (scaled by a factor 4). The figure illustrates the deformed compliant mechanism-like structure with rotation points indicated in red, mirror surface in blue and frame interface in green. Bottom figure shows the first modeshape of the structure. This is the bending mode, as expected relation to the optical performance, in order to investigate the difference with respect to the coupled SDO approach.
The optimization aims to minimize the spot size error due to prescribed boundary conditions of the frame and thermal loads from the propagating beam, while bounded by a the position accuracy, mass, and eigenfrequency constraints. For multi-component systems both position accuracy and spotsize depend on the rigid body motions and SFEs of all components involved. The main target is to investigate whether, and how, the components make use of the capability to interact and compensate for each others optical aberrations and how this may benefit the optical performance.

Problem definition
Consider the schematic structural-thermal-optical system consisting of two flat reflecting surfaces supported by optical mounts with design domains as shown in Fig. 2. An incoming converging beam with a perfect wavefront is reflected by two mirrors before it is focused onto a sensor with a theoretical spot size of zero (if it were unconstrained by the diffraction limit).
The fundamental maximum resolution or minimum spot size of any optical system is limited by the diffraction limit. The diffraction limit for the system as shown in Fig. 2, equals where D is the first minimum in the Airy disk, NA the Numerical Aperture of the system, λ l the wavelength of the light, n r the index of refraction of the medium and θ i the half-angle of the incoming light. Assuming perfect vacuum conditions, i.e. n r = 1, sunlight of λ l = 550 nm (center of sunlight wavelength spectrum) the half-angle of the system of 3.8 • and a NA of 0.0665, the diffraction limit of this system equals 8.3 μm.
Both mirrors are subjected to known rigid body movements from the housing, which is modelled as a rigid interface. The interfaces of the first mirror mount are considered constant at 1 K difference with respect to ambient conditions, whereas the interfaces of the second mirror are constant at 0 K difference. The first mirror is subjected to a decenter rigid body effect of δ x,1 = 200 μm and a the same amount of despace misalignment. The left side of the second mirror is moved out by 200 μm and down the same magnitude. The right side of the second mirror interface is also moved out by 200 μm. This causes the second mirror to initially have a despace and tip/tilt error.
The heat load is modelled as a Gaussian profile over the surface, i.e.
where Q 0,i is the maximum amplitude in domain i in Watts, x is the location on the surface, with x = 0 m in the middle of the mirror surface, μ = 0 m and σ taken equal to 0.1 m. The input heat loads are normalized in relation to the maximum value at x = μ. The first mirror is subjected to a Gaussian heat profile with a maximum input of Q 0,1 = 0.1 W m −1 . Assuming the first mirror absorbs 10% of the heat load, the second mirror is subjected to a heat load with a 90% maximum amplitude in relation to the first mirror. Each design domain is constrained by an individual eigenfrequency constraint, as well as a maximum allowable RMS SFE compared to a perfectly parabolic mirror. Therefore, the MSE SFE response (18), is adopted as a constraint, taking into account proper scaling, i.e.
Additionally, the system is required to keep the position of the image within a certain limit with respect to the optical axis, therefore the positional error, (23) or (24), depending on the optimization problem, has to be adopted as a constraint, such that gε ,i (U a,i ) = log 10 ε 2 i ≤ − log 10 ε 2 i for an uncoupled optimization, and gε(U a,i , ..., U a,N ) = log 10 ε 2 ≤ − log 10 ε 2 (44) for the coupled case. For an uncoupled optimization, the system positional error budget must be split up into the individual components, in this case such that the RSS value equals the total allowed system positional error, with equal weights per mirror, i.e. the positional errorε of the coupled case equals the RSS of the positional tolerancesε 1 andε 2 .
In the uncoupled optimization each mirror is subjected to an individual maximum volume constraint as is typical in industrial projects, however, when performing integrated SDO the full system is subjected to a maximum volume Table 2 Performance of the obtained optomechanical systems, and properties of their individual mirror mounts: RMS spot size, RMS positional error, volume, RMS SFE, mean eigenfrequency rigid body movements and mirror curvature. The values between parentheses indicate the value at the first iteration constraint equal to the sum of the volumes of the individual components. In the example 50% volume is allowed. This allows the optimizer to move material between domains for the benefit of the system performance. The constraint tolerance values, redefined into more relevant measures (e.g. f n instead of ω 2 n ), are given in Table 2. The uncoupled optimization problem of the individual components minimizes the spot size as a function of the misalignments of the component of interest while satisfying the volume, positional error contribution and RMS SFE constraints. The US optimization is a combination of the optimization of mirror one (U-M1), not taking into account the misalignments of the second mirror (U-M2) and vice versa. The problem is stated as for i = 1, ...N, where N = 2 for this system. In the coupled case, the integrated SDO approach simultaneously optimizes all domains, thus objective function, system positional error and volume constraint are a function of the layout of all domains in the system, that is i M φ i = 0 for i = 1, 2, ..., n gε(U a,1 , U a,2 ) = log 10 ε 2 + log 10 ε 2 ≤ 0 for i = 1, 2. The lower bound on the design variables is s = 10 −3 . The density filter radius equals the length of two finite elements and the Heaviside projection parameters are set to β = 1.5 and η = 0.45. The termination criteria are ε s = 0.0015, ε f(s) = 0.002, and ε g(s) = 0.02. Figure 8 shows the resulting RMS spot size diameter as a function of iteration history for both optimization problems. The RMS spot size diameter versus the number of iterations for the coupled case using the SDO approach is illustrated in red. Whereas this approach requires only a single objective function, the uncoupled optimization consists of two separated optimization problems (one for each optical component), of which the RMS spot size diameter versus iteration number are shown in green and purple. In order to compare the system performances the RMS spot size diameter as a function of iteration history of the uncoupled system is calculated afterwards as shown in blue. Note that the optimizer is unaware of the US performance during the optimization process. The resulting topologies for both problems are displayed in Fig. 7, accompanied by the final performances as well as the system and component properties and constraint tolerances shown in Table 2.

Optimization results
The two approaches result in some significant differences, see Table 2 and Fig. 7. Whereas the curvature of the first mirror in both approaches is brought as close as possible to zero, the second mirror in the coupled approach is made even more concave than its original shape. Note that the spot size of US is not simply the average of both mirrors. The positional error constraint of U-M2 is active, although the combined positional error of US is far below the allowable error. Both systems have satisfied the allowable mass, however, the integrated SDO has transferred mass from C-M2 to C-M1.
The optimization process has lowered the RMS spot size diameter from 0.4 mm in the first iteration to 106.5 μm and 4.7 μm for the uncoupled and coupled optimization approaches, respectively. Thus, the integrated SDO approach used for the coupled case proves a ratio of increase in performance of 22.6 times compared to the uncoupled method, equivalent to an improvement of 95.6%. Note that the integrated SDO approach is able to lower the RMS spot size to below the diffraction limit of the system, that means it is not possible to further improve the optical performance measure with a geometry-based performance metric.

Discussion and recommendations
Existing approaches for topology optimization of optomechanical systems focus on the optimization of individual mirror mounts to minimize their surface deformations. Our research extends this by 1. topology optimization of mirror mounts for the systems' optical performance expressed in terms of spot size and position accuracy using a full structuralthermal-optical performance analysis procedure, and 2. simultaneous topology optimization of multiple mirror mounts to exploit the interaction of aberrations between different components due to thermal loads and frame movements and minimize a system optical performance while constrained by dynamic stiffness, weight and optical performance measures, both on a component and system level.
The results of the first case study, shown in Fig. 5, indicate there is a certain bandwidth of minimum eigenfrequency constraint values that influence the ability of the optimizer to minimize the surface deformations. Thus, conflicting requirement tolerance values should be thoroughly   Fig. 2. The system, topologies and deformations are not to scale (deformations scaled by a factor 100). The uncoupled optimization attempts to obtain two flat mirror surfaces, whereas the SDO is able to design both mirrors with respect to each other and utilize the curvature of C-M2 to correct for the defocus error introduced by C-M1 investigated, as their limits can highly influence the topological layout and performance.
The results of the second study indicates that the uncoupled optimization aims to design two perfectly flat mirrors, whereas the layout of the second mirror in the coupled optimization is such that its misalignments (mainly curvature) effectively compensate for the misalignments of the first mirror, resulting in a spotsize improvement of over 95%, without reduction of any other performance aspect considered in the optimization.
The activity of the positional error constraint in the optimal solution of the uncoupled optimization (while the system easily remains within the accuracy limits) indicates the system is unnecessarily overconstrained, see Table 2. The SDO approach makes use of the enlarged feasible domain, which is apparent from the mass transfer between the mirrors in the optimal solution and the significant differences in topologies. Thus, the typical design approach unnecessarily overconstrains problems, whereas the SDO approach enlarges the feasible domain and gives the optimizer more freedom to minimize the objective while still satisfying all constraints.
During optimization the SFE constraints are not always satisfied, which means that optical analysis is based on optical properties that do not accurately describe the deformed surface leading to inaccurate results. However, the converged optimal solutions do satisfy the SFE constraints and hence accurately describe the system's optical performance metrics. A more enhanced surface error determination method may result in faster convergence, and additionally a different optimal solution.
One would expect the SDO problem to be more difficult to solve than the uncoupled problem, due to a larger design space and potentially more complicated cost function. However, the results indicate the opposite as the required number of iterations to solve the problem reduces. The iteration history, Fig. 8, shows that the SDO approach is able to drastically decrease the objective in the first number of iterations, where the uncoupled optimization requires more iterations.
The RMS spot size diameter of the uncoupled system reaches the diffraction limit twice, halfway the iteration history, but the optimization continuous because component responses do not satisfy all constraints and termination criteria. On top of that, when the optimizer is able to decrease the objective function of the second mirror, the system spot size increases again. This indicates that the mirrors exactly counteract each others error during optimization, although the optimizer is unaware of this information. Thus, component interaction should be included to obtain a more optimal solution in terms of system performance. Thus, the case study demonstrates that an integrated structural-thermal-optical design optimization procedure taking into account all system components improves the system optical performance compared to individual component optimization, while subjected to the same (or equivalent) design constraints.
The system designed using the SDO method keeps a considerable margin above the competitor and hence the loads may increase considerably before the system's spot size diameter reaches that of the uncoupled variant. Note that the individual components designed using the SDO method are only applicable to this specific system configuration with these specific loads and boundary conditions. In general, the SDO method will result in designs that are more tailored to a specific case and generally are less robust with respect to other loading conditions not considered in the optimization.
This study opens up further research directions, including: -application towards various different load cases to verify and validate the general applicability of the method, -simultaneous optimization of the housing and optical components and simultaneous optimization when their domains are merged into a single mesh, in order to give the optimizer the freedom to relocate the boundary locations, -extension to multiple, and different types of components, e.g lenses, gratings, prisms and initially curved mirrors, -extension to 3D structures and consideration of manufacturability, -including the uncertainty in both thermal and mechanical loading, i.e. robust design, -extension to multi-material topology optimization to achieve higher performances as there are more possibilities to counter effect thermal expansion, create conductive isolation, as well as damp out external vibrations (Van der Kolk et al. 2017), and -enhancing thermal modeling and control by e.g.
considering design-dependent heat loads affected by the misalignments, including radiation influences, and simultaneously optimizing locations and input of active thermal components (heaters/coolers).
It is expected that the more components the system consist of and the stronger their interaction is, the greater the benefit of the SDO method will be. Allowing the optimizer to distribute unavoidable errors over multiple components in the system, instead of letting the designer impose the error budgets on each component, enlarges the feasible domain and the potential for superior system designs.

Conclusions
The key to satisfy next generation optomechanical system requirements is to not distribute error budgets over components a priori, but to consider and optimize the system as a whole. This allows for focus where it matters, without overconstraining the system unnecessarily. A structuralthermal-optical performance analysis is able to expose the performance metrics that matter for optomechanical systems without relying on intermediate derived performance indicators. For a single component system or multicomponent uncoupled optimized system, only minimizing deformations (and nothing else) leads to better optical systems. However, there is additional room for improvement when multicomponent systems are optimized in a coupled fashion as this allows for error compensation between components. Since the feasible design space of the system level optimization completely encapsulates that of the individual component optimization, the globally optimal performance of the coupled system is always better or equal to the uncoupled optimization approach. This is also shown by the results of the numerical example; coupled optimization based on the full structural-thermal-optical performance analysis is able to reduce the spot size of a two-mirror system with 95% compared to uncoupled component optimization to below the system's diffraction limit. The coupled analysis allowed the two mirrors to compensate for each others errors, which is a mechanism that would be otherwise invisible to the optimizer. Despite the fact that real systems are more complex than the simplified example considered in this study, it shows that optomechanical designers should aim for considering and integrating multiple components and physics simultaneously in the design loop, and thus apply the SDO approach, when requirements seem irreconcilable.