Droplet evolution in expanding flow of warm dense matter
DDroplet evolution in expanding flow of warm dense matter
Julien Armijo ∗ Lawrence Berkeley National Laboratory, Berkeley, California, USA
John J. Barnard
Lawrence Livermore National Laboratory, Livermore, California, USA (Dated: April 21, 2019)We propose a simple, self-consistent kinetic model for the evolution of a mixture of droplets andvapor expanding adiabatically in vacuum after rapid, almost isochoric heating. We study the evo-lution of the two-phase fluid at intermediate times between the molecular and the hydrodynamicscales, focusing on out-of-equilibrium and surface effects. We use the van der Waals equation of stateas a test bed to implement our model and study the phenomenology of the upcoming second neutral-ized drift compression experiment (NDCX-II) at Lawrence Berkeley National Laboratory (LBNL)that uses ion beams for target heating.We find an approximate expression for the temperature dif-ference between the droplets and the expanding gas and we check it with numerical calculations.The formula provides a useful criterion to distinguish the thermalized and nonthermalized regimesof expansion. In the thermalized case, the liquid fraction grows in a proportion that we estimateanalytically, whereas, in case of too rapid expansion, a strict limit for the evaporation of droplets isderived. The range of experimental situations is discussed.
PACS numbers: 64.70.fm, 51.10.+y
I. INTRODUCTION
Warm dense matter (WDM) conditions, correspondingroughly to densites 0 . < ρ/ρ solid <
10 and tempera-tures 0 . < T <
10 eV, can be defined as the regionof thermodynamic space corresponding to the doublecrossover from degenerate to non-degenerate and fromweakly to strongly coupled matter [1], so that the “easy”limiting descriptions in terms of cool plasma and hot con-densed matter meet and have to be somehow connectedto each other. This problem is drawing growing attentionbecause of the serious theoretical challenges involved, andbecause of the occurence of WDM in the contexts of iner-tial fusion energy (IFE), astrophysics (planet cores), andlaser ablation for materials processing, nanoparticles for-mation and film deposition [2–4].Generally, WDM experiments are inertially confinedand explosive. Rapid heating of the material is re-quired, so that the energy deposition (by lasers, ions,neutrons, electrical discharges, etc.) is faster than itsrelease through hydrodynamic expansion. Pressures inthe kbar to Mbar range can be reached before giving riseto supersonic expansion with typical outflow velocities ofseveral km/s.The two-phase region of the phase diagram belongsonly partly to the WDM regime, with its high temper-ature part near the critical point. However, any WDMexperiment will almost inevitably lead to two-phase con-ditions during the expansion and cooling process. Thishappens either from below the critical point (ion heatingexperiments at Gesellschaft f¨ur Schwerionenforschung inDarmstadt [5, 6], or second neutralized drift compres-sion experiment (NDCX-II) at Lawrence Berkeley Na-tional Laboratory (LBNL) [7], low fluence laser ablation, Z machines), or from above it (IFE, high fluence laser ab-lation, upcoming ion heating machines). In the first casean overstretched liquid fragments and evaporates into amixture of droplets and gas, whereas in the second case ahot supersaturated gas nucleates small clusters while ex-panding. In both cases the flow becomes a plume of gasand condensed clusters (most often liquid, so the term”droplet” is appropriate) : the monophasic liquid or gashas undergone phase separation with the creation of sur-faces giving a non-trivial geometry to the fluid, whichmay a priori affect its dynamical properties.Recently, there has been significant progress in the ob-servation of those two-phase flows, from the early abla-tion and plume expansion stages in the the ps and nstimescales [8–11] to the late µ s timescale evolution in-cluding postmortem analysis of the clusters [12, 13].Basic questions arise when considering a two-phaseflow. First, what is the droplets’ size and distribution,and how do they evolve during the expansion? Second,are the droplets and the gas in thermal equilibrium? Theanswer can determine the conditions of validity for hydro-dynamic approaches based on the Maxwell constructionor any two-phase equation of state (EoS) that assumeslocal equilibrium.So far, two main approaches have been used. Onone hand, molecular dynamics (MD) codes [14–20] com-pute the dynamics of each particle separately, and havegiven powerful insight into the processes of fragmenta-tion, phase explosion, and the different mechanisms forablation, but they are inherently limited to only treatthe early times ( < ∼ ). On the other hand, hydrodynamiccodes [13, 22–25] can model experiments completely, butthey deal with mesoscopic fluid cells, and often rely on a r X i v : . [ phy s i c s . f l u - dyn ] M a y crude approximations concerning the molecular and ki-netic processes involved. Complex hydrodynamic codesincluding a treatment of the kinetics of phase change pro-cesses and surface effects in each cell are under develop-ment [13, 26], but providing a complete and reliable de-scription of a whole WDM experiment is still a challenge.Better understanding of two-phase flows should behelpful for the preparation of experiments, including thediagnostics, and for the interpretation of the data. InIFE especially, the problem of debris dynamics is a cru-cial issue due to their impact on the optics and othercomponents of the target chambers [26, 27].In this paper we propose an alternative approach tostudy two-phase flows in the cool or late time WDM sit-uations. Our model was initially conceived to predictthe phenomenology of the upcoming target heating ex-periments with the NDCX II machine at LBNL wherean ion beam will almost isochorically heat a thin metal-lic foil to temperatures of about 1eV. However, the modelshould apply to any two-phase flow.The core of our model is a self-consistent set of kineticrate equations for the particle and energy fluxes betweena droplet and the surrounding gas in an expanding La-grangian cell. This set of equations can be applied to anytwo-phase EoS. The computing cell is considered as partof a larger hydrodynamic code, but in this paper we onlyconsider one cell containing one droplet. We also neglectseveral features that could be added in the near future.To implement the kinetic equations and explore thepatterns of two-phase expansion, we use the van der-Waals fluid model that makes it possible to build a com-plete set of thermodynamic functions.We thus demon-strate the ability of the kinetic model to simulatenonequilibrium two-phase flows in the wide range be-tween themolecular and hydrodynamic scales. In par-ticular, we use it to distinguish the different regimes oftwo-phase expansion: on one side, quasi- or fully ther-malized; on the other side, nonthemalized.We show thatthis distinction depends on the initial target dimensionsand the initial temperature. We then study those regimesanalytically and numerically. II. BACKGROUNDA. Expanding two-phase flows. Supercritical andsubcritical cases
The model that we propose lies at a mesoscopic scalebetween the molecular and hydrodynamic scales, so weneed some preliminar assumptions. Our computing cellis considered as an elementary piece of a larger hydro-dynamic code describing an expanding flow. The lin-ear strain rate η characterizes the expansion of the cell L = L (1+ ηt ), where L is the initial cell size. We definethe hydrodynamic time t hydro = η − . In the following, we assume rapid heating ( t heating < t hydro ) so that en-ergy deposition in the material is almost isochoric. Forsimplicity, we assume instantaneous energy deposition.To get insight into the global flow, it is interesting toreview some analytical and numerical results. The selfsimilar rarefaction wave (SSRW) is the solution [28] de-scribing the one-dimensional (1D) expansion of a perfectgas (semi-infinite at z <
0) of adiabatic coeficient γ af-ter instantaneous uniform heating at the initial temper-ature T . Denoting c s the sound speed in the fluid at T ,the SSRW solution describes an outward expanding fronttravelling at the outflow velocity v = 2 c s / ( γ − c s for a perfect monoatomic gas, while the inward rar-efaction wave propagates at c s [24]. Note that the SSRWcan be computed semi-numerically for any EoS of a non-ideal gas [29] and has been validated by MD simulations[18].As an example of a numerical simulation of expandingflows, Fig. 1 shows a hydrodynamic calculation with thecode DPC using an EoS based on Maxwell construction[23]. Here the liquid and gas are assumed in equilibrium,which is not kinetically justified (see Sec. IV.B), and theoutflow velocity is about 8km/s after 10ns. FIG. 1. Hydrodynamic calculation with DPC code of NDCXII reference case, from [24]. A 3.5 µ m-thick Al foil is heatedwithin 1ns with an ion beam and subsequently cools downduring adiabatic expansion In the following we assume a flow with linear speedprofile and outflow velocity v = 3 c s , but it is worth re-marking that this is quite simplistic. In particular, sev-eral numerical works [22, 24, 25, 29] have reported theoccurrence of ”plateaus”, that is, zones of nearly con-stant density, related to the fluid zones entering into thetwo-phase regime.Let us now present our simple classification of two-phase expanding flows, which is based on equilibrium thermodynamics, and is similar to the one in [30]. Ofcourse, such procedure cannot account for the complex-ity of non-equilibrium situations encountered in experi-ments or simulations [9, 11, 17–19, 30]. The two-phaseregime exists only for temperatures T < T c , and for anaverage density between the value at the liquid and gasbinodals, as seen in Fig. 2.a. Thus, a piece of fluid ex-panding near thermodynamic equilibrium can enter thetwo-phase regime, which amounts to partially undergo-ing the liquid-gas first-order phase transition, in only twoways: by crossing the liquid binodal, which we call the subcritical case , or by crossing the gas binodal, which wecall the supercritical case .In the subcritical case, which corresponds to the calcu-lation of Fig. 1 when it is mapped onto the correspondingphase diagram, an expanding monophasic liquid reachesthe liquid binodal and becomes overstretched. The equi-librium configuration for a piece of fluid in the two-phaseregion is a mixture of liquid droplets (of yet undeterminedsize) and gas whose densities are at the binodals for thesame temperature. We call fragmentation the transfor-mation from the monophasic liquid to the non-connectedcloud of droplets. In the supercritical case, achieved ifthe material is initially heated to higher temperatures,a monophasic gas becomes supersaturated after crossingthe gas binodal, and we call nucleation the process bywhich a certain distribution of liquid droplets is created.Figure 2 represents the two cases and the various ex-perimental situations that they involve. In Fig. 2.a, weshow the Van der Waals phase diagram for Al that weuse in the following, and a schematical representation ofthe sub- and supercritical cases of two-phase expansion(arrow 1 and 2). In Fig. 2.b, reproduced from [17], onesees the particle distribution in a 2D MD simulation oflaser ablation. Due to the inhomogeneous energy depo-sition, different types of thermodynamic evolutions areseen at a same time, and we use this picture to illustratethe various situations of our classification, although notfollowing exactly the terminology of the original work[31]. In zone I, the dense material is still continuous. Inzone II, the expanding liquid has undergone cavitationof gas bubbles. In the upper zone II and in zone III,the liquid is fragmented and the material has entered thetwo-phase regime in the subcritical case. In zone IV, thefully atomized, expanding material is likely to reach thegas binodal in a later stage, thus corresponding to thesupercritical case. In the upper zone III, small clustersare present, which could originate from fragmentation athigh temperature (subcritical case) or recent nucleation(supercritical case). This is why in Fig. 2.a., where thedifferent zones are placed qualitatively , two positions areproposed for zone III. B. Initial droplet size
Both cases lead to droplets formation. In order to ini-tialize the kinetic model that we present further, it isnecessary to know the initial droplet size at the onset ofthe two-phase regime.In the subcritical case, the overstretched liquid startscavitating (see Fig. 2.b, zone II), which we call the bubblesregime and then the bubbles percolate until the liquidphase is not continuous anymore (see Fig. 2.b, zone III),which we call the droplets regime . We assume that thedroplets regime starts when the gas and liquid volumesare equal: V g = V l , which is justified by an argument of (a) (b) IIV
Supercritical fluidGas Liquid2-Phase ! cm " T ! K " III IIIII
FIG. 2. (Color online) (a) Phase diagram of the van derWaals EoS for Al (see Sec. III) showing the liquid and gasbinodals (solid lines), and a schematic representation of thesubcritical (arrow 1) and supercritical (arrow 2) cases of two-phase expansion. (b) 2D MD simulation of laser ablationwith inhomogeneous initial temperature, from [17], showingmaterial in various situations of two-phase expansion, whichwe also locate qualitatively on the schematic classification ofFig. 2.a (roman numbers). surface energy minimization.The mean droplet’s size in a fragmentation scenariocan be obtained by considering a balance between thedisruptive inertial forces and the restoring surface tension[14]. The model proposed initially by Grady [32] hasbeen abundantly validated by MD calculations [14, 17,20] in 2D and 3D, and is in very good agreement withmeasurements on He jets [33]. We note that the scaling ofthe mean radius R of the droplet can be obtained by justsetting to unity the Weber number W e ≡ ρRv /σ [34],where σ is the surface tension ρ the liquid mass density,and v = ηR the typical velocity difference across a pieceof fluid of size R . W e is the ratio of the surface energy tothe inertial energy. In any dimension this criterion yieldsWe ∼ ⇒ R ∼ (cid:18) σρη (cid:19) . (1)Several values of order 1 have been proposed for the pref-actor in this scaling law, either from analytical estimates(prefactor 15 / = 2 .
47 in [25]), or from fits to MD sim-ulation results. In [15], it was shown that both 3D MDresults with a homogeneous strain rate η and data fromhelium free jets experiments from [33] could be fitted toEq. 1 with the same prefactor, thus validating this lawover almost 8 orders of magnitude in cluster mass (theexperimental fragments cover larger sizes than the nu-merical ones).Concerning the size distribution of droplets resultingfrom fragmentation, MD simulations have shown clearlythat it is essentially exponential [14–16, 20], which is con-sistent with simple models of entropy maximization [14].By contrast, it is not so clear how to describe the initialsituation in the supercritical case. This task requires oneto choose a model for nucleation, or to input results fromMD calculations. Nucleation of clusters from a supersat-urated vapor is the situation of nucleation whose kineticsis the easiest to model theoretically [35], but still choiceshave to be made [13], that are beyond the scope of thispaper. Any model for nucleation will depend cruciallyon surface tension, so we make the remark here that es-timating the surface tension for small droplets is delicatebecause of its enhancement at small sizes [36, 37]. III. MODELA. Van der Waals fluid model
The kinetic model that we present in Sec. III.C isapplicable to any EoS. In this paper, for a qualitative in-vestigation of the two-phase expansion regimes, includingkinetic and surface effects, and with emphasis on the latetime and low temperature limits, we use a van der Waals(VdW) fluid model, which for convenience we describefirst. It is important to note, however, that the VdWmodel is not intended to provide a highly accurate de-scription of a fluid, especially in WDM or supercriticalconditions where ionization and radiation effects can bestrong.With only two parameters, the VdW EoS is the sim-plest EoS describing the coexistence of a liquid and a gasphase, and has already been used for theoretical studiesof dynamic two-phase processes [38, 39]. All the thermo-dynamic functions can be derived from the expression forthe mean-field potential energy per particle in such fluid: U = + ∞ if n > /b and U = − an if n < /b , where n is the particle density, b stands for the incompressiblevolume of the particles, and a represents the mean-fieldattractive energy between them.The bulk VdW energy of N particles at temperature T is E = N ( c v T − an ). It can be shown that the spe-cific heat c v is independent of n and can only depend on T [40], so that one has to choose necessarily c v = k B ,where k B is the Boltzmann constant, if one wants theEoS to match the perfect monoatomic gas in the dilutelimit. Writing the partition function, one obtains theother thermodynamic functions. In particular, the pres-sure is P = k B T / ( v − b ) − a/v where v = 1 /n is the vol-ume per particle. This expression implies that the isobars(isotherms) are a cubic relationship between T and v ( P and v ). Hence, below a certain critical temperature T c ,an unstable zone of negative compressibility appears inthe phase diagram, limited by the two spinodals. We ob-tain the equilibrium density of the two stable phases thatcan coexist at certain ( P, T ) by numerically perform-ing the Maxwell construction, which consists of solving P l = P g (i) and µ l = µ g (ii) simultaneously, where µ denotes the chemical potential and the subscripts l and g stand for liquid and gas , respectively, all throughoutthis paper. (ii) is equivalent to (cid:82) gl vdP = 0 and thus (cid:82) gl P ( v ) dv = P l,g ( v g − v l ) [40].Introducing the reduced temperature θ = k B T /l ,where l = a/b is the latent heat at T = 0, and twodimensionless parameters that are small in the low tem-perature limit: v g = b/δ and v l = b (1 + γ ), equations (i)and (ii) become: θγ − γ ) = θ δ − − δ and (2) θ ln (cid:18) δ − γ (cid:19) + δ −
11 + γ = (cid:18) θ δ − − δ (cid:19)(cid:18) δ − (1 + γ ) (cid:19) . (3)It is worth remarking that T c = 8 a/ b , so θ c = 8 / (cid:39) .
3, and therefore one expects that calculations in the”low temperature limit” ( θ (cid:28)
1) should be a good ap-proximation as soon as one is not considering the vicinityof the critical point.Figure 3 gathers the thermodynamic functions of ourVdW model for aluminum. Figure 3.a shows the nu-merical result of the dimensionless Maxwell constructionwhere the VdW parameters a = 9 . × − erg . cm and b = 1 . × − cm , giving l = 3 . eV , have beenadjusted to fit this material. For that, we impose thatthe VdW liquid density matches the aluminum liquiddensity n l ( T m ) = 5 . × cm − at the melting point T m = 933 .
5K (= 0 . l ) [41] and that the VdW latentheat (shown in Fig. 3.b) l = a ( n l − n g ) + P l,g (cid:18) n l − n g (cid:19) (4)coincides with the experimental value l ( T b ) = 4 . × − erg/atom for aluminum at the boiling temperature T b = 2792K (= 0 . l ) [41]. Note that the critical pa-rameters obtained in this way are consistent with the bestestimates to date, although not very precisely [42].In Fig. 3.a are also displayed the simple, useful ap-proximations for n l and n g at lowest orders in θ that oneobtains directly from Eq. 2 and 3 : n l (cid:39) b (1 − θ − θ ) , (5) n g (cid:39) bθ exp (cid:18) − θ (1 + θ ) (cid:19) . (6)Note that Eq. 6 is, at lowest order in θ , equivalent to theClausius-Clapeyron formula applied to a perfect gas. Wealso show in fig 3.c the approximations at first order in θ for the liquid and gas bulk energies per particle e g (cid:39) k B T , e l (cid:39) k B T − ab . (7)For the surface tension, Van der Waals himself had al-ready proposed to model it using density gradients [39],but we have chosen to use a simple formula that is uni-versally verified in simple fluids [43] σ ∝ (1 − T /T c ) r with r = 0 . . (8) oo0 2 4 6 8 100.00.20.40.60.81.0 T ! K " n ! un i t s o f b " ! K " l ! un i t s o f l " ! K " Σ ! un i t s o f a b " (c) (d)(b) n l n g e l e g (a) " " " " " ! K " e ! un i t s o f l " FIG. 3. (Color online) Van der Waals thermodynamic func-tions for Al, in VdW units. The filled circles represent exper-imental data and the empty cirlce the critical point obtainedfrom the fits. (a) Liquid and gas densities (solid lines) withfirst (dashed lines) and second order (dotted lines) low T ap-proximations. (b) Latent heat (solid lines) decomposed in thefirst (dashed lines) and second (dotted lines) term of Eq. 4.(c) Bulk energies per particle, with first order low T approx-imations Eq. 7 (dashed lines). (d) Surface tension (Eq. 8). To model aluminum, we fit this formula to the exper-imental value σ ( T m ) = 1050erg / cm [44], as shown infig. 3.d. Note that in the following, the total liquid en-ergy in the cell is E l = N l e l + σS l , where S l is the surfacearea of the droplet. B. Kinetic equations. Validity condition
Our goal is to compute the evolution of droplets incells. In this paper we limit ourselves to the case of onedroplet in one Lagrangian cell undergoing adiabatic ex-pansion. It is a closed system out of equlibrium, and itscomplete description requires the determination of thefour variables N l , n l , T l , T g . To compute their evolution,we need four rate equations: a liquid-gas particle ex-change rate, an energy exchange rate, a total energy lossrate (work to the outside), and an internal equilibriumcondition to determine the liquid density. As shown infig. 4.a, the particle fluxes between liquid and gas are di-vided between evaporating, condensing, and condensingbut not-sticking particles.The volume expansion V ( t ) shall later be prescribed bya global hydrodynamic code. For our study, we assumecylindrical symmetry and we use a simple model behavior[15, 16] V ( t ) = V (1 + η z t )(1 + η r t ) , (9)where V = L is the inital cell volume and η z =( dv z /dz ) t =0 and η r = ( dv r /dr ) t =0 are the axial and ra- dial strain rates, respectively defined as the initial veloc-ity gradients in the beam direction and the target plane.The evaporation and condensation fluxes are computedusing the standard Hertz-Knudsen formulas [45–47]Φ cond = n g (cid:114) k B T g πm , Φ vap = n ∗ g ( T l , R ) (cid:114) k B T l πm , (10)where m is the particle mass and n ∗ g ( T l , R ) is the equi-librium gas density for a droplet at temperature T l andof radius R . To estimate n ∗ g , we use the Kelvin equa-tion, which describes the increase of the equilibrium va-por pressure surrounding a droplet due to surface tension n ∗ g ( R ) = n ∗ g ( ∞ ) exp (cid:18) σk B T n l R (cid:19) . (11)The Kelvin equation is approximate because its deriva-tion assumes a perfect gas. Also, we use a constant valuefor σ , thus neglecting its increase at small radii [36, 37].Still, this approach is probably not too bad after all [48],and satisfactory enough for our qualitative purpose. ! K " Β (a) (b) FIG. 4. (Color online) The ”droplet-in-cell” kinetic model.a) Sketch of the kinetic fluxes. b) Sticking coefficient β cal-culated with formula from [49] and our aluminum VdW pa-rameters. Considering mass conservation, and combining the twofluxes of Eq. 10, the particle exchange rate equations are d ( N l + N g ) dt = 0 , dN l dt = β ( − Φ vap + Φ cond ) S l , (12)where S l is the surface area of the droplet and 0 < β < sticking coefficient that is usually assumed of order0.5. A recent study [49] has proposed a simple expressionfor β that is in good agreement with MD calculations forseveral simple fluids. This expression depends only onthe ratio of the molecular volumes in the liquid and vaporphase: β = (1 − ( v l /v g ) / ) exp( −
12 1( v l /v g ) / − ) that weplot for our VdW model for Al in Fig. 4.b. Note howeverthat, for the calculations presented in this paper, we usea constant value β = 0 . d ( E l + E g ) dt = − P g dVdt . (13)In this global energy loss rate we have neglected threeterms that could be added in the near future. The firstone is heat conduction between cells. This term mayplay a role, but it cannot be very important as we areconsidering a supersonic flow (see Sec. I). The secondneglected term is radiation. Radiation becomes indeedthe dominant cooling mechanism at long times, as we seelater, but it is negligible for the initial dynamics, so theapproximation is reasonable, because our purpose in thispaper is to study the expansion in a time range where thetwo phases are interacting and the system is not just acollection of isolated clusters flying in vacuum. The thirdneglected term is thermionic emission. One expects elec-trons to be thermally emitted from the droplet, takingaway some energy. Non-neutral effects are totally absentfrom our model, but we expect that the associated cool-ing rates will be small compared to the adiabatic andevaporative cooling rates [50].The energy exchange rate between the liquid and thegas has contributions from the three fluxes of Fig. 4.a.The contribution of the colliding but non-sticking parti-cles can be described with a flux proportional to the tem-perature difference T g − T l , with a relaxation coeficient0 < α < e g . For the evaporating particles,we assume that the energy they individually take awayfrom the liquid depends only on the liquid state. We noteit e ∗ g and define it as the energy of a virtual gas particlethat would be in equilibrium with the droplet of radius R at temperature T l . For the VdW fluid, those energies are e g = c v T g − an g and e ∗ g = c v T l − an ∗ g ( T l , R ), respectively.Note that our definition of e ∗ g is totally analogous to theHertz-Knudsen derivation of the mass evaporation rateof Eq. 12. We finally get the exchange rate dE l dt = (cid:2) β ( − e ∗ g Φ vap + e g Φ cond )+ (1 − β ) αc v ( T g − T l )Φ cond (cid:3) S l . (14)Our set of kinetic equations is fully consistent in thesense that at equilibrium both mass and energy fluxesbetween the droplet and the gas are in balance. In partic-ular, the average energy that a liquid particle takes fromthe rest of the liquid to evaporate is e ∗ g − e l . This term,which for the VdW fluid is − a ( n ∗ g − n l ), corresponds ex-actly to the latent heat per particle for an arbitrary fluid l = ( e g − e l ) − P l,g ( v g − v l ), without the second (work)term, which is expected since the latent heat is an en-thalpy and we are here dealing with energy exchanges atconstant volume.To our knowledge, our set of rate equations is an orig-inal model for the exchanges between a droplet and itsvapor. Other sets of kinetic equations can be found foranalogous systems (see e.g. [52, 53]), but they do not cor-respond to the purely kinetic regime that we are consid- ering, because they deal with larger droplets ( R > µ m)and longer time scales, more relevant to the fields of com-bustion or atmospheric sciences, so they need to combinethe kinetic approach with the more classic hydrodynamictheories of droplet evaporation [54].Our model for WDM situations is simpler, because wedo not distinguish in the gas a Knudsen layer vs a hy-drodynamic layer. We assume that our computing cellsare small enough that the gas density inside them is con-stant. The variations over the whole flow shall insteadbe treated by the global hydrodynamic code that deter-mines the expansion of each kinetic cell. The kinetic andphase change processes in our description are driven bythe hydrodynamic expansion, therefore, the validity con-dition of our model is that the initial cell size should bemuch smaller than the initial dimensions of the expand-ing material : L (cid:28) δr, δz, (15)where δz is the initial foil thickness and δr the heatingbeam diameter. If Eq. 15 is verified, the global hydrody-namic treatment is correct, with gradients properly re-solved. Eq. 15 is verified in the standard situations weconsider, but can break down if the initial droplet or cellsize is not small enough compared to the sample size. C. Equilibrium condition between droplet and gas
In order to get a closed system of equations for particleand energy fluxes, we still need one assumption. Ourmodel is a priori out of thermal equilibrium ( T l (cid:54) = T g ), sothe density of the liquid is not determined yet. It seemsreasonable to assume pressure equilibrium between thedroplet and the gas [13], because we expect that a fewcollisions are sufficient for the droplet to ”experience”the gas pressure, and adjusting the liquid pressure to itrequires only a small density change, due to the very lowcompressibility of the liquid.Due to the droplet curvature, the pressure equilibriumcondition is the Laplace equation P l − P g = 2 σR . (16)An exact numerical implementation of Eq. 16 is difficult,because it requires to solve a non-linear system at eachtimestep in order to determine the liquid density given acertain set of values { V, N l , N g , E l , E g } .To simplify the condition, one can approximate theliquid density by the equilibrium value n ∗ l ( ∞ ) for a flatinterface ( R → ∞ ), but this is wrong for two reasons :first, because due to the fast expansion, the gas pressureis lower than the saturation value corresponding to theliquid temperature, and second, because of the Laplacecompression term of Eq. 16. Running our model, we havechecked that this raw approximation leads to importantinaccuracies in the calculation of the pressure, especiallyat low temperatures where the Laplace term becomesdominant. Still, these errors do not cause important dis-crepancies in the global description of the droplet evolu-tion, which we attribute again to the low compressibilityof the liquid.For more accuracy, we choose to compute the liquiddensity in perturbation from the flat equilibrium value : n l = n ∗ l ( ∞ ) (cid:18) PK l ( T l ) (cid:19) , (17)where K l ( T l ) = n l ( ∂P/∂n l ) T l is the isothermal bulkmodulus of the liquid that we compute directly from theVdW EoS, and ∆ P = 2 σ/R − ( P l ( n ∗ l ( ∞ ) , T l ) − P g ) is thepressure correction that we compute using Eq. 16. As weshow in the next section, this perturbative approach ofthe pressure equilibrium condition is very satisfactory.With Eq. 9-17, we have a complete kinetic model forthe evolution of a droplet and its vapor in a cell. In thefollowing, we combine it with the VdW EoS to study thedifferent regimes of two-phase expansion. IV. RESULTSA. NDCX II reference case (subcritical case)
The reference case envisioned as an upcoming experi-ment on the NDCX II machine at LBNL consists of heat-ing an aluminum foil of thickness δz = 3 . µ m with a shortpulse of ions of duration t heating (cid:46) δr = 1mm.Initial temperatures up to 1eV are predicted for the ex-pected beam fluences [24].In Fig. 5, we present the numerical output of the modelfor a cell containing a droplet and gas initially at equi-librium at T = 8000K with V l = V g = V /
2, becausethis corresponds to the onset of the ”droplets regime”(see Sec. II.B). As we mentioned previously, we makethe crude assumption of a flow with linear speed profileand outward expanding speed v = 3 c s (cid:39) . z > z <
0, where the initial soundspeed c s is estimated roughly as the thermal velocity v th ( T ) = (cid:112) k B T /m . Then, the strain rates in Eq. 9are simply η z = 6 c s /δz and η r = 6 c s /δr . We dis-play a full 3D case (solid lines) and a 1D case (dashedlines) where η r = 0. Here the hydrodynamic time is t hydro = 1 /η z (cid:39) . t D = 1 /η r (cid:39) α = β = 0 . u = ln( t ) tospan a wide temporal range. The result is displayed upto t = 100 µ s because at this time the front has traveled over about 50cm, which is comparable with the size ofan experiment. The initial droplet radius R = 25 . typical droplet . ! ns " N ! " ! ns " T ! K " ! ns " P ! k b a r " ! cm ! " T ! K " (b)(d)(a)(c) t t min t hydro t t hydro N l N g T l T g P l P g g l FIG. 5. (Color online) Droplet and gas evolution in the NDCXII reference case. Initially the droplet of radius R = 25 . T = 8000K. All liquiwd (gas) quantities are labeled with l ( g ). Plotted are the time evolution of (a) the particle numbersand (b) the temperatures, for a 1D (dashed lines) or a full3D expansion (solid lines). (c) Pressure evolution in the 3Dcase. The pressure difference computed (dashed lines) andexpected from Eq. 16 (dotted lines) are indistinguishable. (d)Trajectories in the phase diagram for the 1D (dashed lines)and 3D (solid lines) cases At early times ( t < t min , after which the liquid fraction starts growing slowly. Then, in the purely 1D case (dashed lines), thedroplet continues to grow steadily. In the more realis-tic situation however (solid lines), the droplet evaporatesagain when the 3D regime sets in, at times t > t D .In Fig. 5.b, one sees that, almost instantaneously afterheating ( t < T = T l − T g is established between the gas and the droplet, andremains roughly constant throughout the expansion inthe 1D case. On the contrary, in the 3D case, T g dropsquickly to almost 0 around t D , whereas T l decreasesslowly to a value around 1600K.In the phase diagram trajectories [Fig. 5.d], we see thatin both cases the liquid density remains very close tothe equilibrium value. By contrast, the gas density isclearly below the binodal in the 1D case, and in the 3Dcase it dives deep into non-equilibrium (supersaturated)conditions.In Fig. 5.c, we check the pressure equilibrium conditionin the 3D expansion case. One cannot distinguish thepressure difference in the computed evolution (dashedlines) that uses Eq. 17 from the theoretical value of Eq. 16(dotted lines), as the agreement is better than 2% overthe whole simulation range. The increase of P l at longtimes is due to the Laplace term (Eq. 16).Clearly, from the NDCX II example, two differentregimes can be identified. The first one, where the tem-perature difference is small, and remains constant, isa quasi-thermalized regime . In this regime the dropletgrows. The second one, where the gas becomes muchcolder than the drop, is a non-thermalized regime . Inthis regime the droplet evaporates again, as if it were invacuum. We now discuss the various regimes. B. Thermalization condition, quasi-thermalizedregime
Let us find a thermalization condition. In our equa-tions, the energy is extracted from the system only bythe adiabatic expansion of the gas (Eq. 13) and the gasquenching is then transmitted to the liquid via the liquid-gas energy exchange term (Eq. 14). Therefore, we shouldcompare those two energy fluxes to find the thermaliza-tion condition.Let us assume a small temperature difference ∆
T /T (cid:28)
1, so that we are in the quasi-thermalized regime of ex-pansion, as in the 1D case of Fig. 5. From Fig. 5.a, onesees that N l and N g are almost stationary if T l (cid:39) T g .Hence, let us make the approximation Φ vap = Φ cond (more precisely, | Φ vap − Φ cond | (cid:28) Φ cond ).The ratio between the two energy fluxes can be esti-mated as follows. Let us consider a cell containing fixednumbers N l and N g of liquid and vapor atoms, respec-tively, and let x = N l /N tot be the liquid fraction in thecell, where N tot = N l + N g . In the low T limit, a smallenergy change can be written dE g = N g k B dT g for thegas and dE l = N g k B dT l for the liquid, according toEq. 7. Noting dE tot = dE l + dE g the total energy lost bythe cell, we define ξ = dE l /dE tot . Requiring stationary∆ T (i.e. dT l = dT g ), we obtain ξ = 5 x/ (2 x + 3).In the low T limit, the adiabatic cooling of the gasimplies: dE tot /dt = − P g dV /dt = − n g k B T g ˜ ηV , wherewe define the volume strain rate ˜ η = d ( V /V ) /dt . Notethat, in the 1D expansion regime, ˜ η = η z , whereas inthe 3D expansion regime, for t (cid:29) t D , ˜ η (cid:39) η z η r t anddiverges. The power transferred from the liquid to thegas is: dE l /dt = n g (cid:112) k B T g / πmSχc v ∆ T , where χ = β + (1 − β ) α . To get this expression we have computedthe contributions from the three terms in Eq. 14 andused Φ vap = Φ cond . Expressing S = 4 πR and V =2 × πR , the balance between the fluxes dE l = ξ dE tot finally yields ∆ TT = ξ ˜ ηχ √ π R v th ( T g ) R (18)where we note v th ( T g ) = (cid:112) k B T g /m the thermal speed in the gas. With the approximation R (cid:39) R and neglectingthe prefactors of order unity, we obtain the simple scalingthat is expected to be valid for any EoS :∆ TT ∼ ˜ η R v th ( T g ) . (19) !!!! """" ! ! T " T theo " (a) (b) ! ns " t t t min FIG. 6. (Color online) Test of the thermalization for-mula Eq. 18. (a) Time evolution of the ratio Λ =(∆
T /T ) theo / (∆ T /T ) of the temperature difference computedwith Eq. 18 to the numerical model in the NDCXII refer-ence case, in purely 1D (dashed line) and full 3D (solid line)expansion. (b) Λ in the 1D expansion at t = 1ns (solid sym-bols) and at t = 10 µ s (open symbols) with β =0, 0.2, 0.4,1 (circles), α =0, 0.2, 0.4, 1 (up triangles), δz =1, 2, 4, 8 µ m(squares) and R =5, 10, 20, 40nm (down triangles). In Figure 6 we check the validity of Eq. 18, comput-ing the ratio Λ = (∆
T /T ) theo / (∆ T /T ) of the theoreticaltemperature difference (Eq. 18) to the result of the fullnumerical calculation. To evaluate Eq. 18 we take thevalues of ˜ η , R , T g and x from the result of the numeri-cal simulation. In Fig. 6.a, we show the evolution of theratio Λ in the NDCXII reference case (same calculationas Fig. 5). In the full 3D expansion case (solid line),the prediction becomes bad (error larger than 100%) at t (cid:39) t D , which is expected since the volume expansionrate ˜ η diverges in 3D. In the purely 1D expansion, afterthe first ns, one sees that Eq. 18 is accurate within 20%.More precisely, the error has the same sign as the deriva-tive dN l /dt and vanishes when the droplet is stationary,at t = t min , which is expected since Eq. 18 is obtainedwith the assumption dN l /dt = 0.In Fig. 6.b, we show the ratio Λ at t = 1ns and at t = 1 µ s for the same parameters, but varying one by onethose that are relevant to Eq. 18: β , α , δz (in order tovary η ) and R . The analytic formula overestimates (un-derestimates) ∆ T /T in all cases at t ( t ), and for bothcases the error is larger when the expected (∆ T /T ) theo is larger, which is natural since Eq. 18 holds in the limitof small ∆ T /T . Interestingly, the point β = 0 is sepa-rate from the others at both times, and is closer to 1,which is not surprizing since β = 0 means no particle ex-change, and this again confirms that the main source ofinaccuracy of Eq. 18 is a non-zero value of dN l /dt . Theprediction could thus be refined if this effect was takeninto account, for example using the analytical results ofthe next section. At early times, one can see also thatthe error is larger for droplets of radii smaller than 5nm.This effect is switched off if we set σ = 0, confirming thatit is caused by surface effects.It is also possible to rewrite Eq. 19 in terms of the ini-tial conditions: considering η ∼ c s /δz , one gets the verysimple scaling law: ∆ T /T ∼ R /δz This last expressionis only valid in the case of a 1D linear expansion where ˜ η is constant. In this case, one sees using Eq. 1 that ∆ T /T is expected to decrease slowly when the initial samplesize increases R ∝ η − ⇒ ∆ TT ∝ η ∝ δz − . (20)For larger samples, the thermalization will be better eventhough the droplets are bigger. This justifies again thatin the limit of large samples and slow expansions, anequilibrium hydrodynamic description becomes valid.In summary, Eq. 18 is expected to be always a goodestimate in the quasi-thermalized regime, and Eq. 19 canbe considered as a universal criterion to delimit the quasi-thermalized regime. C. Fully thermalized regime
In the previous section we could distinguish theregimes of quasi-thermalized versus non-thermalized ex-pansion. From Eq. 18 it is clear that the quasi-thermalized regime will become quickly invalid after t D ,because the volume strain rate ˜ η diverges. Nonetheless,in the early times of expansion, or if one is interested insystems of large radial extent, it is worth studying thelimiting case of a fully thermalized flow.In this perspective, let us assume T l = T g = T . Again,we look at the low T regime, which becomes valid veryearly in the expansion process. Using the first order ap-proximations in the VdW model n l (cid:39) (1 − θ ) /b and n g (cid:39) E tot = N tot ( c v T − x (1 − θ ) a/b ),where x = N l /N tot is still the liquid fraction in the cell.The total energy change dE tot = − P g dV becomes, atfirst order in θ : − (1 − x ) θ dVV = ( + x ) dθ − dx . Notingthat dV (cid:39) dV g , we convert θ d (ln( V )) = dθ − dθ/θ usingthe low T approximation Eq. 6, and find finally dx = (cid:18) − − xθ (cid:19) dθ. (21)It is easy to push the approximation to higher orders in θ ,but Eq. 21 already allows one to get good insight into theevolution of the droplet. One sees that the droplet willbe stationary at a temperature satisfying θ min = (1 − x ) corresponding to the time t min already mentioned, itwill evaporate before this point, for temperatures θ >θ min , and grow after it, for θ < θ min . This sequence isin agreement with the NDCX II reference case shown in fig. 5. Note that, at long times, and independentlyfrom the EoS, droplets will always grow in a thermalizedsituation. This is due to the fact that adiabatic expansionof a perfect gas is an algebraic trajectory in phase-space( T ∝ ρ / ), whereas Clausius-Clapeyron law predicts anexponential curve for the gas binodal, meaning that thegas in a two-phase expanding cell will always tend tosaturate and make the liquid fraction grow. (a) (b) ! V x ! V T " K T =7000K8000K9000K10000K (c) " K x FIG. 7. (Color online) Thermalized evolution for a dropletof initial radius R = 25nm and initial temperatures T =7000 K to 10000K (solid lines). Temperature (a) and liquidfraction (b) evolution versus the cell expansion. (c) Full nu-merical calculation compared to solution of Eq. 21 (dashedlines) started at V /V = 10. In Figure 7, we show the exact numerical calculationof the thermalized evolution of a droplet whose initial ra-dius is R = 25nm and for initial temperatures rangingfrom 7000 to 10000K. In the thermalized case, the timeevolution is irrelevant, that is why in fig. 7.a and b thetemperature and liquid fraction are plotted as a func-tions of the volume expansion ratio V /V . An expansionremaining in the thermalized regime over 10 orders ofmagnitude volume expansion is unrealistic in the case ofNDCX II, but for generality it is interesting to study thislimiting trend.In Fig. 7.c, the numerical liquid fraction versus temper-ature is compared to the solution of Eq. 21, starting whenthe volume has expanded by one order of magnitude.The analytic approximation is accurate within 20%. Thisshows that Eq. 21 can be used to make good estimates ofthe asymptotic growth of the liquid fraction in the ther-malized case, which can be, e.g. for T = 9000K, of a bitmore than 10%.This growth of the liquid fraction in the thermalizedregime is a rigorous upper bound for droplet growth. Inthe opposite regime, one can get the reciprocal upperbound for droplet evaporation.0 D. Non-thermalized regime: evaporation in vacuum
If the gas expands too fast for thermalization to oc-cur, one expects the droplet to evaporate as if it were invacuum. The corresponding limit consists of assumingΦ cond = 0. Obviously, in this case the droplet can onlylose particles. Moreover, the evaporation of the dropletis maximal in this regime, because, if there was thermal-ization via collisions with a gas, colder than the drop,it would be a way for the droplet to lose energy with-out losing particles. Also, because the vapor pressuredecreases exponentially with temperature, one expectsthe evaporation in vacuum to slow down fast. However,independently from the kinetics, it is clear that theremust be some upper bound to the evaporation of a drop.Indeed, as every evaporating particle takes away energyfrom the droplet (the latent heat), the droplet gets colderand colder, until the evaporation is ”frozen”, a strict limitbeing that T cannot become negative.Let us find analytic expressions for the maximal evapo-ration of a droplet whose energy is noted E l . Consideringthe evaporating particles and using Eq. 14, the energyloss can be written dE l = e ∗ g dN l . On the other hand,considering the liquid, and neglecting the surface energyterm, one has dE l = e l dN l + N l ( ∂e l /∂T l ) dT l . Equatingthose two expressions, one finds, for the VdW fluid : dN l N l = ( c v − a ( ∂n l /∂T l )) a ( n l − n g ) dT l . (22)Using the development of n l at first order in θ (Eq. 5),one can integrate Eq. 22 from initial T and N l , yieldingat the final T l the remaining fraction N l N l = exp (cid:20) −
52 ( θ − θ l ) −
134 ( θ − θ l ) (cid:21) . (23) (a) (b) ! ns " T l ! K " ! ns " N l N l T =3000K5000K7000K9000K FIG. 8. (Color online) Evaporation in vacuum. Time evolu-tion of the temperature (a) and the non-evaporated fraction(b) for a droplet of initial radius R = 25nm and initial tem-peratures T = 3000 to 9000K (solid lines). The dashed lineson (b) are the numerical integration of Eq. 22 and the dottedlines correspond to Eq. 23 . Figure 8 shows the exact numerical result for the evap-oration in vacuum of a droplet whose initial radius is25nm, and initial temperatures ranging from 3000 to 9000K. The volume expansion is not relevant here, sothe variables are displayed as a function of t only. InFig. 8.a and 8.b, one sees that for all initial tempera-tures, the number evaporation and cooling curves of thedroplet follow a same asymptotic behavior, which is in-creasingly slow at long times. In Fig. 8.b, the numericalintegration of Eq. 22 is shown for each T (dashed lines),with the final T l taken from the full numerical solution.The agreement with the final evaporation ratio is excel-lent, showing that surface effects play a negligible role.We have checked that surface effects cause an overestima-tion of the maximal evaporation of less than 10% even fordroplets of initial radius 1nm. The approximated Eq. 23is also displayed for each T (dotted lines), taking herealso the final T l from the full numerical run. One seesthat it predicts the good limit for evaporation within 5%for initial temperatures up to 7000K. This is very satis-factory because the non-thermalized regime is expectedto be valid only in the late times of expansion so the firstorder low T approximations should be very accurate.Let us now discuss the onset of the radiative coolingregime. At long times and low temperatures, the par-ticle evaporation and the evaporative cooling rates de-crease exponentially (Eq. 6), whereas the radiative cool-ing rate is algebraic ( ∝ T ). Therefore thermal radia-tion becomes the dominant cooling mechanism at longtimes. Within our model it is not difficult to express thetemperature T rad below which radiative cooling becomesdominant over evaporative cooling [55]. Using our VdWparameters for Al, and assuming an emissivity (cid:15) = 0 . T rad (cid:39) T rad (cid:39) − µs and for temperatures below 2000K, although the evapo-rative cooling rates that we obtain within our model aresignificantly larger than the estimates they report. Asan example, for Si at 2000K, our crude model predicts aradiative cooling rate of (cid:39) µ s, while the evaporativecoling rate in vacuum is still (cid:39) µ s. Note howeverthat the rate we compute is a strict upper bound. E. Supercritical case: nucleation and growth ofliquid droplets
The last case to consider is the supercritical case, wherethe material expands first as a supercritical fluid, andenters the two-phase region of the phase diagram crossingthe gas binodal, as a supersaturated gas. In this casenucleation of droplets may occur.We do not propose a model for nucleation but we notethat it has been reported [13] that supersaturation of thevapor does not reach high values and that above a cer-tain threshold value nucleation is very sudden, due to1the exponential dependence of the nucleation frequencyon the supersaturation ratio [35]. Then, our model isexpected to describe correctly the subsequent evolutionof the clusters. In particular, we expect that thermaliza-tion will depend on the droplet radius and the volumestrain rate and that Eq. 18 will be a good estimate inthe quasi-thermalized regime. If clusters and gas are inthermal equilibrium, we expect Eq. 21 to be valid as well.Note that nucleation may also happen in a subcriticalexpansion scenario, if the gas becomes very supersatu-rated at long times, as can be seen for the reference casein Fig. 5.d. In this situation, close to the non-thermalizedlimit, condensation on the existing droplets is too slow,and nucleation of new droplets may happen even withliquid clusters already present in the plume.
V. CONCLUSION
In conclusion, we have studied droplet evolution andthermalization conditions with an original, simple kineticmodel based on a consistent set of rate equations for massand energy exchanges. Using the vdW EoS as a test-bed,we have demonstrated that such a kinetic treatment isable to bridge the gap between the molecular and equi-librium hydrodynamic approaches that have mainly beenused so far. Most of the results of this work are generaland should be extendable to any EoS with which the ki-netic equations are used.In particular, the main output of our study is to iden-tify the different regimes of two-phase expansion. On oneside, the quasi-thermalized case and its limit, the fully-thermalized case, on the other side, the non-thermalizedcase. To distinguish the two situations, we identify a lo-cal thermalization condition (Eq. 18) which depends onthe droplet radius R , the volume expansion rate ˜ η , thegas temperature T g , the liquid fraction x , and the ki-netic parameters α and β , but it can also be traced backto the initial conditions: sample thickness δz and initialtemperature T (Eq. 20). Eq. 19 is a simpler alternativeto Eq. 18 that involves only the initial droplet radius R ,˜ η and T g and gives the scaling expected for any EoS.Due to the crossover from 1D to 3D expansion at thetime t D , the expansion is expected to take place in thequasi-thermalized regime at early times, at least in theNDCXII reference case and for similar parameters, butat long times the non-thermalized regime is almost in-evitable. Eq. 18 shows that this is only a dimensionaleffect, driven by the divergence of ˜ η in 3D.In the quasi-thermalized case, our study shows thatthe relative temperature difference ( T l − T g ) /T l remainsalmost constant throughout the expansion. Eq. 18 is de-rived assuming no net particle exchange, so only kineticenergy terms (but no latent heat) are involved, whichmakes it suitable for generalization to other EoS. Thepredictions of Eq. 18, and, more generally, the validity of the whole droplet evolution model, could also be testedwith MD calculations and experimental measurements.In the fully thermalized case, droplets can grow (mod-erately) at long times and we give an approximate for-mula for the VdW fluid (Eq. 21). In the opposite caseof a fully non-thermalized flow, one can find a strict up-per bound for the evaporated fraction at a given final T l .For the VdW fluid, this limit is Eq. 22. In both limitingcases, simple implementations of the kinetic model, e.g.,with the VdW EoS, can be used to make estimates andprovide upper bounds for the droplets evolution. Notethat, in all the cases we have studied, the droplets nevergrow or evaporate very much from their initial situation.For the moment, the model we have presented is lo-cal, but in the future it could become part of a largerhydrodynamic code that will treat many lagrangian two-phase cells containing droplets and gas. The extensionof our model to several droplets in one cell can be doneeasily. For more realistic simulations, it will probablybe necessary to take into account other effects that arenot two-phase phenomena and that we have left aside,such as radiation, which can be non-blackbody in theearly stages of expansion, but also, thermal conductionbetween cells, and thermionic emission.Before our kinetic model can be used to computedroplets evolutions at the core of a global comprehen-sive code, additional modules are needed to initialize thetwo-phase regime. In the subcritical case, a hydrody-namic code and some model for fragmentation is requiredto determine the droplets distribution at each location,and possibly their velocities, which may differ from thatof the expanding gas. In the supercritical case, a modelfor nucleation is needed. In any situation, the thermal-ization condition (Eq. 18 or 19) can be used to checkwether the two-phase computing cell can be treated asan equilibrium cell or if a non-equilibrium treatment isrequired. The reason being of course that an equilibrium(thermalized) description is much easier to implement.Finaly, we have investigated the role of surface effectsin different cases. Surface tension is expected to play animportant role for droplets of radii R < R σ = σ/k B T n l .This can be seen from Kelvin equation or considering theradius at which the surface energy becomes comparableto the kinetic energy per particle in the liquid. Withour VdW parameters for aluminum, R σ increases from0 . T = 10000K to 64nm at T = 2000K. Surfaceeffects are thus increasingly important in the late stagesof expansion, at low temperature and for the smallestfragments. This is also why a careful treatment of thesupercritical case of in-flight nucleation is more difficultand remains to be done in order to complement this work.2 AKNOWLEDGEMENTS
We wish to thank R. M. More, F. M. Bieniosek, P.A. Seidl, B. G. Logan and I. D. Kaganovich for helpfuldiscussions. One of us (JA) was partially supported byEcole Normale Sup´erieure, France. Work performed un-der the auspices of the U.S. Department of Energy underUniversity of California contract DE-AC02-05CH11231at Lawrence Berkeley National Laboratory and contractDE-AC52-07NA27344 at Lawrence Livermore NationalLaboratory. ∗ Corresponding author: [email protected][1] R. Lee et al. , J. Opt. Soc. Am. B , 770 (2003).[2] S. Eliezer et al. , Phys. Rev. B , 144119 (2004).[3] O. Albert et al. , Appl. Phys. A , 319 (2003).[4] S. Amoruso et al. , Phys. Rev. B , 033406 (2005).[5] N. A. Tahir et al. , Phys. Rev. Lett. , 035001 (2005).[6] N. Tahir et al. , Nucl. Instr. Meth. B , 85 (2006).[7] P. Seidl et al. , Nucl. Instr. Meth. A , 215 (2007).[8] K. Sokolowski-Tinten et al. , Phys. Rev. Lett. , 224(1998).[9] A. Lindenberg et al. , Phys. Rev. Lett. , 135502(2008).[10] N. Zhang et al. , Phys. Rev. Lett. , 167602 (2007).[11] K. Oguri, Y. Okano, T. Nishikawa, and H. Nakano, Phys.Rev. Lett. , 165003 (2007).[12] S. Amoruso et al. , Europhys. Lett. , 404 (2004).[13] E. Lescoute et al. , Phys. Plasmas , 063507 (2008).[14] B. L. Holian and D. E. Grady, Phys. Rev. Lett. , 1355(1988).[15] W. Ashurst and B. Holian, Phys. Rev. E , 6742 (1999).[16] S. Toxvaerd, Phys. Rev. E , 704 (1998).[17] D. Perez and L. Lewis, Phys. Rev. Lett. , 255504(2002).[18] D. Perez and L. Lewis, Phys. Rev. B , 184102 (2003).[19] P. Lorazo, L. Lewis, and M. Meunier, Phys. Rev. B ,134108 (2006).[20] A. Upadhyay and H. Urbassek, Phys. Rev. B , 035421(2006).[21] Typically, when simulating particles of mass m interact-ing via a Lennard-Jones potential of range σ and depth (cid:15) , the natural molecular time scale is τ = ( mσ /(cid:15) ) / ∼ et al. , Phys. Rev. Lett.
86, 12 , 2573 (2001).[23] R. More, H. Yoneda, and H. Morikami, J. Quant. Spec-tro. & Radia. Transfer , 409 (2006).[24] J. J. Barnard et al. , Nucl. Instr. Meth. A , 275 (2007).[25] B. Chimier and V. T. Tikhonchuk, Phys. Rev. B ,184107 (2009).[26] D. Eder et al. , Nucl. Fusion , 709 (2004).[27] A. Koniges et al. , J. Phys. IV France , 587 (2006).[28] L. Landau and E. Lifschitz, Fluid Mechanics (PergamonPress, London, 1959).[29] S. Anisimov et al. , Appl. Phys. A , 617 (1999).[30] C. Cheng and X. Xu, Phys. Rev. B , 165415 (2005).[31] The classification that we use through this paper is notfully compatible with the one in [17–19], where the notionof “fragmentation” is used in supercritical conditions. Such use is absent in our perspective since the conceptof fragmentation of a liquid in several pieces separatedby gas is associated with the first-order liquid-gas phasetransition, which only exist for T < T c . There, the latentheat, and the surface tension, that allows spatial separa-tion of the two phases, are defined. In the works [17–19],dense vs dilute parts are distinguished from a continuoussupercritical medium with density fluctuations, but thecriterium used is not related with the phase transitionupon which our classification is based.[32] D. E. Grady, J. Appl. Phys. , 322 (1982).[33] E. L. Knuth and U. Henne, J. Chem. Phys. , 2664(1999).[34] M. Pilch and C. Erdman, Int. J. Multiphase flow. , 741(1987).[35] S. Balibar, J. Low. Temp. Phys. , 363 (2002).[36] M. P. Moody and P. Attard, Phys. Rev. Lett. , 056104(2003).[37] R. Tolman, J. Chem. Phys. , 333 (1949).[38] A. Onuki, Phys. Rev. Lett. , 054501 (2005).[39] A. Onuki, Phys. Rev. E , 036304 (2007).[40] F. Reif, Fundamentals of statistical and thermal physics (MacGraw-Hill, New York, 1965).[41]
CRC Handbook of Chemistry and Physics , edited by D.Lide (CRC Press, Boca Raton, FL, 1999).[42] We obtain a critical temperature T c = 10545K and acritical density ρ c = 0 . T c = 6250K and ρ c = 0 . Simple dense fluids , edited byH. Frisch and Z. Salsburg (Academic, New York, 1968).[44] V. Sarou-Kanian, F. Millot, and J. Rifflet, Int. J. Ther-mophys. , 277 (2003).[45] H. Hertz, Ann. Phys. , 177 (1882).[46] M. Knudsen, Ann. Phys. , 697 (1915).[47] C. Cercignani, Rarefied gas dynamics (Cambridge Uni-versity Press, UK, 2000).[48] J. Powles, J. Phys. A , 1551 (1985).[49] G. Nagayama, J. Chem. Phys. , 1392 (2003).[50] Using Richardson law, one can estimate that the rate ofelectron emission may be high at temperatures close to T c in the first ns of expansion. In the case where the elec-trons just fly away with no collisions with other droplets,the energy loss for the cell is expected to be smaller thanthe total work term in the early times and smaller thanthe radiative term at long times. In addition the build-upof a charge of the droplet should quickly limit the emis-sion. In the opposite case where the emitted electronswould be recaptured by other droplets, the thermionicemission and recapture terms should only amount to anincrease of the thermalization rate in Eq. 14, withoutcausing a fundamental change to the picture.[51] M. Bond and H. Struchtrup, Phys. Rev. E , 061605(2004).[52] A. Kryukov, V. Levashov, and S. Sazhin, Int. J. HeatMass Transfer , 2541 (2004).[53] S. Sazhin et al. , Int. J. Heat Mass Transfer , 2675(2007).[54] N. Fuchs, Evaporation and droplet growth in gaseous me-dia (Pergamon Press, London, 1959).[55] The evaporative energy flux per unit area is well ap-proximated in the low T limit by: Φ evapE (cid:39) β Φ vap l (cid:39) ( β/bθ ) e − θ ( k B T / πm ) / l , while the radiative energyflux is simply: Φ radE = (cid:15)σ S T where (cid:15) is the emissiv-ity of the liquid and σ S the Stefan constant. Finally thetwo processes are equally important at a temperatureverifying θ − / e − /θ = ( b(cid:15)σ S /β ) l / k − / B (2 πm ) / . Toapply this formula to Si, we obtain the VdW parameters a Si = 10 . × − erg.cm an b Si = 1 . × − cm byfitting the VdW EoS to the experimental values of thelatent heat l ( T b ) = 359kJ/mol at the boiling temperature T b = 3530K and of the liquid density ρ l ( T m ) = 2 . T mm