## 2 Simulation Methods and Evaluation Parameters
### Simulation Methods and Aberration Description
Thus, the amplitude distribution of the aberrated beam behind IO becomes:  
\[\mathrm{E_{GI}(x_{0},y_{0})=E_{GO}(x_{0},y_{0})\exp\bigl{[}-i2\pi W(x_{0},y_{0} )\bigr{]}t_{I}(x_{0},y_{0})} \tag{6}\]  
\(t_{I}\) is the transmittance function of IO (with focal length \(f_{I}\)):  
\[t_{I}(x_{0},y_{0})=\exp\bigl{[}-i\frac{k}{2f_{I}}(x_{0}{{}^{2}}+y_{0}{{}^{2}}) \bigr{]} \tag{7}\]  
Based on the law of Fresnel diffraction, the Gauss beam behind IO at distance \(z\) is:  
\[\mathrm{E_{G}(x^{\prime},y^{\prime},z)} =\frac{\exp(i\mathrm{k}z)}{\mathrm{i}\lambda z}\iint E_{GI}(x_{0}, y_{0})\exp\bigl{\{}i\frac{k}{2z}\bigl{[}(x^{\prime}\ -x_{0})^{2}+(y^{\prime}\ -y_{0})^{2}\bigr{]}\bigr{\}}dx_{0}dy_{0} \tag{8}\] \[=Z\bigl{\{}E_{GI}(x_{0},y_{0}),z\bigr{\}}\]  
Finally, the two-photon excited fluorescence intensity is defined as the fourth power of amplitude:  
\[\mathrm{I_{G}(x^{\prime},y^{\prime},z)=}\bigl{|}E_{G}(x^{\prime},y^{\prime},z )\bigr{|}^{4} \tag{9}\]  
For TPBB, in order to modulate the Bessel beam, an axicon phase pattern is generated by SLM:  
\[\mathrm{\Phi_{B}(u_{0},v_{0})=\frac{2\pi\sqrt{u_{0}{{}^{2}}+v_{0}{{}^{2}}}}{r_ {0}}} \tag{10}\]  
Based on [PERSON]'s research [13], the modulation parameter \(r_{0}{=2\pi/\alpha=2\pi/k\theta_{0}}\) is the width of the concentric ring in this pattern, where \(k=2\pi/\lambda\) and \(\lambda\) is the wavelength of the illumination laser. \(\theta_{0}\) is the oblique angle of equiphase surface shown in Figure 1b, and \(\alpha\) represents the propagation vector of the plane wave components of the angular spectrum to the transverse plane. If SLM is illuminated uniformly with Gaussian beam, this pattern produces a complex amplitude:  
\[\mathrm{x}=\mathrm{rcos}\theta,\,\,\mathrm{y}=\mathrm{rsin}\theta,\,\,\mathrm{ and}\,\,0\leq\mathrm{r}=\sqrt{\mathrm{x}^{2}+\mathrm{y}^{2}}\leq 1 \tag{4}\]0  
where \(w_{B0}\) is the waist radius of input beam and \(A_{B0}\) is the amplitude. Since SLM is conjugated to the sample plane, the complex amplitude at the rear focal plane of IO can be calculated by its 2D Fourier transform [17]:  
\[\mathrm{x}=\mathrm{rcos}\theta,\,\,\mathrm{y}=\mathrm{rsin}\theta,\,\,\mathrm{ and}\,\,0\leq\mathrm{r}=\sqrt{\mathrm{x}^{2}+\mathrm{y}^{2}}\leq 1 \tag{4}\]1 \[\mathrm{x}=\mathrm{rcos}\theta,\,\,\mathrm{y}=\mathrm{rsin}\theta,\,\,\mathrm{ and}\,\,0\leq\mathrm{r}=\sqrt{\mathrm{x}^{2}+\mathrm{y}^{2}}\leq 1 \tag{4}\]2  
Thus, the complex amplitude of the illumination beam through the aberrations can be expressed as  
\[\mathrm{x}=\mathrm{rcos}\theta,\,\,\mathrm{y}=\mathrm{rsin}\theta,\,\,\mathrm{ and}\,\,0\leq\mathrm{r}=\sqrt{\mathrm{x}^{2}+\mathrm{y}^{2}}\leq 1 \tag{4}\]3  
Because \(E_{B2}\) (\(u_{1},v_{1}\)) is located at a distance \(f_{I}\) in front of IO, the complex amplitude immediately behind IO can be written as:  
\[\mathrm{x}=\mathrm{rcos}\theta,\,\,\mathrm{y}=\mathrm{rsin}\theta,\,\,\mathrm{ and}\,\,0\leq\mathrm{r}=\sqrt{\mathrm{x}^{2}+\mathrm{y}^{2}}\leq 1 \tag{4}\]4  
In the same way, to find the distribution behind IO at distance \(z\), Equation (8) is applied:  
\[\mathrm{x}=\mathrm{rcos}\theta,\,\,\mathrm{y}=\mathrm{rsin}\theta,\,\,\mathrm{ and}\,\,0\leq\mathrm{r}=\sqrt{\mathrm{x}^{2}+\mathrm{y}^{2}}\leq 1 \tag{4}\]5Finally, the two-photon excited fluorescence intensity is defined as the fourth power of amplitude:  
\[\mathrm{x}=\mathrm{rcos}\theta,\,\,\mathrm{y}=\mathrm{rsin}\theta,\,\,\mathrm{ and}\,\,0\leq\mathrm{r}=\sqrt{\mathrm{x}^{2}+\mathrm{y}^{2}}\leq 1 \tag{4}\]6  
It is worth noting that the simulated aberration should be set at the plane, which is conjugated to WFS, and the diameter of the aberration should be equal to the diameter of the illumination beam at this plane.