2.1 Standard PALS laser setup

The high-power iodine photodissociation laser system at the PALS facility operates at the fundamental wavelength $\lambda =1315.2\;\mathrm{nm}$ with bandwidth 170 pm and delivers laser pulses with energies up to approximately $700\;\mathrm{J}$ and duration ${\tau}_{\mathrm{FWHM}}\sim 300\;\mathrm{ps}$ ^[ Reference Jungwirth, Cejnarova, Juha, Kralikova, Krasa, Krousky, Krupickova, Laska, Masek, Mocek, Pfeifer, Präg, Renner, Rohlena, Rus, Skala, Straka and Ullschmied ^{8 ]}.

The laser pulse is extracted from a pulse train sequence generated by a mode-locked oscillator. It is then amplified by five flashlamp-pumped laser amplifiers, resulting in a flat-top fluence distribution at the output. To eliminate high spatial frequencies in the fluence distribution, six vacuum spatial filters are used; these filters transfer the beam between each pair of adjacent laser amplifiers, with the final spatial filter transferring the beam from the last amplifier to a focusing lens in the target chamber. The laser system layout is shown in Figure 1.

Figure 1

Schematic layout of the iodine photodissociation laser at the Prague Asterix Laser System facility. The oscillator and preamplifier part is simplified. The image is adapted with permission from Ref. [Reference Dostal, Dudzak, Pisarczyk, Pfeifer, Huynh, Chodukowski, Kalinowska, Krousky, Skala, Hrebicek, Medrik, Golasowski, Juha and Ullschmied26].

The output laser beam with a diameter of 290 mm is focused by an $f/2.2$ aspherical focusing lens with focal length $f=627\;\mathrm{mm}$ at $\lambda =1315.2\;\mathrm{nm}$ . The lens is designed to create a diffraction-limited focus with focal length $f=600\;\mathrm{mm}$ at $\lambda =438\;\mathrm{nm}$ (third harmonic) in the geometrical optics approximation. Therefore, we expect aberrations at $\lambda =1315.2\;\mathrm{nm}$ , specifically the spherical aberration. During the experiment, we did not have a lens corrected for 1315.2 nm at our disposal.

We determine the focus position experimentally by inserting a metal Hartmann plate into the collimated beam immediately upstream of the lens. The Hartmann plate cuts out five beams ( $\mathrm{\unicode{xf8}}=20\;\mathrm{mm}$ ) located on a circle perimeter ( $\mathrm{\unicode{xf8}}=200\;\mathrm{mm}$ ) around the center of the main beam. We conduct a series of four laser shots on an aluminum target at the target chamber center (TCC). For each shot, the focusing lens is positioned at different distances ( $-$ 6, $-$ 3, 3 and 6 mm) relative to the expected focus position. We capture microscope images of the resulting ablation patterns, where each pattern contains five holes corresponding to the five beams, and we measure the distances between the centers of the leftmost and rightmost craters. We plot the measured distances (negative when the ablated pattern is inverted through the center) against the focusing lens positions and perform a linear fit. Its intersection point with the $z$ -position axis determines the focus position.

Beam parameters are monitored online using an attenuated part of the laser beam that is transmitted through a high-reflectivity mirror after the last laser amplifier. The attenuated laser pulse is split and measured by diagnostic devices. A calibrated calorimeter measures the pulse energy. Four photodiodes with different trigger timings and different attenuation filters in front of them measure the laser pulse power distribution in the prepulse before the main laser peak. Together, the photodiodes measure the power distribution in the range of 10 ns before and 0.8 ns after the main laser peak. The highest power peak of the laser pulse is measured by a streak camera. A slit of width 30 μm cuts out a rectangle from the laser beam, which irradiates a photocathode and the emitted electrons are streaked in time. We join the normalized power distributions measured by the photodiodes and the streak camera at $0.05\times$ the peak power to characterize the entire power distribution.

2.2 Measuring the fluence distribution at and around the focus position

We measure the fluence distribution in multiple planes at and around the focus position in the target chamber by far-field imaging (i.e., projecting a magnified image of the fluence distribution onto a camera). Firstly, we build the imaging setup (Figure 2) and determine its magnification factor. We place a calibration target (grid with 10 μm spacing) at the TCC, illuminate it with a defocused alignment neodymium-doped yttrium lithium fluoride (Nd:YLF) laser beam (CrystalLaser IRCL-1W-1313, $\lambda =1313\;\mathrm{nm}$ ) and project the grid’s magnified image onto a near-infrared charge-coupled device (CCD) camera (Goldeye G-032 SWIR TEC, 25 μm pixel size, number of pixels $636\times 508$ ) using a microscope objective without a tube lens ( $50$ × Mitutoyo Plan Apochromatic Infinity Corrected Long Working Distance Objective, $f=4\;\mathrm{mm}$ , $\mathrm{NA}=0.55,\ f\text{-}\mathrm{number}\approxeq 0.9$ , corrected for visible wavelengths $\approx 400$ –700 nm). The microscope objective is corrected for imaging of visible light into infinity. Although we used the near-infrared wavelength of ( $1315.2\pm 0.1$ ) nm and the image distance was 340 mm, the calibration grid was imaged with high quality without significant aberrations. We measured the edge-spread function of the imaging system and determined that the modulation transfer function at 10% contrast is reached for spatial frequency at 500 line pairs per mm ( $\mathrm{MTF}10=500\;\mathrm{lp}/\mathrm{mm}$ )^[ Reference Boreman ^{27 ,} Reference Smith ^{28 ]}. This was made possible by the narrow bandwidth that minimized the potential chromatic aberration, good alignment of the objective and small beam filling factor of the entrance lens that reduced the spherical and other aberrations. After determining the magnification factor of 84 and the scale of $0.3\;\mu \mathrm{m}/\mathrm{px}$ , we focus the alignment laser beam at the TCC and insert large-diameter ( $\mathrm{\unicode{xf8}}=340\;\mathrm{mm}$ ) neutral-density optical filters into its collimated-beam path (upstream of the focusing lens) to check that they do not visibly change the far-field. The alignment laser beam diameter at these optical filters is 290 mm, the same as for the iodine laser. When all seven large optical filters are inserted, they attenuate the fully amplified iodine laser pulse from $700\;\mathrm{J}$ down to approximately $1\;\mathrm{mJ}$ . The pulse then passes through the focusing lens, and the fluence distribution at the TCC is magnified and imaged onto the camera, all in air. We decrease the laser pulse energy so we do not ignite a spark in air with the focused laser beam. We insert small neutral-density $50$ mm $\times$   $50$ mm square filters (Hoya) in front of the camera to further attenuate the beam and stay below the camera saturation limit. We translate the focusing lens along the $z$ -axis (beam propagation direction) between the shots to create a $z$ -scan of fluence distributions along the beam caustic.

Figure 2

Diagram of the setup for fluence distribution measurements along the beam caustic. We shift the focusing lens within a range of $\pm$ 300 μm around the expected best focus position to perform the $z$ -scan. Distances among the object plane, the front principal plane of the microscope objective and the camera are calculated based on the magnification of the calibration target (84×).

3.1 Retrieving the complex amplitude from the measured fluence distributions

We retrieved complete information about the complex amplitude of the electric field from several fluence distribution measurements at different $z$ positions along the beam propagation direction. Subsequently, we propagated the beam throughout its caustic.

Phase retrieval approaches can extract phase from a set of fluence distribution measurements. We used a multi-plane phase retrieval algorithm with dynamic input–output mixing (PhaRe)^[ Reference Chalupsky, Bohacek, Burian, Hajkova, Hau-Riege, Heimann, Juha, Messerschmidt, Moeller, Nagler, Rowen, Schlotter, Swiggers, Turner and Krzywinski ^{19 ]} using functions adapted from the LightPipes library for Python^{[ 29 ]}. The main idea is illustrated in Figure 3.

Figure 3

Illustrative flowchart of the jth iteration loop of PhaRe applied to three fluence distribution measurements. The illustrative iteration starts at the position ${z}_1$ , but the final retrieved complex amplitude is independent of the starting position. Green circles represent the measured electric field amplitudes $\left|{E}_n^{\mathrm{M}}\right|$ , yellow rectangles represent the complex amplitudes $\left|{E}_n^{\mathrm{P}}\right|\exp \left(i{\varphi}_n\right)$ propagated from the previous $z$ position and blue rectangles represent the temporary retrieved complex amplitudes $\left|{E}_n\right|\exp \left(i{\varphi}_n\right)$ . The temporary retrieved electric field amplitude is calculated as a linear combination of the corresponding measured and propagated electric field amplitudes ${\alpha}_j\left|{E}_n^{\mathrm{M}}\right|+{\beta}_j\left|{E}_n^{\mathrm{P}}\right|$ . The forward and backward propagation steps are highlighted in blue and red colors, respectively. The black arrows represent copying of the complex amplitude. The image is adapted with permission from Ref. [Reference Chalupsky, Bohacek, Burian, Hajkova, Hau-Riege, Heimann, Juha, Messerschmidt, Moeller, Nagler, Rowen, Schlotter, Swiggers, Turner and Krzywinski19], copyright 2015, American Physical Society.

Firstly, we preprocess each of the total number of $N$ input fluence distributions measured at ${z}_n$ positions ( $n\in \left\{0,1,\dots, N-1\right\}$ ). We use ImageJ software to subtract the dark image (the laser, amplifier flashlamps and lights are off) and then the mean of the background. We set any negative-value pixels to zero since those pixels do not carry information about the signal. The residual background remaining after the dark image subtraction is caused by broadband radiation from the flashlamps that pump the last laser amplifier. A small part of this incoherent radiation propagates into the experimental chamber and onto the camera, where it is measured as shot-to-shot-variable background. After background subtraction, we manually select the non-irradiated region, where we set the signal to zero. In Python, we smooth each image using a median filter $5\times 5$ and correct the measured nonlinear signal using a calibration function (described in Section 3.5). We insert the image into a larger array filled with zeros and shift the beam centroid to the center to ensure that beam momentum is conserved. To meet the condition of pulse energy conservation when propagating the beam in free space, we normalize each fluence distribution by dividing by pulse energy. We calculate the pulse energy using Equation (6) by multiplying the peak value of the fluence distribution with the effective area. Another possibility is to use the measured pulse energy from the on-shot diagnostic. We address the different normalizations in Section 5.1. Finally, we take square roots of the normalized fluence distributions, which are proportional to the measured (superscript M) electric field amplitudes $\left|{E}_n^{\mathrm{M}}\right|$ at positions ${z}_n$ .

Secondly, we iteratively propagate the complex amplitude ${E}_n\left(x,y\right)$ consisting of amplitude $\left|{E}_n^{\mathrm{M}}\left(x,y\right)\right|$ and phase ${\varphi}_n\left(x,y\right)$ ,

(1) $$\begin{align}{E}_n\left(x,y\right)=\left|{E}_n^{\mathrm{M}}\left(x,y\right)\right|\exp \left[i{\varphi}_n\left(x,y\right)\right],\end{align}$$

to a neighboring position $m=n+1$ (forward steps) or $m=n-1$ (backward steps). The iteration cycle starts at an arbitrary position ${z}_{n_0}$ . The complex amplitude ${E}_{n_0}\left(x,y\right)$ is propagated using Equation (2) from the position ${z}_{n_0}$ forward into the position ${z}_m$ , where $m={n}_0+1$ . The next forward step starts at the position ${z}_{n_0+1}$ and ends at the position ${z}_{n_0+2}$ . The forward steps continue until the rightmost position ${z}_{N-1}$ is reached. This is followed by the backward steps where the complex amplitude is propagated through the positions ${z}_{N-2},{z}_{N-3},\dots, {z}_0$ until the leftmost position is reached. Then, the last forward steps propagate the complex amplitude through the positions ${z}_0,{z}_1,\dots, {z}_{n_0-1}$ to the starting position ${z}_{n_0}$ . This marks the end of one iteration cycle and the following cycle starts at the same position ${z}_{n_0}$ . We propagate ${E}_n\left(x,y\right)$ using the angular spectrum method:

(2) $$\begin{align}{E}_m^{\mathrm{P}}\left(x,y\right)={\mathrm{\mathcal{F}}}^{-1}\left\{\mathrm{\mathcal{F}}\left[{E}_n\left(x,y\right)\right]{H}_m\left({f}_x,{f}_y\right)\right\},\end{align}$$

where the superscript P denotes that the resulting complex amplitude is a result of numerical propagation, $\mathrm{\mathcal{F}},{\mathrm{\mathcal{F}}}^{-1}$ are Fourier and inverse Fourier transforms, respectively, ${H}_m\left({f}_x,{f}_y\right)$ is a transfer function of free space and ${f}_x,\ {f}_y$ are spatial frequencies. We used the Rayleigh–Sommerfeld transfer function in Equation (2):

(3) $$\begin{align}{H}_m\left({f}_x,{f}_y\right)=\left\{\begin{array}{@{}r@{}}\exp \left\{ ik{\zeta}_m\left[\sqrt{1-{\left(\lambda {f}_x\right)}^2-{\left(\lambda {f}_y\right)}^2}\right]\right\}\\[10pt]{}\kern1em \mathrm{for}\sqrt{f_x^2+{f}_y^2}<\frac{1}{\lambda},\\[6pt]{}0\kern9.12em \mathrm{otherwise},\end{array}\right.\end{align}$$

where $k=\left|\left({k}_x,{k}_y,{k}_z\right)\right|=2\pi /\lambda$ is the wavenumber (wavevector amplitude) and ${\zeta}_m={z}_m-{z}_n$ is the propagation distance. The relation between spatial frequencies and the wavevector is $2\pi \left({f}_x,{f}_y\right)=\left({k}_x,{k}_y\right)$ . Spatial frequencies that are too high do not propagate in free space, so their transfer function is set to zero.

After each propagation step, we replace the propagated electric field amplitude by its linear combination with the measured electric field amplitude while keeping the propagated phase, as shown in the following equation:

(4) $$\begin{align}{E}_m=\left({\alpha}_j\left|{E}_m^{\mathrm{M}}\right|+{\beta}_j\left|{E}_m^{\mathrm{P}}\right|\right)\exp \left(i{\varphi}_m\right).\end{align}$$

Coefficients ${\alpha}_j,{\beta}_j$ are dynamic scaling factors dependent on the iteration number $j$ . They always obey the condition ${\alpha}_j+{\beta}_j=1$ . This dynamic input–output mixing is specific for the phase retrieval approach in Ref. [Reference Chalupsky, Bohacek, Burian, Hajkova, Hau-Riege, Heimann, Juha, Messerschmidt, Moeller, Nagler, Rowen, Schlotter, Swiggers, Turner and Krzywinski19].

Throughout the iterative process, we evaluate the goodness of phase retrieval (GoPR) by averaging normalized cross-correlations (represented by $\star$ ) of the propagated and measured electric field amplitudes when they fully overlap:

(5) $$\begin{align}{\mathrm{GoPR}=\frac{1}{N}\sum \limits_{m=1}^N\left(|{E}_m^{\mathrm{P}}|\star |{E}_m^{\mathrm{M}}|\right)\left(0,0\right)\nonumber\\ {}=\frac{1}{N}\sum \limits_{m=1}^N\frac{\iint_{{\mathrm{\mathbb{R}}}^2}\left|{E}_m^{\mathrm{P}}\left(x,y\right)\right|\left|{E}_m^{\mathrm{M}}\left(x,y\right)\right|\mathrm{d}x\;\mathrm{d}y}{\sqrt{\iint_{{\mathrm{\mathbb{R}}}^2}{\left|{E}_m^{\mathrm{P}}\right|}^2 xy\iint {\left|{E}_m^{\mathrm{M}}\right|}^2\mathrm{d}x\;\mathrm{d}y}}.}\end{align}$$

The GoPR equals 1 when the propagated and measured electric field amplitudes are equal ( $|{E}_m^{\mathrm{P}}| = |{E}_m^{\mathrm{M}}|$ ) at all ${z}_m$ positions.

In the first few iterations, ${\alpha}_j=1$ and the measured $\left|{E}_m^{\mathrm{M}}\right|$ fully replaced the propagated $\left|{E}_m^{\mathrm{P}}\right|$ until the GoPR stopped improving (10 iterations in our case). In the next 90 iterations, dynamic input–output mixing took place; ${\alpha}_j$ linearly decreased and the result was weighted more towards the propagated $\left|{E}_m^{\mathrm{P}}\right|$ . In the last five iterations ${\alpha}_j=0$ and ${\beta}_j=1$ , and the complex amplitudes ${E}_m$ did not change between iteration cycles. The result became a self-consistent solution of the Helmholtz equation. This self-consistency means that the complex amplitude ${E}_m$ at position ${z}_m$ in the last iterations ( ${\beta}_j=1$ ) is the same as the complex amplitude we calculate by propagating the final retrieved complex amplitude to this position from any of the other positions.

Care had to be taken when retrieving the phase of highly convergent/divergent beams. A parabolic term in the wavefront (causing focusing/divergence) is zero at the focus position and the wavefront is roughly flat. However, the parabolic term grew significantly as the beam moved away from the focus position. When the wavefront curvature became too large, its high spatial frequency components were insufficiently sampled and the simulated phase exhibited spatial aliasing. The insufficiently sampled complex amplitudes could not be propagated correctly. This applied to fluence distributions measured more than 400 μm away from the focus position ( $|z| \ge 400\;\mu \mathrm{m}$ ). We did not include those measurements in PhaRe.

In addition, we tested a basic phase retrieval algorithm that propagates the complex amplitude between neighboring positions without the dynamic input–output mixing^[ Reference Allen and Oxley ^{15 ]}. The GoPR stopped improving after 10 iterations and the retrieved electric field amplitudes were visibly different from the measured ones. This was because a few imperfect images in the dataset were used as an input in each iteration and their influence was not suppressed in the algorithm. This was mitigated by the dynamic input–output mixing of the retrieved and measured electric field amplitudes^[ Reference Chalupsky, Bohacek, Burian, Hajkova, Hau-Riege, Heimann, Juha, Messerschmidt, Moeller, Nagler, Rowen, Schlotter, Swiggers, Turner and Krzywinski ^{19 ]}. In addition, we tested another algorithmic scheme^[ Reference Hansen ^{18 ,} Reference Jamal and Hansen ^{30 ]} in which a reference position is chosen. The complex amplitudes are propagated from the reference position to all $N$ positions simultaneously and the corresponding measured electric field amplitudes are substituted. The obtained complex amplitudes are propagated back to the reference position and averaged. A new iteration then starts. However, we retrieved significantly different complex amplitudes based on the specific choice of the reference position. This represented a major limitation, and we did not pursue the method further. However, there is a similar approach, which we did not test, described by Brady and Fienup^[ Reference Brady and Fienup ^{17 ]}. The complex amplitude is propagated from the reference position to all the other non-reference positions. Then, it minimizes a weighted sum-of-squared-errors metric calculated between the propagated and measured electric field amplitudes in the non-reference positions.

3.2 Calculating the intensity distribution from the fluence and power distributions

We indirectly determined the intensity distribution in the beam caustic by combining a time-dependent power measurement and the spatial fluence distribution of the complex amplitude output from PhaRe.

One way to describe the widths of laser beam distributions (fluence $F\left(x,y,z\right)$ , power $P(t)$ and intensity $I\left(x,y,z,t\right)$ ) with a single number is by using effective quantities, such as the effective area ${A}_{\mathrm{eff}}(z)$ , effective pulse duration ${\tau}_{\mathrm{eff}}$ and effective time-area product ${\mathrm{TAP}}_{\mathrm{eff}}(z)$ , respectively^[ Reference Pirozhkov, Fukuda, Nishiuchi, Kiriyama, Sagisaka, Ogura, Mori, Kishimoto, Sakaki, Dover, Kondo, Nakanii, Huang, Kanasaki, Kondo and Kando ^{6 ,} Reference Chalupský, Burian, Hájková, Juha, Polcar, Gaudin, Nagasono, Sobierajski, Yabashi and Krzywinski ^{31 ]}. Effective quantities relate integrals of distributions (all equal to the pulse energy ${E}_{\mathrm{pulse}}$ ) to the distributions’ peaks ( ${F}_{\mathrm{peak}}(z)$ , ${P}_{\mathrm{peak}}$ , ${I}_{\mathrm{peak}}(z)$ ). The equation relating the effective area and peak of the fluence distribution is as follows:

(6) $$\begin{align}{{E}_{\mathrm{pulse}}=\underset{{\mathrm{\mathbb{R}}}^2}{\iint }F\left(x,y,z\right)\mathrm{d}x\;\mathrm{d}y\nonumber\\ {}={F}_{\mathrm{peak}}(z)\underset{{\mathrm{\mathbb{R}}}^2}{\iint }f\left(x,y,z\right)\mathrm{d}x\;\mathrm{d}y\nonumber\\ {}={F}_{\mathrm{peak}}(z){A}_{\mathrm{eff}}(z),}\end{align}$$

where $F\left(x,y,z\right)$ and $f\left(x,y,z\right)$ are the non-normalized and normalized fluence distributions at position $z$ , respectively.

The equation for the power distribution is as follows:

(7) $$\begin{align}{E}_{\mathrm{pulse}}={\int}_{\mathrm{\mathbb{R}}}P(t)\mathrm{d}t={P}_{\mathrm{peak}}{\int}_{\mathrm{\mathbb{R}}}p(t)\mathrm{d}t={P}_{\mathrm{peak}}{\tau}_{\mathrm{eff}},\end{align}$$

where $P(t)$ and $p(t)$ are the non-normalized and normalized power distributions, respectively.

Finally, the equation for the intensity distribution is as follows:

(8) $$\begin{align}{{E}_{\mathrm{pulse}}=\underset{{\mathrm{\mathbb{R}}}^3}{\iiint }I\left(x,y,z,t\right)\mathrm{d}x\;\mathrm{d}y\;\mathrm{d}t\nonumber\\ {}={I}_{\mathrm{peak}}(z)\underset{{\mathrm{\mathbb{R}}}^3}{\iiint }i\left(x,y,z,t\right)\mathrm{d}x\;\mathrm{d}y\;\mathrm{d}t\nonumber\\ {}={I}_{\mathrm{peak}}(z)\underset{{\mathrm{\mathbb{R}}}^2}{\iint }f\left(x,y,z\right)\mathrm{d}x\;\mathrm{d}y{\int}_{\mathrm{\mathbb{R}}}p(t)\mathrm{d}t\nonumber\\ {}={I}_{\mathrm{peak}}(z){A}_{\mathrm{eff}}(z){\tau}_{\mathrm{eff}}={I}_{\mathrm{peak}}(z){\mathrm{TAP}}_{\mathrm{eff}}(z),}\end{align}$$

where $I\left(x,y,z,t\right)$ and $i\left(x,y,z,t\right)$ are the non-normalized and normalized intensity distributions, respectively.

The spatial and temporal pulse energy distributions are generally not independent functions. The spatio-temporal coupling is a phenomenon that cannot be neglected for ultrashort laser pulses when using dispersive optics. We estimate its effect for our 350 ps laser pulse with a bandwidth of 170 pm (30 GHz). Pulse chirp from group velocity dispersion is negligible (chirp of 0.2 fs in 1 m of BK7 glass). Much larger influence comes from a radial group delay that quantifies the time delay at the focus position between a marginal ray pulse traveling through the lens edge and a chief ray pulse traveling through the lens center^[ Reference Bor ^{32 –} Reference Bahk, Bromage and Zuegel ^{34 ]}. The estimated radial group delay for our BK7 lens with focal length 627 mm at $\lambda =1315.2\;\mathrm{nm}$ and lens diameter 290 mm is 2 ps, which is still negligible compared to the pulse length of 350 ps. Therefore, we neglect the spatio-temporal coupling and we approximate the intensity distribution as a separable function that can be factorized into spatial and temporal functions. The intensity distribution is calculated from Equations (6)–(8) as follows:

(9) $$\begin{align}{I\left(x,y,z,t\right)=\frac{E_{\mathrm{pulse}}}{{\mathrm{TAP}}_{\mathrm{eff}}}i\left(x,y,z,t\right)\nonumber\\ =\frac{E_{\mathrm{pulse}}}{A_{\mathrm{eff}}{\tau}_{\mathrm{eff}}}f\left(x,y,z\right)p(t).}\end{align}$$

From the intensity distribution, we can calculate the fraction of pulse energy ${E}_{\mathrm{frac}}$ that corresponds to the part of the beam exceeding a given threshold intensity ${I}_{\mathrm{th}}$ at position $z$ :

(10) $$\begin{align}{E}_{\mathrm{frac}}\left({I}_{\mathrm{th}},z\right)=\underset{I>{I}_{\mathrm{th}}}{\iiint }I\left(x,y,z,t\right)\mathrm{d}x\;\mathrm{d}y\;\mathrm{d}t.\end{align}$$

3.3 Calculating the complex amplitude at the focusing lens

We propagated the retrieved complex amplitude from the focus position to the position of the focusing lens. The resulting fluence distribution represents a product of the lens transmission profile and the collimated-beam fluence distribution upstream of the lens. The obtained wavefront is a sum of the collimated-beam wavefront and the wavefront aberrations of the focusing lens, that is, its deviations from the ideal-lens parabolic wavefront.

The relation between the complex amplitude at the position of an ideal lens and the complex amplitude at its focal plane is given by the Fourier transform^[ Reference Goodman ^{35 ]}:

(11) $$\begin{align}{{E}_{\mathrm{f}}\left({x}_{\mathrm{f}},{y}_{\mathrm{f}}\right)=\frac{\exp \left\{j\left[\frac{k}{2f}\left({x}_{\mathrm{f}}^2+{y}_{\mathrm{f}}^2\right)+ kf\right]\right\}}{j\lambda f} \nonumber\\\quad\times \mathrm{\mathcal{F}}\left[{E}_{\mathrm{l}}\left({x}_{\mathrm{l}},{y}_{\mathrm{l}}\right)\right]\left(\frac{x_{\mathrm{f}}}{\lambda f},\frac{y_{\mathrm{f}}}{\lambda f}\right),}\end{align}$$

where ${E}_{\mathrm{f}}$ and ${E}_{\mathrm{l}}$ are complex amplitudes in the focal plane and at the position just downstream of the lens with the focus term subtracted, respectively, and $f$ is the focal length of the ideal lens. We calculated ${E}_{\mathrm{l}}$ from ${E}_{\mathrm{f}}$ by inverting Equation (11):

(12) $$\begin{align}{{E}_{\mathrm{l}}\left({x}_{\mathrm{l}},{y}_{\mathrm{l}}\right)={\mathrm{\mathcal{F}}}^{-1}\left\{\exp \left\{-j\left[\frac{k}{2f}\left({x}_{\mathrm{f}}^2+{y}_{\mathrm{f}}^2\right)+ kf\right]\right\} j\lambda f \right.\nonumber\\\quad\times \left.\left[{E}_{\mathrm{f}}\left({x}_{\mathrm{f}},{y}_{\mathrm{f}}\right)\right]\vphantom{\frac{k}{2f}}\right\}\left(\frac{x_{\mathrm{l}}}{\lambda f},\frac{y_{\mathrm{l}}}{\lambda f}\right).}\end{align}$$

3.4 Strehl ratio

We calculated the Strehl ratio of the focusing system by comparing the retrieved or measured effective area ${A}_{\mathrm{eff}}$ with the diffraction-limited effective area ${A}_{\mathrm{eff},\mathrm{DL}}$ , all at the focus position.

The Strehl ratio $\mathrm{SR}$ is a measure of the quality of an optical system. $\mathrm{SR}=1$ for a theoretically ideal system and $0< \mathrm{SR}<1$ for real optical systems. The Strehl ratio can be calculated as the ratio of peak fluences (or effective areas) of the real and diffraction-limited fluence distributions at the focus position:

(13) $$\begin{align}\mathrm{SR}=\frac{F_{\mathrm{peak}}}{F_{\mathrm{peak},\mathrm{DL}}}=\frac{A_{\mathrm{eff},\mathrm{DL}}}{A_{\mathrm{eff}}},\end{align}$$

where the DL subscript denotes the diffraction-limited beam^[ Reference Pirozhkov, Fukuda, Nishiuchi, Kiriyama, Sagisaka, Ogura, Mori, Kishimoto, Sakaki, Dover, Kondo, Nakanii, Huang, Kanasaki, Kondo and Kando ^{6 ]}. The diffraction-limited beam was simulated as the retrieved complex amplitude at the focusing lens with the wavefront corrected to be completely flat. We then propagated the beam back to the focus using Equation (11).

3.5 Determining the camera nonlinearity calibration function

We calibrated the near-infrared camera (Goldeye G-032 SWIR TEC) to ensure a linear response of the signal to the incident fluence. The camera nonlinearity was characterized by a calibration function that mapped the measured signal onto the incident radiant energy by applying a NoReFry algorithm.

We directed the continuous-wave alignment Nd:YLF laser onto a piece of laser alignment burn paper used as a scatter screen, and the illuminated spot was imaged onto the camera through an objective lens. We placed a near-infrared bandpass filter and neutral-density optical filters in front of the camera to bring the measured signal close to the pixel saturation limit. Subsequently, we inserted neutral-density filters into the laser beam path, measured the attenuated laser power using a power meter and acquired 20 images of the illuminated spot. The laser beam was stable while acquiring the set of images. We selected the attenuation levels to cover the entire camera dynamic range with the exposure time set to $5\;\mathrm{ms}$ . We monitored the sensor temperature stability throughout the measurement.

The NoReFry algorithm^{[ 36 ]} was developed to determine the nonlinear response of material desorption when irradiated by multiple laser pulses. The algorithm simultaneously calculates the fluence profile of the incident laser beam and the so-called response function, that is, the dependence of the local crater depth on the local accumulated radiant energy. Its complete description can be found in Ref. [Reference Vozda, Burian, Hájková, Juha, Enkisch, Faatz, Hermann, Jacyna, Jurek, Keitel, Klinger, Loch, Louis, Makhotkin, Plönjes, Saksl, Siewert, Sobierajski, Strobel, Tiedtke, Toleikis, de Vries, Zelinger and Chalupský24]. The algorithm can be applied to other nonlinear phenomena if the response function is strictly monotonically increasing with increasing accumulated energy. The camera signal response function is monotonically non-decreasing rather than strictly increasing because continuous accumulated energy values are discretized into digital signal levels. However, given the dense sampling resolution of our camera, we can approximate the response function as strictly monotonically increasing.

Before we input the measured camera data into the NoReFry algorithm, we process the acquired images. We average 20 images taken for each incident laser power and subtract a master dark image (the average of 20 dark images). We then smooth the images with a median filter $5\times 5$ . Each image is assigned radiant energy ${E}_{\mathrm{r}}= Pt$ incident on laser alignment burn paper, where $P$ is the incident laser power and $t=5\;\mathrm{ms}$ is the camera exposure time. The paired images and energies are input into the NoReFry algorithm. The resulting nonlinear response function was fitted with the following function:

(14) $$\begin{align}S\left({E}_{\mathrm{r}}\right)={aE}_{\mathrm{r}}-b\left[\exp \left(-{cE}_{\mathrm{r}}\right)-1\right],\end{align}$$

where $S$ is the measured camera signal, ${E}_{\mathrm{r}}$ is the radiant energy and $a,b,c$ are fitting parameters. The function is approximately linear with slope $a+ bc$ in the low ${E}_{\mathrm{r}}$ limit and follows an increasing asymptotic trend ${aE}_{\mathrm{r}}+b$ in the high ${E}_{\mathrm{r}}$ limit where the exponential approaches zero.

The numerical inverse ${E}_{\mathrm{r}}(S)$ of the response function $S\left({E}_{\mathrm{r}}\right)$ is proportional to the calibration function of the camera, which assigns the true signal $\tilde{S}\propto {E}_{\mathrm{r}}$ to the measured signal $S$ . It is practical to express the true signal in units of the measured signal because it makes the signals identical for a camera with an ideal linear response function. We regard the response of a real physical camera as linear at least in the low-signal regime, since pixel saturation is negligible in this region. Therefore, we introduce the true signal calibration function

(15) $$\begin{align}\tilde{S}(S)=\left(a+ bc\right)\cdot {E}_{\mathrm{r}}(S)\end{align}$$

that fulfills the condition $\tilde{S}=S$ in the low $S$ limit and is proportional to ${E}_{\mathrm{r}}$ .

We used the calibration function $\tilde{S}(S)$ to process the camera images and ensure that the signal used in PhaRe is linearly proportional to the incident fluence.