## 2 Materials and Methods
### Numerical Model of Beam Propagation in a Transparent Medium  
In order to reveal the transformation of the high-intensity pulses within a thin water layer, we performed some simplified numerical simulations, considering the paraxial beam propagation along the z-axis of a linearly polarized laser beam. We model the linearly polarized beam with cylindrical symmetry around the propagation axis \(z\) by the envelope \(A\) of the electric field, written as \(A=Re[Aexp(-ik_{L}z-i\omega_{0}t)]\), where \(k_{L}=n_{0}\omega_{0}/c\) and \(\omega_{0}\) are the wave number and frequency of the carrier wave, and \(n_{0}\) denotes the refraction index of water at \(\omega_{0}\). The energy \(E_{in}\) pulse had a Gaussian amplitude profile in the temporal and spatial domains and a temporal half width \(\tau_{P}\) (FWHM) duration \(\tau_{FWHM}=t_{P}\) (2 log 2)\({}^{1/2}\) and \(\rho_{FWHM}=\rho_{0}\) (2 log 2)\({}^{1/2}\).  
\[A=A_{0}\exp(-\frac{t^{2}}{\tau_{0}^{2}}-\frac{r^{2}}{\rho_{0}^{2}}). \tag{1}\]  
Under the conditions of the experiments regarding the propagation of the beam within the water layer, the effects related to the pulse group velocity and dispersion were ruled out, since the  
\begin{table}
\begin{tabular}{c c} \hline \hline
**Parameters** & **Value Range** \\ \hline Position of focus relative to the upper glass interface with water & \(-\)1.5\(-\)1.5 (mm) \\ Scanning speed along the surface & 20\(-\)1500 (mm/s) \\ Repetition rate of the laser & 33\(-\)200 kHz \\ Thickness of the water layer & 0.3\(-\)1.5 mm \\ \hline \hline \end{tabular}
\end{table}
Table 1: Parameters varied during femtosecond laser micromachining experiments. The average power was kept constant at 10 W.  
Figure 1: Micromachining setup. The laser beam is directed by a scanner and is focused using an f-theta lens (f = 100 mm) onto the sample that is covered with a thin layer of water.  
spectral width of 260 fs is below 10 nm and water having a group velocity dispersion coefficient at 1030 nm of 29.96 fs2/mm does not induce pulse dispersion over the thickness of the water layer. Therefore, the goal of the research was to predict the alterations to the intensity distribution within the water layer which impinged on the surface of the sample, since it is believed that the change in the intensity distribution is the main driving factor for the increase in the micromachining throughput for ablation in water. The energy distributions striking the sample were acquired by solving the modified non-linear Schrodinger equation (NLS) in 2 + 1 dimensions. The equation includes diffraction (first term), self-focusing (second term), and non-linear losses (third term) due to multi-photon absorption. For the scalar envelope \(A(r,\,z,\,t)\) evolving along the propagation axis \(z\) according to the non-linear envelope, the equation is:  
\[\frac{\partial A}{\partial z}=\left(\frac{i}{2k_{L}}\Delta_{r}+ik_{L}n_{2}|A|^ {2}-\frac{\beta^{(K)}}{2}|A|^{2K-2}\right)A, \tag{2}\]  
where, \(r\) and \(\Delta_{r}\) stands for the radial coordinate and Laplace operator, respectively, \(\beta^{(\mathrm{K})}[\mathrm{cm^{2K\text{-}3}W^{1\text{-}K}}]\) is the multi-photon absorption (MPA) coefficient and \(n_{2}|A|^{2}\) responds to the intensity-dependent changes of the refractive index. We evaluate the amplitude of the beam at different distances \(L\) along the propagation direction for the incident beam power (at the Gaussian maximum). The input power is computed from the energy and pulse duration \(P_{in}=E_{in}/\tau\). Accordingly, the initial pump amplitude becomes \(A_{0}=\sqrt{2P_{0}/\pi\,\,\rho_{0}^{2}}\), where \(\rho_{0}\) stands for the beam waist.