Power-law scaling of plasma pressure on laser-ablated tin microdroplets

The measurement of the propulsion of metallic microdroplets exposed to nanosecond laser pulses provides an elegant method for probing the ablation pressure in dense laser-produced plasma. We present the measurements of the propulsion velocity over three decades in the driving Nd:YAG laser pulse energy, and observe a near-perfect power law dependence. Simulations performed with the RALEF-2D radiation-hydrodynamic code are shown to be in good agreement with the power law above a specific threshold energy. The simulations highlight the importance of radiative losses which significantly modify the power of the pressure scaling. Having found a good agreement between the experiment and the simulations, we investigate the analytic origins of the obtained power law and conclude that none of the available analytic theories is directly applicable for explaining our power exponent.


I. INTRODUCTION
High-density laser-produced plasmas find many applications, ranging from inertial confinement fusion 1-3 , over the propulsion of small spacecrafts 4,5 , to sources of extreme ultraviolet (EUV) light for nanolithography [6][7][8][9][10] . The thermodynamic and radiation transport properties, particularly of high-Z laser-produced plasmas (LPPs), are extremely challenging to measure because of the transient nature of these plasmas, combined with complex equations of state and atomic plasma processes. One thermodynamic variable -the pressure -can however be elegantly obtained by measuring the propulsion velocity of metallic liquid microdroplets as a result of a laserpulse impact 11,12 . In an industrially relevant setting for EUV light production such droplets are irradiated by relatively long (∼ 10-100 ns) laser pulses at modest intensities (∼ 10 9 -10 12 W/cm 2 ), where the laser absorption takes place mostly through the inverse bremsstrahlung mechanism.
If the pulse length is large compared to the hydrodynamic time scale of the ablation flow, a quasi-stationary regime sets in, where the structure of the ablation front only slowly varies in time. The structure of such quasistationary ablation fronts has been extensively studied under various simplifying assumptions for more than 40 years 1, [13][14][15][16][17][18][19] . However, none of these theoretical works is directly applicable to our system. One of the reasons is a) Electronic mail: o.versolato@arcnl.nl the treatment of energy transport by thermal radiation. Another reason is departure from the ideal-gas equation of state (EOS) due to multiple temperature-dependent ionization of the target material. These two effects are of major importance for tin (Z = 50) targets at the here considered irradiation intensities 20 . A significant further issue is the non-trivial geometry of the laser-target configuration in our experiments, where a spherical target is irradiated from only one side and an essentially twodimensional (2D) ablation flow develops. It is likely to alter the scaling laws obtained within one-dimensional (1D) models.
Here, we present measurements of the propulsion velocity of free-falling microdroplets of liquid tin and two of its alloys over three decades in the driving Nd:YAG laser pulse energy, operating at its fundamental wavelength of 1064 nm. The propulsion velocity is obtained by means of high-resolution stroboscopic shadowgraphy techniques. Our data exhibit a remarkable, near-perfect power law dependence of the propulsion velocity on the laser pulse energy, when allowing for a certain threshold energy below which no propulsion occurs. Further, we provide results of simulations performed with the RALEF-2D 21-23 radiation-hydrodynamic code and compare these critically to the experimental data. We find a very good agreement between the simulations and the experimental power law in cases well above the threshold energy, but establish a significant disagreement regarding the threshold behavior itself.
Next, we investigate whether the obtained power law can be derived within the conventional approach based on the approximation of a steady-state planar ablation flow but corrected for the strong radiative loss. Interestingly, we conclude that none of the analytic theories available in the literature is directly applicable for explaining the power exponent observed in our experiments. We interpret this as evidence that our scaling belongs to a more complex class of scalable phenomena. Two-or three-dimensional effects, possibly combined with an essentially non-steady-state behavior, are crucial. Inevitably, the respective power-law exponents can only be calculated by solving numerically an appropriate system of partial differential equations.

A. Experimental setup
The experimental setup is described in detail in Ref. 11 and is summarized in the following. Droplets of pure liquid tin (99.995%), or one of its alloys with indium (50%) or antimony (5%), are dispensed from a piezodriven droplet generator at a repetition rate 10 kHz with a flight speed of 12 m/s in a vacuum environment ( 10 −7 mbar). The droplets relax to a spherical shape with a fixed initial diameter D 0 , which slightly varied between different experimental campaigns but stayed in the range D 0 = 2R 0 ≈ 45-47 µm, where R 0 is the droplet radius.
The produced droplets pass through the focus of an auxiliary He-Ne laser beam, whose scattered light triggers an injection-seeded Nd:YAG drive laser, operating at a 10-Hz repetition rate. The drive laser pulse, emitted at a λ = 1064 nm wavelength, is circularly polarized and has a Gaussian temporal shape with a duration t p = 10.0 ns, defined as the full-width at half-maximum (FWHM). By using an appropriate plano-concave lens, the laser beam is focused down to a circular Gaussian spot. The experiments were performed for three different focusing conditions with spot sizes of d f oc = 50, 100, and 115 µm (FWHM). Note that, due to a finite geometrical overlap, the droplets in all cases capture only a fraction of the full laser pulse energy. The pulse energy is varied over three decades, spanning the range 0.15-300 mJ as measured by using calibrated energy meters, in a manner that does not affect the transversal mode profile of the laser beam.
The position of the laser-impacted droplet is obtained from shadowgraphs generated by pulsed backlight in combination with long-distance microscopes and CCD cameras. This system provides front-view (at 30 • with respect to the drive-laser light propagation direction) and side-view (at 90 • ) images. By varying the time delay of the backlight pulse with respect to the drive laser pulse, stroboscopic images of consequent droplets are obtained (see Fig. 1). The analysis of the images is realized by a code that recognizes the center-of-pixels of the propelled and deformed droplet. Knowing the time delay between the backlight shots with a nanosecond accuracy, the droplet propulsion velocity is obtained from the slope of a linear fit to the time-dependent position of the center-of-pixels.

B. Experimental results
The measured values of the propulsion velocity U are plotted in Fig. 2 versus the energy-on-droplet E od that is defined as the fraction of the incident laser energy E given by the geometrical overlap of the spatial beam profile in focus and the droplet; in particular, for a Gaussian beam and a spherical droplet we have The thus defined energy-on-droplet appears to be a very convenient parameter, characterizing the effective portion of the laser pulse energy that gives rise to a given value of the propulsion velocity U . It also enables the comparison of the results of measurements for different focal spot sizes. As seen from Fig. 2, using this energy parametrization all data fall on a single curve. Figure 2 further demonstrates that, above a certain threshold region of E od,a ≈ 0.1 − 0.2 mJ, the dependence U (E od ) is well represented by a power law with constant values of the proportionality factor K U (m s −1 mJ −α ) and the exponent α. A fit of a power law to the full concatenated data set, using the energy range E od ≥ E od,a , yields α =0.60 (1). Fitting separately to the individual experimental data sets yields a weighted value of 0.60(1), an identical number, that is bounded by a minimum obtained value of 0.56 and a maximum of 0.63. We note that fitting only the data with a 50-µm focus size gives a slightly larger power, at 0.62(1). This value, however, is still consistent with the aforementioned result of the fit of the full concatenated data set. Similarly considering only the data from the 100-and 115-µm size focus cases, yields a power of 0.59(1), consistent with the average of 0.60(1) which is the number used in the comparisons in the following. The value obtained for K U is, in all cases, consistent with 34(3) m s −1 mJ −α , where the quoted uncertainty is the error in obtaining the absolute magnification of the imaging system. For E od < E od,a , the U (E od ) curve deviates downward from the simple power law described by Eq. (2), with a threshold at E od = E od,0 . The parameter range E od,0 < E od < E od,a corresponds to a transition regime between the onset of the ablation flow at E od = E od,0 , and the fully ablative stage at E od > E od,a . To incorporate the threshold behavior, the entire set of the experimental points in Fig. 2 is fitted by a single shifted power law, defined as The value of the offset energy E od,0 is obtained by fitting Eq. (3) to the experimental data with K U and α fixed to the values determined above, i.e. 34 m s −1 mJ −α and 0.60 respectively. The result is shown in Fig. 2 and yields a value of E od,0 = 0.04(1) mJ. Remarkably, the naive form of Eq. (3) is able to capture all the data to excellent accuracy.
These values are consistent with, and in fact nearly identical to, the values found in our previous work (α = 0.59(3), K U = 35(5) m s −1 mJ −α , E od,0 = 0.05(1) mJ), dealing with a much smaller data set for solely indiumtin droplets 11 . Consequently, the here demonstrated excellent reproducibility of the data strongly improves the statistical significance of our findings and the broad applicability of the power law. It presents a solid basis for drawing conclusions on the underlying physics.
As is explained in more detail in Section III, the energy E od,a marks the lower boundary of a distinct pattern of laser ablation. In such conditions, the hot plasma with T 5-10 eV envelopes the entire front-illuminated (laser-facing) hemisphere of the droplet, the velocity field across the laser absorption zone approaches that of a quasi-spherical flow, and all the laser flux contributing to E od is efficiently absorbed in the ablated plasma cloud by the inverse bremsstrahlung mechanism. Accordingly, we designate the regime above E od,a as the fully ablative regime. In this regime the peak laser intensity on target spans the range 10 9 W/cm 2 < I l < 3 × 10 11 W/cm 2 .

A. RALEF-2D code
The simulations reported in this work have been performed with the two-dimensional (2D) radiationhydrodynamics code RALEF 21,22 , which has lately been extensively used to simulate laser-driven, droplet-based EUV sources for nanolithography applications 20,23,24 . The hydrodynamics module of RALEF is based on the upgraded version of the CAVEAT package 25 , where the second-order Godunov-type algorithm on an adaptive quadrilateral grid is used. The thermal conduction and the spectral radiation transfer (in the quasi-static approximation) are treated within a unified symmetric semi-implicit scheme 21,26 with respect to time discretization. To describe the spatial dependence of the spectral radiation intensity, the classical S n method is used, combined with the method of short characteristics 27 to integrate the radiative transfer equation.
The equation of state (EOS) of tin is constructed by using the FEOS model 28 that provides, within a unified model, an adequate and thermodynamically consistent description of high-temperature plasma states together with the low-temperature liquid-gas phase coexistence region. The model for thermal conductivity is based on a semi-empirical expression for the transport cross section of the electron-ion collisions 29 , which enables a smooth matching of the Spitzer plasma conductivity to that of metals near normal conditions. All the simulations are performed for a spherical droplet of pure tin with initial radius R 0 = 25 µm and initial density ρ 0 = 6.9 g/cm 3 , assuming that slight differences between the physical properties of pure Sn and its two alloys used in the experiments, are insignificant. The adaptive numerical mesh has a topological structure as displayed in Fig. 3. It extends with 360 zones over the π interval of the polar angle θ, and with 350 radial zones over the interval 20 µm ≤ r ≤ 1 mm. This totals to 142 200 mesh cells over the simulated half-circle in the rz plane. The mesh is progressively refined in the radial direction towards the droplet surface to resolve the skin layer of the liquid tin. The minimum cell thickness of this layer is 4.5 nm. The outer region 25 µm ≤ r ≤ 1 mm is initially filled with a tenuous tin vapor at a density of ρ v0 = 10 −10 g/cm 3 .
In all the simulation runs, the same Gaussian temporal power profile of the 1064 nm laser pulses is used, with the pulse duration t p = 10 ns (FWHM), peaking at t = 1.5t p = 15 ns. The spatial laser profile is also Gaussian, with two values of the focal spot diameter (FWHM): d f oc = 115 µm (series A) and d f oc = 50 µm (series B). Propagation and absorption of the laser light is calculated within a hybrid model 30 , which accounts for refraction in the tenuous corona. In addition, it ensures a physically correct description of reflection from the critical surface, including the Fresnel reflection from the metal-vacuum interface. Lastly, the incident light is assumed to be unpolarized.
For all cases in the fully ablative regime, radiative energy transport is important. Radiation generation and transport is treated with the same opacity model as in Ref. 24, where the conversion efficiency into the 13.5nm EUV emission is investigated for a CO 2 -laser-driven plasma. The angular dependence of the radiation inten-sity is modeled with the S 6 quadrature, while the spectral dependence is simulated with 28 discrete spectral groups of variable width. Two spectral groups belong to the 2% band at 13.5 nm, where the strongest emission from the Sn plasma is expected at sufficiently high laser intensities.
B. Simulation results

Droplet propulsion
The calculated propulsion velocity U for various E od values is plotted in Fig. 4. In the RALEF code it is computed as the velocity of the center of mass, comprising all the material with the density in excess of 1 % of its maximum value at the time t = t f = 200 ns. Similarly to the experimental results, for E od > 0.1 − 0.2 mJ, the dependence U (E od ) is almost a perfect power law: the deviations of the calculated points from Eq. (2) with the best-fit values of calculated for the combined set of points from series A and B in the range E od ≥ 0.2 mJ, do not exceed ±2.5% -which is practically the intrinsic accuracy of the simulations. Fig. 4 confirms that within the same ±2.5% accuracy the energy-on-droplet E od proves indeed to be an adequate universal parameter, which unites the d f oc = 115 µm and d f oc = 50 µm points into virtually a single curve. For the variation of the coefficient K U with the droplet size R 0 and the laser pulse duration t p , we refer to the Appendix. Judging from Fig. 4, the agreement between the calculated and the measured U values Dependence of the propulsion velocity U on E od calculated with the RALEF-2D code. The focus diameter d f oc (µm) and the droplet diameter D0 (µm) for different simulation series are indicated in the legend as d f oc /D0.The black curve represents the best fit to the experimental points (see Fig. 2). The vertical line at E od = 0.04 mJ corresponds to the threshold for droplet propulsion as inferred from that fit.
in the fully ablative regime could hardly be better: the deviations from the best experimental fit do not exceed 11%, which lies within the experimental errors. However, the droplet diameter D 0 = 50 µm, used in the simulations, slightly exceeds the actual values of D 0 ≈ 45-47 µm. For instance, the correction to a smaller value D 0 = 46 µm would raise the calculated U values in the fully ablative regime in Fig. 4 by some 20%, leaving the power α unchanged. The fact that the model tends to slightly overestimate the propulsion velocity can, on the one hand, be attributed to a systematic experimental uncertainty, combining possible measurement errors in the spatial beam profile and the droplet diameter. Alternatively, the RALEF simulations may, for example, systematically underestimate the radiation energy losses, whose modeling could still noticeably be improved.
All in all, a very good agreement between the simulation and the experiment is found in the fully ablative regime. Particularly, concerning the scaling exponent α, the best-fit experimental value α = 0.60(1) is practically the same as the theoretical value in given Eq. (4). This provides a strong evidence that the RALEF code sufficiently accurately accounts for the key physical processes governing the Sn plasma dynamics in this regime. Therefore, it can be used to extract additional information about the relative role of these processes.
At the low energies E od < 0.1 mJ, the simulation results begin to significantly deviate from the experimental values. Here we have to deal with the initial phase of the onset of ablation, which is controlled by physical processes that are quite distinct from those governing the fully ablative regime. The key role in this initial phase should belong to an adequate modelling of laser-optical properties and propagation of a non-steady thermal wave across a thin surface layer of tin. In such conditions, this layer is driven into a non-trivial thermodynamic state of superheated metastable liquid, followed by a phase transition into a state of dense hot vapor. We leave the full investigation of this regime for future work.

Plasma characterization in the fully ablative regime
A general perception of the plasma dynamics in the fully ablative regime can be obtained from Fig. 5, which displays the 2D density and temperature distributions for the two cases of E od = 0.2 mJ and 30 mJ at time t = 15 ns, coinciding with peak laser power. As is seen in Figs. 5(b) and (c), a characteristic feature of the fully ablative regime is a stabilized geometry of the plasma flow across the laser absorption zone. The latter manifests itself in Figs. 5(e) and (f) as the region with highest plasma temperatures. Note that the peak temperature in the ablative regime varies with E od over a wide range of 5 eV T 100 eV. In all cases with E od ≥ E od,a , by the middle of the pulse, the plasma plume attains a size of several R 0 and occupies the entire 2π of the solid angle above the illuminated droplet hemisphere; the ve- locity field stabilizes to a quasi-steady, quasi-spherically diverging pattern; the laser-absorption zone itself reaches its maximum size, which becomes practically independent of E od .
Intuitively it is clear that, once the 2D (or 3D) geometry of the plasma flow and laser absorption settles down to a stable pattern, the principal ablation parameters (like the characteristic pressure, temperature, ablation velocity, etc.) can be expected to become scalable. On the other hand, in the low-energy cases with E od < E od,a (see Figs. 5(a) and (d)), intense laser absorption takes place in a narrow plasma plume near the target pole while a large portion of the incident flux contributing to E od is reflected from a cooler and sharper liquid-vapor boundary at θ 40 • -50 • . Therefore, the ablation parameters from these low-energy cases cannot be expected to be scalable in the same way as those in the fully ablative regime. Fig. 5 also demonstrates that the ablation flow is subject to hydrodynamic instabilities, the most salient of which appears to be the self-focusing instability due to laser refraction in the underdense plasma, inherent in the laser deposition model 30 . The resulting irregular fluctuations of the plasma parameters in space and time manifest themselves as "spotty" temperature distributions and "wavy" n e isocontours in Fig. 5. The temporal variation of the ablation pressure at a fixed location, illustrated in Fig. 6 for the target pole, becomes especially violent for low E od values.
Although the self-focusing instability has a clear physical origin, the amplitude of the ensuing fluctuations tends to be overestimated in the present RALEF simulations (especially on length scales comparable to, or smaller than the laser wavelength λ) due to the absence of diffraction effects in the laser propagation model 30 . However, when averaged over space and time, the impact of this "noise" on the calculated U values turns out to be negligible, i.e. on the level of ±1% as ascertained by dedicated computer runs. Having verified it in 2D, we expect no more than only a moderate, by about a factor of 1.5, increase of this effect in the full 3D approach. This is similar to what has firmly been established for the nonlinear stage of the Rayleigh-Taylor instability 31 .

Ablation pressure
The ablation-plasma parameter most directly related to the propulsion velocity U is the ablation pressure. More specifically, the velocity U can be determined from the relationship where M is the total mass, and P is the total momentum of the liquid tin at a certain moment t f t p . As the entire simulated configuration is axisymmetric, the total momentum vector P lies along the z axis. In our case, the results become insensitive to t f for t f 100 ns, thus we present results for t f = 200 ns. From the simulations we learn that the ablated mass fraction δ M , defined as the relative fraction of the total tin mass with ρ < 0.1 g/cm 3 , does not exceed 10% for the entire range of E od ≤ 40 mJ (see Tables I and II). The subsequent deformation of the ablated surface is not significant (see Fig. 5). Then, the propulsion momentum P can be evaluated as where p a (t, θ) is the ablation pressure at the spherical droplet surface as a function of time t and polar angle θ, and j p (θ) is the local impulse of the ablation pressure. Note that θ is measured with respect to the negative direction of the rotation axis z, as is shown in Fig. 3.
Equations (5) and (6) can be used to relate the established scaling of U with E od in Fig. 2 to existing analytic scaling laws for the ablation pressure p a . However, all the previous analytic results on the scaling of p a with the incident laser flux I l have been obtained under a few assumptions. It is assumed that the ablation flow either (i) has a 1D planar geometry (p a is constant in space), or (ii) is in a steady state (p a is independent of time), or both 18 . Unfortunately, neither of these assumptions can be considered as adequate for our situation. Nonetheless, the effects of the spatial, along the droplet surface, and the temporal variations of the ablation pressure p a (t, θ) can be separated as follows. One can rewrite Eq. (6) as (7) where j pθ = 2 π 0j p (θ) sin θ cos θ dθ,j p (θ) ≡ j p (θ)/j p0 . (8) Our simulations demonstrate that in the fully ablative regime the dimensionless spatial form-factor j pθ of the pressure impulse barely depends on the incident laser flux when the focal spot is fixed (see Tables I and II). For d f oc = 115 µm, for instance, it fluctuates in the range j pθ ≈ 0.57-0.59, remaining virtually constant within our simulation accuracy. Hence, as long as we can neglect small variations of the mass M and size R 0 of the irradiated droplet, the problem of the analytic derivation of the scaling of U with E od is reduced to the derivation of the analogous scaling for the local (at the pole) pressure impulse j p0 . Before tackling this issue, we provide some additional information on the angular dependence of the ablation pressure that might be helpful for a general analysis of the hydrodynamic response of liquid droplets to laser pulses 11,32 . Fig. 7 shows several angular profiles of the normalized pressure impulsej p (θ), calculated with the RALEF code. Despite the fact that thej p (θ) curve for the highestenergy case E od = 30 mJ is clearly broader than those for lower pulse energies, its integral (see Eq. (8)) remains practically the same because of the negative contribution from the backward hemisphere θ > 90 • . A salient local rise ofj p (θ) at θ 150 • for the 2-mJ case is explained by the plasma flowing around the droplet and accumulating on its horizontal axis. It leaves a local cloud of relatively dense and hot vapor, which exerts a noticeable backward pressure onto the droplet for some 30-50 ns after the laser has already been off. We further note that, for the same E od = 2 mJ, a tighter laser focus (the d f oc = 50 µm curves) produces an only slightly narrower pressure profilej p (θ).

IV. ANALYTIC SCALING LAWS
Having found an excellent agreement between the experiment and simulations, we will attempt to derive the obtained scaling law analytically on the basis of an appropriately simplified model. Additional information, available from the simulations, provides guidance for working out such a model.
Analytic scaling laws are usually derived for the ablation pressure p a as a function of the hydrodynamically absorbed flux I lh (W/cm 2 ), assumed to be constant in time and fully converted into the kinetic and internal energy of the ablated material 18 . To simplify the argumentation, we focus our attention on the simulations (series A) with a fixed spot size d f oc = 115 µm. Then, because all the pulses have the same temporal profile, the polar incident flux I l,0 (t), the incident laser energy E, and the energy-on-droplet E od are all directly proportional to one another, as well as to the polar energy fluence F l,0 = I l,0 (t) dt. Consequently, an approximate analytic scaling of U with E od could be obtained by (i) relating the incident laser fluence F l,0 to the hydrodynamically absorbed one F lh,0 , and (ii) making an assumption that the time-integrated quantities j p0 and F lh,0 = I lh,0 (t) dt scale with one another in the same way as p a and I lh in a steady-state planar 1D ablation front, for which analytic results are available. Here we assume that the droplet mass M and the 2D form-factor j pθ in Eqs. (5) and (7) are constant. Note that assumption (ii) is by no means obvious, and might, in fact, be rather inaccurate.

A. Laser absorption and radiative losses
There are two main loss mechanisms that reduce the incident laser energy fluence F l,0 to the hydrodynamically absorbed one F lh,0 , namely, partial reflection of the laser light and radiative losses. Accordingly, since F l,0 is directly proportional to E od , we can, following our logic, introduce a hydrodynamically absorbed energy-on-droplet In Eq. (9) f la is the laser energy absorption fraction, and φ r is the fraction of the absorbed laser energy which escapes from the plasma by thermal emission. Having introduced effective corrections for the laser reflection and radiative losses by means of Eq. (9), we take the next step and relate the resulting scaling of j p0 with E od,h to an analytic scaling of p a with I lh predicted by an appropriate 1D model. If a close agreement were found, we could accept the invoked 1D model as an appropriate one for the interpretation of our experiments.
Strictly speaking, both factors f la and (1 − φ r ) in Eq. (9) must be calculated at the target pole. But even a simplest analytic model for evaluating f la and φ r would be too cumbersome for the present work 20 . Instead, we take their values from the RALEF simulations. The problem, however, is that the local polar value of φ r cannot be extracted from the simulations. Moreover, it is an ill-defined quantity because of the non-local nature of radiation transport. Thus, we are forced to use the integral values of φ r , calculated for the whole plasma volume and listed in Tables I and II. For the laser absorption, whose impact on the scaling is considerably less important (∆α ≈ 0.03), we also use the integral values of f la = f la,od , calculated for the laser energy fluence over the cross-section πR 2 0 of the droplet. These values are consistent with the integral values of φ r and exhibit weaker instability variations than the local polar values f la,0 .
First of all we note that the calculated values of φ r , ranging from 20% to 70% as E od increases from 0.2 mJ to 40 mJ, provide clear evidence of the important role played by radiative losses in our situation. For the scaling exponent it is important that the coefficient (1 − φ r ) changes by about a factor of 2.5-3 over the considered range of E od , which implies an exponent shift by ∆α ≈ 0.17. Fig. 8 shows the dependence of the calculated pressure impulse j p0 on the incident, E od , and hydrodynamically absorbed, E od,h , energy-on-droplet. Solid lines represent the respective power-law fits, that yield the following exponents.
The results of the fits significantly differ from one another. This difference of ∆α ≈ 0.14 provides a quantitative measure of the influence of radiative losses on the discussed scaling law. In fact, this influence is even stronger (∆α ≈ 0.17) since the two factors f la and (1 − φ r ) in Eq. (9) change in opposite directions (see Tables I and  II). Clearly, it is the second exponent α = 0.724 (14) that should be compared with the known analytic scalings for p a (I lh ). A noticeably larger statistical uncertainty in this exponent (±0.014 versus ±0.005, thus comparable to the experimental error), related to the goodness of fit, is apparently caused by using the integral values of φ r and f la , which "feel" the 2D ablation geometry of a spherical droplet. Note that the exponent α = 0.583(5) for the j p0 (E od ) dependence differs slightly from the previously quoted value of α = 0.610(5) for the U (E od ) scaling (see Section III B 1). This difference of ∆α ≈ 0.03 arises from the fact that the remaining liquid mass M in Eq. (5) decreases by about 9% as E od increases from 0.2 mJ to 30 mJ, and less impulse is needed to attain a given velocity U .

B. Effects of the equation of state
Well-known theoretical models of 1D quasi-stationary ablation fronts, based on the ideal-gas equation of state (EOS) with the adiabatic index γ = 5/3, yield two limiting scaling laws for the ablation pressure. Namely, the one for the case where laser absorption occurs in an infinitely thin layer at the critical surface 13,16,17 (case I), and the other one for the case where laser light is absorbed in an extended region by the inverse bremsstrahlung mechanism before reaching the critical surface 17,33,34 (case II), case I (ideal-gas EOS), lh L −1/9 , case II (ideal-gas EOS). (11) In case II, an additional relevant parameter enters the scaling, which is the density-gradient length L in the absorption zone. For quasi-spherical (or cylindrical) diverging flows, where a steady-state solution with a sonic point exists, L should be set equal to the radius of the sonic point 33 . In the planar geometry, where no steady-state solution is possible 33 , one can assume the laser to be absorbed in a non-steady rarefaction wave in expanding plasma, where L ∝ c s t, and c s is the characteristic sound velocity. In this way one arrives at yet another well-known analytic scaling p a ∝ I 3/4 lh t −1/8 , applicable to non-steady planar ablation flows with the ideal-gas EOS 13,17,33,35 .
All the above analytic scalings with rational-number exponents, based on the ideal-gas EOS, can definitely be applied to interpretation of experiments on low-Z targets (like plastic foils) that are fully ionized by a sufficiently high laser energy flux. None of them, however, can be employed in our case, where a temperature-dependent ionization of tin (Z = 50) changes the appropriate planar analytic scalings in Eq. (11) The experimental situation analyzed here lies between these two cases but closer to case II. We compare the exponent α = 0.724 (14) in Eq. (10) with 0.56 α 0.64 in Eq. (12). The effect of variation of the density-gradient scale L with the laser intensity I lh for case II is small and only enhances the discrepancy because L can only grow with I lh . From comparison between Figs. 5(b) and (c) one infers that the radius of the absorption zone increases by no more than a factor of 1.7 as E od increases from 0.2 mJ to 30 mJ, implying an effective reduction of the scaling exponent by ∆α ≈ −0.02. Thus, a good agreement with the appropriate analytical scaling could have been claimed if Fig. 8 demonstrated j p0 ∝ E α od,h with 0.56 α 0.62 -which is obviously not the case. A superficial observation that the scaling (10) of j p0 with E od,h is very close to the theoretical result p a ∝ I 3/4 lh (with t ≈ t p being fixed) should be qualified as incidental. Summarizing, we conclude that the scaling (2), (4) of the propulsion velocity U with the energy-on-droplet E od , established in this work, cannot be derived from the previously published 1D analytic models of the laser ablation fronts.

V. CONCLUSION
Having performed an extensive series of experiments with Nd:YAG laser pulses at different focusing conditions, we have found that within a certain range of laserpulse energies, covering more than three decades in magnitude, the propulsion velocity of tin droplets scales as a power law U ∝ E α od of the energy-on-droplet E od (the incident laser energy intercepted by the cross-section of the droplet). The theoretical analysis, based on 2D simulations with the radiation-hydrodynamics RALEF code, has revealed that the scalability range corresponds to a fully developed regime of laser ablation, where the zone of laser absorption (by inverse bremsstrahlung) in the ablated plasma settles to a stable configuration. For droplets with radii R 0 ≈ 25 µm it starts at E od 0.1-0.2 mJ. The scaling exponent α = 0.610(5), obtained from the RALEF results, agrees perfectly with the experimental value of α = 0.60(1). The performed analysis demonstrates how the propulsion of metallic microdroplets by a laser-pulse impact can be a good probe for the plasma ablation pressure.
It should be noted that our study was done under a rather unique combination of conditions. A spherical target composed of a high-Z material was irradiated from one side and propelled by an essentially 2D ablation flow. Since the vast majority of previous measurements of the laser ablation pressure were done on low-Z planar targets or on pellets with spherically symmetric irradiation geometry (see, e.g., [36][37][38][39][40]), we chose to avoid direct comparison of our results to those obtained in these other works, as spurious coincidence of two numbers from different experiments could obfuscate the underlying physics. Instead, we focused our efforts on analyzing the main physical effects that determine our scaling power.
A thorough examination, facilitated by additional information from the RALEF simulations, of the physical processes governing the fully ablative regime in our series of experiments has revealed that the scaling law cannot be directly derived from any of the existing analytic models of quasi-steady 1D ablation fronts. Moreover, this cannot be done even after the effects of radiation energy losses and realistic EOS of tin have been accounted for. The cause must be a complex, essentially 2D (or even 3D) structure of the ablation plasma flow, where the non-local energy transport by thermal radiation in both lateral and radial directions plays an important role. An additional complication comes from the finite pulse length t p = 10 ns. It is difficult to justify the steady-state approximation, usually implied by analytic evaluation of the scaling exponent, when t p remains fixed. While the timescale of flow relaxation 20 to a quasi-steady state is comparable with t p at E od = 0.2 mJ, it decreases by about a factor of 3-4 at the upper end E od = 30-50 mJ of the explored range.
In conclusion, the established scaling of the plasmapropulsion velocity U of tin microdroplets with laser en-ergy E od belongs to a class of scaling laws where theoretical evaluation of the scaling exponent requires the numerical solution of partial differential equations that capture the relevant physical effects in two or three dimensions.

ACKNOWLEDGMENTS
Part of this work has been carried out at the Advanced Research Center for Nanolithography (ARCNL), a public-private partnership between the University of Amsterdam (UvA), the Vrije Universiteit Amsterdam (VU), the Netherlands Organization for Scientific Research (NWO) and the semiconductor equipment manufacturer ASML. We acknowledge the assistance of the mechanical workshop and the design, electronics, and software departments of AMOLF as well as the direct technical support at ARCNL. We also thank SURFsara (www.surfsara.nl) for the support in using the Lisa Compute Cluster. This work was partially (the contributions by M. M. Basko and D. A. Kim) funded by the Russian Science Foundation through grant No. 14-11-00699-Π.
Appendix: Dependence of the propulsion velocity on the droplet size and laser pulse duration Having established the scaling Eqs. (2) and (4) of the propulsion velocity U with the energy-on-droplet E od , one can, following the logic of Section III B 3 and making some reasonable assumptions, evaluate the dependence of U on the droplet radius R 0 and the laser pulse duration t p . This might be useful for practical applications.
First of all, we suppose that the exponent α in Eq. (2) does not vary with R 0 and t p , and only the dimensional coefficient K U changes. If, when varying R 0 , we keep the values of the polar energy fluence F l,0 = I l,0 (t) dt and of the ratio R 0 /d f oc fixed, both the polar pressure impulse j p0 and the form-factor j pθ should remain practically unchanged. Then, having noted that in Eq. (5) M ∝ R 3 0 and, as it follows from Eq. (7), P ∝ R 2 0 , we obtain U = K U E α od ∝ R −1 0 . Finally, because for fixed F l,0 and R 0 /d f oc one has E od ∝ R 2 0 , we arrive at Similarly, we can deduce the scaling with the pulse duration t p by assuming that the Gaussian pulse profile is simply stretched in time by a factor a (t p → at p ), with the peak laser intensity kept fixed. Then, because the local (polar) ablation pressure p a (t, 0) depends primarily on the local laser intensity, one can surmise that the corresponding pressure pulse will also be simply stretched in time by the same factor a. As a result, the propulsion velocity would scale as U → aU . Since E od in Eq. (2) is directly proportional to t p , the factor K U should scale as