## 2 Materials and Methods
### Numerical Simulations of Crater Formation when Multiple Successive Pulses Imping on the Surface  
In order to establish how multiple successive pulses form the ablation crater, a semi-empirical model was constructed. The model takes the core basis that if a pulse is focused on the surface of the sample, the energy density (fluence, F) is not constant at different places on the surface since the intensity distribution of the pulse is non-uniform. Since the energy density varies at different positions in the impingement zone, the removed mass (volume) must also be different. In order to estimate how much volume is removed from the surface, a simple experiment was performed by firing pulses of different energy at the sample and the shape was registered with an optical profiler (see Section 2.1). Knowing the actual intensity distribution (which was measured previously with a CCD camera), the shape of the formed crater is put side-by-side with the impinging fluence distribution and a graph depicting how much the different portions of the sample have increased in depth in relation to different fluence values. Figure 2 demonstrates this approach. The left side (a) of the figure shows a crater formed by a Gaussian pulse. Only half of the crater and Gaussian distribution is displayed since the other half is symmetrical and does not carry any significant information. This graph represents only experimental data of a crater that was formed after a single pulse (having energy of 160 uJ, fluence = 29 J/cm\({}^{2}\)) struck the surface of the sample. The red curve shows how the fluence was distributed in infinitesimally small regions inside the impingement zone before the first pulse struck the sample, i.e., the sample was flat before the first pulse arrived and the energy contained in the region \(dx\) produced a change in the \(Z\) coordinate of the surface \(Z(x)\). From this, with rearranging, the left side (b) of Figure 2 is constructed showing how much the depth changes given a specific energy density (or fluence). We observed that the dependency takes on a logarithmic form (for greater values than the ablation threshold) and two parameters are sufficient to fit the curve. For the second and all consecutive pulses, the energy stored in the area \(dx\) was projected onto the area \(dl\) (see Figure 2). Note that although in real-world experiments dl is an area, for ease of computational load, a two-dimensional case was considered, and therefore instead of an area, \(dx\) and \(dl\) are line segments. We did not take into account the defocusing of the beam, since throughout our experiments an f = 100 mm focusing lens was used (the input beam was 3.5 mm at 1/e\({}^{2}\)) which focused the beam down to 37.4 \(\mu\)m (1/e\({}^{2}\)). From this, the resulting Rayleigh range (\(z_{R}=\pi\rho_{0}^{2}/\lambda\)) of approximately 1.1 mm was calculated, and therefore defocusing was neglected. By taking into account the above-mentioned considerations, a semi-empirical model (since the coefficients of the logarithmic function are obtained from the experiment) is constructed that predicts the evolution of the craters when multiple pulses hit the surface. The changing area \(dl\) is calculated after each pulse and the ratio of energy projected onto this area is taken and plugged into the empirical logarithmic function. This produces a change (\(dZ\)) in the surface \(Z(x)\) (see Equation (3)). The sum of these changes produces a new function \(Z(x)_{n}\) which denotes the shape of the surface after each n-th pulse (see Equation (4)). In the equations, \(F(x)=F_{0}\exp\left(-4\log(2)\left(x/\rho_{0}^{2}\right)\right)\) takes the form of a Gaussian shape (\(F_{0}\) is the amplitude of the Gaussian distribution, \(\rho_{0}\) is the width of the beam), however different distributions are also valid as will be displayed later.