## 2 Theory and Method  
The wave equation is the fundamental equation, which describes electromagnetic waves. In cylindrical coordinates, the solution to the wave equation is  
\[u(r,\varnothing,z)=\frac{ke^{\mathrm{i}\mathrm{k}z}}{i2\pi z}\int_{0}^{r}\int_{ 0}^{2\pi}u_{0}(\rho,\varnothing,\ 0)e^{\mathrm{i}\hbar\frac{r^{2}+p^{2}-2\mathrm{ terms}(\varnothing-\theta)}{2z}}\rho d\rho d\theta. \tag{1}\]  
The wavenumber is defined as \(k=2\pi/\lambda\), where \(\lambda\) represents the wavelength. \(\rho\) represents the radial distance in the polar coordinate system defined on the plane perpendicular to the optical axis. The \(i\) represents the imaginary unit. And the radial \(\mathbf{r}\) and azimuthal \(\varnothing\) coordinates, along with the longitudinal distance \(z\), form the spatial coordinate system. The initial wavefront at \(z=0\) is expressed as \(u_{0}(\rho,\varnothing,\ 0)=A(\rho)exp\left[ikC(\frac{\rho}{\omega_{\rho}})^{ \gamma}\right]\), where \(A(\rho)\) is amplitude function, \(\rho\) is the polar diameter of the polar coordinate system perpendicular to the optical axis direction, \(\omega_{\rho}\) is the phase normalization coefficient, \(C\) (dimensionless) is modulation depth, and \(\gamma\) is shape factor. Subsequent studies show that the modulation depth \(C\) determines the magnitude of the wavefront phase perturbations, while the shape factor \(\gamma\) controls the rate at which the beam forms. Together, these parameters influence the propagation characteristics of the optical beam through turbulent media. Using the Gaussian amplitude \(e^{-\frac{\rho^{2}}{\epsilon_{0}}}\) as the initial \(A(\rho)\), substituting the polar coordinate amplitude function \(A(\rho)\) into \(\upsilon_{0}\) and evaluating the Fresnel integral (1) along the optical axis (\(r=0\)) gives the on-axis light field distribution in Equation (2):  
\[u(0,z)=\frac{k}{iz}e^{ikz}\int_{0}^{\infty}e^{-\frac{\rho^{2}}{\epsilon_{0}}}e^{ ik(|C\frac{\rho}{\epsilon_{0}}\gamma+\frac{\rho^{2}}{\epsilon_{0}}|)}pd\rho. \tag{2}\]  
Due to the axial symmetry of the input optical field, the integral is evaluated in one-dimensional polar coordinates (\(\omega_{\rho}=1\)) using the phase stabilization method. The phase term of the integrand is  
\[\phi(z)=C(\rho)^{\gamma}+\frac{\rho^{2}}{2z}. \tag{3}\]  
An equilibrium point satisfies  
\[\frac{d\phi}{d\rho_{0}}=0. \tag{4}\]  
This yields the equilibrium point:  
\[\rho_{0}=(\,|\,C\,|\,\gamma z)^{1/(2-\gamma)}. \tag{5}\]  
The second derivative at \(\rho_{0}\) is  
\[\phi^{\prime\prime}(\rho_{0})=\frac{2-\gamma}{z}. \tag{6}\]  
Since \(\phi^{\prime\prime}(\rho_{0})>0\), we introduce the phase factor \(e^{i\frac{\pi}{4}}\). Applying the phase stabilization formula gives  
\[\int_{0}^{\infty}e^{-\frac{\rho^{2}}{\epsilon_{0}}}e^{ik\phi(\rho)}\rho d\rho \approx e^{-\frac{\rho^{2}}{\epsilon_{0}}}\rho_{0}\sqrt{\frac{2\pi z}{k(2- \gamma)}}e^{ik\phi(\rho_{0})+\frac{i\pi}{4}}. \tag{7}\]  
Here, \(\phi(\rho_{0})\) is defined as  
\[\phi(\rho_{0})=\frac{|C|\gamma}{2-\gamma}(|C|\gamma z)^{\frac{\gamma}{2- \gamma}}+\frac{\rho_{0}^{2}}{2z}. \tag{8}\]  
It should be noted that this formula serves as a simplified representation to facilitate focal position calculations and does not accurately reflect the actual optical field distribution. The complex amplitude distribution on the optical axis is  
\[u(0,z)=\frac{k}{iz}e^{ikz}e^{-\frac{\rho^{2}}{\omega_{0}}}\rho_{0}\sqrt{\frac {2\pi z}{k(2-\gamma)}}e^{ik\phi(\rho_{0})+\frac{i\pi}{4}}. \tag{9}\]  
This holds for \(0<\gamma<2,C>0,\ z>0\).