Transient dynamical-thermal-optical system modeling and simulation

In modern high resolution optical systems like astronomical telescopes or lithographic objectives, performance degradations can be caused by various disturbances. Holistic optical system simulation is required to predict the performance or the high precision systems. In this paper a method for transient dynamical-thermal-optical system modeling and simulation is introduced. Thereby, elastic deformation, rigid body motion, and mechanical stresses due to dynamical excitation are calculated using elastic multibody system simulation and temperature changes are determined using thermal finite element analysis. The deformation, the motion, and the mechanically and thermally induced stress index changes are then considered in a gradient-index ray tracing. Finally, the presented method is applied in a dynamical-thermal single lens system.


Introduction
Nowadays, the developed high precision optical instruments are very sensitive to external influences. Holistic system simulation is used to predict the system performance accurately and to prove that the design meets the demanded requirements. In lithographic objectives even small mechanical disturbances like deformations or rigid body motions of the mirrors or lenses resulting from external excitation, e.g., from the waver motion system lead to performance degradations. Besides that, temperature changes in the optical elements due to the heating from the laser beams influence the imaging performance. In laser-based additive manufacturing the laser beam has alternating on and off times and the laser beam moves during the exposure time. Thus, transient simulations are necessary to simulate the optical system behavior. All together, a holistic optical system simulation must be able to consider transient mechanical and thermal disturbances. In this paper a method for transient dynamical-thermal-optical system modeling and simulation is introduced and an exemplary system simulation is presented for a single lens system.
The dynamical-thermal behavior of optical systems can be calculated by a combination of rigid body motion, small deformations, the related mechanical stresses, and the temperature distribution. For the mechanical analysis various literature exists, where structural finite element (FE) analysis results are used to describe the deformation of the optical elements [1][2][3]. In order to consider mechanical stress in the optical simulation three dimensional shape-function interpolation is used to map the changes of the refraction index with the mechanical stress resulting from the FE analysis [4]. In contrast to this, the introduced method of this paper uses more efficient elastic multibody system (EMBS) simulations to calculate the transient elastic deformations and the rigid body motions of the system. Due to a model order reduction and the possibility to combine rigid bodies with flexible bodies, the EMBS simulation is usually faster than the FE analysis and still sufficiently accurate. Instead of the stresses from FE analysis we get the stresses from the EMBS simulation and transfer them to cylindrical polynomials which can be evaluated efficiently during the optical simulation. The thermal effects on the optical system are taken into account in the presented method by considering the changes of the refraction index due to temperature changes. A thermal FE analysis is executed in various papers and the refraction index changes are transferred to the optical system using three dimensional shape function interpolation [4,5] and optical path difference (OPD) maps [6,7]. In the presented method cylindrical polynomials are used to describe the three dimensional refraction index changes. With this method, arbitrary three dimensional refraction index fields can be considered in the optical simulation.
Additionally to the dynamical and the thermal analysis, and the data transformation from the mechanical and the thermal results to the optical system the presented method includes a gradient index (GRIN) ray tracing algorithm which considers the dynamical and thermal disturbances.
There are several existing tools for mechanical-thermaloptical system simulation. An overview is given in [8]. The novel contribution of the presented method is the usage of transient EMBS simulations instead of structural FE analysis and the usage of cylindrical polynomials combined with GRIN ray tracing instead of OPD maps or shape-function interpolation for the stress and temperature dependent changes of the refraction index. A workflow for the presented method is given in Fig. 1.
In the following, the mechanical, thermal, and optical methods are introduced briefly. Additionally, a numerical example of a holistic simulation of a dynamical-thermaloptical single lens system is given.

Methods and numerical example
The holistic simulation of transient dynamical-thermaloptical systems combines mechanical simulations, thermal system analysis and optical simulations. The method and an exemplary numerical example are described in more detail in the following.

Mechanical modeling of dynamical-optical systems
In the mechanical part of the holistic simulation the optical system, consisting of optical elements and their holdings, is modeled as EMBS [9]. In a first step, each body is modeled separately. Optical sensitive components such as optical lenses are considered to be elastic bodies while the remaining less important elements such as the holdings are considered as being rigid. For each elastic body, a finite element model with N degrees of freedom (DOF) is built. The deformations of each modeled body are described by its shape functions and its nodal displacements vector q(t) ∈ R N which varies dependent on the time t. The equations of motion for each elastic body result from the linear FE model and read ( 1 ) Thereby, the system matrices are the mass matrix M ss , the damping matrix D ss , and the stiffness matrix K ss . The system input is the vector of the exciting forces f . The normally large number of DOF of FE models leads to a high computational effort for the analysis of their dynamical behavior. By means of a model order reduction (MOR) method, the number of DOF can be significantly reduced after the determination of the system matrices and the force vector. In the MOR scheme, a projection matrix V ∈ R N×n with n N is calculated. This matrix enables the approximation of the nodal displacements by q ≈ Vq (2) with the reduced coordinatesq(t) ∈ R n . The projection matrix is applied to the system matrices and the system input vector which leads to a reduction of the whole equation system size of Eq. (1). The reduction allows for an efficient analysis of the dynamical behavior of each modeled body. The detailed procedure description for a MOR and different methods for the calculation of V are given in [10,11]. The application of MOR methods on optical systems is discussed in [12,13]. Subsequently to the MOR, the rigid and elastic bodies are assembled to the EMBS. Then, the dynamic simulation of the assembled system, i.e. a time integration, yields the time-dependent kinematics of the complete system and the deformations of the elastic bodies. Apart from the usage in the MOR, the projection matrix is also provided to the FE model of each elastic body if stresses need to be considered. There, it is used for obtaining the deformation state of the elastic body. By the use of this deformation state and material laws the matrix of the reduced stress functions¯ σ ,P ∈ R 6×n at a discrete point P within the elastic body is calculated. The matrix¯ σ ,P and the reduced nodal displacements vectorq allow for the calculation of the vector of the time-dependent stresseŝ at the point P. The normal stresses σ x , σ y , and σ z in the direction of the reference system in P and the shear stresses τ xy , τ yz , and τ xz in the associated planes can be rearranged to the stress tensor P . By the use of an eigen analysis of P the principal stresses σ I,P , σ II,P , and σ III,P are computed. A more detailed explanation on the calculation of the stresses can be found in [14,15]. In summary, the EMBS dynamical simulation allows for the computationally efficient calculation of the timedependent deformations and principal stresses in the elastic bodies and the time-dependent rigid body motion of all bodies.

Thermal finite element analysis
In the thermal analysis part of the holistic simulation of the dynamical-thermal-optical system, an FE model is built for all optical lenses of the system. Normally, the FE mesh for the thermal analysis is coarser than the mechanical mesh. This saves computing time and nevertheless leads to sufficiently precise results. The calculation of the time-dependent nodal temperatures ϑ(t) in each optical lens is performed using the dynamical equation for the heat conduction The matrix D tt for the specific heat, the matrix K tt for the thermal conductivity, and the thermal input vector φ which represents the thermal heat flux at the nodes, are computed with FE methods. While thermally induced changes of the refraction index are considered in the presented method, thermally induced deformations are not considered here.

Optical modeling of dynamical-thermal-optical systems
The dynamical and thermal effects in optical systems can be analyzed using ray tracing. Thereby, the propagation characteristics of light rays in different media and especially media with an inhomogeneous refraction index must be taken into account. The used ray tracing approach is described in the following. Ray tracing is a part of the geometrical optics where light rays determine the imaging performance of an optical system and wave-optical effects are neglected. In homogeneous media the rays propagate along straight lines and change their direction of travel only at the surfaces. In sequential ray tracing through homogeneous media, the ray positions r and the ray directions d are calculated in a relative surface formulation from surface to surface, from the object plane to the image plane. Thereby, the calculated ray position r i at a surface S i is the intersection point of a ray and the surface. The corresponding ray direction d i results from the previous ray direction d i−1 , the intersection point r i , the surface normal n i , the refraction index of the previous media n i−1 , and the current one n i [16]. Light propagates as straight line through homogeneous media. An inhomogeneous refraction index within a medium influences the path of a light ray and makes it curved. The resulting bending of the ray depends on the gradient of the refraction index. In order to trace a ray through a GRIN medium, the partial differential ray equation must be solved numerically. A derivation of this equation can be found in [17]. The path parameter s can be used to approximate the differential geometrical path length ds by ds = dsn. In combination with (5) one gets which can be solved using a fourth order Runge-Kutta method [15,18]. With this method a ray can be traced efficiently through a medium with an arbitrary inhomogeneous refraction index in small spatial steps till the back surface of the inhomogeneous medium is reached. In every spatial step the changes in the refraction index are calculated and considered in the calculation of a new ray position r and a new ray direction d. The only requirement is that the refraction index distribution within the inhomogeneous medium is known during the ray tracing.

Dynamical-thermal-optical system modeling and simulation
The analysis of transient dynamical-thermal-optical systems requires the combination of the mechanical, the thermal, and the optical system analysis. Therefore, the results of the dynamical EMBS simulation and the results of the thermal FE analysis must be transformed in such a way that the results can be considered in the ray tracing.
From the EMBS simulation we get the transient kinematic behavior of the system. Thereby, the motion is described relative to an inertial system frame. For the sequential ray tracing, all rigid body motions and elastic deformations must be transferred in a relative surface description for each time step and each surface. The time-dependent deformations of the surfaces must additionally be transformed in a way that the surfaces are described smooth enough to avoid discretization effects during the calculation of intersection points of the rays and the surfaces. Hence, the deformations must be approximated using polynomial functions and in this case we use Zernike polynomials. These polynomials can be evaluated very efficiently during the ray tracing and they can reproduce random deformation states of circular lenses. Deformations of free-form surfaces can be approximated smoothly using B-splines. The third result from the mechanical simulation, are the time-dependent principal stresses which lead to a spatial variation of the refraction index within deformed optical elements. The change of the refraction index n σ ,j due to stress at a position r j in the medium or on the surface of an optical medium result from the stress optic law n σ ,j = C σ I,j + σ II,j + σ III,j with the material dependent stress-optical coefficient C. Thereby, polarization effects are neglected. The stress dependent refraction index changes are approximated using cylindrical polynomials, to avoid discretization effects. Therby The height and the radius are chosen so that all nodes of the deformed lens are inside the cylinder. The refraction index changes at a random point p = (x p , y p , z p ) within the lens can be approximated using with the cylindrical polynomial coefficients c. These coefficients can be calculated by the evaluation of Eq. (8) with the known refraction index changes at all nodes of the thermal FE mesh. The Eq. (8) is then evaluated during the numerical solution of Eq. (6) and thus, the changes of the refraction index due to stress can be considered efficiently in the ray tracing. A more detailed description of the application of cylinder polynomials in the GRIN ray tracing is given in [19]. The result of the thermal FE analysis of the optical elements are the time-dependent nodal temperatures. The temperature changes ϑ lead to a spatial variation of the refraction index. The changes of the refraction index in a position r j due to temperature changes can be described by Thereby, D 0 , D 1 , D 2 , E 0 , E 1 , and λ tk are material dependent constants and n 0 is the refraction index at the reference temperature of 20 • C. The spatially disturbed refraction index changes are then approximated with cylinder polynomials and evaluated during the ray tracing, such as previously described for the mechanical stresses. The refraction index changes due to stress and due to temperature changes are added during the ray tracing. Thus, mechanical stress effects and thermal effects are considerable in the ray tracing.
With the transformed results from the dynamical and the thermal system analysis the ray tracing can be performed. Thus, a holistic simulation of optical systems is possible with this method. Thereby, transient rigid body motions, dynamical deformations, refraction index changes due to mechanical stresses and thermally induced refraction index changes are considered computational efficiently during the ray tracing.

Implementation
To investigate rigid body motion, deformation, mechanical stress and thermal effects in a dynamical-thermaloptical system, the introduced simulation procedure is implemented. The holistic simulation is performed using ANSYS for the building of the FE models, MatMorembs 1 for the MOR, Neweul-M 22 for the EMBS simulation, TheelaFEA for the thermal analysis, and OM-Sim 3 for the data transformation and the optical simulation. The latter four are MATLAB based software tools which are developed at the Institute of Engineering and Computational Mechanics at the University of Stuttgart. A holistic workflow diagram is shown in Fig. 1.

Numerical example
The presented holistic dynamical-thermal-optical system simulation method is applied to a numerical example in the following. The used system and its excitation are shown in Fig. 2. It consists of a single lens with a diameter of 100 mm and a homogeneous refraction index of n 0 = 1.6 as long as no mechanical stresses and no temperature changes occur. The lens is fixed by a holding, which has six support beams evenly spaced around the circumference. The lens and the holding beams are modeled as elastic bodies, while the holding itself is considered to be rigid. As shown in the figure, the system is excitated mechanically with a predefined movement over the time in all spatial directions translational a t,x , a t,y , a t,z and rotational a r,x , a r,y , a r,z . Therefore, a time-dependent random signal based on a Gaussian distribution with a maximum translation of 1 × 10 −3 mm, a maximum rotation of 1 × 10 −2 rad, and a maximum frequency of 100 Hz is applied on the center of gravity of the optical system. A time-integration of the EMBS system leads to the deformations and rigid body motions of the lens and to the corresponding mechanical stresses. The thermal excitation are seven heat fluxes applied on the lens as shown in Fig. 2. They represent the heat transfer from seven light rays which come from infinity, propagate through the optical lens, and focus on the optical axis on the image plane if no disturbances exist. A time-integration for the thermal FE model leads to the temperatures over the time. For the analysis of the dynamical-thermal-optical system behavior, the results of the ray tracing algorithm are investigated next for the previously mentioned seven rays.

Fig. 2 Single lens optical system with thermal (red) and mechanical (black) excitation
In order to investigate the dynamical-thermal-optical system for a large number of rays, another thermal boundary condition is applied. A heat flux is applied on all nodes of the thermal FE lens model which are within a cylinder with a radius of 20 mm, with the center at x = y = 11 mm and with the height over the whole lens thickness. This applied thermal boundary condition is completely unrelated to the ray paths. It is assumed to visualize thermal influence on images only. The dynamical excitation is the same as described before. In this case the light propagates from one point uniformly disturbed through the lens and forms a grid on the image plane. So, a geometrical image simulation is performed. The unexcited system and its image are shown in Fig. 3.

Results and discussion
The results of the ray tracing for the dynamical-thermaloptical system with the first thermal boundary condition are shown in Fig. 4. Thereby the thermal effects, the surface deformations and the mechanical stresses are scaled so that the effects get visible. The figures show the ray tracing results for different exciations all at the same time instant. In each subfigure the blue dashed lines show the ray paths for the ideal, undisturbed system. There, the rays focus on the optical axis on the image plane.
In Fig. 4a only the changes of the refraction index due to the heating of the lens are considered in the ray tracing. Deformations or stresses are not taken into account. The temperature distribution within the lens is the same for each ray. Therefore, each ray is subject to the same changes in the refraction index which results in rotational symmetric changes of the ray paths around the optical axis. As it can be seen, the thermal influence leads to a shorter focal length.
If only the surface deformations and the rigid body motion of the lens due to the dynamical excitation are considered in the ray tracing, the ray paths shown in Fig. 4b result. No refraction index changes are taken into account. The refraction angle and the intersection position of each ray with the lens surface differs due to the non rotational symmetric deformation. This leads to the visible deviation of the ray positions on the image plane from the ideal focal point.
In Fig. 4c the surface deformations and rigid body motions due to the dynamical excitation as well as the resulting stresses within the lens are considered during the ray tracing. All three disturbances are not rotational symmetric around the optical axis. Thus, not only the changes in the surface intersections and refraction angles but also different refraction index changes occur for each ray. Each considered effect leads to different ray path changes for each ray. The additional considered refraction index changes, compared to the ray tracing where only the deformation is considered, intensify the individual path If the deformation due to the dynamical excitation, the resulting mechanical stresses, and the thermally induced refraction index changes are considered during the ray tracing, the presented holistic dynamical-thermal-optical system simulation is executed. The effects shown in Fig. 4a and c summarize and the ray paths shown in Fig. 4d result. The reduction of the focal length due to the thermally induced refraction index changes as well as the deviation of the ray positions on the image plane from the ideal focus due to the dynamical excitation are present and clearly visible in the figure. Thus, in the holistic simulation of optical systems the presented tool can consider changes in the refraction index due to temperature changes and due to mechanical stresses and it is able to consider rigid body motion and elastic deformations due to a dynamical excitation.
In order to visualize thermal effects within a lens on an image, the cylindrical thermal excitation described in the previous section is applied on the single lens. The resulting temperature distribution in the lens after three seconds is shown in Fig. 5 on the left. In the same figure, results from the geometrical image simulation are shown for three different time steps. The center of the heated cylinder almost coincides with the center of the quarter top right of the image and is marked in the image for t = 0s with the red dot. The changes of the refraction index due to temperature changes for each ray are visualized in the image plots by the color. The higher the changes in the refraction index, the more orange are the The temperature distribution on the lens after three seconds (left) and the images for t = 0s, t = 2s and t = 3s. The higher the changes in the refraction index due to temperature is, the more orange are the ray positions on the image plane. The red dot shows the center of the cylinder within which the heat is applied on the lens shown ray intersection points with the image plane. As it can be seen in the three images of Fig. 5, the heating of the lens leads to a distortion of the image. The grid expands in the heated area due to the thermally induced refraction index changes. In the parts of the grid where no or nearly no temperature changes occur, the grid does not change.
In order to visualize dynamical and thermal effects within a lens on an image, the dynamical exciation and the cylindrical thermal excitation described in the previous section are applied on the single lens. The resulting geometrical images at t = 3s are shown in Fig. 6. Thereby, the geometrical images are shown for different combinations of the disturbances.
In Fig. 6a the geometrical images of the optical system with only dynamical excitation are shown. The geometrical image shown in light grey result from the Fig. 6 Geometrical image after three seconds excitation time with different excitations undisturbed system. If only the rigid body motions and deformations are considered in the ray tracing, the image shown in orange results. As it can be seen, the grid is irregulary aberrated due to the irregular disturbances. If the mechanical stresses and the resulting changes of the refraction index n σ are considered additionally, the blue image results. Due to the irregular stress effects additionally to the irregular deformation and motion disturbances the aberrations in the blue image are larger than in the orange image.
In Fig. 6b the undisturbed image in light grey and the orange image resulting from only dynamical excitation are shown again. The pink image results if the disturbances due to lens motions and deformations and additionally the thermal refraction index changes n ϑ due to the heating of the lens are considered in the ray tracing. Therein, it becomes visible that the distortions due to the lens deformations (orange) and the stretching of the grid due to the heating (Fig. 5) add up in the pink image.
If the dynamical excitation, the resulting mechanical stresses with the associated refraction index changes n σ , and the cylindrical thermal excitation with the resulting thermal refraction index changes n ϑ are considered we get the image shown in Fig. 6c in violet. The geometrical image is now strongly affected. The stretching of the grid at the heated areas is still visible but due to the rigid body motion, the deformation, and the mechanical stress effects it is supperimposed by other disturbances.
Altogether, Fig. 6 shows that the thermal distortions and the shifting of the ray intersections with the image plane due to the dynamical excitation superpose and visibly influence the imaging performance of the system. Thus, it is important to consider time dependent thermally and mechanically induced refraction index changes, rigid body motions, and elastic surface deformations in the geometrical image simulation of optical systems if high precision images are required.
In order to evaluate the efficiency of the implemented method, the calculation time of the geometrical image simulation was tracked. The thermal and structural FE model building, the system matrices extraction, the model order reduction, the dynamical time simulation with the rigid body motion, deformation, and stress calculation, the thermal simulation, and the data transformation to the optical system are calculated in less than five minutes for 30 calculated time steps in a time period of three seconds. The image simulation for one time step, so the GRIN ray tracing of 2356 rays, takes less than 30 seconds. Thereby, for each ray about 230 positions and new ray directions within the lens were calculated. At each of these positions the refraction index changes due to temperature changes and mechanical stresses are considered. At the front and the back surface of the lens, the surface deformations, the rigid body motion and the refraction index changes are considered. The given calculation times show, that the holistic simulation is efficiently applicable for the simulation of dynamical-thermal-optical systems.
The described method enables a qualitative characterization of the influence of dynamical and thermal distortions on optical systems. A quantitative characterization of the aberrations would be possible if the optical path difference and the line of sight error are calculated additionally. For materials with a homogeneous refraction index this is possible in the software OM-Sim as shown in [16,20]. For materials with a gradient in the refraction index the implementation is not yet realized.

Conclusion
In order to efficiently carry out holistic time simulations of dynamical-thermal-optical systems, a numerical ray tracing procedure has been presented and applied in this paper. In the system simulation, finite element thermal analysis, model order reduction methods, elastic multibody system simulation, and optical ray tracing with gradients in the refraction index are combined. The presented procedure is able to consider time-dependent dynamical disturbances, such as rigid body motions and elastic deformations of optical elements and their holdings. Moreover, it can calculate the time-dependent stresses resulting from the deformations and transfer this information into refraction index changes which are considered during the ray tracing. Besides the dynamical disturbance consideration, thermal FE analyses are implemented and the resulting time-dependent temperatures are used to consider the thermally induced refraction index changes in the ray tracing. In a numerical single-lens example the presented method was applied. The example shows the importance of the consideration of the dynamical and thermal effects in the simulation of high precision dynamical-thermal-optical systems. Besides that, it shows that the presented method enables the numerically efficient calculation of time-dependent rigid body motions, elastic deformations, mechanical stresses, and temperature changes of optical systems and it allows for the consideration of the resulting disturbances in the gradient-index ray tracing.