## 2 The Method of Ray Tracing
### Ray Tracing in the Non-Uniform Refractive Index Medium  
In the process of ray tracing, we use Equation (1) to map the discrete density field data from the fluid simulation to a refractive index field. Since the refractive index field is not regular, an exact analytical solution to Equation (1) cannot be obtained. Thus, based on the concept of differentiation, we divide the beam into small segments along the propagation path for the ray tracing process. We assume that within each segment, the refractive index of the medium remains approximately unchanged. The refractive index within each microelement can be approximated as the index at the centre point of that microelement. Thus, the following formula can be obtained:  
\[OPL=\int_{0}^{L}n(s)ds=\sum_{i=1}^{N}OPL_{i}=\sum_{i=1}^{N}l\times n_{i}, \tag{5}\]  
where \(OPL_{i}\), \(l\), \(n(s)\), and \(n_{i}\) denote the optical path length, geometric length, refractive index profile, and mean index of refraction of the \(i\)-th microelement respectively, and \(N\) is the total number of microelements. According to the differential hypothesis, a refractive index gradient exists at both the incident and exit surfaces of each microelement. Therefore, according to [PERSON]'s law:  
\[n_{1}\sin\theta_{1}=n_{1}\sin\theta_{2}, \tag{6}\]  
When a beam passes through each interface of the element body, its propagation direction will be deflected. Since the isosurface of the refractive index field in the beam expanding system is not a plane but an irregular surface and the incident plane of a microelement cannot be either coincident with or parallel to a constant refractive index surface, we cannot simply take the surface normal of each element body as the normal to the incident plane for a ray. Therefore, we use the Barron gradient operator to solve for the gradient of the refractive index [19; 20]:  
\[\begin{array}{rl}\frac{\partial n(i,j,k)}{\partial x}=\frac{1}{12}n(i-2,j,k) -\frac{8}{12}(i-1,j,k)+\frac{8}{12}n(i+1,j,k)-\frac{1}{12}n(i+2,j,k)\\ \frac{\partial n(i,j,k)}{\partial y}=\frac{1}{12}n(i,j-2,k)-\frac{8}{12}(i,j-1, k)+\frac{8}{12}n(i,j+1,k)-\frac{1}{12}n(i,j+2,k)\\ \frac{\partial n(i,j,k)}{\partial z}=\frac{1}{12}n(i,j,k-2)-\frac{8}{12}(i,j, k-1)+\frac{8}{12}n(i,j,k+1)-\frac{1}{12}n(i,j,k+2)\end{array} \tag{7}\]  
where \(\left(\frac{\partial n(i,j,k)}{\partial x},\frac{\partial n(i,j,k)}{\partial y },\frac{\partial n(i,j,k)}{\partial z}\right)\) is the direction of the refractive index gradient at the point \((i,j,k)\). In this way, we can calculate the phase change and deflection of the light in the beam expanding system using the simultaneous Equations (5), (6) and (7).  
The beam quality of a beam expanding system can be evaluated by calculating the wavefront aberration of a beam passing through the system. Moreover, Zernike annular polynomials can be used to accurately fit the annular wavefront aberration to effectively extract the components of the wavefront aberration, which can serve as an important basis for wavefront aberration correction. Therefore, in this paper, the wavefront aberration obtained through simulation is fitted with Zernike annular polynomials, and the compositions of the OPD profiles under different inlet conditions are compared to determine the optimal intake layout [21; 22].