Ray tracing in thermally loaded optical systems

High‐power laser beams are usually imaged with the help of suitable optics. The propagation of such beams, through optical systems, often result in a focus shift. This is mainly due to the deposited power in the beam acting as a heat source in the optical elements. Finally, the resulting change in temperature leads to both, a deformation of the lens surface, and a change of the refractive index in the lens. Ray tracing and the finite element method are needed in order to analyze the above effects. The accuracy and the computational time of focus shift calculations, in both simulation techniques, are strongly influenced by the interpolation of the data. This paper will present, a newly proposed interpolation method, allowing a fast and accurate interpolation of data using an unstructured finite element grid. This method will specifically focus on the accurate calculation of laser beam thermal load inside an optical lens.


Introduction
High power laser systems are needed in multiple applications in the fields of material processing, medical treatment, or fundamental research. One of the most crucial points, in engineering such lasers systems, is stability in a high power regime. The temperature of a lens usually changes due to the induced thermal load, deposited by the beam, and leads to a deformation of the lens and a change in its refractive index. This causes the thermal lens effect, which induces a focal shift, as shown in Fig. 1. It is therefore of high interest, to simulate optical systems under the thermal influence, in order to obtain a thermally stable operating state and investigate compensation methods.
A new method of interpolating continuous data on the unstructured grid and efficiently modeling the heat source is presented in two parts, using the proposed framework. Continuous data interpolation and mapping of the heat source are the main focus of the computational effort in the presented simulation.
Interpolation is usually applied with help of an additional regular grid, on which one can define e.g. splines or polynomials, for high accuracy, or use trilinear interpolation on regular grids, for fast results. We propose a new interpolation method, which uses a smart trilinear interpolation directly on unstructured grids, allowing both fast and accurate simulation results. Additionally, a new method for mapping the heat source is presented, by only evaluating the surface intersection points of propagated rays. This contrasts with regular heat source mapping methods, as a significantly lower number of rays can be used, thereby speeding up simulation time.
The combined use of both proposed methods enable extensive parameter studies allowing for significantly more new design opportunities and further compensation methods.

Heat source modeling
An optical beam can be modeled with help of rays where each ray is associated with a certain power. All rays can then be propagated independently, starting at the vertices of the initial intensity distribution (see Fig. 2). It is first necessary to obtain the intensity in normal direction I n to the element surface in order to represent the optical beam as a heat source. Also, the power flow inside the element, in propagation direction s, is required.
In the geometrical optics regime, the intensity normal to a surface can be accurately approximated with rays. The intensity variation in a ray tube, formed by closely spaced rays, can then be written in the following manner, as in [1] I n (ξ, η, σ f ront ) = I n (ξ, η, σ init ) dA init dA f ront , where dA is the ray tube cross section at the corresponding planes, σ the physical arc length of the ray, n is the surface normal, and s is the ray direction in the optical element. 2k 0 n ′′ = α abs is the local absorption coefficient, attenuating the intensity,  where k 0 is the wave number and n ′′ is the imaginary part of the refractive index. c abs is the surface absorption coefficient. The relation of intensity and cross section is also known as the intensity law of geometrical optics [2]. Therefore, Eq. (1) can be used to approximate the intensity along optical beams, only by calculating the corresponding area increment dA. The area is defined orthogonally to the surface normal n on an optical element or a stop plane along a ray, and each ray has its own area increment. The area increment dA can be calculated by summing one-third of the adjacent triangle areas, 1 3 A T (marked in red), see Fig. 2. This respectively allows calculation of intensity on arbitrary continuous surfaces, such as lenses or mirrors.
The total heat source can now be split into two main components: the surface absorption ∂Q and the bulk absorption Q • . From the intensity on the surface I n , the absorbed power in the medium Q • , as well as the absorbed intensity on the surface ∂Q, can be calculated. On the surface of the optical elements, the absorption ∂Q can directly be obtained by multiplying the intensity on the surface I f ront where n is the surface normal on Ω and s the direction of the energy flow. σ f ront and σ back are the points of the ray-surface intersection. Assuming that the ray can be approximated by a straight line, σ is the corresponding trajectory in between. The ray trajectory is not a straight line in general, however, it is very close to one, such that the resulting error is negligible. In comparison to the regular ray tracing method, where each ray is associated with a certain power, this method has the benefit of significantly reducing the number of rays necessary to obtain an accurate heat source. Additionally, this method utilizes power density rather than traditional ray tracing methods, which provide power per ray. Calculation of heat source requires power density rather than pure power.

Heat equation
The weak formulation of time dynamic heat equation is given by where k 0 is the thermal conductivity, H the heat transfer coefficient, and c p ρ the volumetric heat capacity. The heat source terms ∂Q and Q • are modeled by Eq. (2). The boundary ∂Ω can be split in actively cooled surfaces ∂Ω a and passively cooled surfaces ∂Ω p . Actively cooled surfaces have contact with metal, which efficiently transfers heat, where passively cooled surfaces transfer heat only via air or thermal radiation. In case of strong cooling, −κ 0 ∂T ∂n ≪ H (T − T c ) at the boundary Ω a , the Dirichlet boundary condition can be used, such that T = T c on ∂Ω a .

Ray tracing through in homogeneous media
Both, change in temperature and deformation, influence the propagation of rays through media. Temperature change induces a refractive index gradient, which deflects rays inside the medium, whereas deformation mainly alters refraction direction on the surface of optical elements. This article will not discuss deformation and will instead focus on refractive index.
The refractive index depends in good approximation linear on the temperature, such that As the temperature T is not constant, this effect induces refractive index gradients. In order to propagate rays in gradient media, it is necessary to solve the ray equation in the following manner where r is the position of a ray, and p its directional vector. This equation has to be discretized properly, e.g. by using Runge-Kutta [3] or Euler discretization. For all discretization schemes, it is crucial to obtain fast and accurate continuous data on the domain Ω, in order to solve Eq. (5). Additional to the refractive index n, also the gradient ∇n is necessary.
In order to solve the ray equation Eq. (5), and additionally consider deformation, it is necessary to obtain continuous data on discretized domain Ω h . Partial derivatives ∂ ∂x , ∂ ∂y and ∂ ∂z are also needed accordingly. Simulation accuracy thus not only depends on the finite element simulation but also on data interpolation.
We additionally introduce a new normal projection (NP) interpolation method in order to allow fast and accurate approximation of the trilinear Eq. (6). Using the data locality of ray tracing, this method is almost as fast as using the O(1) trilinear interpolation on regular grids, while at the same time increasing accuracy.
For interpolating a given point x Int in a hexahedron, the trilinear Eq. (6) has to be solved. Using the NP interpolation method, this can be approximated by calculating the shortest distance of the point x Int to each of the six surfaces of the cell, such that where d i is the distance to the surface of the cell, c i is the central point of face i, and n i the corresponding normal vector. As an example, on the top side of Fig. 3, c i and n i are given by www.gamm-proceedings.com In general, the hexahedron H(r i ), in which a point r i is contained, is unknown and has to be found first. This can be very expensive, especially for randomly distributed points. In the case of ray tracing, however, which accesses data at subsequent nearby points, the hexahedron finding can speed up significantly. Consider the points r n−1...n+2 of Fig. 3. If the step size ∆σ is small enough, H(r n+2 ) is always directly adjacent to H(r n+1 ). The new hexahedron can be found by evaluating Eq. (7) with the new point in the previous hexahedron. One or more of the coordinates α, β, γ will not be in [0, 1], which can be used to directly find the new hexahedron. For example, if Eq. (7) gives α = 0.2, β = 0.3, γ = 1.1, then the correct hexahedron will be found by incrementing the hexahedron index, corresponding to the direction of γ, by 1. Similarly, the same method can be applied to all rays, as they are propagated in close vicinity to each other. As a result, an expensive time consuming full hexahedron search is only once required for each simulation.

Simulation results
A comparison of errors, for different interpolation methods, are shown in Fig. 4. Different functions, similar to Lee Et.al. [4], are used. These figures show the interpolation error and time for evaluating a single point in a discretized lens domain Ω h , shown in Fig. 3, for different interpolation algorithms. Trilinear interpolation on an auxiliary regular grid, direct interpolation on the unstructured grid with both, NP and Newton's iteration, as well as the Multilevel B-Splines (MBA) interpolation [4] are compared accordingly. The auxiliary regular grid is selected, as it has a finer mesh size than an unstructured grid, and results in approximately 20× more grid points. Fig. 4 shows the major benefits of using the NP interpolation method. NP is almost as fast as using interpolation on a regular, trilinear grid, with much greater accuracy. This is especially true for the boundary of the domain Ω h , which is why the maximal error difference is even more significant. Also, the NP method is much faster than solving Eq. (6) with Newton's method, while giving no significant difference in interpolation error. Finally, the direct interpolation benefits from significant memory savings, as the additional regular auxiliary grid does not have to be initialized and stored.
For comparison, the MBA interpolation method, described in [4], is also shown. It has an even lower average interpolation error, for the price of being ∼ 10 3 × times slower than the auxiliary grid interpolation. Interestingly, the MBA has a better average interpolation error, however lacks in the maximums norm. Thus, the newly developed NP interpolation method combines both fast interpolation results and good accuracy.

Conclusion
A new interpolation and heat source mapping method was developed for simulation of thermal effects in arbitrary optical systems. This new method allows for fast, accurate, and data saving interpolation while significantly reducing the number of rays required for heat source mapping.