High-energy ions from Nd:YAG laser ablation of tin microdroplets: comparison between experiment and a single-fluid hydrodynamic model

We present the results of a joint experimental and theoretical study of plasma expansion arising from Nd:YAG laser ablation (laser wavelength λ = 1.064 μm) of tin microdroplets in the context of extreme ultraviolet lithography. Measurements of the ion energy distribution reveal a near-plateau in the distribution for kinetic energies in the range 0.03–1 keV and a peak near 2 keV followed by a sharp fall-off in the distribution for energies above 2 keV. Charge-state resolved measurements attribute this peak to the existence of peaks centered near 2 keV in the Sn3+–Sn8+ ion energy distributions. To better understand the physical processes governing the shape of the ion energy distribution, we have modelled the laser-droplet interaction and subsequent plasma expansion using two-dimensional radiation hydrodynamic simulations. We find excellent agreement between the simulated ion energy distribution and the measurements both in terms of the shape of the distribution and the absolute number of detected ions. We attribute a peak in the distribution near 2 keV to a quasi-spherical expanding shell formed at early times in the expansion.


Introduction
Laser-produced plasmas (LPPs) formed on tin microdroplets are now established as the light source of choice in newgeneration lithography machines for high-volume manufacturing of integrated circuits below the 10 nm node [1][2][3]. Their incorporation in modern-day lithography machines relies on their ability to provide sufficiently high fluxes of shortwavelength radiation to enable the patterning of nanometrescale features on integrated circuits.
The superposition of millions of lines arising from transitions between complex configurations in n = 4 shell Sn 11+ -Sn 15+ ions are the atomic origins of this light [9][10][11][12][13]. Importantly, this feature overlaps with the 2% reflective bandwidth (13.5 ± 0.135 nm-the so-called 'in-band' region) of molybdenum/silicon multilayer mirrors (MLMs) [14]. Such mirrors are an integral component of EUV lithography tools, transporting EUV photons from the light source to their final destination at the wafer stage.
One key aspect of industrial EUV light source development has focussed on optimising the photon output of LPP EUV light sources. To-date, efforts have concentrated on increasing (i) EUV power and (ii) the so-called conversion efficiency (CE-the ratio of in-band EUV energy emitted into the half sphere back towards the laser to input laser energy) of the light source [1,7]. To meet the high power levels required for industrial applications, a dual-pulse irradiation scheme is employed [15]. First, a low-intensity prepulse is used to deform the droplet into an elongated disk-like target. This target is then irradiated by a second, high-energy CO 2 laser pulse (λ = 10.6 μm) which generates a hot, EUV-emitting plasma. This EUV light is then focussed by an MLM (known as the collector mirror) to an exit port of the light source vessel whereupon it enters the scanner tool for use in the lithographic process.
A second and no-less crucial aspect of EUV light source development has focussed on the design and implementation of so-called debris mitigation schemes. In the context of the current application, plasma expansion will lead to the bombardment by tin ions on the plasma-facing collector mirror. The combined effects of sputtering and ion implantation will, over time, degrade the performance of the collector mirror and reduce EUV throughput. In an industrial setting, the light source vessel is typically filled with a background hydrogen gas to stop energetic ions from reaching the collector mirror [1,16,17]. One can also introduce a strong magnetic field in the region surrounding the droplet to deflect plasma ions away from the collector mirror [18][19][20]. A comprehensive understanding of the characteristics of the plasma expansion (distribution of ions over kinetic energy, angular distribution of ions, etc.) can greatly assist the design of effective debris mitigation schemes.
A number of studies examining tin plasma expansion have been performed over the past two decades. In the mid-2000's, Murakami et al [21] and Fujioka et al [22] demonstrated that ion energy distributions recorded from minimum-mass plasmas driven by 10 ns-long Nd:YAG (λ = 1.064 μm) pulses can be described using a model of isothermal plasma expansion. In this model, the distribution of the number of ions N as a function of kinetic energy E is written where E 0 = 2Zk B T e ln[R(t)/R 0 ] is a characteristic energy scale dependent on the charge state Z, electron temperature T e , initial plasma size R 0 and a time-dependent characteristic system size R(t) (k B denotes the Boltzmann constant).
In the above equation α is the dimensionality of the expansion (α = 1, 2 and 3 correspond to planar, cylindrical and spherical geometries, respectively), Γ is the gamma function and N 0 = ( √ πR 0 ) α n i00 is the total number of ions and n i00 is the initial ion number density at the origin of the model, i.e., n i00 = n i (r = 0, t = 0). Experiments performed by Bayerle et al [23] also qualitatively demonstrate the applicability of Murakami's model (equation (1)) to Nd:YAG-irradiated (6 ns pulse duration) planar tin targets.
Plasmas driven by shorter, ps-duration pulses [23,24] exhibit ion energy distributions whose shapes are better described by the planar isothermal expansion model of Mora [25]. In this case the ion energy distribution reads where N * 0 is the total number of ions and the characteristic energy scale E * 0 = Zk B T e . The main difference between the two models lies in the choice of the plasma density profile. The model of Murakami et al [21] employs a Gaussian profile for the plasma density ρ ∝ exp(−[r/R(t)] 2 ) whereas an exponential density profile is used in the work of Mora [25]. As noted by Murakami et al [21], this latter form for the density profile is better suited for plasmas generated by short pulse lasers or those formed on thick targets.
Other work on the topic of tin plasma expansion has explored, for example, the role of laser pulse duration and laser wavelength on the ion energy distribution [26][27][28], the angular distribution of ion kinetic energies [29][30][31], the suppression of fast ions using a low-energy prepulse [32,33], and the role of electron-ion recombination during the expanding phase of the plasma [34]. It is important to note that the vast majority of these studies have investigated plasma expansion from laser-irradiated planar tin targets rather than from industriallyrelevant droplet targets. Much work still remains to be done on this latter topic.
The goal of the present study is to investigate plasma expansion in the form of emission of energetic charged particles from Nd:YAG-irradiated tin microdroplet targets. This study serves to complement recent work on photon emission from solid-state laser-driven EUV light source plasmas [7,35,36]. In contrast to the current industrial implementation, solid-state laser-driven plasmas may not require the use of a pre-pulse for efficient EUV production [7]. As such, they are a promising candidate for future laser-driven EUV light source plasmas. First we present measurements of the ion energy distributions using an electrostatic analyser (ESA). These measurements reveal the existence of peaks near 2 keV in the high-energy tails of the Sn 3+ -Sn 8+ ion energy distributions. These features combine to yield a peak near 2 keV in the charge-state summed ion energy distribution. To elucidate the origin of this peak, we have performed two-dimensional radiation-hydrodynamic simulations of the plasma formation and its subsequent expansion using the radiative arbitrary Lagrange-Eulerian fluid dynamics in two dimensions (RALEF-2D) code. The ion energy distribution obtained from the simulations compares favourably to the measurements both in terms of the shape of the distribution and the absolute number of detected ions. We attribute the peak in the ion energy distribution to a high-velocity, quasi-spherical expanding shell formed at early times in the plasma expansion. The current work advances on the work presented in references [8,23] to provide a quantitative understanding of absolutelycalibrated measurements via radiation-hydrodynamic modelling of the expanding plasma, beyond the aforementioned idealized plasma expansion models.
The layout of this paper is as follows: in section 2 we discuss the experimental setup and provide details of the ion energy distribution measurements. This is followed by a description of the single-fluid, single-temperature model implemented in the RALEF-2D code and a brief discussion of the simulation parameters. In section 4 we discuss the results of the simulations, focussing on the temporal and spatial evolution of the speed and ion number density profiles in the expansion. In section 5 we compare the ion energy distribution obtained from the simulations with our experimental measurements. Comparisons are drawn with the predictions of wellknown analytical models of plasma expansion into vacuum. Finally, we summarise this work in section 6.

Experimental setup, method and results
In the experiments, tin droplets were dispensed from a droplet generator mounted at the top of a vacuum chamber (backing pressure ∼10 −7 mbar). The droplet generator consists of a heated (260 • C) molten tin reservoir connected to a nozzle. The diameter of the tin droplets used in the experiment was 28 μm. Upon crossing the centre of the chamber, the droplets pass through a light sheet created by a He:Ne laser. Light scattered by the droplets was detected by a photomultiplier tube which triggered the plasma-generating laser pulse and the acquisition apparatus. Plasmas were generated by focussing the output of a commercial Nd:YAG laser system onto the tin droplets. The laser pulses exhibited Gaussian-like temporal and (focussed) spatial laser profiles. The temporal full-width at half-maximum (FWHM) was 10 ns and the FWHM of the focussed pulses was approximately 60 μm. Employing a laser pulse energy of 60 mJ resulted in a laser power density on the targets of I L ≈ 2 × 10 11 W cm −2 . This choice of power density is known to yield optimum CE's for Nd:YAG-driven tin plasmas [7]. We note that this particular choice of laser parameters, combined with the given droplet diameter, will not lead to the full ablation of the tin droplet.
Charge-state resolved ion energy distributions were measured using an ESA. The opening aperture of this device was located 1.12 m away from the droplet targets and was positioned at 60 • with respect to the incident laser axis. The ESA consists of a radial electric field deflection region followed by a calibrated channeltron detector. The radial electric field between the two electrodes of the ESA selects charge states based on the ratio of their kinetic energy E to charge state Z according to E/Z = 5 × U ESA where U ESA is the voltage across the ESA electrodes (measured in volts). A time-offlight (ToF) analysis is used to obtain charge-state resolved ion counts for a given E/Z. By scanning the ESA voltage U ESA over a desired range, one can obtain charge-state resolved ion energy distributions. The ESA-ToF measurements have been benchmarked against charge-state integrated measurements made using a Faraday cup (FC) which was positioned at an angle of −60 • with respect to the incident laser axis.
The total ion energy distribution d 2 N/dE dΩ was derived from the measurements via where ΔE (E) is the energy-dependent absolute energy resolution, N m Z (E) is the number of tin ions of charge state Z having kinetic energy in the range E ± (ΔE(E)/2) detected by the channeltron, η det is the detection efficiency of the channeltron, η CR ( 1) is a correction for undetected counts due to high count-rate effects and ΔΩ is the solid angle of the input aperture of the ESA device. The absolute energy resolution scales linearly with E according to ΔE/E 10 −2 . Further details regarding the ESA calibration can be found in [37].
In figure 1 we present the results of our charge-state resolved ion energy distribution measurements. Examining this figure, we first note that the distributions associated with Sn 1+ and Sn 2+ ions are rather broad, spanning energies in the range 0.03-2 keV. The kinetic energy for which the ion energy distribution peaks, E peak , clearly increases with increasing charge state. Both distributions also exhibit a nearexponential fall-off for E > E peak . While the aforementioned trend in E peak continues for Sn 3+ and Sn 4+ , we note the emergence of a second, high-energy peak located just below 2 keV (this peak is also present in the Sn 2+ distribution although it is less pronounced than in the Sn 3+ and Sn 4+ distributions). With increasing charge state this peak grows in intensity until E peak ≈ 2 keV in the Sn 5+ and Sn 6+ ion energy distributions (we also make a tentative observation of two peaked features in the Sn 7+ ion energy distribution near 1.6 and 2.5 keV, respectively). While we do detect Sn 8+ ions in the experiments, the amplitude of the ion energy distribution is an order-of-magnitude lower than the Sn 7+ distribution. No traces of higher charge states could be reliably detected. Importantly, the kinetic energy associated with this high-energy feature is independent of charge state. Shown in red in figure 1 is the total ion energy distribution obtained by summing the individual Sn 1+ -Sn 8+ ion energy distributions. Three distinct regions emerge: (i) a near-plateau in the ion energy distribution between 0.03-1 keV (ii) a peak near 2 keV followed by (iii) a sharp fall-off for E > 2 keV. Finally, we note that the EUV-generating tin charge states Sn 11+ -Sn 15+ , whilst generated in the hot, dense region of the plasma, are not detected in the measurements. The absence of Sn 11+ -Sn 15+ charge states may in part be attributed to the process of recombination, whereby free electrons in the expanding plasma recombine with these ions through processes such as three-body or radiative recombination [34].

Radiation hydrodynamic simulations and the RALEF-2D code
In order to elucidate the dynamics of the plasma expansion and its influence on the ion energy distribution, we have undertaken numerical modelling of the plasma formation, growth and subsequent expansion using radiation-hydrodynamic simulations. In the following, we discuss the underlying assumptions of the single-fluid, single-temperature approach adopted in the present work and we provide details of the simulations we have performed with the RALEF-2D code.

Single-fluid single-temperature radiation hydrodynamics
We have chosen to model the plasma formation and its subsequent expansion using a single-fluid, single-temperature hydrodynamic model including the effects of radiation transport and thermal conduction. In this approach, the free electrons and ions are treated as a single fluid having a single temperature T e = T ion = T. Although more complex approaches such as the two-fluid, two-temperature [38] or single-fluid, two-temperature models [39][40][41][42] have been pursued, the single-fluid, single-temperature description should be adequate for the current purposes. For one, simulations of Nd:YAG-irradiated lithium, plastic and gold targets performed by Sunahara and Tanaka [39] indicate a rather small difference (less than 20%) between T e and T ion in the plasma. This behaviour has also been observed in simulations of laser-driven tin plasmas [43]. Second, the moderate ionisation degrees (Z ≈ [11][12][13][14][15] found in EUV source plasmas implies that the free-electron contribution to the pressure p e = Zn ion kT e (in the ideal gas approximation) far outweighs the ion contribution to the pressure p i = n ion kT e (n ion is the ion number density). As such, the ion temperature will play a near-negligible role in the context of the current study.
The equations of single-fluid, single-temperature hydrodynamics take the form In the above equations, ρ is the fluid mass density, v is the fluid velocity, p = p e + p i is the pressure, E = e int + |v| 2 /2 is the mass-specific total energy (sum of the internal and kinetic energy contributions), S T represents thermal conduction, S R is the volume-specific heating rate provided by thermal radiation and S ext represents any external energy sources, e.g. energy deposition from a laser beam, ion beam, etc.
In the single-fluid approach, the plasma is treated as a quasineutral fluid, i.e., the electron number density n e is related to the ion number density n ion through n e = Zn ion . It is gradients in the plasma pressure which drive plasma expansion in a quasi-neutral hydrodynamic framework. The principal criterion justifying this approach is that the local Debye length λ D = [T e /(4πn e e 2 )] 1/2 must be significantly smaller than the scale length of the flow variation L = n ion /|∇n ion |. While this condition is typically met in the hot, dense region of the plasma (λ D < 10 nm), departures from it may exist as the plasma expands and rarefies. In such situations, a kinetic description of the plasma is often employed, where particle-in-cell methods are used to evolve the ion and electron particle distribution functions [15,[44][45][46]. We will discuss the validity of using the hydrodynamic description of a plasma in the current work in section 4.

RALEF-2D
We have performed radiation-hydrodynamic simulations using the RALEF-2D code. This code was originally developed to provide theoretical support for laser-plasma experiments at GSI Darmstadt under moderate laser intensities 10 13 -10 14 W cm −2 [47,48]. More recently, the code has found application in modelling laser-driven plasma sources of EUV light [13,[49][50][51][52]. The hydrodynamic component of the code is based on an upgraded version of the fully-explicit CAVEAT code for ideal hydrodynamics [53,54]. RALEF-2D solves the single-fluid, single-temperature hydrodynamic equations (equations (4)-(6) in two spatial dimensions on a structured quadrilateral mesh in either Cartesian (x, y) or axisymmetric (z, r) coordinates using a second-order Godunov-type method [54]. The axisymmetric coordinate system has been used in the present simulations. In the RALEF-2D code, radiation transfer and heat conduction are coupled into the fluid energy equation using a symmetric semi-implicit method with respect to time discretisation [55]. The code solves the LTE radiation transfer equation in the quasi-static approximation using pre-tabulated absorption coefficients generated with the THERMOS code [56,57]. The equation of state (EOS) of tin was built using the Frankfurt equation of state (FEOS) model which can treat both lowtemperature liquid-gas phase coexistence regions as well as high-temperature plasma states [58]. The FEOS model supplies the RALEF-2D code with key thermodynamic quantities such as the pressure, mass-specific internal energy as well as the average charge state of the plasma.
The simulations were performed on a computational mesh shaped in a half-disk consisting of multiple blocks with initially distinct properties. A simplified representation of the mesh is shown in figure 2. Centered in the origin of the (z, r) coordinate system we define a tin 'droplet' having a mass density of 6.9 g cm −3 and a temperature of 592 K (0.051 eV). As in the experiments, the droplet diameter is set to 28 μm. The bulk of the droplet is defined in a mesh section constructed as a rectangle stretched to a half-disk with dimensions of 45 × 90 mesh cells. As shown in figure 2, the outer region of the droplet exhibits a more refined mesh structure. Here, the length of each successive mesh cell along the radial direction decreases with increasing r. The mesh cell length on the outer droplet boundary is approximately 10 nm. Outside the droplet the mesh is filled with a tin vapour having a mass density of 10 −12 g cm −3 . This section is divided into quadrilateral cells by approximate concentric and radial edges and extends to a radial distance of 10 mm.
In essence, the experimental laser parameters have been replicated in the simulations. The laser beam is circular and coaxial to the positive z-axis. The simulations have employed unpolarized laser light. The laser absorption coefficient was calculated from the complex dielectric permittivity of the plasma [59].
The ion energy distribution is extracted from the fluid simulation by considering mass flow through the computational mesh. The main fluid variables density and velocity are assigned to each mesh cell throughout the simulation and are converted to the quantities mass and speed (velocity magnitude). These variables are cell-centered and form the basis of keeping track of the fluid throughout the simulation. As the curved boundary of the computational mesh is defined as a free-outflow boundary, matter flowing out of the mesh leaves the computational domain; it leaves the simulated area in space. This is closely related to the treatment of the ion energy distribution by RALEF-2D. Mass flowing out of the mesh is recorded ('binned') in the (i) energy bin corresponding to its speed and (ii) the angular bin corresponding to the angle between the laser axis and the escape velocity vector. This module is called at every hydrodynamic time step, summing the number of particles equivalent to the outflowing mass. This procedure explicitly constructs the distribution of the number of particles into 360 predefined discrete energy bins in the range [1, 20 × 10 3 ] eV. The bin width increases exponentially with increasing energy. The angular domain is divided into 36 bins (over 180 • ). In the current simulation we consider mass flow into two angular bins extending over the range [55 • , 65 • ]. The duration of the simulation is 1 μs which allows accounting for ions with energies down to ∼ 70 eV leaving the computational domain in this time window.

Plasma formation and expansion
In figure 3 we present the evolution of the plasma expansion through the variables speed and ion number density, where the pseudocolour indicates the magnitude of these variables. At distances larger than 0.5 mm the velocity vector effectively points radially outwards. The ion number density n ion is obtained from n ion = ρN A /A, where ρ is the fluid mass density, N A is Avogadro's constant and A is the atomic weight of tin. For visibility, we reflect the ion number density information into the lower plane (this is possible as the simulations were performed using the axisymmetric (z, r) coordinate system). We define t = 0 ns as the time that the laser pulse is switched on in the simulations. The left column shows times t = {11, 15} ns, the middle column t = {25, 35} ns and the right column the late times t = {60, 120} ns. The laser propagates along the positive z axis (laser axis) and its (local) intensity is represented by the black shading seen in the t = {11, 15} ns frames. Frames grouped in the same column, e.g. t = {11, 15} ns share the same axial and radial domains. In order to follow the plasma expansion in space, we increase the axial and radial coordinate domains in the t = {25, 35} and {60, 120} ns frames.
The overall dynamics of plasma formation and expansion, as displayed in figure 3, can be qualitatively described as a succession of two distinct bursts of laser-induced ablation from the droplet surface. These two bursts are clearly identified in the t = 15 ns and t = 25 ns frames as two concentric red regions exhibiting high speed. In the following two subsections we describe the formation and evolution of these two ablation bursts.

Initial burst of laser-induced ablation
The initial burst of laser-induced ablation forms in the first 2-3 ns after the laser pulse is turned on. Initially, the laser pulse has an intensity I ≈ 3 × 10 8 W cm −2 which lies only moderately above the ablation threshold of liquid tin. In figure 4 we plot one-dimensional (1D) profiles of the mass density ρ (black), temperature T (red), fluid speed |v| (orange dashed curve) and mass-specific heating rate of the laser q (green) along the negative laser axis (starting from the droplet center) at times t = (a) 0.5, (b) 1 and (c) 1.4 ns, respectively. These data are extracted along a so-called 'lineout' taken at θ = 0 • with respect to the axial coordinate (laser) axis.
As is evident from figure 4(b), the temperature of the ablated plasma T ≈ 1-2 eV exceeds the critical temperature of tin (T critical ≈ 0.5 eV) by only a moderate factor. In addition, we notice the emergence of a 'hump' in the speed profile at d 0 • ≈ 15 μm which is coincident with the location of the peak value of q, the mass-specific laser heating rate. Material associated with this hump accelerates and eventually 'catches up' with the initially-ablated material. By t ≈ 2.5 ns, the speed profile exhibits a near-linear dependence on distance d, i.e., |v| ∝ d and the speed of the front edge of the plasma cloud stabilizes to |v front | ≈ 33 km s −1 . This expanding plasma cloud drives a shock into the low-density (ρ = 10 −12 g cm −3 ) ambient gas which, however, has a negligible effect on the overall expansion dynamics.

Second burst of laser-induced ablation
The increase in laser intensity (and subsequent increase in the mass-specific laser-heating rate q) as t → 15 ns generates a second burst of ablation visible in the upper halves of the t = 11, 15 and 25 ns frames in figure 3. This second burst is characterized by a significantly higher ablation rate, density and temperature (T > 30 eV) of the ejected plasma whose leading edge quickly reaches a speed |v front | ≈ 60 km s −1 .
To better understand the dynamics of this ablation burst along the line-of-sight of the ESA and FC devices, we show in figure 5 the (a) speed and (b) ion number density profiles along a lineout taken at θ = 60 • with respect to the laser axis at times t = 5 (black), 10 (red), 15 (blue), and 20 ns (green), respectively. This lineout is shown as a black dashed line in the t = 11 ns frame in figure 3. The shaded vertical bars indicate the location of local maxima in the number density lineouts and serve to guide the eye between both variables. From figure 5 we see that as the second and more powerful ablation burst rams into the background plasma (left behind by the initial ablation burst), it rakes up material into a quasi-spherical expanding shell. This shell is evident as a hump in the t = 15 and 20 ns n ion profiles in figure 5(b). The deposition of energy (via laser absorption) and the expansion that follows can be likened to the effects of a fast piston pushing a gas; the shockwave launched by the piston effectively sweeps up material in front of it, driving a compression wave.
In figures 6(a) and (b) we plot the speed and ion number density profiles along the θ = 60 • lineout at times t = 30 (orange), 50 (light blue), 80 (purple) and 120 ns (dark green). As the plasma expands, we note that the two speed profiles merge to form a single profile exhibiting a nearlinear dependence on distance. A linear dependence of speed on distance, i.e., |v| = d/t is exactly the asymptotic t → ∞ behaviour of any 'explosion-like' expansion [60]. Our simulations therefore recover the late-time behaviour expected from such an expansion. The evolution of the ion number density profiles over time and space is less dramatic. First, we note that the quasi-spherical expanding shell observed in figure 5(b) persists throughout the expansion. As time evolves the shell is observed to broaden, a direct consequence of the non-constant speed profile across the shell. It is interesting to note that the spatial variation of the ion number density up to the quasi-spherical shell appears to follow a power law of the form n ion ∼ d −n where 2 < n < 3. Unlike the speed profile, a universal analytic form for n ion as t → ∞ does not exist for an explosion-like expansion. As discussed by Zel'dovich and Raizer [60], the asymptotic solution |v| = d/t is satisfied for n ion = φ(d/t)/t 3 where φ(d/t) is an arbitrary function of d/t. This function can only be evaluated through numerical simulations of the system and is case-specific.
Before proceeding with the comparison of experimental and simulated ion kinetic energy distributions, we wish to make two remarks. First, recall that the validity of the single-fluid description of a plasma relies on the condition λ D L where λ D is the Debye length and L is the characteristic flow length. We have calculated λ D and L for the late time case t = 120 ns and have found that λ D ≈ 1-10 μm and L ≈ 500 μm in the vicinity of the density hump. These results validate the use of the quasi-neutral hydrodynamic approach in the current context. It is important to mention that the mechanism of ion acceleration in the other extreme case, i.e., plasmas in which λ D L is often referred to as 'Coulomb explosion' [38]. Second, recall that in the simulations the region outside the droplet is filled with a low-density tin vapour having ρ = 10 −12 g cm −3 . The experiments, on the other hand, have been performed in a near-vacuum environment. This choice of density is sufficiently low as not to distort the vacuum-like expansion we wish to emulate in our simulations. This is evident from  figure 6 where we do not observe a decrease of the peak velocity as the fluid propagates through the low-density background gas.

Ion kinetic energy distributions: experiment and simulation
We now wish to compare our measured ion energy distribution with that obtained from the simulations. These two quantities are compared in figure 7. The experimental data, shown in red, corresponds to the total ion energy distribution (also shown in red in figure 1) and the solid black curve is the ion energy distribution obtained from the RALEF-2D simulations. As described in section 3.2, the RALEF-2D ion energy distribution is obtained by recording mass flow into two angular bins subtending an angle 55 • < θ < 65 • with respect to the laser axis. From the two-dimensional computational mesh we have calculated the three-dimensional solid angle Ω S ≈ 0.95 sr by revolving the arc length in the mesh around the laser axis. The RALEF-2D ion energy distribution, once corrected for this solid angle, is then convolved with bin-specific Gaussian functions having full width at half-maxima equal to 5% of the lower boundary of the energy bin. The purpose of this convolution is to account for processes which may broaden the distribution, e.g. mass distribution of tin isotopes.
It is clear from figure 7 that the RALEF-2D ion energy distribution closely resembles the experimental measurements. First, the simulation reproduces the high-energy peak observed in the experimental data near 2 keV. This Figure 7. The distribution of the number of ions over ion kinetic energy. The experimental ion energy distribution is shown in red (solid curve) and the RALEF-2D ion energy distribution is shown in black (solid curve). Also illustrated are the predictions of the analytic models of Murakami (green dash-dot) [21], Mora (green dash) [25] and the Riemann wave (green solid) [60]. high-energy feature originates from fast-moving material associated with the quasi-spherical expanding shell (a tin ion with |v| = 57 km s −1 has a kinetic energy E ≈ 2 keV). Second, the simulations reproduce the near-constant behaviour of the experimental ion energy distribution in the 0.07-1 keV range. Fluctuations observed in the RALEF-2D ion energy distribution most likely arise from spatial fluctuations in the density during the expansion (visible in figure 6 for t = 120 ns). Above 2 keV, both the simulations and experimental data exhibit a sharp fall-off with increasing kinetic energy. This fall-off is sharper in the case of the simulations, which do not predict any ions having kinetic energies above 3 keV. It is also interesting to note that, within the limits of the experimental uncertainties, the simulations provide a reliable prediction for the absolute number of ions detected in the experiments.
We show in figure 7 the predictions of the models of Murakami (equation (1)) and Mora (equation (2)). Guided by the work of Torretti et al [13], we have taken Z = 12, k B T = 35 eV and N 0 = 2 × 10 12 . For the model of Murakami we have chosen a spherical expansion (α = 3) and a value ln[R(t)/R 0 ] = 4. We also provide in figure 7 the ion energy distribution arising from a planar isentropic expansion, better known as the Riemann wave [60]: Here γ is the adiabatic index (we have taken γ = 4/3 [61]) and E max = 2Zk B T e γ/(γ − 1) 2 is the maximum ion kinetic energy. Both the planar isothermal expansion model of Mora and the Riemann wave solution predict a similar monotonic decrease in dN/dE with increasing E which is not observed in the experimental measurements. The shape of the experimental data is qualitatively better described by the spherical isothermal expansion model of Murakami, which exhibits a slow rise in dN/dE up to a peak at E peak = E 0 /2 ≈ 1.7 keV. The fall-off in dN/dE at energies above E peak is far less steep compared to the experimental data. The reason why the models of Murakami and Mora and the Riemann wave solution cannot be expected to reproduce the current experimental distribution ultimately lies in the plasma density profiles adopted in these analytic models. The function φ(d/t) obtained with the RALEF-2D simulations is significantly more complex than the Gaussian and exponential density profiles assumed in the models of Murakami and Mora, respectively. It is the interaction of many complex processes (2D expansion of a non-uniformly heated, non-planar (radiating) plasma) that ultimately determines the function φ(d/t).
Four possible causes have been identified which may contribute to the observed differences between the experimental and simulated ion energy distributions. The first cause is the effect of numerical diffusion in the large mesh cells at larger mesh radii. This has been partly tackled by increasing the radial detail in the mesh at larger distances. The three other causes are inherent to the ansatz of the simulation code. As mentioned in section 3.2, RALEF-2D uses Godunov's method for the Lagrangian phase of each hydrodynamic cycle. As the internal energy component of the total energy determines the pressure, rounding errors can propagate especially if the internal energy is small. A third possible cause is related to the EOS model employed in RALEF-2D. The EOS model adopted in this work assumes LTE ionization throughout the entire simulation. This assumption breaks down at late times in the expansion when ionization and recombination processes cease to exist, leading to the well-known 'freezing' of charge states [34].
We note that the simulated domain (10 mm) is much smaller than the experimental flight path (∼1 m) and the assumption is made that neither the experimental nor simulated ion energy distributions change significantly between the two distances. This assumption is supported by two arguments: (1) within the simulated spatial scale the velocity profile attains its asymptotic 'triangle-like' shape |v| = d/t before leaving the mesh, having converged to the late-time behaviour; (2) on the experimental side the aforementioned freezing of charge states will occur on a length scale similar to that of the simulation spatial scale [34] and, thus, no significant changes over the remaining flight path will occur in our high-vacuum environment. Keeping these remarks in mind, the results presented in this paper demonstrate that the single-fluid single-temperature approach implemented in RALEF-2D can (i) reproduce the general shape of the experimental ion energy distribution and (ii) provide a reliable prediction for the absolute number of ions detected in the experiments.

Conclusion
We have undertaken a joint experimental and theoretical study of plasma expansion arising from Nd:YAG laser irradiation of tin microdroplets. The experimentally-recorded ion energy distribution is found to exhibit a complex, non-monotonic dependence on ion kinetic energy. Charge-state resolved measurements of the ion energy spectra reveal the existence of peaks centered near 2 keV in the Sn 3+ -Sn 8+ distributions. Two-dimensional radiation-hydrodynamic simulations performed using a single-fluid single-temperature approach are shown to reproduce the overall shape of the experimentallyrecorded ion energy distribution and provide a reliable prediction for the absolute number of ions detected in the experiments. The existence of a peak in the experimental ion energy distribution near 2 keV is attributed to the formation of a quasi-spherical expanding shell at early times in the plasma expansion. Our interpretation of the plasma dynamics in terms of two distinct bursts of laser-induced ablation indicates that the observed ion energy distribution would in general be sensitive to the temporal profile of the laser pulse. The results of the present work are therefore specific for a Gaussian temporal profile with a laser intensity varying on the timescale of ∼10 ns.