## 2 Theoretical Background
In this work, the values of \(x\) and \(y\) are considered relatively small as compared to the radii of the two spherical surfaces of lens, so that the following approximations hold true with sufficient accuracy [15],  
\[\sqrt{1-\frac{x^{2}+y^{2}}{{R_{1}}^{2}}}\approx 1-\frac{x^{2}+y^{2}}{2{R_{1}}^{2}}, \tag{9}\]  
for \(i=1\), \(2\), respectively. The central component part is a flat glass with constant thickness, \(\Delta_{3}(x,\ y)=\Delta_{03}\). Finally, we get the total thickness function given by  
\[\Delta(x,\ y)=\Delta_{0}-\frac{x^{2}+y}{2}^{2}(\frac{1}{R_{1}}-\frac{1}{R_{2}}), \tag{10}\]  
where \(\Delta_{0}=\Delta_{01}+\Delta_{02}+\Delta_{03}\). It is noted that the nominal focal length of optical lens, \(f_{0}\), is defined by  
\[\varepsilon(\omega)=\varepsilon_{0}\Big{[}1+\chi^{(1)}+\frac{3}{4}\chi^{(3)} \big{|}\mathbf{E}(\omega)\big{|}^{2}\Big{]}, \tag{2}\]0  
Then, by combining Equations (4), (10), and (11), we obtain the transfer function of lens  
\[\varepsilon(\omega)=\varepsilon_{0}\Big{[}1+\chi^{(1)}+\frac{3}{4}\chi^{(3)} \big{|}\mathbf{E}(\omega)\big{|}^{2}\Big{]}, \tag{2}\]1  
where the spatial function \(\beta(x,\ y)=\alpha[\Delta_{0}-\frac{x^{2}+y^{2}}{2f_{0}(n_{0}-1)}]\) is proportional to the parameter \(\alpha\). We further assume that the coordinates of the rear nominal focal plane of lens with \(z=f_{0}\) are \((u,\ v)\), and use the Fresnel diffraction integral formula to calculate the focal spots formed by the optical lens. In fact, the field pattern of focal spots in the rear nominal focal plane is given by  
\[\varepsilon(\omega)=\varepsilon_{0}\Big{[}1+\chi^{(1)}+\frac{3}{4}\chi^{(3)} \big{|}\mathbf{E}(\omega)\big{|}^{2}\Big{]}, \tag{2}\]2 The three component parts of an optical lens and their structure parameters. \[\varepsilon(\omega)=\varepsilon_{0}\Big{[}1+\chi^{(1)}+\frac{3}{4}\chi^{(3)} \big{|}\mathbf{E}(\omega)\big{|}^{2}\Big{]}, \tag{2}\]3  
where \(E(x,\;y)\) is the complex light field distribution of an incident laser beam in the front plane of the lens. On substituting Equation (12) into Equation (13), the quadratic phase factor within the integrand are canceled, and one obtains  
\[ \begin{array}{ll}E(u,\;v)&=\frac{\exp[\mathrm{i}\frac{k}{2f_{0}}(u^{2}+v^{2} )]}{\mathrm{j}\lambda f_{0}}\exp[\mathrm{j}kn_{0}\Delta_{0}]\\ &\times\iint E(x,\;y)\exp\left[\mathrm{j}k\beta(x,\;y)\left[E(x,\;y)\right]^{2 }\right]\exp[-\frac{\mathrm{i}2\pi}{\lambda f_{0}}(xu+yv)]\mathrm{d}x\mathrm{d} y.\end{array}  \tag{14}\]  
Based on the light field distribution in the nominal focal plane given by Equation (14), we can further calculate the light field distribution in several transverse planes adjacent to the nominal focal plane according to the angular spectrum algorithm. As a matter of fact, we can define the transfer function of angular spectrum diffraction as  
\[H_{\mathrm{AS}}(f_{x},f_{y})=\exp\left(\mathrm{i}\frac{2\pi\Delta z}{\lambda} \sqrt{1-\lambda^{2}(f_{x}^{2}+f_{y}^{2})}\right), \tag{15}\]  
where \(f_{x}\) and \(f_{y}\) are the spatial frequency in the \(x\) and \(y\) directions, respectively. Then the laser transmission in free space from the transverse plane \(U(x_{0},\;y_{0},\;z_{0})\) to the transverse plane \(U(x,\;y,\;z)\) with interval distance \(\Delta z\) can be calculated by the formula of angular spectrum,  
\[U(x,\;y,\;z)=F^{-1}\Big{\{}F[U(x_{0},\;y_{0},\;z_{0})]H_{\mathrm{AS}}(f_{x},\; f_{y})\Big{\}}, \tag{16}\]  
where \(F\) and \(F^{-1}\) represent the two-dimensional Fourier transform and inverse Fourier transform, respectively [16; 17]. Then the light field distribution of focal spots in the transverse plane of interest with the specified interval distance \(\Delta z\) away from the nominal focal plane can be numerically calculated by applying the 2D-FFT and 2D-IFFT algorithms on Equation (16).