Dynamics of heterogeneous clusters under intense laser fields
DDynamics of heterogeneous clustersunder intense laser fields
A dissertation submitted to the Technical University of Dresden,Faculty of Science, Department of Physicsin partial fulfillment of the requirements for the academic degree ofDoctor rerum naturalium (Dr. rer. nat.)by
Pierfrancesco Di Cintio
Thesis supervisor: Prof. Dr. Jan-Michael RostSecond Thesis supervisor: Prof. Dr. Ulf SaalmannDresden, 7 August 2014 supported by a r X i v : . [ phy s i c s . a t m - c l u s ] A ug ontents N − body methods for systems of charged particles . . . . . . . . . 855.2 The N − body codes used for the simulations . . . . . . . . . . . . . 935.3 Including quantum processes in classical dynamical simulations: aMonte Carlo approach . . . . . . . . . . . . . . . . . . . . . . . . . 985.4 The analysis of the end products . . . . . . . . . . . . . . . . . . . 102 Contents
B Dynamical friction 111C Some remarks on kinetic theory 115Bibliography 119 hapter 1
Introduction
Understanding the response of matter when exposed to intense radiation is rele-vant to many areas of modern research, such as for example particle acceleration[110], astrophysics of active galactic nuclei [132] and biomolecular imaging withx-rays pulses [16]. The latter has been considered in recent years, as a promis-ing application for the modern x-ray free electron lasers (XFEL) like the LinearCoherent Light Source (LCLS) at Stanford [167], The Japanese XFEL SACLAat Kouto [168], or the European XFEL currently under construction in Hamburg[169].Since the pioneering paper by Neutze et al. [126], where single-molecule diffrac-tive imaging was suggested for the first time, many theoretical and experimentalstudies have been aimed at finding an optimal parameter range for the imaging oflarge biological molecules using x-ray pulses (see e.g. Refs [31], [18]). In a nutshell,as sketched in Fig. 1.1, a jet of “different copies” of the molecule to be imagedis irradiated by a pulsated x-ray laser field and the diffraction pattern resultingfrom the interaction of the radiation with matter is collected for every illuminatedsample [85]. Due to the fact that every molecule of the jet comes with a differentorientation with respect to the laser’s direction of propagation, a complicated re-construction procedure is then needed to extract the three dimensional structureof the molecule from many different two dimensional patterns with unknown ori-entation, see [71] and [165].However, this method is additionally plagued by the fact that at photon en-ergies of a few keV (soft x-rays), the cross section for K − shell photoionization ofcarbon and oxygen (that are the most common elements in biological moleculesthat one is interested in imaging) are considerably larger -roughly one order ofmagnitude- than the electronic elastic scattering cross sections. Due to this, thesample is likely to be destroyed by radiation damage within the duration of thepulse [72]. Obtaining pulse lengths of a few femtoseconds (1fs = 10 − s) and keep-ing large intensities (of the order of 10 W / cm ) in order to have enough photonshitting the targets has encountered several theoretical and technical difficultiesthrough the years. Modern x-ray laser sources, such as the aforementioned LCLS Chapter 1. IntroductionFigure 1.1: The steps of x-ray biomolecular imaging are illustrated.[78], are now able to reach pulse lengths of the order of one femtosecond for photonenergies of up to 10-12 keV (see e.g. [21], [23], [103]) reaching peak brightness upto 10 photons per s mrad mm , as shown in Fig. 1.2. Therefore, such machinesalso represent perfect candidates to attempt single-molecule imaging, as well asto probe the properties of matter under extreme conditions of irradiation [125].As optimal “test systems” with respect to the molecular imaging, clusters (i.e.droplets of atoms or molecules containing from a few up to ∼ particles andwith typical number densities of 10 − ˚A − ) have been suggested. They have par-ticularly simple structure and can be tailored in size from 1 up to 10 nm to matchthat of more complex organic molecules that one currently aims at imaging.It must be pointed out that clusters, when exposed to laser pulses with focus-ing of the order of few micrometers, feel the same laser intensity on their wholevolume contrary to what happens for larger solid-state targets. Thus clusters canabsorb a larger fraction of the energy “pumped in” by the laser field with respectto solids with similar particle density. The interaction of strong lasers with clus-ters leads to different fragmentation products (as sketched in Fig. 1.3), such asenergetic electrons [148], ions [46], [49], [151], as well as photons [111], [48], [157],emitted in the decay of inner shell vacancies or due to bremsstrahlung .The strong x-ray pulses generated by contemporary machines, with their highintensities (up to 10 W / cm ) and lengths from 1 up to 100 femtoseconds, ef-ficiently depose large amounts of energy in the target. Once charged, due tomultiple almost simultaneous photoionzation events, the cluster becomes what wewill hereafter refer as nanoplasma (i.e. a plasma of ions and electrons with typical Figure 1.2: Peak brightness as function of the photon energy for different modernlight sources. Note the enormous differences in the brightness. Figure taken fromRef. [78].size of a few up to thousand nm, [63]). Note that, with respect to long wave lengthlasers (for instance infrared or VUV light), the charging of the target happens viadifferent processes when exposed to x-ray pulses. In particular, photons with en-ergies of the order of one keV ionize mainly the inner electronic shells of elementssuch as carbon, nitrogen and oxygen that are among the principal constituents ofbiological molecules. The photoelectrons released in this way have typical kineticenergies of the order of 500-700 eV allowing them to leave the charged clusteron a time scale of few femtoseconds, and for electrons absorbing 10 keV photons,such times are even of the order of attoseconds (1as = 10 − s). The inner shellvacancies, with life times comparable if not shorter than the pulse lengths con-sidered here, decay via Auger processes thus forming a secondary population of Chapter 1. Introduction
Intense Femtosecond X-ray pulse
Clusters
Multicharged ions Electrons Photons (UV, X-ray)
Figure 1.3: The different products of laser irradiation of clusters.electrons with kinetic energies of roughly 200 eV, implying that the absorption ofone photon results in the emission of two electrons. We stress the fact that thispicture is radically different from that which one has when longer wavelength areemployed. In that case, the cluster is charged by the laser removing the electronson time scales, typically of the order of one picosecond (1ps = 10 − s), and in themeantime ions have started moving under the combined action of their repulsiveCoulomb forces and the laser electric field. The fast and at the same time mas-sive charging, obtainable with the contemporary x-ray sources, was until a decadeago unreachable when employing that time’s conventional laser machines. Theseextreme ionization conditions (i.e. high ion charge states produced in short time)thus open the door on previously unexperimented regimes of non-neutral plasmas,deserving therefore the interest of the theorist.When photoelectrons (and the faster Auger electrons) leave the cluster havingkinetic energies K e larger than their potential energy in the cluster’s electrostaticpotential, the system is rapidly torn apart by mutual repulsive Coulomb forcesamong the ions. The latter have suffered initially little to no displacement due tothe pulse. This regime of expansion is called Coulomb explosion, and for a clusterof initial radius R and number density n , is obtained when K e | e ¯ q πn R / | > q is the average charge state of the ions and e the charge of the electron.Figure 1.4 sketches the formation of the cluster potential and the system’s expan-sion. At the beginning of the pulse ( t ) the first photoelectrons leave. Then, asthe overall charge increases the cluster potential deepens reaching its maximum“depth” ( t ). At this stage, some of the less energetic electrons, produced via sec- Figure 1.4: Upper row: schematic view of the explosion of a charged cluster withinthe laser pulse. Lower row: evolution of the cluster potential as the photoelectronsleave the system (red arrows). Figure taken from Ref. [71].ondary processes such as Auger decay of inner shell holes or ionization of valenceelectrons due to the clusters electrostatic field, are trapped. The cluster expansionhowever rapidly, makes the potential well shallower ( t ) letting them evaporate.Albeit being structurally simpler than biological molecules, clusters still presentseveral “degrees of freedom” that make their dynamics after strong ionization nontrivial [55], and therefore interesting to study. To mention only a few, in molecularclusters the different species of ions can be accelerated differently in the clusterpotential leading to dynamics with different time scales. Moreover, the kineticenergy spectrum of the ions is strongly influenced by the initial structure of thecluster and its shape. Systems starting with flattened or elongated geometries areexpected to behave qualitatively very differently. In addition, charge migrationdue to rapid electron motion after an almost homogeneous spatial charging, also,influences the energetics of the cluster fragments and it is therefore important tohave a clear picture also of the electronic component.In this thesis, with the focus of shedding some light on the points mentionedabove, in a regime of laser parameters relevant for molecular imaging with x-raypulses, we have carried out a quantitative study of the dynamics of clusters irradi-ated by femtosecond x-ray pulses modelling those produced by contemporary lasersources such as LCLS or the European X-FEL. Our study is based on classical N − body simulations for the dynamics of particles, coupled with rate equationsand Monte Carlo samplings to treat photoionization processes.To prepare the field, the pure Coulomb explosion (i.e. no electrons are con-sidered) has been treated for different systems. We started by reviewing the Chapter 1. Introductionsimple continuum model idealizing ionized clusters as uniformly charged sphericalsystems. Since clusters are in reality constituted by particles, we discussed thediscrepancy with the approach based on particles. As mentioned before, the initialshape of the cluster, or in general of the ionized target, determines the energies ofthe products of its explosion.To this purpose, we have studied the expansion of non-neutral cluster plasmaswith non spherical symmetry for different families of ellipsoidal system, both bymeans of a (semi-)analytic continuum model and numerical simulations. We havealso analyzed the effects of initial conditions characterized by non uniform charg-ing and non negligible ion temperature, as they are of some relevance with respectto regimes of laser irradiation where some ion motion is possible within the laserpulse. Finally, we have implemented a detailed model of laser cluster interactionincorporating Auger transitions and electron recombination and used it to studythe response of molecular hydride clusters (i.e. H O, NH , CH clusters) to fem-tosecond x-ray pulses. This was motivated by the puzzling experimental findingsfor methane clusters exposed to short and intense x-ray irradiation [86].The thesis is structured as follows: First of all, in Chapter 2, the physics of theCoulomb explosion is introduced and discussed for the case of spherical systemswith homogeneous and heterogeneous composition. In the latter case, particularattention is devoted to the effect of the charge to mass ratio on the dynamics ofthe explosion, the results here reported will be finally contained in a forthcomingpublication.In Chapter 3 we extend the discussion to (homogeneous) systems significantlydeparting from the spherical symmetry (i.e. ellipsoidal geometry), focusing on thestructure of the final energy spectrum. Part of the work contained here, appearsin publication [65].Chapter 4 is devoted to the theoretical study of molecular hydrides clustersirradiated by extreme XFEL pulses. Its content has been published in [86] and[45].The numerical methods used throughout this work are discussed in Chapter5. In Chapter 6 the main results of this study are summarized and the futureapplication and aims are discussed.Finally, the thesis is completed by three appendixes treating respectively thestructure and production of clusters, the effect of multiple binary collisions ona charge travelling through a background of particles, and Coulomb explosionmodelled with kinetic theory. hapter 2 Coulomb explosion of sphericalcluster plasmas
In this chapter we discuss the physics of the Coulomb explosion of small sphericallysymmetric targets (i.e. atomic or molecular clusters) irradiated by intense laserpulses. We start by reviewing the simple case of single-component homogeneouslydense systems in the non-relativistic regime, where exact analytical results areavailable. We compare these results to calculations, where the spherical clusteris composed of particles. The characteristic differences between the two casesare discussed in detail. We then treat the case of non uniform density profilesand the formation of shock shells. Finally, we discuss the case of systems withheterogeneous composition.
As outlined before in the introduction, when an initially neutral cluster is irra-diated by an intense laser pulse it gets charged and therefore starts to expand.According to the duration of the pulse, its intensity, photon energy and the struc-tural and atomic properties of the target itself, different types of expansion cantake place.When almost all the electrons stripped from the atoms in the cluster remaintrapped by the space charge of the overall cluster, we speak of quasi-neutral orhydrodynamical expansion (see Refs. [36], [120] and [147]). On the contrary,when all the electrons are taken away, leaving a non-neutral plasma (see e.g. [38],cfr. also Eq. 1.1) whose dynamics is governed only by the repulsive inter-particleCoulomb forces, one has instead the so called Coulomb explosion (see Refs. [94],[128]). Note that intermediate regimes are also possible, for instance, due to nonhomogeneous charging of the cluster (see e.g. [2]).In this thesis, we concentrate on the case of Coulomb explosion, since we dealwith x-ray pulses with intensities and photon energies for which the majority ofthe electrons escape the system. Chapter 2. Coulomb explosion of spherical cluster plasmasThe dynamics of an expanding non-neutral plasma can be modelled using dif-ferent approaches, depending on what quantity or observable one is interested in.First we present the main analytical results in the continuum model and then wediscuss the results of numerical calculations in a particle-based approach.
Continuum model
Let us consider an isolated single component cluster that has been stripped of allits electrons composed by N ions of the same charge and mass q and m . In theso called continuum approximation the cluster is replaced with a smooth numberdensity n ( r ), so that its radial mass and charge densities are given by ρ m ( r ) = mn ( r ) and ρ c = qn ( r ), where m and q are the unit mass and charge respectivelyand r is the radial coordinate. We assume here no angular dependence. Therefore,the spherical symmetry of the system is preserved during the explosion.Under the assumption of incompressibility, one can treat this model with theequations of (non relativistic) fluid dynamics, see e.g. [95], namely the continuityand momentum equations, that are respectively ∂n ( r, t ) ∂t + 1 r (cid:2) r n ( r, t ) v ( r, t ) (cid:3) = 0 (2.1)and n ( r, t ) (cid:20) ∂v ( r, t ) ∂t + v ( r, t ) ∂v ( r, t ) ∂r (cid:21) = QM ∇ Φ( r, t ) , (2.2)where Q and M are the total mass and charge of the system, whose ratio equalsthe ratio q/m between unit mass and charge, and v ( r, t ) is the velocity field. Theelectrostatic potential Φ at position r is given byΦ( r ) = 4 πr (cid:90) r ρ c ( r (cid:48) ) r (cid:48) d r (cid:48) + 4 π (cid:90) + ∞ r ρ c ( r (cid:48) ) r (cid:48) d r (cid:48) . (2.3)Equation (2.3) is obtained from the radial Poisson equation∆Φ( r ) = 1 r ∂∂r (cid:18) r ∂ Φ( r ) ∂r (cid:19) = 4 πρ ( r ) . (2.4)Here ∆ is the radial part of the Laplace operator which we give explicitly. Notethat we are using here the atomic units for which the constant κ = 1 / π(cid:15) is setto 1. It must be pointed out that continuum approximation or continuum model is not the sameconcept as continuum limit (i.e. N → ∞ ; q i → N of the original system. For an extensive discussion see Refs. [87]and [88] and references therein .1. Mono-component systems The derivative of the potential Φ( r ) needed in Eq. (2.2) can be obtained fromEq. (2.3) and reads ∇ Φ( r ) = 4 πr (cid:90) r ρ c ( r (cid:48) ) r (cid:48) d r (cid:48) . (2.5)The system of partial differential equations (PDEs) given by Eqs. (2.1), (2.2)and (2.5) formally describes in closed form the Coulomb explosion of a mono-component cluster plasma and could be easily generalized to heterogeneous sys-tems, as well as to cases where an electron density is present (see Refs. [121] and[10]).We will now discuss a case which represents the continuum version of a homo-geneously charged cluster of total charge Q and mass M , initial radius R and withinitial kinetic energy K = 0. Therefore the initial charge density and velocitydistribution read ρ ( r,
0) = ρ n ( r,
0) = ρ θ ( R − r ) v ( r,
0) = 0 (2.6)where θ ( x ) is the Heaviside’s step function and ρ = 3 Q/ πR . It turns out thatin this case the full dynamics can be modelled simply by a second order ordinarydifferential equation (ODE) with a class of self similar solutions, instead of asystem of PDEs.For such initial conditions the electrostatic field at a generic r < R is givenby Eq. (2.5) and corresponds to the linear function of the radial coordinate r itself E ( r ) = ∇ Φ( r ) = 4 π ρ r . (2.7)In an infinitesimal time δt , an infinitesimal volume element of mass δm and charge δq (so that δq/δm = Q/M ) placed initially at radius r will reach the velocity δv ( r , δt ) = δqδm E ( r ) δt = δqδm π ρ r δt. (2.8)Due to the linearity with r , matter initially “sitting” at a given radius can not overtake matter initially placed at larger radii, always attaining larger velocities.The differential equation for the dynamics of such element of volume readsd r ( t )d t = δqδm π (cid:82) r ( t )0 ρ t r (cid:48) d r (cid:48) r ( t ) = δqδm Q ( r, t ) r ( t ) . (2.9)Since no overtaking is taking place, Q ( r, t ) = Q ( r ) = 4 πρ r /
3, and the equationabove can be rewritten as d r ( t )d t = Cr ( t ) , (2.10)where C = δq πρ r / δm , and its first integral reads (cid:20) d r ( t )d t (cid:21) = 2 C (cid:20) r − r ( t ) (cid:21) . (2.11) Chapter 2. Coulomb explosion of spherical cluster plasmas r / R ! t 00.511.52 0 0.2 0.4 0.6 0.8 1 n ( " * ) " * Figure 2.1: Trajectories of different volume elements initially placed at differentradii expressed in units of the initial cluster radius R . The red line marks theexpansion of the surface.By further integrating the latter equation, one gets the time t that takes toincrease from radius r to radius r ( t ) as t = (cid:114) ω − (cid:104)(cid:112) ξ ( ξ −
1) + ln (cid:16)(cid:112) ξ + (cid:112) ξ − (cid:17)(cid:105) , (2.12)where ξ = r ( t ) /r and ω = (cid:112) πρ δq/δm = Q (cid:112) /M R .When ω t (cid:28) r is given asymp-totically by r ( t ) (cid:39) r (cid:18) ω t (cid:19) , (2.13)while in the limit of ω t (cid:29)
1, one has r ( t ) (cid:39) (cid:114) r ω t. (2.14)Finally, the asymptotic velocity reads v ∞ = (cid:114) ω r . (2.15)In Fig. 2.1, the trajectories r ( t ) of different volume elements are shown. Theyshow a quadratic increase for early times cfr. Eq. (2.13), and become linear atlarge times, cfr. Eq. (2.14). It is clearly evident that they do not intersect.The absence of overtaking and the fact that every choice of r in the interval.1. Mono-component systems (0; R ) initial condition of Eq. (2.12) leads to identical (rescaled) dynamics, meansthat an initially cold homogeneously charged sphere expands retaining as spatiallyuniform density profile . In other words Eq. (2.12) represents a self-similar solu-tion.Knowing that a homogeneously charged sphere expands self-similarly it re-mains to determine its asymptotic number energy distribution n ( E ) = d P d E (2.16)defined as the fraction of the system with energy E . Whereas the total energy E tot is initially given by the potential energy U if the initial kinetic energy K equals0, in the limit t → ∞ it is E tot = K .For the homogeneously charged sphere considered here, one has E tot = U = 2 π (cid:90) + ∞ ρ ( r )Φ( r ) r d r = 35 Q R . (2.17)Using the fact that shells of matter starting from different radii in an uniformlycharged sphere do not overtake each other and the first Newton theorem (seee.g. [127], [81] and [146]), one can extract the asymptotic n ( E ), from the system’sinitial configuration as n ( E ) = d P d E = d P d r / d E d r . (2.18)Here we have defined the probability at t = 0 to find an infinitesimal element ofvolume at radius r as d P d r = 3 r R θ ( R − r ) . (2.19)The asymptotic kinetic energy (per unit charge) of a volume element is propor-tional to its initial radius through Eq. (2.15) and reads E ( r ) = Qr R , (2.20)therefore d E d r = 2 Qr R . (2.21)Substituting Eqs. (2.21) and (2.19) into Eq. (2.18) and expressing r = (cid:112) E R /Q ,which follows from Eq. (2.20), leads to a square root distribution n ( E ) = 32 (cid:115) E R Q θ ( E ( R ) − E ) = 32 (cid:115) EE θ ( E max − E ) , (2.22)where E max = Q/R , is the maximal energy per unit of charge, reached by the A uniformly charged shell of charge q , exerts at its exterior the same force as due to apoint-like particle of the same charge sitting at its centre, cfr. Ref. [91]. See also Eq. (2.5) Chapter 2. Coulomb explosion of spherical cluster plasmas n ( ε * ) ε * Figure 2.2: Asymptotic number energy distribution given by Eq. (2.22). Thescaled energy is defined as E ∗ = E / E max .fraction of the system initially sitting at its surface, see Fig. 2.2. Alternatively,the energy distribution can be also derived by substituting Eq. (2.20) into n ( E ) = d P d E = 3 R (cid:90) R δ ( E − E ( r )) r d r , (2.23)and performing the integration in r .The total (conserved) energy E tot is recovered by E tot = (cid:90) E max n ( E ) E d E . (2.24)Remarkably, the expression for the asymptotic n ( E ) given in Eq. (2.22) holds truealso in case of relativistic velocities, as it has been proved in [26].The simple case discussed here, albeit bearing a high level of abstraction, servesas a reference system for problems involving Coulomb explosion. Numerical simulations using particles
If one takes into account its particle nature, the system discussed in Sect. 2.1, andin general every mono-component plasma, can be described by the Hamiltonian H = N (cid:88) i =1 p i m + q N (cid:88) j (cid:54) = i =1 | r j − r i | (2.25).1. Mono-component systems n ( ε * ) ε * Analyticvariable δ fixed δ Figure 2.3: Normalized differential energy distributions at t = τ so that E = 99% U , for two simulated systems of N = 4 × particles starting withhomogeneous density profiles and K = 0. The red curve corresponds to the casewhere a minimum nearest neighbour distance δ is enforced in the initial conditionat every radius. The green curve instead, to the case of initially randomly dis-placed particles (i.e. an arbitrary small δ may occur). The two systems have thesame initial radius and total charge, and therefore the same averaged density . Ascomparison, the black line marks the theoretical expression given by Eq. (2.22).Energies are normalized in units of the final cutoff energy E max .where r i and p i are position and momentum of particle i and m and q its chargeand mass. Particle’s trajectories are obtained by integrating the system of 6 N Hamilton equationsd r i d t = ∂H∂ p i ; d p i d t = − ∂H∂ r i ; i = 1 , N. (2.26)It is possible to study the dynamics of a charged particles system, accounting forits particulate nature, only with the aid of N -body numerical simulations, alsoknown as molecular dynamics simulations (MD). Several numerical approaches doexist in order to compute the forces between particles and we redirect the readerto Chap. 5 for a description of the most widely used ones.Here we discuss the simulations of homogeneous systems performed with directforce calculations (see Sect. 5.1) and the origin of discrepancies with the contin-uum model.Figure 2.3 shows the number energy distribution n ( E ), when the initial po-tential energy U has been essentially completely converted into kinetic energy, Chapter 2. Coulomb explosion of spherical cluster plasmasfor two systems of N = 4 × identically charged particles homogeneously dis-tributed at rest in a spherical volume of initial radius R . The only differencebetween the two, is the way particles have been distributed in the initial condi-tion. In one case (red curve) a minimum inter-particle distance δ is fixed, whilein the other (green curve) the positions are randomly generated with no limit onthe minimum distance between nearest neighbours.One notices immediately that both curves reproduce the theoretical squareroot trend up to E ∗ ∼ .
9. For larger energies they are instead characterized bya sharp peak which makes the numerical curve departing considerably from itsanalytical counterpart (black curve in Fig. 2.3). For the system with arbitrarilyclose particle at t = 0, the peak is “milder” and broader.The origin of such feature in the case of a homogeneous, albeit made by parti-cles, density profile is not straightforward. It has been shown recently (see [146]),that due to its granular nature, a homogeneously charged sphere of radius madeby particles exerts on a test particle placed close to its surface a force that isconsiderably lower than that due to a continuous distribution with the same totalcharge.For an ion placed at position r inside the cluster, the probability to find an-other ion at radius r (cid:48) vanishes if | r − r (cid:48) | < δ , where δ is the radius of the so calledcorrelation hole only in one of the two cases shown in Fig. 2.3.It is possible to account analytically for the effect of a correlation hole in adistribution of charge. Let us first consider the electrostatic field at position r dueto an infinitesimally thin shell of radius r ∗ and charge q that is given by angularlyintegrating E r ∗ ( r ) = q π (cid:90) π (cid:90) π ∂∂ r | r − r ∗ | sin θ d θ d φ, (2.27)where r ∗ = r ∗ (sin θ cos φ, sin φ cos θ ) is a vector on the shell.Without loss of generality, due to the spherical symmetry of the problem, wetake into account only the the radial component of E r ∗ ( r ). Setting τ = cos θ andperforming the integration over φ , the latter reads E r ∗ ,δτ ( r ) = q (cid:90) +1 − δτ − dd r (cid:112) r + r ∗ − rr ∗ τ d τ = qr (cid:34)
12 + (1 − δτ ) r − r ∗ (cid:112) r ∗ + r − − δτ ) r ∗ r (cid:35) , (2.28)where restricting the upper boundary of integration to +1 − δτ accounts for thepresence of the correlation hole.Note that in the limit of δτ → q/r outside, cfr. Ref. [91].In Fig. 2.4 the field produced by a shell with a hole is shown for two values of δτ . Note how the field is non zero at the interior of the shell. It is expected that,.1. Mono-component systems E r * , δ τ (r) / ( q / r *2 ) r/r * pure Coulomb δτ =0.01 δτ =0.1 Figure 2.4: Electric field E ( r ) produced by a shell with a hole with δτ = 0 . δ placed closeto the surface (particle radius r in the sketch in Fig. 2.5), due to an extendeddistribution with homogeneous charge density. Setting δτ = ( δ − ( r − r ∗ ) ) / (2 r ∗ r )for the shells intersecting the hole (i.e. | r − r ∗ | ≤ δ ) and zero otherwise in Eq.(2.28), and integrating over r ∗ gives E R ,δ ( r ) = (cid:40) Qr/R ; for r < R − δA R ,δ ( r ) Qr/R ; for R − δ < r < R , (2.29)where the “attenuation factor” A R ,δ ( r ) is defined by A R ,δ ( r ) = ( r + R − δ ) (2( r + R ) δ + δ − r − R ) )16 δr . (2.30)As seen in Fig. 2.5 (bottom panel), the radial electric field E R ,δ ( r ), does notincrease linearly with r up to the surface, but starts abruptly to increase markedlysublinearly for r > R − δ .Since the asymptotic energy distribution is entirely determined by the initial In other words, the closer a particle is to the surface, the smaller is the compensation on theunderestimated radial electrostatic field generated by particles at lower radii, due to the spurious(with respect to the continuum picture) internal field of charged discrete shell placed at largerradii. Chapter 2. Coulomb explosion of spherical cluster plasmas with r ij the distance between the particles i and j and ~ g ð r Þ ¼ ð r=R $ Þ ð r=R þ Þ the distribution of distances r in a sphere with radius R [20]. The normalization with ~ g ð r Þ is necessary in order to remove the trivial dependenceon r due to the finite size of the system.We can assess how the correlation hole affects the forceacting on an ion at position r analytically. To this end,we determine the Coulomb repulsion of the ion from aninfinitesimally thin spherical shell of radius r and charge q by integrating over all angles f r ð r Þ ¼ q ! Z ! d " Z ! d sin @@ r j r $ r j ; (5)with r & r ð sin cos " ; sin sin " ; cos Þ a vector on theshell. The radial component of this force can be determinedwithout loss of generality by choosing r ¼ r ^ z along the z axis. Integrating over " and using $ & cos it reads f r ; %$ ð r Þ ¼ q Z þ $ %$ $ d $ ddr ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r þ r $ rr $ p ; (6a) ¼ qr " þ ð $ %$ Þ r $ r ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r þ r $ ð $ %$ Þ r r p : (6b)Hereby, the integration was restricted to the upper limit $ %$ < in order to account for a correlation hole, cf. the sketch in Fig. 4 and its inset, which shows %$ explicitly.Equation (5) corresponds to %$ ¼ . In this case, perform-ing the integration over $ yields—as expected—Gauss’law, i.e., f r ; ¼ for r
1, the heavy blue line marks the averagedradial electric field (cid:104) E ( r ) (cid:105) in the particle distribution, note how in both cases itbecomes markedly sublinear as r approaches the surface.the initial state. The asymptotic kinetic energy as a function of r reads E ( r ) = r (cid:90) + ∞ E ηR ,ηδ ( ηr )d η = rE R ,δ ( r ) . (2.32)By plugging the latter into Eq. (2.23), and having assumed a finite energy reso-lution δ E so that the Dirac δ ( x ) is replaced by D δ E ( x ) = exp(( − x/δ E ) ) / √ πδ E , (2.33)one obtains an asymptotic n ( E ) following the square root behavior for a broadinterval of energies and peaking close to the cutoff energy in a similar fashion ofthe numerical curves shown in Fig. 2.3.The high energy behavior of n ( E ), when a correlation hole of radius δ is present, Chapter 2. Coulomb explosion of spherical cluster plasmasis to be interpreted as the fact that the kinetic energies reached by the particlesinitially placed at a distance from the surface smaller than δ are lower than thoseattained by an element of volume in an idealized continuous distribution startingfrom the same radius. This causes the “bunching” of such energies that is seen asa peak in the differential energy distribution.We shall call this feature a discreteness peak to distinguish it from the similarlylooking feature of systems with non-uniform initial density profile that we willdiscuss in the next section.It must be stressed out that the discrete nature of the system implies that, thediscrepancy between E ( r ) computed in the continuum model and for particles,depends strongly on how the positions of the particles in the initial state areselected.In Fig. 2.3 the discreteness peak is sharper for the system where particles areinitially placed enforcing a minimal nearest neighbours distance (blue curve). Inthe case without a correlation hole, i.e. fully random positions, (green curve), thisfeature is broadened by the combined effect of a more randomized contribution toforce due to the near neighbours and position specific size of the correlation hole.However, the discreteness peak is not entirely removed since the local arrangementof the particles in proximity of the surface acts like a correlation hole.In Fig. 2.6 the radial component of the electrostatic field E ( r ) acting onthe particles is shown for the initial states relative to Fig. 2.3. In both cases,where a minimum inter-particle distance is or is not imposed, E ( r ) presents largedeviations from its theoretical value (indicated by the thin black line) due to therandomization of the contribution of the nearby particles. The averaged radialelectric field (cid:104) E ( r ) (cid:105) (blue lines in right panels of same figure) clearly drops closeto the surface with respect to the linear trend of E ( r ) predicted by the continuummodel, in the same fashion of Fig. 2.5 (lower panel). For the model with enforcedminimum inter-particle distance δ (red dots), (cid:104) E ( r ) (cid:105) departs more from the lineartrend although overall E ( r ) is less noisy than in the model with arbitrarily small δ (green dots). Nevertheless, in the latter case the noisiness of E ( r ) partiallycompensates the effect of discreteness related to the correlation hole which resultsalso in a less pronounced peak in the final n ( E ). Smoothed Coulomb interaction
Curiously, in direct N -body simulations, the modification of the pair Coulombpotential and force introduced in order to avoid its divergence for vanishing sep-aration (see Chap. 5), introduces an effect that lowers in a similar way the elec-trostatic field at the surface of an homogeneous sphere made by particles.If the Coulombian 1 /r potential is substituted with V ( r ) = 1 √ r + (cid:15) , (2.34)where (cid:15) is the so called softening length, the electrostatic field produced by a.1. Mono-component systems E ss o f t (r) / ( q / r *2 ) Analyticpure Coulomb ε =r * /100 ε =r * /50 ε =r * /2000.20.40.60.81 0 0.5 1 1.5 2 E ss o f t (r) / ( q / r *2 ) r/r * Numerical ε =r * /100 ε =r * /50 ε =r * /20 Figure 2.7: Top panel: for different values of (cid:15) in units of r ∗ , softened field E softs ( r )at radius r evaluated analytically for a shell of density σ . Bottom panel: same asabove but computed for a discrete distribution of N = 10 particles placed on ashell, computed with direct summation and averaged over 20 realizations.point charge is also non Coulombian.Since 1 / √ r + (cid:15) is not a valid Green function for the Laplace operator, also thefield produced by an extended distribution of particles (or continuum) is expectedto differ from that calculated using the real Coulomb interaction.Let us consider an infinitesimally thin shell of radius r ∗ and surface density σ = q/ (4 πr ∗ ), where q is its total charge. With the modification of the Coulombinteraction given in Eq. (2.34), the radial component of the electric field exertedby the surface element δs on a point placed at distance r from the centre of theshell, reads δE softs ( r ) = σδs ( r ∗ cos θ − r )( a + (cid:15) ) / , (2.35)where a = (cid:112) r ∗ + r − r ∗ r cos θ is the distance of the shell element from the pointwhere the field is evaluated and θ the angle between the vectors of length a and r , with origin set in that point. Setting δs = r ∗ d ψ sin θ d θ in Eq. (2.35) andperforming the integration over the two angular variables ψ and θ , the field atdistance r due to the whole shell is Chapter 2. Coulomb explosion of spherical cluster plasmas E R , ε (r) / ( Q / R ) r/R ε =0 ε =R /30 ε =R /20 ε =R /10 Figure 2.8: For different values of (cid:15) in units of the system’s R , softened electricfield at E R ,(cid:15) ( r ) inside an homogeneous sphere with softened interaction. E softs ( r ) = 2 πσr ∗ (cid:34) r ∗ r + r ∗ + (cid:15) r (cid:112) (cid:15) + ( r ∗ + r ) + r ∗ r − r ∗ − (cid:15) r (cid:112) (cid:15) + ( r ∗ − r ) (cid:35) , (2.36)from which in the limit of (cid:15) → E softs ( r ) = q/r for r ≥ r ∗ and E softs ( r ) = 0 for r < r ∗ .In Fig. 2.7 (top panel) we show the analytical estimation of the softened E softs ( r ) for different choices of (cid:15) , and its numerical computation averaged over20 realizations, with a direct summation code for a system of N = 10 particlesdistributed on a spherical shell (bottom panel).It is clearly evident that with this modification of the Coulomb interaction,a test particle placed inside a spherical charged shell would be prone to a radialforce, directed towards the centre if the particle and the shell have charges of op-posite sign, or directed outside in the opposite case. For r > r ∗ instead, the field E softs ( r ) results underestimated respect to the field E s ( r ) calculated with the realcoulomb interaction. The effect is more pronounced the larger (cid:15) is in units of theshell’s radius r ∗ .For a homogeneous charge distribution of density ρ and radius R , integratingthe contribution of every infinitesimal shell given by Eq. (2.36), gives a trend ofthe softened electrostatic field E R,(cid:15) ( r ) similar to that obtained in Eq. (2.29) whichtakes into account the effect of the correlation hole. Figure 2.8 shows the softenedelectric field produced by a uniformly charged sphere acting on a particle of unit.2. Shock shells in non-uniform density profiles charge placed at its interior.It is evident how E R ,(cid:15) ( r ) substantially departs from the linear trend of thereal force as r approaches the system’s edge R . The larger (cid:15) is the more impor-tant such deviation is.In conclusion, this means that using a softened interaction in numerical calcu-lations may spuriously increase the effect of the correlation hole between particles.This is overcome by choosing (cid:15) < δ . To this point we have studied the Coulomb explosion of systems having a uniforminitial density profile. In reality, clusters ionized by strong laser pulses may haveradially non uniform charge densities, due to their intrinsic initial structure or tothe details of the ionization process.In this case, the dynamics of the expansion might be highly non trivial (seee.g. [90] and [94]). For instance, when the radial component of the electric field E ( r ) is non monotonic with r and has its maximum inside the cluster, inner ionsmay eventually move faster than those initially placed at larger radii. This impliesthat sooner or later they will reach and then overtake them. Hereafter, we willrefer to overtaking as “shell crossing”. While inner ions rapidly catch up the outerones, the charge density depletes in the central region of the cluster and startsincreasing near its edge leading to the formation of what is called a shock.In the continuum or fluid picture, the shock wave, or shock front, is defined asa propagating discontinuity in the systems properties (i.e. velocity field, density,pressure), [115]. In our case it corresponds to a divergence in the cluster’s densityprofile.The Coulomb explosion of non-uniform systems and the control of the shockwave’s dynamics have received much attention (see [129], [131]) due to their possi-ble application in the context of intra-cluster fusion, [98], [5]. We are now going tosketch the problem in the continuum picture and then we will discuss the resultsof our particle based numerical calculations. Continuum model
Under the assumption of non collisionality (i.e. the effects of particles encoun-ters are negligible), when passing to the continuum approximation, the Coulombexplosion can be treated in principle with the equations of pressureless ( P = 0) hy-drodynamics written in Sect. 2.1. Note that, in correspondence of a shell crossing,the solution of the latter becomes multivalued since at the same radius elementsmatter with different velocities are found.The problem of an exploding cluster with non uniform initial density profilehas been studied in Ref. [90], for the one parameter family of initial density profilesgiven by Chapter 2. Coulomb explosion of spherical cluster plasmas ρ (r / R s ) β =0.5 β =1.0 β =10.00.20.40.60.81 0 0.5 1 1.5 2 2.5 3 3.5 4 E (r / R s ) r /R s Figure 2.9: Top panel: density profiles given by Eq. (2.37) for β = 0 . , R s . ρ ( r ) = 3 Q πR s (cid:34) (cid:18) r R s (cid:19) β (cid:35) − − /β ; 1 / < β < ∞ , (2.37)where Q is the total charge and R s a scale radius. The charge enclosed by radius r then reads Q ( r ) = (cid:90) r ρ ( r (cid:48) ) r (cid:48) d r (cid:48) = Qr R s (cid:20) (cid:16) r R s (cid:17) β (cid:21) /β . (2.38)Note that for finite β in Eq. (2.37), the density falls off to zero at infinity, whilein the limit of β → ∞ , ρ ( r ) tends to the uniform profile of Eq. (2.6), where R s = R . In Fig. 2.9 we show for some values of β the density profile ρ ( r ) andthe radial component of the electric field it generates.Formally, the dynamics of the expansion is given by the solutions of the systemof ODEs ρ ( r ) = 14 πr d Q ( r )d r ;d r d t = QM Q ( r ) r , (2.39).2. Shock shells in non-uniform density profiles that are nothing but the characteristics of the equations of hydrodynamics.For reasons of clarity, let us now use normalized variables with respect to thesystem’s initial parameters, s = r/R s and τ = t/t s where t s = (cid:112) M R s /Q . (2.40)With this choice, the normalized density is ρ ( s ) = ρ ( r ) R s /Q (2.41)and the system (2.39) becomes ρ ( s ) = 14 πs d Q ( s )d s ;d s d τ = Q ( s ) s . (2.42)Up to the critical time τ c when the shell crossing generates the shock, the firstintegral of the second equation is simply12 (cid:18) d s d τ (cid:19) = Q ( s ) (cid:18) s − s (cid:19) , (2.43)where Q ( s ) is the charge enclosed by the radius s at τ = 0. Trajectories s ( τ ) ofelements of volume initially placed at s are the implicit solutions of (cid:112) ξ ( ξ −
1) + ln( (cid:112) ξ + (cid:112) ξ −
1) = τ (cid:113) Q ( s ) /s ; ξ = s/s , (2.44)in the same fashion of the case of a homogeneous sphere, cfr. Eq. (2.12). Thecharge density then reads ρ ( s ) = (cid:16) s s (cid:17) d s d s . (2.45)For τ ≥ τ cr , the solutions above are invalid since instead of Eq. (2.43), one hasnow 12 (cid:18) d s d τ (cid:19) = (cid:90) Q ( s ) s d s. (2.46)At τ = τ cr the derivative of the velocity profiled v ( s )d s = dd s (cid:18) d s d τ (cid:19) (2.47)diverges at the critical coordinate s cr where the crossing takes place. Such diver-gence in the velocity profile corresponds to diverging density in s cr , as one canstill obtain from Eq. (2.45).At later times, v ( s ) stays multivalued in the region (shock shell) s < s < s Once ρ ( s ) becomes multivalued, Eq. (2.45) is invalid and it is substituted by ρ ( s ) =(4 πs ) − d Q ( s ) / d s . Chapter 2. Coulomb explosion of spherical cluster plasmas -4 -3 -2 -1
0 1 2 3 4 5 6 7 8 ρ (r) / ρ (r = ) r/R s τ =0.0 τ =1.0 τ =3.5 τ =6.0 τ =8.0 00.20.40.60.81 0 1 2 3 4 5 6 7 8 v (r / r s ) r/R s Figure 2.10: Left panel: formation and evolution of the shock shell in a clusterstarting with density profile given by Eq. (2.37) where β = 1. Right panel: velocityfield at the same times. To improve visualization, the densities are normalized inunits of the initial central density ρ ( r = 0) and radii in units of the initial scaleradius r s . Times τ are given in units of (cid:112) M r s /Q .where s and s are the radii of the so called leading and trailing shocks respec-tively.Numerical integration of (2.42) for various initial density profiles given by Eq.(2.37), revealed that the initial width of the shock shell is narrower in the limitof large β . In each case, however, it always broadens with time, spanning almostthe entire system.Interestingly, even in the limit of large β , the coordinate and the critical timeof the shock formation do not tend respectively to the systems’ edge and to 0 asone would expect. Instead one has s cr (cid:39) .
635 and τ cr (cid:39) .
237 as roots of (cid:112) s ( s −
1) + ln (cid:0) √ s + √ s − (cid:1) = 2 s / √ s − (cid:112) s ( s −
1) + ln (cid:0) √ s + √ s − (cid:1) = τ √ , (2.49)see e.g. [90].This implies in other words, that even an infinitesimal initial perturbation nearthe edge of the homogeneous sphere treated in Sect. 2.1, is enough to give rise toshocks, thus breaking the self similarity of the expansion.Remarkably, if one considers truncated initial density profiles, defined as ρ ( r ) = (cid:37) ( r ) θ ( R − r ) (2.50) Actually this is the case for all densities smoothly going to 0 at infinity. .2. Shock shells in non-uniform density profiles where (cid:37) ( r ) is a monotonically decreasing function of r and R the cutoff radius,the leading shock front disappears as it rapidly meets the discontinuity of thedensity gradient d ρ/ d r at the edge of the system. This is qualitatively similar towhat happens to a shock front moving through a non homogeneous medium, (seee.g. [28], [27]).If a singularity arises in the density profile, associated with a singularity in thevelocity profile, one may expect that also the kinetic energy distribution presentsa divergence in correspondence of the kinetic energy of the shock front.As example, in Fig. 2.10, we show at different times the density for a clusterwith initial profile given by Equation 2 .
37 for β = 1, as well as the associatedvelocity profile. The curves are obtained by numerically integrating the hydrody-namics equations. The formation of a “singularity” (due to the finite resolutionit is a sharp peak) in the density profile, happens already at early stages of theexplosion ( τ ∼ .
5) when the velocity profile has become multivalued (see rightpanel and analogous plots in Ref. [90]).Unfortunately, to derive the asymptotic n ( E ) from the potential energy of theinitial state is in principle not possible for an arbitrary initial density profile, usingthe method described in Sect. 2.1. One has to use instead an approach based onkinetic theory, we redirect the interested reader to Appendix C and the literaturetherein referenced. Numerical simulations using particles
Numerical simulations based on particles, aiming at studying intra-cluster nuclearfusion, have been performed in [129], [131] and [130] with more realistic initial con-ditions where a non negligible electron density were considered. This has revealedthat radial shocks in the ion density profile still occur, even when the residualelectron density screens part of the ions.In this work we are interested mainly in the structure of the asymptotic numberenergy distribution n ( E ), since it is an experimentally deducible quantity (fromthe time of flight spectrum), and how it is affected by an initial velocity distri-bution. To this scope, we have performed N -body simulations of pure Coulombexplosion (i.e. no electron contribution) for different non uniform initial densityprofiles, with a single kind of particles of mass m and charge q .We discuss here the properties of the final states of two families of initial den-sity profiles. In the first case, to model systems characterized by a flat core (i.e.almost homogeneous central region) and a smoothly decaying density in the outerlayer, we use the expression for ρ ( r ) given by Eq. (2.37) where β controls thesteepness of the system’s edge. In the second case, where we model instead sys-tems with a highly dense central region and an outer layer decaying with a fixedslope, we use the family of γ − models ρ ( r ) = Q (3 − γ )4 π R s r γ ( R s + r ) − γ ; 0 ≤ γ < , (2.51) Chapter 2. Coulomb explosion of spherical cluster plasmaswhere R s is a scale radius and γ = dlog( ρ ) / d s is the so called logarithmic densityslope. The charge enclosed by radius r is given by Q ( r ) = Q (cid:18) r r + R s (cid:19) − γ . (2.52)Note that, all the profiles given by Eq. (2.51) fall of as r − for r (cid:29) R s and for γ (cid:54) = 0 diverge for r →
0; the latter is not really an issue since we consider discreteparticles in the simulations. The expression (2.51) was introduced by Dehnen [39]in the context of galactic dynamics to model spherical galaxies with prominentcentral density cusp. We use it here to generate our initial conditions only for thefact that varying γ allows one to control the “importance” of the central densitycusp, from a flat core ( γ = 0) to an extremely steep cusp ( γ → ρ falls to zero. In particles simulations they are obviously intrinsically truncateddue to the finite number of particles, however the initial conditions are not char-acterized by a sharp edge.The initial ion velocities v ,i are sampled from a position independent Maxwelliandistribution n ( v ) = 4 π (cid:18) π (cid:19) / v exp( − v / . (2.53)and then renormalized to obtain the wanted value of the total initial kinetic energy K = N (cid:88) i =1 mv ,i . (2.54)For a system of N particles sampled from a given density profile, the initial po-tential energy is given by U = q N (cid:88) j (cid:54) = i =1 | r ,i − r ,j | , (2.55)where r ,i are the particle positions at time 0. Hereafter, we refer to the ratio η = K /U as “initial ion temperature”.We define the dynamical time scale of the simulated system as t dyn ≡ π (cid:114) mq ˜ n , (2.56)where ˜ n is the average particle density inside the radius r h enclosing half of theparticles at initial time. As a rule, the calculations are run up to the time whenthe system’s kinetic energy equals 99% of the total energy E tot = K + U .Figure 2.11 shows n ( E ) at final time for systems starting from the smoothed-step density profile of Eq. (2.37). In order to compare the curves on the samescale, energies are rescaled with respect to the maximal energy E max attained by.2. Shock shells in non-uniform density profiles n ( ε * ) ε * β =1 0 0.2 0.4 0.6 0.8 1 ε * β =2 0 0.2 0.4 0.6 0.8 1 ε * β =10 η =0.00 η =0.01 η =0.05 η =0.10 Figure 2.11: From left to right, scaled final differential energy distribution for β = 1 , η . -8 -7
40 80 120 ρ / ρ ( R s ) r/R s β =1 40 80 120r/R s β =2 40 80 120r/R s β =10 Figure 2.12: Density profiles for the systems of Fig. 2.11.the particles (averaged over the 20 fastest ones). We observe that when η = 0 (redcurves), the final n ( E ) is characterized by a peak close to E max that gets sharperas β increases (i.e. the initial density profile tends to the uniform), and a long tailat low energies.If the particles have already some kinetic energy at t = 0, as for instance insystem where some energy has been transferred from the laser pulse also to theions, the peak in the final energy is “smeared” and n ( E ) has its maximum value Chapter 2. Coulomb explosion of spherical cluster plasmasat lower values of E . n ( ε * ) ε * γ =0 0 0.2 0.4 0.6 0.8 1 ε * γ =1 0 0.2 0.4 0.6 0.8 1 ε * γ =2 η =0.00 η =0.01 η =0.05 η =0.10 Figure 2.13: From left to right, scaled final differential energy distribution, for γ = 0 , η . -9 -8 -7
40 80 120 160 ρ / ρ ( R s ) r/R s γ =0 40 80 120 160r/R s γ =1 40 80 120 160r/R s γ =2 Figure 2.14: Density profiles for the systems of Fig. 2.13.The corresponding density profiles ρ ( r ) are shown in Fig. 2.12. A density peakin proximity of the cluster’s edge is evident for all the systems starting form coldinitial conditions (red curves). Such feature is the “remnant” of the trailing frontof the shock shell and appears to be narrower for larger values of β , consistentlywith what predicted in [90].When ions are starting with a non zero temperature, the density peak becomes.2. Shock shells in non-uniform density profiles less prominent and disappears entirely for η > .
02. Generally, different initial ε P β γ η =0.00 η =0.01 η =0.05 η =0.10 η =0.20 Figure 2.15: Left panel: as a function of β , position of the maximum of the final E for different values of the initial energy ratio η in the interval (0; 0.2). Right panel:same as left but in this case the initial distribution of the simulation particles isextracted from a truncated γ − model.temperature produce qualitatively different final density profiles with several slopechanges. This is less evident for large values of β , see the case β = 10 in Fig. 2.12.Simulations of clusters with initial density profiles with central cusp, show acompletely different picture. In Figure 2.13 the final differential energy distribu-tion is shown for γ = 0 (no central density cusp), γ = 1 (mild cusp) and γ = 2(extreme cusp). The cases starting with η = 0 have, as expected, three differentbehaviors, with n ( E ) peaking at high energy for γ = 0, and at low energy for γ = 2. For γ = 1 the distribution is instead flat for a broad range of energyvalues.From the values of γ presented here, it appears that introducing a non zeroinitial temperature has no effect on the final n ( E ) for large values of γ . In fact, (seeFig. 2.13), in the largest γ = 2 case, the normalized n ( E ) for different initial val-ues of η are practically indistinguishable. For γ = 0 the situation is qualitativelysimilar to that of the above discussed flat cored models, with the peak energysmoothly shifting towards lower values. Curiously, the final n ( E ) for systems withintermediate values of γ (in this case γ = 1) has a much more abrupt transitionfrom the perfectly cold case ( η = 0) to gradually initially hotter systems, loosingits “plateau structure” even for very small values of η .The density profiles shown in Fig. 2.14 display essentially the same picture,with no significant effect (except at low radii) due to the initial temperature forsystems with large γ , and a similar behavior to that of the flat-cored systems forthe low γ cases.As a general remark, we observe from both n ( E ) and ρ ( r ), that even a smallamount of initial kinetic energy can significantly affect the explosion dynamics,leading to end states that considerably depart from those of the initially coldmodel, whether the initial density belongs to one or the other family. Chapter 2. Coulomb explosion of spherical cluster plasmasFigure 2.15 summarizes for the two families of initial density profiles, the effectof the initial temperature in shifting the maximum of the energy distribution. Forflat cored systems, Eq. (2.37) (left panel), the maximum of the energy distributionof the final states at fixed η falls to roughly the same energy E p (in units of thecutoff energy), for β ≥ .
5. This implies that the energy spectrum tends to differless and less as the initial density profile approaches the perfect step.For systems with central density cusp, Eq. (2.51) (right panel), it is foundthat for γ ≥ .
4, the normalized distributions peak at the same E p independentlyof the ratio η , while at lower slopes of the initial density cusp, it is shifted to lowervalues the greater η is. To study the Coulomb explosion of systems composed of two or more differentspecies of particles is particularly relevant with respect to the modelling of ionizedclusters of heteronuclear molecules, as well as core-shell clusters, [154], [108] oratomic clusters embedded in helium droplets [97].From a theoretical point of view, the problem of multi-component (also calledmulti-species) Coulomb explosion has been studied in [3] by means of kinetic the-ory, see also Appendix C, and in [123] and [2] using a simple two fluid model. Theeffect of a residual electron density has been treated in [135].In this part of the work we have performed simulations of pure Coulomb explo-sions of clusters containing two species of particles, aiming at studying their finaldifferential energy distributions. More detailed calculations involving ionizationend treating dynamics of the electrons are presented in the next chapter.
Continuum model
As we have seen, the dynamics of mono-component clusters is already quite com-plex if the initial density profile is non-uniform. Adding another degree of freedom,by making the system heterogeneous in composition, makes the problem even morecomplicated. However, assuming a continuum picture and restricting ourselves tothe uniform initial density still allows, under certain assumptions, to derive theasymptotic n ( E ) for one of the two components analytically.Following the approach of [123], let us consider a heterogeneous cluster ofinitial radius R , whose initial total charge density is given by ρ ( r ) = ρ θ ( R − r ) + ρ θ ( R − r ) , (2.57)as the sum of the individual (constant) densities of the two components ρ and ρ . The total charge Q is given in this case by Q = Q + Q = 4 πρ R πρ R , (2.58).3. Multi-component systems and the radial component of the electric field is inside R reads E ( r ) = 4 π ρ + ρ ) r. (2.59)The infinitesimal element of volume occupied by the component with ρ has unitmass δm and unit charge δq , for the other component we have instead δm and δq . Hereafter we assume δm /δm ≥
1, regardless of the ratio of their unitcharges. Therefore we refer to component 1 as heavy component and to component2 as light component .In order to characterize the system, we define now the two parameters a and b as a = δq /δm δq /δm ; b = ρ ρ + ρ . (2.60)Note that in the limit of a → ∞ the heavy component does not move, while for b = 1 (or b = 0) one retrieves the single component case treated in Sect 2.1.For large values of a , during the explosion the light component overtakes en-tirely the heavy one. Assuming that no shell crossing happens among elements ofthe same component, one can write the asymptotic energy (per unit charge) of anelement of light component as function of the initial coordinate r as E ( r ) = U , + U , = Q R (3 R − r ) + Q R r . (2.61)From Eq. (2.58) and the definition of b , Equation (2.61) can be rewritten as E ( r ) = Q R (cid:2) − b ) R + (3 b − r (cid:3) . (2.62)With these assumptions and making the same steps as in the case of a single com-ponent uniform system (cfr. Sect. 2.1 and Ref. [146]), the asymptotic differentialenergy distribution for the light component is obtained as n ( E ) = 3(3 b − / (cid:115) EE − − b ) E θ ( E max − E ) , (2.63)where E max = Q/R . For b = 1 from the formula above one obtains Eq. (2.22),while for b = 1 / a → ∞ ) one has n ( E ) = δ ( E − E max ) , (2.64)as E ( r ) is independent on r , crf. Eq. (2.62).Note also, that independently on b , in case of large a , the asymptotic distribu-tion for the heavy component n ( E ) is expected to resemble qualitatively that forthe one component model, since the dynamics of the latter is barely influenced bythe fast explosion of the lighter component. Chapter 2. Coulomb explosion of spherical cluster plasmas
Numerical simulations using particles
We have investigated the effect of a heterogeneous composition in the Coulombexplosion of a charged cluster by means of molecular dynamics calculations. Theinitial conditions are implemented as follows: N particles of mass m and charge q and N particles of mass m and charge q , are homogeneously distributedinside a spherical volume of radius R . No minimum inter-particle distance isenforced and, in all cases, the initial velocities v ,i are set to 0, therefore K = 0.The total energy of the system is given by its initial potential energy U = 12 N (cid:88) i (cid:54) = j =1 q i q j | r ,i − r ,j | ; N = N + N . (2.65)As usual, we define the final state of the system when the kinetic energy K equals99% of U .In presence of two species, the scale time in the simulations is redefined as t dyn ≡ π (cid:114) ¯ m ¯ q ˜ n , (2.66)where ¯ q and ¯ m are the average charge and mass respectively. Hereafter, we char-acterize each system with its fractional charge to mass ratio (cfr. also Eq. (2.60))given by a = q /m q /m , (2.67)and for reasons of convenience, we will label here with 2 the species with thehigher charge to mass ratio, so that a ≥ i − th particle placed at position r i , due to thecluster’s electrostatic field E , is given by a i = q i E ( r i ) /m i , (2.68)therefore, as expected, particles of the component with higher charge to mass ratio q/m are accelerated to higher velocities with respect to particles of the other com-ponent. Due to that the cluster expansion is not uniform and even if the initialdensity profile is homogeneous, one of the two components overtakes the other.This is clearly evident in Figure 2.16 (upper row) where we show at three dif-ferent times, for a cluster with a = 2, the projection of the particle’s positions. Anouter shell containing the faster particles is already present at early stages of theexplosion. From the phase-space sections radial velocity v r versus radial coordi-nate r (same figure, lower row), it can be noticed how once the faster componenthas overtaken the other, its velocity profile is not linear with r (showing a littlemultivalued region). However, (as expected, e.g. [90] and [123]) at later times(here τ = 10) it has again a linear trend.Figure 2.17 shows the final n ( E ) for the two components of three clusters with.3. Multi-component systems Figure 2.16: Upper row: particles positions in the x − y plane for a multi-component cluster ( N = N = 10 ) with a = 2 at τ = 2 t dyn , 5 t dyn and 10 t dyn .Red dots mark the particles of the component with higher charge to mass ratio andblue dots those with the lower. Lower row: for the same times as above, phase-space pairs r and v r . Coordinates are normalized in units of the initial clusterradius R while the radial velocities in units of the scale velocity v typ = R /t dyn .different values of the parameter a . In all cases N = N = 7500. The differen-tial energy distribution of the fast component peaks in all cases at high energies,containing roughly the 70% of the particles between 0 . E max and E max . For abroad range of energies is remarkably well fitted by a power law ( E ) k , see solidlines in Fig. 2.17. For increasing values of a , the exponent k increases making n ( E ) steeper. For the case shown here, k (cid:39) .
98 for a = 2, k (cid:39) .
89 for a = 4 and k = 9 .
33 for a = 10. In the limit of infinite a , the component 1 does not move,and the asymptotic distribution for component 2 is expected to have a delta-likestructure, see e.g. Refs. [102], [123] and [135].The qualitative explanation of such behavior is, if one has that for the initialcharge density of the two components ρ , (cid:28) ρ , , the cluster’s electrostatic fieldis always dominated by the effect of component 1 that moves on a considerablylarger time scale. Due to that, when a particle of component 2 coming from theinner regions of the cluster reaches the radius R , it has already a certain amountof kinetic energy and therefore will reach a larger asymptotic energy than a par-ticle of the same species initially sitting at R . The more the charge density ofcomponent 2 is small, the less the kinetic energies of particles coming from differ-ent inner radii once reaching R , differ from each other, thus steepening n ( E ). Chapter 2. Coulomb explosion of spherical cluster plasmas -2 -1
0 0.2 0.4 0.6 0.8 1 n ( ε * ) ε * a=2a=4a=10 00.511.522.5 0 0.2 0.4 0.6 0.8 1 n ( ε * ) ε * Figure 2.17: Left panel: asymptotic number energy distribution for the “fast”component of three initially cold cluster with N = N = 7500 and a = 2, 4 and10 (points), the best fit curves are also shown (lines). Right panel: same quantityfor the “slow” component. For comparison, the thin black line line marks thetheoretical asymptotic n ( E ) of the single species homogeneous sphere. Each curveis normalized to its individual cutoff energy.The number energy distribution n ( E ) for the slow component is, as expected,reminiscent of n ( E ) observed for a single species initially uniform system. Re-markably, the cusp at high energies, is milder than that observed in n ( E ) forsingle-component systems.Let us now consider the a = 1 case. Intuitively, independently of the densityprofile, if the combinations of charge and mass are such that q /m = q /m (forinstance the situation that one would have in a mixture of deuterons D + andcarbon ions C ) the two normalized asymptotic energy distributions n ( E ) and n ( E ) should in principle coincide. However, as seen in Fig. 2.18, for the endproducts of 3 direct N − body simulations this is not exactly true. In fact, theenergy distribution for the component with larger mass (and charge), in this case m = 2 m , shows a more pronounced high energy peak than the other. Thishappens regardless of the ratio N /N and persists using different binning for thenumerical energy distribution.By contrast, this is not the case for the final states of particle-mesh simulations(see Chap. 5) shown in Fig. 2.19 where the electric field acting on particles isnot computed by direct summation (as in MD simulations) but solving Poissonequation on a cartesian grid. In this case the two curves perfectly coincide.The high energy peak is to be interpreted as an effect of the discrete gridbased electric field, that close to the system’s surface is slightly underestimated,leading to an energy bunching effect, analogous to that induced by the correlationhole treated in Sect. 2.1.To interpret this discrepancy, let us consider now a test particle of mass m k and charge q k moving through a homogeneous spherical background of N par-ticles with charges and mass q and q , and N particles with charge and mass.3. Multi-component systems n ( ε * ) ε * N /N =2/3m m =2m
0 0.2 0.4 0.6 0.8 1 ε * N /N =1 0 0.2 0.4 0.6 0.8 1 ε * N /N =3/2 Figure 2.18: Final normalized energy distributions of the two components. Forthe three combinations of N /N indicated above, for systems with total numberof particles N + N = N = 2 × , a = 1 and m = 2 m . m and m , so that q /m = q /m . For consistency we assume m > m (andobviously q > q ).The average charge of the background particles is given by¯ q = ( q N + q N ) / ( N + N ) , (2.69)while their total number density is simply the sum of the individual numberdensities of the two species ¯ n = n + n .The Langevin-type equations for the radial motion of the test particle, (seee.g. [123], see also [161]) read˙ v = q k Q ( r ) m k r − ω coll v ; v = ˙ r, (2.70)where Q ( r ) = (4 π/ r ¯ n ¯ q is simply the charge inside radius r and ω coll = 8 π ¯ nq k ¯ q ln Λ m k ¯ µ k v (2.71)is the collision frequency with the background ions, where ln Λ is the Coulomblogarithm , (see [150]) and ¯ µ k = m k ¯ m ( m k + ¯ m ) (2.72)is the species averaged reduced mass.If the term depending on ω coll in (2.70) is non negligible, the test particle ex-periences an effective drag force along the radial direction due to the two body This quantity is defined as the natural logarithm of the ratio between the maximum andminimum impact parameter in the collision experienced by the test particle. See also AppendixB Chapter 2. Coulomb explosion of spherical cluster plasmas n ( ε * ) ε * lightheavy Figure 2.19: Final normalized n ( E ) for the two components of an N + N = N =8 . × multi-species spherical and homogeneous cluster, from a collisionlesssimulation with a particle-mesh code. Each curve is normalized respect to itscutoff E . For E > . √E ∗ trendshowing a dip followed sharp cusp, signatures of an energy bunching, due to theworse representation of ρ c and the potential for a spherical system on a cartesiangrid. See also Chap. 5.encounters with the background particles.Note that in such case, due to the different dependence of ˙ v on q k and m k , atest particle with m k = m and q k = q is prone to a different deceleration than atest particle with m k = m and q k = q starting with the same velocity from thesame initial position. In particular, for m = 2 m the more massive (and morecharged) test particle with mass m suffers a deceleration of a factor roughly 1.96larger than a particle with mass m .The Coulomb explosion of a spherical heterogeneous system with one charge tomass is still to a large extent a self similar expansion without particles propagatingthrough the others. However, at early stages of the explosion, the contribution tothe force due to near neighbours dominates over the mean field for most of theparticles, depending also on the local structure of the system. Particles of bothspecies do not have pure radial trajectories as for t → ∞ , and some energy cantransferred between radial and tangential motion via few “two body collisions”.Since, cfr. Eq. (2.70), such energy exchange depends on the (averaged) reducedmass, the heaviest component is expected to lose more kinetic energy in favour ofthe light ones. Thus, in addition to the energy bunching due to the discretenessdiscussed before, the heavier component suffers an additional shift towards lowerenergies causing the height of the peak in the normalized n ( E ) to increase.Note that this is also influenced by the percentages of the two species N and N since they enter the definition of ω coll twice, through ¯ µ and ¯ q .In order to characterize the effects of different mass ratios and different per-.3. Multi-component systems Figure 2.20: Left panel: for four different number ratios 0.5, 1, 1.5 and 2, ratioof the final kinetic energies of the light and heavy particles. Right panel: as afunction of 0 . < m /m < . < N /N <
1, the ratio of the averagekinetic energies for the two species is shown.centages of the two species, we now consider test systems of fixed total charge Q where q = q = Q/ ( N + N ), and we span the two intervals of m /m and N /N .We assume, m < m . In this way, if in the initial conditions K = 0 and thenumber density ¯ n is homogeneous inside the spherical volume of radius R , thetotal energy of the system, given by Eq. (2.65) is always the same in all the casesconsidered.For the final states of the numerical simulations we define the quantity R = (cid:104)E(cid:105) (cid:104)E(cid:105) (2.73)as a function of m /m and N /N , where (cid:104)E(cid:105) i is the final average energies of thetwo components.In Fig. 2.20 we show R as a function of m /m for some some relevant numberratios as well as the quantity as function of both N /N and m /m . Remarkably,see right panel, for mass ratios larger than roughly 0.75, the final energy ratio is notsignificantly influenced by different percentages of the two species in the cluster.For small values of m /m , R grows indefinitely for larger values of N /N whilegrows sublinearly for N /N < Chapter 2. Coulomb explosion of spherical cluster plasmas
We have reviewed here the basic aspects of the Coulomb explosion of spherical sys-tem in order to set the stage for our study of laser irradiated atomic and molecularclusters. By means of numerical simulations we have investigated the explosionof systems starting with different density profiles and ion temperature as well aswith heterogeneous composition.We have observed the discrepancy between the continuum model and sys-tems where the particulate structure is accounted for the homogeneously chargedsphere. This comes in the form of a spike in the number energy distribution n ( E )caused by the presence of a “correlation” hole in the initial particle distribution.The effect appears to be considerably reduced if no minimum inter-particle dis-tance is enforced in the initial condition (i.e. the initial state is characterized bya non lattice-like structure). In addition, we showed that the regularization onshort distances of the Coulomb (or Newton) force induces a spurious external fieldeffect in the interior of perfectly spherical systems of the same entity of that dueof the discreteness of the system itself.The formation of shock shells in systems with non uniform initial charge den-sity profile is found to be inhibited when an ion temperature η > . η .Finally, from the simulations of two component clusters we find that such sys-tems are characterized by a mean field segregation effect induced by the differentcharge to mass ratios of the different species. The effect of different percentageson the redistribution of the total energy on the two species is strongly influencedby the mass ratio and tend to vanish already at m /m ∼ . hapter 3 Coulomb explosion ofellipsoidal systems
Initially spherically symmetric cluster, once irradiated by strong laser pulses, mayassume non spherical (i.e. ellipsoidal) charge distributions, either due to the laserspatial polarization or to resonances between the laser frequency and the initialplasma frequency (see e.g. Refs [117], [118], [149] and [97]). Moreover, electron orion beams confined by electromagnetic fields in accelerators are known to adjustto ellipsoidal shapes depending on the parameters of the confining field (see e.g.[8], [88], [57], [67] and references therein).It is therefore useful to extend our discussion of Coulomb explosion carried onin the previous chapter, to the case of ellipsoidal systems. limiting ourselves tosingle-component systems, we study the cases of homogeneously charged triaxialand axisymmetric ellipsoids (i.e. spheroids) and finally axisymmetric systems withnon uniform initial density. In the same line of Chap. 2, we treat the problemfirst in the continuum picture, and then we discuss the analysis of our N − bodycalculations. To tackle the problem of non-spherical Coulomb explosion in the continuum ap-proach we first need the expressions for the potential due to ellipsoidal chargedistributions.Let us consider an infinitesimally thin ellipsoidal shell of semi-axes a (cid:48) , b (cid:48) and c (cid:48) , total charge q and uniform surface density. We shall call such object a homeoid .The electrostatic potential exerted by the homeoid at a point r = ( x, y, x ) placedon its exterior is given (see [91], [30]) byΦ hom ( r ) = q (cid:90) + ∞ λ d u (cid:112) ( a (cid:48) + u )( b (cid:48) + u )( c (cid:48) + u ) , (3.1) Chapter 3. Coulomb explosion of ellipsoidal systemswhere the lower boundary of integration λ is obtained as the largest root of thealgebraic equation x a (cid:48) + λ + y b (cid:48) + λ + z c (cid:48) + λ = 1 , (3.2)and its equipotential surfaces are ellipsoids confocal to it. The potential insidethe homeoid is constant, therefore it exerts no force at its interior. Such results itis known as third Newton’s theorem, see e.g. [91]. We now consider a continuumellipsoidal charge distribution with with semi-axes a, b, c and density ρ ( s ) strat-ified on concentric similar homeoids, for which we have introduced the so calledellipsoidal radius s = (cid:112) x /a + y /b + z /b . (3.3)Using Eq. (3.1) and the third Newton theorem, one can formally construct theexpression for the potential at a point r = ( x, y, z ) generated by such a genericellipsoidal charge distribution integrating the contribution of each homeoidal shell.We define the integration variable u from the solution of x a + u + y b + u + z c + u = s , (3.4)and the auxiliary integral function ψ ( s ) = (cid:90) s ( u )1 ρ ( s (cid:48) )d s (cid:48) . (3.5)The potential due to the charged ellipsoid then readsΦ( r ) = πabc (cid:90) + ∞ l (cid:2) ψ (1) − ψ ( s ( u )) (cid:3)(cid:112) ( a + u )( b + u )( c + u ) d u, (3.6)where l = 0 if r is inside the charge distribution and l = λ , with λ from Eq. (3.2)otherwise.In the special case when the density is equal everywhere to the constant value ρ c , ψ ( s ) = s and Eq. (3.6) becomesΦ( r ) = πρ c abc (cid:90) + ∞ l (cid:18) − x a + u − y b + u − z c + u (cid:19) d u (cid:112) ( a + u )( b + u )( c + u ) . (3.7)Note that the equations of motion for a test particle moving under a potential ofthe form (3.6), are in general not separable in cartesian coordinates. For particular forms of ρ ( m ), the equations are indeed separable in polar ellipsoidal coordi-nates ( ξ, η, ζ ) such that x = ( B + ζ ) cos ξ cos η ; y = ( B + ζ ) cos ξ cos η ; z = (cid:2) (1 − e ) B + ξ (cid:3) sin η ,where e = ( a − b ) /a and B = a/ (cid:112) − e sin ξ , if the potential can be expressed, see [42], asΦ( ξ, η, ζ ) = F ( ξ )( ξ − η )( ξ − ζ ) + F ( η )( η − ζ )( η − ξ ) + F ( ζ )( ζ − ξ )( ζ − η ) , where F ( x ) is a generic smooth function. A potential of such form is called a St¨ackel potential[152]. .2. Coulomb explosion of uniform axisymmetric systems The first cases of interest are those of homogeneously charged rotational ellipsoids(also known as spheroids), In our discussions of both continuum and particleapproaches we set a = b = a ⊥ and c = a (cid:107) . In addition, we define the quantity α ≡ a ⊥ a (cid:107) , (3.8)that we will call hereafter aspect ratio . When α <
1, the spheroid has an elongatedshape and it is called prolate, when α > α = 1 case corresponds to the sphere.Given the axial symmetry of such systems, it is convenient to introduce cylin-drical coordinates: r = (cid:112) x + y ϕ = , for x = y = 0arcsin( y/r ) , for x ≥ π − arcsin( y/r ) , for x < z = z, (3.9)where we have assumed conventionally that a (cid:107) is oriented along the z axis. Continuum model
The electrostatic potential inside a uniformly charged spheroid can be written in amore compact form using cylindrical coordinates and setting a = b = a ⊥ , c = a (cid:107) ,and u ∗ = u/a (cid:107) in Eq. (3.8). With a little algebra it readsΦ int ( r, z ) = 2 πρ c (cid:2) a ⊥ ζ ( α ) − z ζ (cid:107) ( α ) − r ζ ⊥ ( α ) (cid:3) , (3.10)where the three “shape functions” ζ , ζ ⊥ and ζ (cid:107) , are ζ ( α ) = 12 (cid:90) + ∞ d u ∗ ( α + u ∗ ) √ u ∗ = sec − α √ α − , (3.11) ζ (cid:107) ( α ) = α (cid:90) + ∞ d u ∗ ( α + u ∗ )(1 + u ∗ ) / = α (cid:16) √ α − − sec − α (cid:17) ( α − / , (3.12) ζ ⊥ ( α ) = α (cid:90) + ∞ d u ∗ ( α + u ∗ ) √ u ∗ = √ α − α sec − α α − / . (3.13)Note that, once made explicit they depend only on α . Alternative expressionsare given in term of hyperbolic functions and logarithms in Ref. [92], while their Chapter 3. Coulomb explosion of ellipsoidal systemsFigure 3.1: Upper panels: Potential in the coordinate planes x, z (left) and x, y (right) for a homogeneous oblate spheroid with α = 10. Lower panels: Potentialin the coordinate planes x, z (left) and x, y (right) for a homogeneous prolatespheroid with α = 0 .
1. The potential is normalized to its value at the origin Φ ∗ and the x, y, z coordinates to s = a (cid:107) for the prolate case and s = a ⊥ for the oblatecase.asymptotic forms for large and vanishing values of α are listed in Tab. (3.2)Once written in form (3.10), it is evident that the potential inside the homoge-neous spheroid is a harmonic function of the coordinates ( r, z ) only, of which ζ (cid:107) and ζ ⊥ are its two eigenfrequencies.Note that, any regular ellipsoidal distribution of charge or mass stratified onsimilar concentric homeoids, generates at its exterior a family of confocal ellip-.2. Coulomb explosion of uniform axisymmetric systems Case α → α = 1 α → ∞ ζ ( α ) ln(2 /α ) 1 π/ (2 α ) ζ (cid:107) ( α ) α [ln(2 /α ) −
1] 1/3 1 ζ ⊥ ( α ) 1/2 1/3 π/ (4 α )Table 3.1: Asymptotic trends of ζ , ζ (cid:107) and ζ ⊥ in the limit of large and vanishing α , as well as α = 1 (sphere).soidal equipotential surfaces, see Refs. [30] and [91], therefore equipotential andequidense surfaces do not coincide. For a discussion of the analogous problem forthe gravitational collapse of spheroidal initially cold distributions of matter seee.g. [107], [53], [14], [17] and references therein.The internal electrostatic field, given by E int ( r, z ) = −∇ Φ int ( r, z ) is then,written by components (cid:40) E int z ( z ) = − ∂ z Φ int ( r, z ) = 4 πρ c ζ (cid:107) ( α ) zE int r ( r ) = − ∂ r Φ int ( r, z ) = 4 πρ c ζ ⊥ ( α ) r. (3.14)In other words, this means that the perpendicular and parallel components of E int are linear functions of their coordinates r and z respectively.In Fig. 3.1, as an example, we show a section of the equipotential surfaces forthe potential generated by a prolate ( α = 10) and an oblate ( α = 0 .
1) spheroid.Having established an expression for the electric inside the spheroid, the sep-arated nonrelativistic equations of motion, for an infinitesimal element of volumeof unit charge δq and mass δm readd d t ˜ z = 4 π δqδm ρ c ζ (cid:107) ( α )˜ z, d d t ˜ r = 4 π δqδm ρ c ζ ⊥ ( α )˜ r (3.15)where ˜ z = z/z and ˜ r = r/r and r and z are the element’s initial coordinates.Due to the linearity in r and z of the components of E int , and to the symmetryof the problem, using the same arguments of Sect. 2.1 we find that differentelements of volume coming from two nested ellipsoidal shell can not cross eachother. Thus a homogeneously charged spheroid starting with initial conditionscharacterized by vanishing kinetic energy expands retaining a homogeneous chargedensity, see also Ref. [8]. However, since Eqs. (3.15) are essentially those for aparticle in a 2d harmonic repulsor with two different eigenfrequencies, we expectthat the aspect ratio α changes as a function of time since the ellipsoid expandsat different rates along the z and r directions.Since no overtaking can happen among different ellipsoidal shells, the Coulomb Potential and electric field on the exterior of the ellipsoid are obtained simply by putting λ as inferior extreme of integration in the shape functions. Chapter 3. Coulomb explosion of ellipsoidal systemsexplosion of a uniform charged spheroid of total charge Q and mass M and initialaspect ratio α = a ⊥ , /a (cid:107) , , is fully determined by the time evolution of its twosemi-axes a ⊥ and a (cid:107) . Essentially this is done by integrating the equations ofmotion (cfr. also Eqs. 3.15) for two infinitesimal elements of volume initiallyplaced at one of the two poles ( z = a (cid:107) , ) , r = 0 and on the equatorial plane( z = 0 , r = a ⊥ , ), that readd d t a (cid:107) = ω a ⊥ , a (cid:107) , a ⊥ ζ (cid:107) (cid:18) a ⊥ a (cid:107) (cid:19) , d d t a ⊥ = ω a ⊥ , a (cid:107) , a ⊥ a (cid:107) ζ ⊥ (cid:18) a ⊥ a (cid:107) (cid:19) . (3.16)In the equations above ω = Q (cid:113) / ( M a (cid:107) , a ⊥ , ) is the initial plasma frequency. Inorder to simplify the notation we use the normalized quantities τ = ω t , ˜ a ⊥ = a ⊥ /a ⊥ , and ˜ a (cid:107) = a (cid:107) /a (cid:107) , , so that Eqs. (3.16) becomed d τ ˜ a (cid:107) = ζ (cid:107) (cid:16) α a ⊥ ˜ a (cid:107) (cid:17) ˜ a ⊥ , d d τ ˜ a ⊥ = ζ ⊥ (cid:16) α a ⊥ ˜ a (cid:107) (cid:17) ˜ a ⊥ ˜ a (cid:107) . (3.17)Assuming initial kinetic energy K = 0, the set of (normalized) initial conditionsreads ˜ a (cid:107) ( τ = 0) = 1 , ˙˜ a (cid:107) ( τ = 0) = 0 , ˜ a ⊥ ( τ = 0) = 1 , ˙˜ a ⊥ ( τ = 0) = 0 . (3.18)The first integral of (3.17), accounting for the conversion of potential energy intokinetic energy is formally12 (cid:18) d˜ a (cid:107) d τ (cid:19) = (cid:82) ˜ a (cid:107) ( τ )1 ζ (cid:107) (cid:18) α a ⊥ ˜ a (cid:107) (cid:19) ˜ a ⊥ d˜ a (cid:107) , (cid:18) d˜ a ⊥ d τ (cid:19) = (cid:82) ˜ a ⊥ ( τ )1 ζ ⊥ (cid:18) α a ⊥ ˜ a (cid:107) (cid:19) ˜ a ⊥ ˜ a (cid:107) d˜ a ⊥ . (3.19)Noting that at every time ˜ a (cid:107) and ˜ a ⊥ can be written as function of the normalizedtime dependent aspect ratio ˜ α , the latter become12 (cid:18) d˜ a (cid:107) d τ (cid:19) = (cid:82) ˜ a (cid:107) ( τ )1 ζ (cid:107) ( α ˜ α )˜ α ˜ a (cid:107) d˜ a (cid:107) , (cid:18) d˜ a ⊥ d τ (cid:19) = (cid:82) ˜ a ⊥ ( τ )1 ζ ⊥ ( α ˜ α )˜ α ˜ a ⊥ d˜ a ⊥ . (3.20).2. Coulomb explosion of uniform axisymmetric systems In the limit of τ → ∞ , at the right hand sides of the equations above one has the -3 -2 -1 α ∞ α Figure 3.2: Asymptotic aspect ratio α ∞ as a function of α . The thin black lineshows the limit ∼ .
86 for α → ∞ .asymptotic kinetic energies E (cid:107) , ∞ and E ⊥ , ∞ for the two elements of volume startingfrom one of the poles, and from the equator. Such energies are also the maximalkinetic energies that elements of volume can attain along the r and z coordinates,we can therefore assume that α ∞ ∝ (cid:113) E ⊥ , ∞ / E (cid:107) , ∞ . (3.21)Unfortunately, due to the dependence of ζ (cid:107) and ζ ⊥ on the time dependent as-pect ratio α , no analytical solution is available for the system of Equations (3.17)in terms of simple functions for a general value of α , and one is forced to solve itnumerically for instance with explicit finite difference schemes.However, in the two special cases α → ∞ and α → α ∞ making use of (3.20) and (3.21) and the asymptotic forms of ζ ⊥ and ζ (cid:107) , ob-taining α ∞ ∼ .
86 and ∞ respectively. In Fig. 3.2, we show the asymptotic finalaspect ratio for a large range of α . Note that the convergence to the limit valuefor large α is particularly fast.Figure (3.3) shows the time evolution of α for some values of α between 0.1and 10 obtained integrating numerically Eqs. (3.17). It is evident how initiallyprolate spheroids become oblate as they Coulomb-explode and initially oblatespheroids become instead prolate. This means that in the first case the expansionoccurs mainly in the transverse direction, while in the latter along the symmetryaxis. Remarkably, while for extremely prolate configuration the α ∞ is unbound,for extremely initially oblate systems, the limit aspect ratios are bounded by a Chapter 3. Coulomb explosion of ellipsoidal systemsfinite value (i.e. one can not obtain final arbitrarily prolate configuration). -1 -1 a ( t ) t a =0.1 a =0.2 a =2.0 a =5.0 a =10. Figure 3.3: Evolution of the aspect ratio starting from different values in theinterval (0.1-10) obtained via the finite difference integration of Eqs. (3.17).In Fig. 3.4, for the same systems introduced above, the time dependent max-imal energies along the transverse and parallel directions are plotted. In the caseof a prolate spheroid ( α < α >
1) spheroid.It remains to derive an expression for the asymptotic differential energy dis-tribution n ( E ). Due to the linearity at every time of the components of the electricfield in the coordinates z and r (cfr. Eq. (3.14)), the two components v ⊥ ( r ) and v (cid:107) ( z ) of the velocity of each infinitesimal element of volume are at any time thelinear functions of its coordinates z, rv (cid:107) ( z ) = ( z/a (cid:107) ) v (cid:107) , max , (3.22) v ⊥ ( r ) = ( r/a ⊥ ) v ⊥ , max , (3.23)where v (cid:107) , max and v ⊥ , max are the maximum values of the velocity along the symme-try axis of the spheroid and in the equatorial plane respectively. As a consequenceof that, the initial potential energy of the system maps into its asymptotic kineticenergy, and thus we can use the same arguments of Sect. 2.1 to extract its finaldistribution.Since in a homogeneous spheroid the charge contained in a thin cone directed.2. Coulomb explosion of uniform axisymmetric systems Figure 3.4: Upper panel: Maximum energy along the symmetry axis in units ofits asymptotic value for the same values of α of Fig. 3.3. Lower panel: Maximumenergy in the equatorial plane in units of its asymptotic value.along a ⊥ , of infinitesimal opening angle δϑ increases with z , and analogously thecharge in an infinitesimal cylinder at the equatorial plane increases with r , thetwo partial energy distributions n ( E (cid:107) ) and n ( E ⊥ ), along the symmetry axis and inthe equatorial read n ( E (cid:107) ) = = 32 (cid:115) E (cid:107) E (cid:107) , max θ ( E (cid:107) , max − E (cid:107) ) , (3.24) n ( E ⊥ ) = = 32 (cid:115) E ⊥ E ⊥ , max θ ( E ⊥ , max − E ⊥ ) , (3.25)where E (cid:107) , max = δmv (cid:107) , max / E ⊥ , max = δmv ⊥ , max / t , it is evident from Eqs. (3.22), that all the equivelocity surfaces(i.e. surfaces where the modulus of the expansion velocity is constant) have thesame aspect ratio α v ( t ) = a ⊥ ( t ) a (cid:107) ( t ) v (cid:107) ,max ( t ) v ⊥ ,max ( t ) . (3.26)We stress the fact that such aspect ratio is different from that of the charge dis-tribution (and it is the reason why the latter changes), see the sketch in Fig. 3.5. Chapter 3. Coulomb explosion of ellipsoidal systemsFor the reason that no overtaking happens, the asymptotic velocity on each ellip- zz zr r
Figure 3.5: Section along the symmetry axis of the equivelocity surfaces (bluedashed lines) for a prolate (left) and an oblate (right) spheroid (gray shaded area).soidal shell of aspect ratio α v ( t ) (and therefore the corresponding kinetic energy)depends on its initial potential energy. From the charge enclosed inside everyequipotential surface one obtains n ( E ) α < = 32 E ⊥ , max (cid:113) EE (cid:107) , max , for E < E (cid:107) , max (cid:113) E ⊥ , max −EE ⊥ , max −E (cid:107) , max , for E (cid:107) , max < E < E ⊥ , max (3.27)for the prolate case, and n ( E ) α > = 32 E ⊥ , max (cid:113) EE (cid:107) , max , for E < E ⊥ , max (cid:113) EE (cid:107) , max − (cid:113) E−E ⊥ , max E (cid:107) , max −E ⊥ , max , for E ⊥ , max < E < E (cid:107) , max (3.28)for the oblate case. In Fig.3.6, left panel, the ratio χ ∞ = ( E ⊥ , max / E (cid:107) , max ) α < ;= ( E (cid:107) , max / E ⊥ , max ) α > is shown. In the right panel the asymptotic differentialenergy distributions are shown for different values of χ ∞ .Note that if E is normalized to its maximum value E max = E (cid:107) , max for the oblatecase and E ⊥ , max for the prolate case, the latter equations consent one parameterfits with the end products of numerical simulations of the Coulomb explosion ofhomogeneous spheroids, as we will see in the next Section. Note also, that for α (cid:54) = 1 the scaled energy distribution vanishes for E ∗ = 1. This is not surprisingas the latter energy is in all cases that of a subset of measure 0 of the whole chargedistribution (i.e. the poles for the initially prolate systems and the equatorial ring,for the initially oblate ones)..2. Coulomb explosion of uniform axisymmetric systems -2 -1 χ ∞ α -2 -1 n ( ε * ) ∞ ε * Figure 3.6: Left panel: ratio χ ∞ of the energy corresponding to the max of n ( E ) tothe maximal energy. Right panel: differential energy distribution n ( E ) for differentvalues of χ ∞ , the thick green curve marks the α = 1 case. Numerical simulations using particles
We have run N − body simulations of single component spheroidal Coulomb explo-sion with both direct molecular dynamics and with particle-in-cell (PIC) schemes.For an extended treatment of the numerical methods here involved, again we di-rect the reader to Chap.(5) and the references cited therein.In both sets of simulations we have spanned the range of initial aspect ratio α going from 0.1 to 10 and considered only systems starting with no kinetic en-ergy. The dynamical time t dyn used to normalize times is again defined as in Eq.(2.56) and, as usual, we run the calculations up to τ .As discussed in Ref.[65], the products of our numerical simulations are in goodagreement with the analytic predictions discussed above. An example, Figure 3.7shows at different times the positions of the particles from two MD simulations ofthe Coulomb explosion of initially cold spheroids with α = 0 . α = 10, andthe same number density and particle number N = 10 . Note how the expansionis faster in the transverse direction for the first case, while is faster along the z axis in the latter.As one would expect from the theoretical time evolution of E ⊥ , max and E (cid:107) , max (see Fig. 3.4), for fixed initial particle density and dynamical time, for the α = 1case the conversion of the initial potential energy U in kinetic energy happensfaster than for other values of α . This is shown in the right panel of Fig. 3.8.However, (see right panel of the same figure), for larger and smaller values of α , higher energies (in units of E max for the sphere with same particle density) canbe obtained along the symmetry axis (for initially oblate spheroids), and in themeridional plane (for initially prolate spheroids).Knowing its analytical expression in the continuum picture, it comes naturalto compare the numerical n ( E ) of the simulations’ end products with its analyti-cal counterpart. For direct MD simulations, in Fig. 3.9 we show the differentialenergy distribution at τ for the indicated values of α as well as the fitted Chapter 3. Coulomb explosion of ellipsoidal systemsFigure 3.7: Projection at 3 different times of the positions in the x, z (upper row)and x, y (bottom row) planes, for of an initially cold prolate spheroid ( α = 0 . α = 10, bottom panel) withthe same initial density. Particle’s coordinates are given in units of a ⊥ , , in bothcases N = 10 ..2. Coulomb explosion of uniform axisymmetric systems curves given by Eqs. (3.27) and (3.28). The analytic expression derived for theFigure 3.8: Left panel: maximal energies along the symmetry axis (circles) andin the equatorial plane (crosses) expressed in units of the maximal asymptoticenergy for the sphere for N − body simulations with N = 2 . × . Right panel:Time dependence of the kinetic energy K in units of the total initial electrostaticpotential energy U for α = 0 . , . , . , , α = 1 case for which the energy conversion happens in the shortest time inunits of t dyn .continuum model fits the numerical n ( E ) for a broad range of energies (on averageroughly the 75% of the whole spectrum). However, non negligible discrepanciesare found at low and high values of E ∗ , where the residuals (see bottom panel)exceed 0.3. The reasons why one finds such deviations from the theoretical curvenear the normalized cutoff energies E ∗ = 0 and E ∗ = 1 are essentially different.Particles occupying the low energy region of the spectrum are those sitting closeto the centre of the system in the initial condition. In homogeneous spheroidalsystems the number of particles increases linearly with the ellipsoidal radial co-ordinate s . This means that in a homogeneous distribution sampled with pointparticles at low s (i.e. close to the centroid) there is a relatively small number ofparticles, implying that the numerical n ( E ) is depopulated close to the minimumenergy.If the deviation in the lower energy part of the spectrum has essentially a triv-ial origin, what happens near the cutoff energy is different and has the same originas in the case of the homogeneous sphere discussed in Chap. 2 and Ref.[146].Note that for a generic homogeneous spheroidal system with α (cid:54) = 1, the the-oretical n ( E ) vanishes for E = E max , see Equations (3.27) and (3.28), the α = 1case stands out as the the only one having the maximum value of n ( E ) right at E max . Discreteness effects leading to energy bunching, affecting essentially the con-tribution to the differential energy distribution at E max , are for spheroidal systemsmitigated by the fact than particles reaching E max do no come from a surface, asin the case of a sphere, but from the poles for an initially oblate spheroid, andfrom the equatorial ring for an initially prolate spheroid, hence, the less the frac- Chapter 3. Coulomb explosion of ellipsoidal systems -1 n ( e * ) a =0.53 a =0.84 a =2.00 a =10.0-0.4-0.200.20.4 0 0.2 0.4 0.6 0.8 1 D n ( e * ) / n ( e * ) e * Figure 3.9: Upper panel: differential energy distribution n ( E ∗ ) of the end productsat τ (points) for different direct N − body simulations ( N = 2 . × ), with theindicated values of the initial aspect ratio α . Fits with their analytic expressions(lines) are also shown. Each curve is normalized to its maximum energy, so that E ∗ = E / E max . Lower panel: residuals of the fits ∆ n ( E ∗ ) /n ( E ∗ ).tion of the system contributing to E max , the less the numerical n ( E ) is prone todiscreteness effects.Spanning a broad range of values of α , we found that the peak in the nu-merical n ( E ) disappears for α < .
25 and α >
2, for initial conditions whereno minimal inter-particle separation is fixed, while it is always present (in partic-ular for initially oblate systems) if a fixed inter-particle distance δ between nearneighbours is enforced, albeit not as in the spherical case. Fig. 3.10 shows thefinal energy distributions for systems with and without fixed minimum distance δ ,for initially prolate and oblate spheroids; the systems with fixed minimum δ (redcurves) are clearly departing from the analytic function.The final number energy distributions from the particle-in-cell simulationsmade with the code calder (see Refs. [134] and [104] for its details, see alsoChapter 5 for a more general description of Particle-mesh codes), starting fromanalogous initial conditions also show the high energy peak as shown in Fig. 3.11.Since in this type of numerical approach, the potential and force are based on thesolution of PDEs rather than on the direct sums of individual particles contribu-tions, one would expect it to better reproduce the quantities obtained analyticallyfrom continuum models , having “erased” discreteness effects.Nevertheless, the fact that charge densities and electromagnetic fields are dis-.3. Uniform triaxial systems n ( ε * ) ε * α =0.1variable δ fixed δ n ( ε * ) ε * α =10 Figure 3.10: Differential energy distributions for the case with imposed initialminimal inter-particle distribution (red curves) and random distribution (greencurves) for initially prolate ( α = 0 .
1) and initially oblate ( α = 10) spheroids.The thin black lines marks the analytic expressions. All cases have the sameaverage particle number density and number of particles N = 10 .cretized on a mesh implies, contrary to what happens in direct simulations, thatnear the surface of the system the density is spuriously lower, even if the parti-cles are distributed homogeneously in the continuum space. This resulting in aneffectively decreasing density profile, and thus prone to generate shocks of whichthe high energy peaks in n ( E ) are a clear signature.For the PIC simulations presented here and in Ref. [65], we found that n ( E ) isreproduced fairly well for the α = 1 case (for more than the 80% of the normalizedenergy interval), except for the usual prominent peak near E max , while the finaldifferential energy distribution for the spheroids with α < α > n ( E ) is well reproduced for all the systemssimulated with the PIC code, since the combined effect of a larger number ofparticles and the way the force is computed do not allow for its depopulation. The imaging of complex biomolecules with x-ray lasers (see Refs. [31] and [58], seealso [126]) presents the problem of recovering the tridimensional structure from2-dimensional images (projections) of different samples whose orientation withrespect to the laser focus is unknown. Chapter 3. Coulomb explosion of ellipsoidal systems -1
0 0.2 0.4 0.6 0.8 1 n ( ε * ) ε * α =0.1 α =1.0 α =10. Figure 3.11: Normalized differential energy distributions n ( E ∗ ) for the cases with α = 0 .
1, 1 and 10 for the end products of PIC simulations. Note that the energiesare rescaled to each case’s maximal energy so that E ∗ = E / E max . Note also, thatin the lower energy region the numerical points do not fall considerably underthe analytic curve as for the direct simulations. All systems have the same initialnumber density and number of particles N = 1 . × .If one assumes as a first crude approximation that such molecules are triaxialellipsoids, the energies of their fragments may in principle give us informationon their orientation. This makes us some interest to investigate the Coulombexplosion of triaxial systems. Continuum model
The potential at a given point r = ( x, y, z ) generated by a triaxial ellipsoid filledwith homogeneous charge density is a harmonic function of the ( x, y, z ) coordinatesof the form Φ x,y,z = Φ − (cid:0) ω a x + ω b y + ω c z (cid:1) , (3.29)where Φ is a constant and the three (time dependent) eigenfrequencies ω a , ω b , ω c depend on the three semi-axes a, b and c through incomplete elliptic integrals (seeRef.[30]).In analogy to what has been done in cylindrical coordinates for a rotationalellipsoid, in Cartesian coordinates the Coulomb explosion of a triaxial ellipsoidof initial charge density ρ c, and initial semi-axes a , b and c oriented along x, y and z , is fully described by three coupled ODE (envelope equations) for thethree infinitesimal volume elements of mass δm and charge δq initially placed.3. Uniform triaxial systems at ( a , , , b ,
0) and (0 , , c ). Setting C = 2 πρ c, a b c δq/δm the threeFigure 3.12: Time evolution of the ratios b/a (top panel) and c/a (bottom panel)for different triaxial uniform ellipsoids with intermediate semi-axis a , obtained byintegrating numerically the three envelope equations. Times are normalized to theinverse of the initial plasma frequency ω .equations read d a d t = Ca (cid:90) + ∞ ψ ( a, b, c, u ) d u ( a + u ) , (3.30)d b d t = Cb (cid:90) + ∞ ψ ( a, b, c, u ) d u ( b + u ) , (3.31)d c d t = Cc (cid:90) + ∞ ψ ( a, b, c, u ) d u ( c + u ) , (3.32)where ψ ( a, b, c, u ) is defined by ψ ( a, b, c, u ) ≡ (cid:112) ( a + u )( b + u )( c + u ) . (3.33)Unfortunately, in Cartesian coordinates it is not possible to derive simple limitexpressions for the asymptotic axial ratios b ∞ /a ∞ and c ∞ /a ∞ due to their im-plicit dependence on elliptic functions in the auxiliary variable u . The analogous Chapter 3. Coulomb explosion of ellipsoidal systems /a b / a ε c , ∞ / ε a , ∞ b / a ε b , ∞ / ε a , ∞ /a b / a c ∞ / a ∞ b / a b ∞ / a ∞ Figure 3.13: Left panels: final axial ratios c ∞ /a ∞ and b ∞ /a ∞ for combinations ofinitial ratios c /a and b /a in the interval (0 . ,
1) and (1, 10) respectively. Rightpanels: maximal kinetic energies ratios E b, ∞ / E a, ∞ and E c, ∞ / E a, ∞ along the samedirections.problem arises also in polar ellipsoidal coordinates, and there is in general lessanalytic work to do.Intuitively, it must be pointed out that triaxial ellipsoids that depart onlyslightly from prolate or oblate spheroids, are expected to behave not very differ-ently from their regular counterparts. By integrating numerically Eqs. (3.30),(3.31) and (3.32) for the time evolution of a, b and c , one obtains the asymptoticfinal axial ratios b ∞ /a ∞ and c ∞ /a ∞ as well as the maximal kinetic energies ener-gies reached along the three axes during the the explosion.Figure 3.12 shows for some initial values of the three semi-axes the time evo-lution of the ratios b/a and c/a . It is clearly evident how the intermediate semi-axis (in these cases a ) remains the intermediate as the system expands, whilethe initially longer is asymptotically the shorter and vice versa. As a generaltrend, for fixed initial charge density, charge and mass, the more the initial con-figuration is markedly triaxial, the slower (in units of the dynamical time scale ω − = (cid:112) M a b c / /Q ) is the convergence to the limit axial ratios.Obviously, the largest among the asymptotic maximal energies along the threesemi-axes, E a, ∞ , E b, ∞ and E c, ∞ , is that along the direction of the initially short-est semi-axis, as it is evident from the left panels of Fig. 3.13 where the ratios E b, ∞ / E a, ∞ and E c, ∞ / E a, ∞ are shown. The ratios of the corresponding semi-axis.3. Uniform triaxial systems r a ti o b/ac/a00.20.40.60.81 0.1 1 10 ε / ε ∞ , m a x τ ε a ε b ε c Figure 3.14: Top panel: evolution of the axial ratios b/a and c/a (top panel)for an N − body simulation starting from a markedly triaxial initial configuration( b /a = 0 . c /a = 0 . a, b and c in units of the asymptotic maximal value of the three E ∞ , max = E c, ∞ .are also shown (right panels), and present the same trend with the initial combi-nations of a , b and c .In contrast to the spheroidal systems, to derive the asymptotic expression for n ( E ) is not possible in terms of simple functions since an integration over a triax-ial ellipsoidal volume is involved, therefore we can in principle derive them only“empirically” from the end products of our N − body simulations. Numerical simulations using particles
Following the same criterion for axisymmetric systems, we performed direct N − bodysimulations of Coulomb explosion of initially cold triaxial ellipsoids.We spanned the range of ratios (1, 10) in b /a and (0.1, 10) in c /a keep-ing fixed the system’s number density n = 3 N/ πa b c and total charge andmass. In none of the runs we impose a minimal inter-particle distance. As for thespheroidal systems, the dynamical time scale t dyn of the simulation is defined byEq. (2.56).The evolution of the maximal energies along the three semi-axis as well asthat of the ratios of the latter reproduces nicely what found in the continuum Chapter 3. Coulomb explosion of ellipsoidal systems n ( ε * ) ε * b /a =0.2; c /a =0.1b /a =0.5; c /a =0.1b /a =0.9; c /a =0.1 Figure 3.15: Normalized differential energy distribution n ( E ∗ ) at τ with respectto E max for three different combinations of b /a and c /a that are slightly prolate,markedly triaxial and slightly oblate. In all cases N = 2 . × .model. In Fig. 3.14 we show the time evolution of such quantities for an initiallymarkedly triaxial ellipsoid.In Fig. 3.15 the n ( E ∗ ) is plotted for the end products of three N − body simu-lations starting from three particular initial conditions quasi-prolate ( b /a = 0 . c /a = 0 . b /a = 0 . c /a = 0 .
1) and quasi-oblate( b /a = 0 . c /a = 0 . a (cid:39) b > c or a (cid:39) b < c while for particularly triaxial initial configurations,the n ( E ) appears to be a combination of those of the two (regular) types, oblateand prolate.It must be noted, that none of the analyzed systems in the ranges of initialaxial ratios presents a spike at high energies in n ( E ) as observed for spherical andspheroidal systems. To conclude, in the same line of Chapter 2, we now briefly treat the Coulombexplosion of spheroids with nonuniform density profile. Such problem is of somerelevance for instance in the field of particle acceleration from the interaction ofstrong lasers with nanostructured targets with non uniform densities (see e.g. [64]and references therein), and as also pointed out in Ref. [90].As it happens to spherical systems with non uniform initial density profiles,.4. Non uniform axisymmetric systems also spheroidal systems with non uniform density are expected to undergo shell α ∞ α homogeneous γ =0.0 γ =0.5 γ =1.0 γ =1.5 γ =2.0 Figure 3.16: Aspect ratio at τ as function of α in the range (0.1, 10). Differentsymbols refer to different values of the logarithmic density slope γ . The solid blackline marks as in Fig. 3.2 the theoretical relation for the case of a homogeneousdensity profile.crossing as they Coulomb explode. The additional complication due to their, nonpreserved axial ratios makes the problem even more prohibitive and, in principle,complete analytic treatment is impossible and one is forced to rely on numericsIn our MD simulations, we extracted the initial positions of the particles, fromthe family of triaxial models with density given by ρ ( s ) = Q (3 − γ )4 πabc s − γ (1 + s ) γ − ; 0 ≤ γ < , (3.34)where the “ellipsoidal” radius s is defined by s ≡ (cid:112) ( x/a ) + ( y/b ) + ( z/c ) (3.35)and γ = dlog( ρ ) / d s is the (angularly averaged) logarithmic density slope. Asusual, Q is the total charge and a, b and c are as usual the three semi-axes.The density given by Eq. (3.34) is nothing but the triaxial generalization of ofDenhen’s one parameter γ − models (see e.g. Ref. [113]), also used to in the fieldof beam physics to model charged particles bunches in accelerators, see Refs. [87]and [89]. We made this choice since it allows one to control (even in this case) theimportance of the density cusp at the system’s centroid, from flat cored systems( γ = 0) to highly concentrated systems ( γ = 3). Chapter 3. Coulomb explosion of ellipsoidal systems n ( ε * ) ε * γ =0 0 0.2 0.4 0.6 0.8 1 ε * γ =1 0 0.2 0.4 0.6 0.8 1 ε * γ =2 α =0.2 α =1.0 α =10. Figure 3.17: Differential energy distribution n ( E ∗ ) for γ = 0 , α .For reasons of simplicity, and having established that homogeneous triaxialellipsoids do no show a considerably different behavior from that of the spheroids,here we limit ourselves to the axially symmetric a = b = a ⊥ , c = a (cid:107) case. Notethat in principle the density given by Eq. (3.34) falls to 0 only at s → + ∞ , herewe have defined truncated models with a surface defined by s t = (cid:112) ( x/a ) + ( y/b ) + ( z/c ) = 1 . (3.36)Figure (3.16) shows for different values of γ the aspect ratio of the end products ofthe simulations as a function of the aspect ratio of their initial conditions. In theinterval of initial aspect ratio 0 . ≤ α ≤
10 shown here, the aspect ratio of theend products of the expanding initially prolate γ − models takes for every γ largervalues with respect to those reached by spheroids of uniform density starting fromthe same α , by contrast if the model is initially oblate, α ∞ falls always under thevalue attained by the correspondent homogeneous spheroid independently of γ .As a general trend, at fixed α , α ∞ is smaller for larger values of γ if α > α <
1. For a very prominent initial density cusp at thecentroid 2 ≤ γ <
3, the relation α ∞ vs α seems to flatten to limit values both athigh and low values of α (not shown here). On the other hand, the differentialenergy distributions of the end products, shown in Fig. 3.17 are qualitativelysimilar to those of the correspondent homogeneous spheroids when γ ∼ < γ <
2, have quite complex finaldifferential energy distributions with several changes of slope, both for initiallyprolate and oblate initial conditions due to the interplay between shell crossingand time dependent aspect ratio. In the spherical cases, n ( E ) are indeed very.5. Summary similar to those shown in Chapter 2 for other non homogeneous radial densityprofiles. For large central density cusps (2 ≤ γ < α , n ( E )peaks at low energies and falls to high energies with similar power law decays alsofor α = 1.Remarkably, see Fig. 3.18, for initially very flattened ( α < .
5) cuspy systems,the spheroidal symmetry is not retained during the expansion, instead, after few t dyn , the systems assume a toroidal shape, with most of the particles on a thickequatorial ring.Figure 3.18: Projections of the particles positions in the x, y (upper row) and x, z (bottom row) planes for a system of N = 15000 particles starting with α = 0 . a ⊥ . Note how an initialdense core results in particles bunching at a dense equatorial ring. The main results presented in this chapter, in part published in Ref. [65], canbe summarized as follows: by means of a simple semi-analytical model for theCoulomb explosion of a uniformly charged spheroid in the limit of nonrelativisticregime we have studied the expansion of non-neutral aspherical nanoplasmas. Themodel allows one to express quantities such as the maximum energy a particle canreach at a given time, the time-dependent particle energy distributions, as a func- Chapter 3. Coulomb explosion of ellipsoidal systemstion of the initial spheroid aspect ratio α , charge density ρ c, , and total conservedcharge Q .Our theoretical predictions for aspect ratio evolution and final number energydistribution are found to be in good agreement with numerical simulations withboth direct molecular dynamics and particle-mesh approaches. The ”discreteness“peak in the number energy distribution observed in MD simulations of sphericaland homogeneous system is also present, albeit of a lesser magnitude, also forspheroidal systems and in a certain measure for triaxial ellipsoid.In addition we performed simulations of spheroidal systems with non uniformcharge densities finding a similar behavior to that of their spherical counterpartstreated in the previous chapter. However, in all cases the high energy peaks indue to the particle overtaking are less prominent the more the symmetry departsfrom the spherical one.We stress the fact that, albeit highly idealized, a pure Coulomb explosionmodel is indeed useful and to some extent “realistic” in the vertical ionizationregime, where electrons are expelled from the target on a time much shorter thanthe characteristic time of ion motion, attainable using ultra-intense lasers or x-raypulses.Moreover, with the recent progress in nanotechnology, Coulomb explosion ofspecifically designed nanostructured targets can be considered. Our results maythus give us simple design guidelines how to optimize target properties; for exam-ple, for inertial fusion applications or to maximize ion collision events for neutronproduction.Finally, our results can also be helpful to model laser-solid-target interactionfor ion acceleration, which is characterized by the emission of short, compact,and highly charged ion bunches. Propagation of these bunches (e.g., through avacuum) is strongly affected by space charge effects. By approximating the ac-celerated ion bunches as uniformly charged spheroids, the results presented heremay allow us to derive the conditions required for limited energy and angulardispersions. hapter 4 Dynamics of molecular hydrideclusters irradiated by intenseXFEL pulses at LCLS
Atomic clusters have received much attention in recent years as they are tunabletargets for intense laser-matter interaction. This applies to “conventional” laserpulses as well as to VUV and X-ray pulses available from new and upcoming free-electron laser machines.Molecular clusters add another degree of freedom and may thus be a tool toapproach the radiation damage processes of large organic molecules exposed toX-ray radiation and they are as well interesting with respect to the possibility todrive fusion of deuterium via Coulomb explosion of deuterated samples, see [24][166], [6] and [5].In this chapter we apply the results on the heterogeneous clusters presented inthe previous chapters to study the real case of hydrides clusters (i.e. clusters ofmolecules containing an element of the first row and one or more hydrogens) irradi-ated by short and intense X-ray pulses causing multiple K − shell photoionizations.In particular, we discuss, in the light of our model, the experimental findings atSLAC’s Linac Coherent Light Source (LCLS) on methane clusters exposed to softX-ray (¯ hω ∼ Recent experiments [86] performed with the x-ray Free Electron Laser (xFEL)of the LINAC Coherent Light Source in Stanford, where jets of methane clustershave been irradiated with short pulses ( ∼ −
100 fs) of soft x-rays, showed a largeproton signal in contrast to an almost vanishing yield of carbon ions or molecularfragments containing carbon in the time of flight spectrum, as it is evident in Fig4.1 where the different curves show the fragments yield for different cluster sizesin the jet. Chapter 4. Dynamics of molecular hydride clusters irradiated by intense XFELpulses at LCLSSuch effect appears to be independent on the cluster size, in the range observedFigure 4.1: Time of flight spectrum in microseconds for the fragments of (CH ) N clusters with N = 10 , , and 10 . irradiated by short X-ray pulses (roughly20 femtoseconds) with ¯ hω = 1 keV. Figure taken from Ref. [86].here, and hints towards a charge migration from the two species in the molecules.Photons are prevalently absorbed via K − shell photoionization in carbon, that atsuch photon energies has a cross section of the order of 250 kilobarn at photonenergy of 1 keV, almost two orders of magnitude larger than the cross section forphotoionization of the valence electrons of the molecule, (see Fig. 4.2, see also[12]).Although the ionization mechanism is different, what is observed here is aspecies segregation like that seen for the same cluster irradiated with VUV sources[47] and [160], where the light protons overtake the heavy ions once reaching highvelocities. What is puzzling is the low signal of the carbon ions implying that alarge fraction of the carbon component emerges as neutrals. Prompted by these findings, we performed a parametric numerical study of theinteraction with short X-ray pulses of molecular clusters, in order to investigatethe effect in a more general frame. We do not restrict ourselves to the case ofmethane clusters, but we consider the series of first row hydride molecules H O,NH and CH . As an “atomic limit” of our system we use rare gas clusters of neon.2. The model Figure 4.2: Total photoionization cross section for a carbon atom (black curve)and for hydrogen (red curve) as function of the photon energy in eV. The Carboncurve has always larger values than that for hydrogen and it is dominated by the K − shell contribution for ¯ hω larger than roughly 300 eV. In a CH molecule, thecross section of the K − shell of carbon does not significantly differ from that ofthe atomic species, while it is reasonable to assume that the cross section of thevalence shell (molecular orbitals) is of the same order as that of atomic hydrogen.(Ne) since such species has the same electronic configuration of the aforementionedhydrides and comparable K − shell ionization cross section. We report here the keyassumptions and summarize the numerical model detailed in Chapter 5. Setting the stage
In the same spirit of Refs. [59], [62] and [60], we use a hybrid quantum-classical model where ions and electrons are considered as classical particles and propagatedwith simple N − body schemes, while quantum processes, such as photoionizationor Auger decays, are treated with the Monte Carlo method based on rates com-puted from photoionization cross sections, laser intensity and energy amplitudesof the transition.Such approach, contrary for instance to that used in Ref. [13] where electronswere approximated by a continuum fluid, allows us to obtain a better descriptionof the electron component for processes such as Auger decay or recombination as Chapter 4. Dynamics of molecular hydride clusters irradiated by intense XFELpulses at LCLSwell as to book-keep the electronic shell occupation of the ions, and in addition,to stick to easier to implement propagation algorithms.In this study, we assume theoretical X-ray pulses with Gaussian time envelopeand uniform space envelope. The latter choice is motivated by the fact that in theexperiments discussed above, the laser focus diameter is ∼ µ m and the typicalcluster radius is of the order of a few nm. We fix the photon energy ¯ hω = 1keVwhile we span a range of full-width-at-half-maximum for the time envelope, from1 to 50 fs. Peak intensities I span the interval 10 − W / cm . The K − shellphotoionization cross sections σ and other atomic parameters of the elementshere considered are summarized in Tab. 4.2.For such laser parameters, the interaction is always in the so called weak fieldregime , corresponding to values larger than unit attained by the Keldysh param-eter [50] γ K ≡ (cid:115) ∆ E E pond , (4.1)where ∆ E is the ionization energy for the K − shell electron and the ponderomotiveenergy E pond is given as E pond = e I m e cω , (4.2)where e and m e are the electron charge and mass, c is the speed of light and ω is thelaser carrier frequency. For the range of intensities explored here γ K (cid:39) .
5, thatis strictly larger than 1, implying that we are in the regime of non-perturbativesingle photon tunneling ionization, (see, e.g. Ref. [145] and references therein).From the outcome of the experiments with methane clusters, it is evident thatthe occurrence of a strong hydrogen peak in the spectra is not influenced by thecluster size. We have considered therefore only clusters containing a number ofmolecules (or atoms) in the range N ∗ = 50 − t and t + ∆ t for the j -th shell of the i -th atom or moleculein the cluster, as well as the probability of an Auger transition, are computed withthe schemes discussed in Sect. 5.3, while, due to the low size interval consideredhere the dynamics is resolved with a direct scheme where the Coulomb interactionpotential is smoothed at distances r ≤ a with a linear spline, see Sect. 5.2. The choice of the initial conditions
The photon energies considered here are much larger than the typical binding en-ergy of the electrons in the molecular orbitals, which are of the order of 4.5 eV,as well as of that the intermolecular van der Walls bonds in the cluster (of theorder of 0.5 - 2 eV). In addition, the time scales of the exposure to the laser pulse( ∼
10 fs) and therefore of the charging process, are shorter than any chemistryrelated process in the cluster, we thus neglect completely the latter. Therefore,.3. Numerical simulations and results all simulation particles are initialized with v = 0.The “molecules” are simply implemented by placing their constituting ele-ments at equilibrium distance. To include the effect of field ionization (i.e. atomsor molecules are further ionized by the global cluster electrostatic field), each atomhas associated an electron represented by a real simulation particle that is initiallyplaced at the same position. Note that smoothing the Coulomb potential impliesthat the latter electron is confined by an harmonic-like potential when placed at r ≤ a from the mother atom. With respect to the Auger processes, the effectiveelectron (if still present) has always the priority on the “virtual electrons” (i.e.those not yet initialized but just accounted in the electronic configuration of theion) for being emitted as Auger electron, whenever such event takes places.Atoms or molecules are displaced in the clusters according to the typical trun-cated icosahedral shell structure, see also [60] and [59] and references therein, sothat on each shell the charge density is constant.Specie E bind [eV] σ [kbarn] τ moleculeauger [fs] τ atomauger [fs]CH O 543.1 121.9 4.46 5.07Ne 870.2 248.2 2.90Table 4.1: K shell parameters of the heavy atom in the molecular and atomicspecies used in the simulations. Cross sections σ and binding energies E bind areassumed to be identical in the two cases while the average K shell hole lifetimessensibly differ, therefore in the last column we put as a reference the value of τ auger for the pure atomic species. All values listed here refer to neutral systems. Previous numerical studies have addressed the problem of atomic or compositeclusters exposed to strong laser pulses, aiming either at the description of ionsdynamics (such as for instance in [13], see also the simplified models introducedin Chapter 2) or at a detailed description of the electron spectra, at expensesof a more simplified treatment of the ion motion, that in some studies are evenapproximated by a smooth charge background or jellium model , see e.g. [62]. Thereason of this being that ions move significantly on time scales considerably largerthan that of the electrons. The typical kinetic energy of a K − shell photoelectronextracted from a carbon atom by a 1 keV photon is of the order of 700 eV,implying that it crosses a distance of 50 atomic units of distance (radius of atypical nanocluster) in roughly 0.2 femtoseconds, while the typical ion motionscale is of the order of ten femtoseconds.In this thesis we tackle the problem with a rather simplified approach, in a Chapter 4. Dynamics of molecular hydride clusters irradiated by intense XFELpulses at LCLSway that the interplay of the ionic and electronic components can be treated asfunction of few tunable parameters. Moreover, using an approach based on rateequations and Monte Carlo sampling, allows one to model the coupling to thelaser field in a computationally “cheaper” way than a full quantum treatment.We have studied both the time evolution and the dependence on the lasersintensity of global quantities such as the total charge of the system and the averageand maximal energies of the different species of particles, as well as spectra suchas the charge states of the ions or their energy distribution n ( E ). The interplay between photoionization and Auger decay, qualita-tive picture
Photoionization and consequent Auger decay produce two families of electrons.Their initial energies may in principle have a complicated spectrum since differentshells can be ionized and many channels of Auger decay are possible with differentenergies. However, photons of a few keV energy ionize mainly the K shell (1sorbitals) for the elements of the first row under consideration. In our test simula-tions for CH clusters, less than 3% of the photoelectrons come from the valenceshell at peak intensities I > W / cm , making reasonable to neglect its contri-bution. Hence, for the sake of simplicity, all results presented have been obtainedwith photoionizing exclusively K − shell electrons, this means that photoelectronscan be produced with two possible energies only, one relative to the case of thefull K shell and the other for the case of the one fold ionized K shell.On the other hand, the typical Auger life-times of K shell vacancies, see Tab.4.2, are always shorter than the pulse length assumed in this study, so that inprinciple autoionization may refill the inner shell within the time the laser pulseis effective, in a way that the two processes are strictly entwined.Assuming for the pulse a Gaussian time envelope and constant space profile,means that in principle the number of absorbed photons per atom n γ theoreticallyscales with time with an error function n γ ( t ) = c ∗ I σ ¯ hω (cid:114) π (cid:20) erf (cid:18) t √ σ (cid:19) + 1 (cid:21) , (4.3)where σ = T / √ c ∗ is to be inter-preted as a constant containing the information relative to the number of atoms(or molecules) in the cluster N ∗ , the photoabsorption cross section and the Augerlifetime τ . Note that, this is valid only for combinations of cluster sizes and inten-sity such that no saturation takes place (i.e. the total number of absorbed photonsis smaller than N ∗ ). Equation (4.3) can be used in principle for a one parameter( c ∗ ) fit of the correspondent numerical curve as it is shown in Fig. 4.3 where thequantity n γ extracted from the realization average of 20 runs with identical initialconditions is perfectly fitted by its theoretical counterpart.It is evident in this specific case, how Auger decay sets in early within thepulse, so that at t = 0, the time when the peak in intensity is attained, half of.3. Numerical simulations and results fr ac ti on o f e v e n t s time [fs] n g (fit)n g (simulation)n auger (simulation) Figure 4.3: Number of absorbed photons (and released photoelectrons) as functionof time (thin red line) fitted with Eq. (4.3) (heavy black line) and number ofAuger electron produced by the decay of the inner shell vacancies (blue curve) ina (CH ) cluster irradiated by a pulse with I = 10 W / cm , T = 10 fs andphoton energy of 1 keV so that c ∗ I σ/ ¯ hω ∼ .
22. The gray shaded area representsthe time envelope of the pulse.the molecules that have been photoionized also produced an Auger electron. Forthe species of interest, once released such electrons have typical kinetic energies K auger corresponding to the ∆ E between unstable and decayed configurations, ofthe order of 250 eV, allowing them to escape the cluster if K auger > eQ t R t , (4.4)where Q t is the total charge of the ionized cluster when its radius increased to thevalue R t .The K − shell photoabsorption cross sections σ differ little between the atomicand molecular species (here we have assumed them to be identical). Therefore,the different Auger lifetimes of the molecular and atomic species marks the mostimportant difference between molecular and atomic clusters exposed to the sameX-ray pulses, since this implies that the inner shells are “refilled” on different timescales in the two cases. However, for given photon energy and pulse length, theeffect becomes appreciable only at large intensities I , as it is clearly evident fromFig. 4.4.At fixed pulse length, the higher the laser intensity is, massive photoionization Chapter 4. Dynamics of molecular hydride clusters irradiated by intense XFELpulses at LCLS -1 n γ I [W/cm ] (C) N * (CH ) N * (N) N * (NH ) N * (O) N * (H O) N * Figure 4.4: Number of absorbed photons per unit n γ for molecular systems (heavylines) and parent atomic systems (thin lines) for different laser peak intensities.For values of I smaller than roughly 8 × W / cm , no difference between thetwo cases can be noticed. Each point of the curves is averaged over 30 statisticalrealizations with identical initial conditions.in the cluster and subsequent decay of the K − shell holes, cooperate in building inboth atomic and molecular systems, a deeper cluster potential. Having assumeda constant particle density on every shell, and since the ionization occurs isotropi-cally (photoionization cross sections σ of the elements here considered (cfr. Tab.4.2)), we can idealize it being harmonic for r < R t and Coulombian for r > R t .The cluster electrostatic field induces the so called field ionization, whenever ithappens to be strong enough to strip the ions of their external electrons. Theeffect is in general more effective on the ions sitting at the cluster’s surface, wherethe electric field is stronger.Electrons produced this way are typically confined by the space charge of theions and have kinetic energies lower than eQ/R . Molecular cluster versus pristine atomic clusters
Figure 4.5 illustrates the time evolution of characteristic parameters total chargeper unit
Q/N ∗ , radius R and “pseudo temperature” of the confined electronnanoplasma (i.e. the average electron kinetic energy E el ), for an atomic clus-ter (C) and a molecular cluster (CH ) under the influence of the laser pulse(gray shaded area). It is immediately clear that the dynamics of the pristinecarbon cluster (dashed line) and the methane cluster (solid line) is completely.3. Numerical simulations and results Figure 4.5: Top: time dependent total charge Q per unit enclosed by the surfaceof the ion core of a (CH ) N cluster (solid line) and total charge per unit enclosedby the surface of an atomic carbon cluster (dashed line) both irradiated by agaussian pulse with I = 10 W / cm , T = 10fs and photon energy of 1 keV. Thegray shaded area shows the time envelope of the pulse. Middle: radii of the ioncore (blue solid line) and proton shell (red solid line) and radius of the pure carboncluster in units of the initial cluster radius R . Bottom: average kinetic energy ofthe trapped electrons in the cluster in the atomic and molecular case.different.We observe that, as one would expect, the carbon atoms get successivelycharged through photoionization leading to more than 90% carbon ions (toppanel). The cluster ions create a deep binding potential from which most Augerelectrons can not escape, forming a nanoplasma together with the field ionizedelectrons.The maximum kinetic energy of the trapped electrons is limited by the depthof the cluster potential and the average kinetic energy (bottom panel, dashed line)is indicative of the nanoplasma temperature . This is the normal behavior as is Chapter 4. Dynamics of molecular hydride clusters irradiated by intense XFELpulses at LCLS n ( q ) / N * qno focus averaging(C) (CH )
0 1 2 3 4 5 6qfocus averaging
Figure 4.6: Left panel: charge states at t = 500 fs for the same systems and laserparameters of Fig. 4.5 normalized respect to the total number of heavy ions in thecluster N ∗ . Right panel: same quantity as left but spatially averaged over focus.well known from rare-gas clusters exposed to X-ray pulses, see e.g. Refs. [145]and [61].The molecular cluster, however, does not follow this scheme: While initiallysimilarly charged as in the pristine cluster, the total charge Q ( R ) enclosed by theradius of the core experiences a dramatic drop at around 9 fs, so that this regionis in the end almost neutral on average (Fig. 4.5 middle panel, solid line). At thesame time, the kinetic energy of the trapped electrons and, hence, the temperatureof the nanoplasma remains comparatively low (Fig. 4.5, bottom panel, solid line).Both phenomena originate in the ejection of fast protons from the molecularcluster (see upper (red) line in middle panel of Fig. 4.5). Although the carbon K − shells are initially photoionized, the charge distribution of the doubly chargedmethane after Auger decay is such that the carbon ion is screened and the posi-tive charge is dominantly localized on two hydrogen atoms which are likely to beejected from the entire cluster as protons. These protons take away the excesspositive charge created by photoionization which is, of course, not possible in the It shall be pointed out that technically speaking, one can not in principle speak of electron“temperature”, as the time scales at which we are looking at the systems might be considerablyshorter respect to the typical time scales at which the electron plasma attains a Maxwellian energydistribution due to Coulomb collisions [142], for our combinations of density, total cluster chargeand heavy ions mass. In addition, it is important to mention that confined Coulomb systemsand in general systems of particles interacting with long range forces undergo other processes ofenergy relaxation due to collective oscillations and other non collisional phenomena, see e.g. [87],[105] and [43] .3. Numerical simulations and results pristine carbon cluster. The remaining positive charge in the cluster is small giv- -3 -2 -1 K m a x / ε t o t n γ NeH ONH CH Figure 4.7: Kinetic energy of the fastest ion K max , 0.5 ps after the peak of the pulse( T = 10 fs), versus the average number of photons absorbed per atom/molecule n γ . Cluster size is N ∗ = 689. The solid line marks the “atomic limit” of the neoncluster. Note that the kinetic energies are normalized in units of the total energyabsorbed by the cluster E tot .ing rise to a weak potential which can only trap low-energy electrons. Therefore,the temperature of the nanoplasma is by a factor of 6 smaller than for the purecarbon cluster after t ∼
65 fs; this holds true also for the Coulomb explosion of thecarbon ions with the charging of the carbon ions even more dramatically reduced,almost by an order of magnitude. In Figure 4.6 the charge states spectrum definedas n ( q ) = N ∗ (cid:88) i =1 δ qq i , (4.5)are presented for the methane and pure carbon clusters at 0.5 ps after the peak ofthe pulse. In the left panel the spectra are averaged over 50 realizations assumingconstant space envelope, while in the right panel a Lorentzian space envelope isassumed at the pulse’s focal spot (see Eq. 5.47), and the spectra are averagedover 200 realizations distributed in space with normal distribution. In both cases,for the molecular cluster almost 85% of the carbon emerges as neutrals and therest is singly charged. In the spectrum of the pure atomic system, up to six fold Chapter 4. Dynamics of molecular hydride clusters irradiated by intense XFELpulses at LCLScharged ions are present when the space envelope is assumed constant, while only n γ H ONH CH < K > H / < K > X I [W/cm ] Figure 4.8: Top panel: fraction of absorbed photons per molecule n γ a functionof the x-ray intensity I for the same systems shown in Fig. 4.7. Bottom panel:Ratio of the average kinetic energy of protons and heavy atoms.up to five fold charged ions are observed in the other case, but still it appears tobe qualitatively different than the spectrum of the molecular cluster.The radical difference between the final charge states of carbon ions in molec-ular and atomic clusters is not imputable to a different photoabsorption (see Fig.4.4), but instead to a more efficient recombination of the nanoplasma electronsfor the molecular system due to the presence in the molecular cluster of a largernumber of plasma electrons per ion with lower average kinetic energy.As a general remark, for the other cases of NH and H O the effect is qualita-tively similar to what we have discussed for the methane clusters; as we will seein the following, the difference in behavior of the atomic and molecular clustersbecomes more marked for intermediate intensities. A time scale of a few hundreds of femtosecond is not sufficient to equilibrate electrons andions, the latter are characterized by a non-thermal energy distribution. The recombination hasto be intended here as electrons being classically bound to the neighbouring ion. .3. Numerical simulations and results Intensity and pulse length dependence of the ion dynamics
In the limit of short pulses and extreme intensities the dynamics of the ions andprotons is obviously expected to be qualitatively similar to the idealized case ofmulti-component Coulomb explosion treated in Sect. 2.3. For a given pulse length T , one may expect that the absolute difference in velocity of heavy and light ionsgets larger with increasing intensity and, therefore, higher charging of the cluster.Figure 4.7, however, reveals that the ratio of the kinetic energy for the fastestheavy ion K max , in units of the total energy absorbed by the cluster E tot = ¯ hωN ∗ n γ , (4.6)exhibits a non monotonic behavior as a function of the fraction of photons absorbedper atom/molecule n γ with a dip at a certain critical number.The latter depends moderately on the species considered, but is otherwise auniversal feature of hydride clusters in obvious contrast to the isoelectronic atomicneon cluster (black solid line in Fig 4.7). Such feature in the curves for the hydrideclusters is a signature of the dynamical segregation of protons and heavy ions ascan be seen in bottom panel of Fig. 4.8, where a more global quantity, namely,the ratio of the average energy of all protons versus that of all heavy ions R = (cid:104) K (cid:105) H (cid:104) K (cid:105) X , X = O; N; C (4.7)is shown as a function of peak laser intensity I . Again, one sees a qualitativelysimilar non monotonic behavior of all three hydrides, even though the photo-absorption goes monotonically in all three cases (top panel same figure). Thereis a maximal segregation at a well-defined intensity for each of the three hydridessince they differ in their respective ionization energy for the 1s electrons, Augerrates, and photoionization cross sections, as listed in Tab. 4.2. Hereby, the shiftof the peak positions is in principle due to the different photoionization crosssections.The final ion charge states n ( q ) taken at t = 500 ps are shown in Fig. 4.9 forthe three hydrides cluster as well as their atomic counterparts, for N ∗ = 689 and I = 2 × , × and 5 × W / cm . As expected, a much lower chargingof heavy ions as compared to the pristine cluster of the heavy-atom species at thelowest and intermediate intensities.At I = 2 × , most heavy atoms remain neutral or singly charged in thepristine as well as in the hydride cluster. This changes drastically for intermediateintensities of about 10 W / cm , where the fraction of neutral atoms survivingthe light pulse illumination is small to vanishing in the pristine cluster. In thehydride cluster (CH and NH cases), on the other hand, about 80% neutral heavyatoms result from recombination with the cold electrons after proton segregationin the surface layer, which has been fully charged due to efficient field ionization.The water cluster has an almost negligible fraction of neutrals ( ∼ Chapter 4. Dynamics of molecular hydride clusters irradiated by intense XFELpulses at LCLS n ( q ) / N * q 0 1 2 3 4 5 6q 0 1 2 3 4 5 6q n ( q ) / N * q 0 1 2 3 4 5 6 7 8q 0 1 2 3 4 5 6 7 8q n ( q ) / N * q 0 1 2 3 4 5 6 7q 0 1 2 3 4 5 6 7q I = 2 × W / cm I = 1 × W / cm I = 5 × W / cm Figure 4.9: Normalized charge-state distribution n ( q ) /N ∗ for carbon ions C q + from(CH ) (top row), for nitrogen ions from (NH ) (middle row) and oxygen ionsfrom (H O) (bottom row). The black bars mark the case of the correspondingpristine atomic systems, and the x-ray intensities I are indicated on top.spectrum is shifted towards lower charge states than that of pure oxygen cluster..For higher intensities, we expect the proton segregation to cease (as seen inFig. 4.7) and, as a consequence, similar charge spectra for the pristine and thehydride cluster. This is indeed true with respect to a vanishing yield of neutral.4. A synthetic model atoms. However, the form of the charge distribution is still somewhat different. InFigure 4.10: From left to right for (H O) N ∗ , (NH ) N ∗ and (CH ) N ∗ , the ratio R of the average kinetic energy of protons and heavy atoms is shown as function ofthe laser’s peak intensity I and length T is shown. In this case N ∗ = 689.particular, NH and H O cluster present even larger percentages of higher chargestates than their parent atomic clusters, in contrast to what observed for CH and C clusters. The latter is due to the fact that at I = 5 × W / cm forthese species, the number of absorbed photons in the atomic and molecular caseis considerably different (cfr Fig. 4.4), and in addition less protons are “available”to carry away the excess positive charge.Figure 4.10 shows that the dynamical segregation of protons and heavy ionsdoes not happen only for the pulse length considered here, but takes place fora broad interval of pulse lengths. The peak in the ratio R broadens and shiftstowards lower I for increasing T , this is a signature of the fact that the emergenceof such segregation effect is principally due to the amount of charge created in thesystem. From the analysis carried out to this point it is clear that the proton segregation isstrongly influenced by the laser’s intensity for fixed pulse length, as it is revealedby the total charge yields and particle-averaged kinetic energies. However, it isstill needed to clarify whether this segregation is a local effect due to the heavy-light character of the hydride molecules or whether is a consequence of the clusternature of the entire system. In order to clarify this, we have set up a simple modelwhere N = 10 singly charged ions are distributed homogeneously in a sphere ofinitial radius R , neither supporting any molecular substructure nor allowing anyintra-atomic or intramolecular electronic processes. Three fourths of the ions havethe mass of the proton, while the other one fourth has a 20 times higher mass(roughly the mass of Neon). N − Q electrons are placed at randomly selected Chapter 4. Dynamics of molecular hydride clusters irradiated by intense XFELpulses at LCLSions. With this initial configuration, ions and electrons, interacting via Coulombforces regularized at short distance, are propagated up to 1 picosecond. < e > l / < e > h Q/NA B C
Figure 4.11: Ratio of the average energies of light (cid:104)E(cid:105) l and heavy (cid:104)E(cid:105) h ions withscenarios A, B and C discussed in the text.Because of its net positive charge, the system undergoes expansion. The ratioof the average final energies of light and heavy ions, exhibits one central maximumsee Fig. 4.11 (region B, 0 . < Q/N < .
25) similarly, as for the fully microscopiccalculations for hydride clusters discussed in the previous Sections. In Figure 4.12we show the final energy of the light and heavy ions as a function of their initialradial position r/R . It is clear how the dependence of the final energy of an ionon its initial radial position in the cluster influences its final energy and inducesthe heavy-light ion segregation.In region A ( Q/N ≤ . Q/N ≥ . Q is sufficiently strong totrigger field ionization of surface ions, as it has been discovered for homogeneousclusters [60]. As a consequence, one expects the cluster core to be screened by thefield ionized electrons up to the radius indicated by the green line in left and rightpanels of Fig. 4.12.However, the protons in the hydride cluster core are light enough (or more pre-cise: have a sufficiently large charge-to-mass ratio) to have started moving beforethe repelling forces are compensated by the screening electrons. Hence, protonsleave the cluster core and, as a result, the surplus of screening electrons preventsheavy ions in the surface layer even beyond the screening radius of a homogeneous.4. A synthetic model Figure 4.12: Normalized ion kinetic energies E ∗ of a heterogeneous cluster modelwith N = 10 particles in a sphere of initial radius R as a function of the totalcharge Q = N of the system and their initial position in units of R . The energiesof the heavy ions (bottom panel) are normalized to Q/R while those of the lightions (middle panel) to Q/ R . The thin green line marks the part of the clusterwhich is supposed to explode if the charge would concentrate in a shell at thesurface. The top panel reproduces the ratio R also shown in Fig. 4.11. Figureadapted from Ref. [45] Chapter 4. Dynamics of molecular hydride clusters irradiated by intense XFELpulses at LCLScluster from exploding. In contrast, protons escape with high final energy fromthe surface layer.In region C, the initial charge Q further increases which weakens the screeningeffect through field ionized electrons for two reasons. First, the fraction of screen-ing electrons available versus initial charge Q decreases. Second, the temperatureof the screening electrons is higher than in B due to the deepening of the trappingpotential with increasing Q . Hence, the surface layer no longer forms efficientlyand a Coulomb explosion, as in region A, although more violently, results.We may conclude that heavy-light segregation is not a local effect of heavy-light molecules. Rather, it happens in a surface layer of the heterogeneous clustertriggered by field ionization.In contrast, “dynamical acceleration” in such clusters see e.g. [76] and [98], doesnot require electron screening and relies exclusively on ion-ion repulsion through-out the cluster. We now summarize this chapter, whose results are also published in Ref. [45].First of all, it must be pointed out that the behavior of molecular hydride clustersis qualitatively very different from pristine atomic clusters as the first have anadditional channel of energy loss, namely the light protons, to release the energyabsorbed from the laser pulse.As a general qualitative fact, there is a higher percentage of neutrals and lowcharge states in the final n ( q ) of the heavy ions in the molecular case than in thepure atomic case.The dynamics of two species of ions is strongly dependent on the laser intensity I .Three regimes can be distinguished where, the weakly charged cluster expands as awhole as a quasi-neutral plasma (hydrodynamic expansion), the cluster is chargedenough do induce significant field ionization at the surface, electrons collapse tothe core neutralizing it and a charged shell of protons explodes faster carryingaway most of the energy absorbed by the system, and finally the cluster almost isentirely charged and behaves as in the limit of multi-species pure Coulomb explo-sion.The intensity dependent heavy ion-proton segregation in the intermediateregime is found to be a universal feature of hydride clusters and can be classi-fied as a cluster effect.In conclusion, we speculate that, due to the experimentally relevant intensitywindow where such segregation effect happens, the results discussed here couldturn useful for the study of radiation damage of large bio-molecules as it may oc-cur during coherent diffractive imaging with intense X-ray pulses, see Refs. [126],[164], and [70, 31]. hapter 5 Numerical methods
In this chapter we present the different numerical methods used in this thesis tostudy the dynamics of finite size plasmas generated by ionization of cluster targets.In addition we discuss the techniques based on Monte Carlo sampling and rateequations, used to treat photoionization, and the Auger decay and recombinationprocesses.It must be also pointed out, that many of the concepts concerning numericalsimulations using particles discussed here, also apply in the field of gravitational N − body simulations, as the Coulomb electrostatic force and Newtonian gravityshare the same 1 /r nature.In particular, the codes used for the simulations discussed in this work, anddeveloped by the Author during its Ph.D., are presented in detail with the resultsof some of their performance tests. N − body methods for systems of charged particles As previously introduced in Chapter 2, in the simple electrostatic picture wherethe contribution of self-induced magnetic fields and radiative losses are neglectedand the velocities are always non relativistic, a plasma can be thought as anensemble on N particles of charge q i , mass m i with positions and velocities x i and v i , interacting via Coulomb forces and it fully described by the Hamiltonian H = N (cid:88) i =1 p i m i + 12 N (cid:88) j (cid:54) = i =1 q i q j || r j − r i || (5.1)where p i = m i v i . Given the initial conditions R = ( r , r , ..., r N ) t =0 and V =( v , v , ..., v N ) t =0 , the configuration of the system at time t is formally obtainedby integrating the system of 6 N coupled ODEs (Hamilton equations)d r i d t = ∂H∂ p i ; d p i d t = − ∂H∂ r i ; i = 1 , N. (5.2) Chapter 5. Numerical methodsIn the limit of large N and vanishing inter-particle correlations, the evolutionof the system is that of its one-particle phase space distribution defined on the6-dimensional one particle phase space f ( r , v , t ), through the Collisionless Boltz-mann Equation (CBE), or Vlasov Equation in the jargon of Plasma Physics, seeAppendix C, D f D t ≡ ∂f∂t + v · ∂f∂ r + ∇ Φ · ∂f∂ v = 0 . (5.3)The electrostatic potential Φ( r ) is related to f by the Poisson equation∆Φ( r ) = 4 πρ ( r ) = 4 πq (cid:90) Ω f ( r , v )d v , (5.4)where Ω is the region of the phase space occupied by the system and q is the unitcharge .In this picture what one has in principle to solve is a system of two PDEs(Vlasov and Poisson equations) and an integral equation for the density calledVlasov-Poisson system. When the effect of (external and self-induced) electro-magnetic fields and the coupling of matter with radiation are taken into account,the term ∇ Φ · ∂f /∂ v in Eq. (5.3) is substituted by the full Lorenz force per unit ofmass F = q/m ( E + v /c × B ) and 3 of the 4 Maxwell’s equations have to be addedto the Vlasov equation. In this case one speaks of a Vlasov-Maxwell system.Numerical models based on Equation (5.3) deal essentially with the discretiza-tion of the DF f ∗ ( r , v , t ) = n (cid:88) i =1 δ [ r − r i ( t )] δ [ v − v i ( t )] (5.5)where n is the effective number of macroparticles used in the simulation that issome orders of magnitude smaller than N , and the electrostatic potential at r i isgiven by Φ( r i , t ) = (cid:90) g ( r i − r )d r (cid:90) f ∗ ( r , v , t )d v = n (cid:88) j (cid:54) = i q j g ( r i − r j ) , (5.6)where now q j is the charge of the j − the macroparticle and g ( r i − r j ) = 1 / | r i − r j | is the Green’s function of the Laplacian operator ∆, and | ... | is the standard Eu-clidean norm. Given Φ( r ) and the acceleration as a ( r ) = ∇ Φ( r ), what a numeri-cal code does, is essentially to propagate in time with finite difference integrationmethods (see e.g. Ref. [137] for a complete description) the n particles, by solvingtheir equations of motion, that is substantially returning to the picture of Eqs.(5.2).What characterizes the different plasma (as well as gravitational) N − bodycodes is the scheme that allows one to compute the electrostatic potential Φ (aswell as the electric fields). The main families of codes based on particles are: Note that this approach can be obviously extended to multi-species plasmas where eachspecies of particle of mass m i and charge q i has its own f i and Φ is computed by solving thePoisson equation due to all the species. .1. N − body methods for systems of charged particles N − body codes This category of codes (see e.g. Ref. [1] for a complete review), also known as particle-particle codes (PP) uses the conceptually simplest approach, the potentialΦ i = Φ( r i ) and the acceleration a i at the position of the i − th simulation particleare computed directly, summing over the contributions of all the other particlesof the system as Φ i = q i N (cid:88) j (cid:54) = i =1 q j || r i − r j || , a i = q i m i N (cid:88) j (cid:54) = i =1 q j ( r i − r j ) || r i − r j || . (5.7)Such operation involves for each step two nested cycles over the number of parti-cles implying that the number of operations (numerical complexity) for the forcecalculation, and hence the computational time, scales with N as O ( N ). Becauseof that, codes based on direct force computation are not suitable to treat morethan ∼ × particles, even on modern processors, since for larger systems thetime taken by the force computation, that is the bottleneck of the scheme, be-comes prohibitive even within a single time step ∆ t .However, if one is interested in simulating systems of small size where actuallythe dynamical collision among particles play, an important role, the direct codeshave in this case the advantage that the forces on the particles can be better es-timated (depending only on the machine precision) than with the other methodsdiscussed next.There is unfortunately an additional drawback, due to the singular nature ofthe Coulomb interaction for vanishing separation, the potential and the force com-puted for a couple of particles placed at a very small distance is large and therefore,the velocity change attained within a single time step can be spuriously high. Inaddition, due to the finiteness of the machine representable reals, the pair forcemay diverge even for no zero separation. As we will see in Sect. (5.2), these com-plications are generally circumvented by modifying the Coulomb potential at smallseparation by the introduction of the softening length (cid:15) , (see e.g. [40] and [41], seealso [114]), so that in Eq. (5.5) q i δ [ r − r i ( t )] is replaced by D [ r − r i ( t )] where D is a distribution with scale length (cid:15) that smoothes q on a finite or infinite support.Consequently, the softened potential that it generates satisfies ∆Φ soft = 4 π D , andthe softened acceleration is given as usual by a soft = −∇ Φ soft . Particle-Mesh codes (PM)
A second category of N − body schemes, sometimes also dubbed Particle in Cell(PIC), used when the number of particles in the simulation larger than 10 andwhen collisional effect are negligible, that is extensively discussed in Ref. [75], isso-called particle-mesh. In this case, Φ is first computed using ∆Φ = 4 πρ on a Chapter 5. Numerical methods3-dimensional cartesian grid superimposed to the system where the density ρ is dis-placed, and then the potential and its gradient are interpolated at the particles po-sitions. To assign the density ρ i,j,k to a given mesh points P i,j,k first a shape func-tion S ( x, y, z ) is introduced so that the contribution W P ( x l − x i , y − y j , z − z k ) ofthe l − th particle to the mesh cell P i,j,k is W P ( x l − x i , y − y j , z − z k ) = (cid:90) x i +∆ x/ x i − ∆ x/ (cid:90) y j +∆ y/ y j − ∆ y/ (cid:90) z k +∆ z/ z k − ∆ z/ S ( x l − x i , y − y j , z − z k )d x l d y l d z l , (5.8)where ∆ x, ∆ y and ∆ z are the mesh spacings along the three axis. The density ρ i,j,k then reads ρ i,j,k = 1∆ x ∆ y ∆ z N (cid:88) l =1 q l W P ( x l − x i , y l − y j , z l − z k ) . (5.9)The two main choices of S that are generally used are the nearest-grid-point (NGP)method that assigns each particle’s charge q l and coordinates ( x l , y l , z l ) to thenearest point of the grid with the shape function S NGP ( x l − x i , y l − y j , z l − z k ) == 1∆ x ∆ y ∆ z θ (cid:18) − | x l − x i | ∆ x (cid:19) θ (cid:18) − | y l − y j | ∆ y (cid:19) θ (cid:18) − | z l − z k | ∆ z (cid:19) , (5.10)and the cloud-in-cell (CIC) method that instead associates to the particle a finitesize and a cubic shape with side of length ∆ p , and assigns to P i,j,k a fraction ofcharge corresponding to the portion of its volume overlapping with it. In this case S ( x l − x i , y l − y j , z l − z k ) CIC = 1(∆ p ) b (cid:18) x l − x i ∆ p (cid:19) b (cid:18) y l − y j ∆ p (cid:19) b (cid:18) z l − z k ∆ p (cid:19) , (5.11)where b is the zeroth order b -spline function defined as b ( ξ ) = (cid:40) | ξ | < /
20 otherwise . (5.12)Higher orders b l are obtained recursively via b l ( ξ ) = (cid:90) + ∞−∞ b ( ξ − ξ (cid:48) ) b l − ( ξ (cid:48) )d ξ (cid:48) . (5.13)As an example, b , b and b are shown in Fig. 5.1. Using a higher order b k in(5.11), instead of b refines the charge deposition algorithm. For l = 1 one has theso called triangular-shape-cloud (TSC)..1. N − body methods for systems of charged particles Explicit Particle Simulations of Space Weather
The subsequent b-splines, b l , are obtained by successive in-tegration via the following generating formula: b l ( ξ ) = ! ∞−∞ dξ b ( ξ − ξ ) b l − ( ξ ) (4.12)Figure 4.2 shows the first three b-splines. − − − − b ( ξ ) − − − − b ( ξ ) − − − − ξ b ( ξ ) Fig. 4.2. First three b-spline functions.
Based on the b-splines, the spatial shape function of PICmethods is chosen as: S x ( x − x p ) = 1∆ x p ∆ y p ∆ z p b l " x − x p ∆ x p b l " y − y p ∆ y p b l " z − z p ∆ z p (4.13)where ∆ x p , ∆ y p and ∆ z p are the lengths of the support ofthe computational particles (i.e. its size) in each spatial di-mension. A few PIC codes use splines of order 1 but the vastmajority uses b-splines of order 0, a choice referred to as Figure 5.1: From top to bottom, first three spline functions b k ( ξ ).Once the charge density is placed on the grid with one of the described meth-ods, the main step (that constitutes actually the core of a PM-code) is representedby the solution of the Poisson equation to obtain Φ i,j,k as a function of ρ i,j,k , thathas then to be derived and interpolated at the particle’s positions. Two differentmethods exist to solve ∆Φ = 4 πρ , the finite difference method and the solution inFourier space. In the first, the Poisson equation is discretized and reorganized inthe formΦ i +1 ,j,k + Φ i − ,j,k + Φ i,j +1 ,k + Φ i,j − ,k + Φ i,j,k +1 + Φ i,j,k − − i,j,k = 4 πρ i,j,k (∆ s ) (5.14)where ∆ x = ∆ y = ∆ z = ∆ s (i.e. a cubic grid with equal mesh spacing in alldirections is assumed); then by rearranging the elements of the resulting 3d arraysfor Φ and ρ on two 1d arrays, with the index relabeling l = iN g + jN g + k , where N g is the number of grid points in each direction, one is left withΦ l +2 N g +1 + Φ l + N g +1 + Φ l − (2 N g +1) + Φ l − ( N g +1) + Φ l +1 + Φ l − − l = 4 πρ l (∆ s ) (5.15)that is solved, for instance with the relaxation technique (see e.g. [137]), as alinear system of the form M · a = b (5.16) Chapter 5. Numerical methodswhere M is the tridiagonal matrix of the coefficients, vector b contains the “knowninformation” (i.e. the density) and a the unknown potential.The second method, is based on the fact that in Fourier space the PoissonFigure 5.2: Sketch of the extended grid for the simulation of an isolated object.For simplicity only two dimensions have been represented.equation reads ˆΦ = ˆ ρ ˆ g (5.17)where the hats denote the Fourier-transformed potential, density and Green func-tion of the Laplace operator. The latter in real space is 1 /r , and it is defined onthe N g × N g × N g mesh, as g i,j,k = 1 (cid:113) [∆ x ( i − + [∆ y ( j − + [∆ z ( k − . (5.18)Note that for i = j = k = 1, g diverges. To overcome this inconvenient usuallyone takes g , , = 1 (cid:112) (∆ x + ∆ y + ∆ z ) (5.19)or, alternatively g , , = 23 (cid:112) (∆ x + ∆ y + ∆ z ) . (5.20)Such regularization introduce an effective softening of the Coulomb interaction ona length scale given by the size of the grid step..1. N − body methods for systems of charged particles The simplest approach is again to take cubic grid with ∆ s = ∆ x = ∆ y = ∆ z ,with this choice the discrete Fourier transforms of ρ i,j,k and g i,j,k areˆ ρ α,β,γ = N g (cid:88) i,j,k =1 ρ i,j,k exp (cid:20) − i 2 πN g ( αi + βj + γk ) (cid:21) , ˆ g α,β,γ = N g (cid:88) i,j,k =1 g i,j,k exp (cid:20) − i 2 πN g ( αi + βj + γk ) (cid:21) . (5.21)It is important to state that ˆ g has a diagonal structure in Fourier space, makingthe solution of Eq. (5.17) computationally faster than in the real-space basedapproach. The potential Φ i,j,k is finally obtained by back-transforming ˆΦ α,β,γ =ˆ ρ α,β,γ ˆ g α,β,γ and readsΦ i,j,k = 8 π ( N g ∆ s ) N g (cid:88) i,j,k =1 ˆ ρ α,β,γ ˆΦ α,β,γ exp (cid:20) i 2 πN g ( αi + βj + γk ) (cid:21) . (5.22)Several numerical algorithms scaling in the best cases with N g , as O ( N g log N g ),are nowadays publicly available to solve this problem, and of course machineoptimized routines do exist. However, a detailed description, even of the mostcommonly used, is far outside from the scope of this chapter, so we redirect thereader to the specialized literature.The above discussed method gives the solution of the Poisson equation for an infinitely extended system with periodically symmetric density on a N g × N g × N g grid were N g is a power of 2. This is not our case and to treat isolated systemsthe most common approach is that used for instance in the code superbox , (seeRef. [54] and references therein), where by doubling grid size in each direction,the domain is extended by 7 ghost boxes, similarly as what depicted in Fig. 5.2for the two-dimensional case, where ρ is imposed to be 0 and the Green functionis written as g N g − i,j,k = g N g − i, N g − j,k = g N g − i,j, N − k = g N g − i, N g − j, N g − k = g i, N g − j,k = g i, N g − j, N g − k = g N g − i,j, N g − k = g i,j,k , (5.23)so that it is periodic on the extended grid. The potential Φ i,j,k is computed es-sentially in the same way as for infinite systems. In both cases, potential and ac-celeration at the particles’ positions are obtained by interpolating from grid basedquantities using standard techniques such as cubic splines or Lagrange polynomi-als, see [137].When self induced and external magnetic fields are present, in some PM codessuch as for example calder [134], that we have used in some of our numericalsimulations, the Poisson equation is not solved at every step but instead, usingthe density decomposition technique introduced by Esirkepov in [52], is solved only Chapter 5. Numerical methodsonce at the beginning of the computation, and the Maxwell Equations are coupledFigure 5.3: Sketch of the domain decomposition in a tree code. Again, the thirddimension has been omitted.instead to the continuity equation enforcing in this way the charge conservation.Finally, in some other PM codes, the co called particle-particle-particle-mesh (PPPM or P M), the interaction between the particles is refined by computingalso the direct force due to neighbours inside each cell.
Tree-codes
Although it has not been directly used in this study, we have to mention thislatter family of codes particularly suited for highly clustered or dishomogeneoussystems. Based on 1986 Barnes and Hut’s scheme (see Refs [7], [74] and [75]),the tree codes divide the force and the potential at a particle’s position in twocontributions, due to the neighbours, computed by direct sum and due to thedistant particles, computed instead with multipole expansion. The domain of thesimulation (i.e. the cubic volume space containing all the system, called tree ) isdivided hierarchically, first in 8 cubic cells with sides half the length of the initialcube side, the branches , and then recursively until each particle has its own cubeor leaf , see Fig. 5.3. Cubes containing no particle are discarded from the count. Typically since the system of 4 Maxwell equations is superdetermined , see [83], only two ofthem need to be solved, usually Amp`ere and Faraday equations. .2. The N − body codes used for the simulations The introduction of the dimensionless parameter θ , the so called opening angle (ingeneral, 0 . < θ < . r distance R from the centre of mass r cm of a cell of side l , lR < θ, (5.24)the contribution of the cell to Φ and a in r is computed via multipole expan-sion, otherwise the cell is partitioned in sub-cells that are analyzed with the samemethod, in case that leaves -cells, containing only one particle are reached, theinteraction is obviously computed by direct sum.The potential at position r due to a cell containing n c particles, of total charge Q tot = (cid:80) n c i =1 q i , is given in multipole expansion truncated at the quadrupole term(see e.g. Ref. [83]), byΦ( r ) = Q tot | r − r cm | + 12 | r − r cm | ( r − r cm ) · Q · ( r − r cm ) . (5.25)The components of traceless quadrupole tensor Q are Q i,j = n c (cid:88) k =1 q k (cid:2) r k,i − r cm ,i )( r k,j − r cm ,j ) − | r k − r cm | δ i,j (cid:3) , (5.26)where r k = ( r k,x , r k,y , r k,z ) are the particles positions inside the cell of centreof mass r cm and δ ij is the Kr¨onecker delta with i, j running over x, y, z . Thecorrespondent electric field E ( r ) = −∇ Φ( r ) is E ( r ) = Q tot | r − r cm | ˆ e − | r − r cm | Q · ˆ e + 52 (ˆ e · Q · ˆ e ) ˆ e | r − r cm | (5.27)where ˆ e = ( r − r cm ) / | r − r cm | .Tree codes, in their most straightforward implementation, have a numericalcomplexity in the force calculation that scales with the number of particles N as O ( N log N ). N − body codes used for the simulations We now discuss the details of the N − body codes used for the simulations presentedin this thesis. Since we had to model both systems with relatively small number ofparticles (i.e. small atomic or molecular clusters with ∼ − ions),two force and potential calculation schemes have been employed among thosediscussed above, the direct sum and the particle-mesh. The first being suited forsmall systems where inter-particle effects can not be neglected, while the latter isinstead better to describe large systems where the force on most of the particlesis dominated by the mean field. Chapter 5. Numerical methods Φ (r * ) g p g spl1 g spl2 g spl3 pure Coulomb00.511.522.53 0 1 2 3 4 5 - d Φ (r * ) / d r * r * Figure 5.4: Potential (top) and acceleration (bottom) for different choices of thesoftening kernel (color lines) compared to the real Coulombic interaction (blackline) as function of r ∗ = r/(cid:15) . Quadratic and cubic splines g spl2 and g spl3 , that havemore complicate polynomial expressions are also shown. The softening of potential and forces in direct codes
As anticipated before, due to the singular nature of the 1 /r interaction fro r = 0,the numerical computation of the potential and force due to a point-like particlepresent problems for small values of r .In the direct codes, where a large number of pair interaction computationis involved, as well as in tree codes for the near neighbours contributions, thepotential at distance r from a charge q is substituted byΦ soft ( r ) = q(cid:15) g ∗ ( r/(cid:15) ) , (5.28)where (cid:15) is the softening length and g ∗ , the so called softening kernel, is a functionthat determines how the interaction is regularized, that has to approach 1 /r forlarge separations. In our direct code we implement two of the most widely adoptedchoices g p = 1 (cid:112) r /(cid:15) + 1 , (5.29).2. The N − body codes used for the simulations that gives a softened form Φ corresponding to that generated by a spherical dis-tribution of total charge q with density ρ = 3 q/ (4 π(cid:15) )(1 + r /(cid:15) ) − / , and g spl1 = (cid:40) − ( r /(cid:15) − / r/(cid:15) < (cid:15)r ; r/(cid:15) ≥ , (5.30)that corresponds to the linear spline. In Fig. 5.4 we show the softened potentialand acceleration regularized using different functions against the real interaction.Higher order spline functions as well as different fictitious densities with finiteor infinite support can be used, (see Ref. [41]). Alternatively, if the number ofparticles in the simulations is not prohibitively large, the problems introducedby the misrepresentation of 1 /r for small r can be partially overcome using anadaptive or particle dependent time step.The value of (cid:15) used in the simulations is chosen each time by bench-markingbetween physical consistency and numerical stability criterions. In the numericalcalculations discussed in this thesis, where each particle of the run represent anion or an electron and the softening kernel g spl1 of Eq. (5.30) is used, such asfor example the simulations of molecular clusters irradiated by X-ray pulses, weadopt (cid:15) = a , where a is the Bohr radius. With this choice the interaction for r > a is the true Coulomb one, while for r ≤ a is substituted by a harmonicoscillator.In the direct simulations where the numerical particles represent only a MonteCarlo sampling of a given macroscopic charge density (for instance our simulationsof expanding spheroidal ion bunches, see Chap. 3), where we adopt instead thesoftening kernel g p that implicitly smooths out the total charge density, which isconstituted by a sum of deltas, and reduces unwanted discreteness effects, we takeas optimal (cid:15) half of the minimum inter-particle distance of the initial condition,which allows also to chose sufficiently large time steps, see also Sect. 5.2.It must be pointed out that the substitution of the Coulomb potential withone of its softened versions based on an infinite kernel may introduce unphysicaleffects even at the level of the mean field felt by a test particle of the simulation. The normalization of the equations of motion and the choice of theoptimal time step
The normalization of particles’ equations of motion is an important point of the N − body simulation, since it affects also the choice of the optimal time step ∆ t .In the majority of the simulations performed for this thesis, where atomic physicsprocesses are considered, the natural choice is to express all the physical quantitiesin atomic units (a.u.) defined so that m e = e = ¯ h = κ = 1, with this choice, theBohr radius for the hydrogen atom a is also 1 and becomes the unit in which thelengths are expressed, while from the fine structure constant α = κ e / ¯ hc we havethe speed of light as c = 1 /α (cid:39)
137 which means that velocities are then expressedin units of αc and energies in terms of the so called Hartree energy E h = α m e c . Chapter 5. Numerical methodsOn the other hand, in simulations of systems dominated by the mean field orin those extended to large times a different approach is used. First a natural timescale t ∗ in which the times will be expressed is set up from physical considerationsor with dimensional arguments, for instance in the case of the Coulomb explosionof a homogeneous sphere of charge Q , mass M and initial radius R a typicalchoice is t ∗ = (cid:34) √ √ (cid:35) (cid:115) M R Q , (5.31)that is the time it takes to the sphere to increase its radius to 2 R , second, byintroducing a scaled length r ∗ to express the lengths one finally has the velocityscale as v ∗ = r ∗ /t ∗ . In general the radius containing half of the particles of thesimulation is a good choice of r ∗ .Choosing the time step ∆ t for the integration of the equations of motion de-pends on the kind of normalization used by the code and also by other factors suchas the intention to neglect or not phenomena happening on time scales consider-ably smaller than t ∗ . In the type of problems treated here, normalizations in a.u.are always associated with ∆ t ∼ ω − B /
10, where ω B is the electron frequency of thehydrogen atom in the fundamental state of the Bohr hydrogen atom. Whenevera macroscopic timescale is chosen, we typically use ∆ t ∼ t ∗ / . The propagation scheme
Several methods do exist in order to integrate the system of ODEs (5.2), In bothour direct and particle mesh codes we use the so called leapfrog scheme, also knownas Verlet algorithm, (see e.g. Ref [66]). This particular second order method isbased on the splitting of Hamiltonian (5.1) in its two terms, kinetic and potential,so that a particle of mass m , momentum p = m v and position r is moved in phasespace at each time step ∆ t , in two half steps alternatively due to K and to U : r n +1 / = ∆ t p n m + r n ; (5.32) p n +1 = ∆ tm a n +1 / + p n ; (5.33) r n +1 = ∆ t p n +1 m + r n +1 / , (5.34)where a n +1 / is the acceleration computed at half time step.Albeit being very simple to implement, the leapfrog scheme has the importantproperties of being time reversible (contrary to some higher order schemes suchas Runge-Kutta or Hermite, see [137]) and, at least theoretically within machine In principle for some particular problems involving electrons oscillations it might be worthtaking instead the inverse Langmuir frequency of the electrons t ∗ = ω − = (cid:112) m e (cid:15) / πn e e asscale time and the Debye length λ D as scale length. .2. The N − body codes used for the simulations precision, to conserve in time average the invariants of the dynamical system iden-tified by Eqs. (5.2) such as the Hamiltonian itself and the angular momentum L .Algorithms of this kind are called symplectic.It must be pointed out that if one uses relativistic mechanics or if more gen-erally a depends on v , the symplecticity properties are not fulfilled anymore andlarge errors on the particles’ orbits can be introduced even with small and adaptive∆ t . In the case of relativistic simulations, the most popular approach is the Borisscheme [19], that needs only another intermediate step to update the relativisticfactor γ n +1 / but unfortunately introduces spurious drifts in particles’ orbits forparticular configurations of electric and magnetic fields for which E + v × B = 0.General velocity dependent accelerations can be treated with the recently intro-duced method by Rein and Tremaine [138], that is symplectic and involves just acoordinate change to a non inertial frame.In our implementation, whenever the problem involves relativistic velocitiesand or non zero magnetic field, we adopt instead the corrected Boris scheme in-troduced for the first time in the context of merging of relativistic beams, see [162].At the cost of a large number of substeps, this method manages to circumventthe problems plaguing the standard relativistic methods and is time reversible.A particle of mass m , charge q , velocity v n , position r n and relativistic factor γ n = (cid:112) − v n /c at step n and undergoing the action of (in general) time depen-dent fields E and B , is advanced at step n + 1 in the following way: r n +1 / = ∆ t v n + r n ; (5.35) u n +1 / = γ n v n + qm ∆ t (cid:0) E n +1 / + v n × B n +1 / (cid:1) ; (5.36) u (cid:48) = u n +1 / + qm ∆ t E n +1 / , γ (cid:48) = (cid:112) u (cid:48) /c ; (5.37) τ = qm ∆ t B n +1 / , u ∗ = u (cid:48) · τc , σ = γ (cid:48) − τ ; (5.38) γ n +1 = (cid:115) σ + (cid:112) σ + 4( τ + u ∗ )2 , t = τγ n +1 , s = 11 + t ; (5.39) u n +1 = s (cid:2) u (cid:48) + ( u (cid:48) · t ) t + u (cid:48) × t (cid:3) , v n +1 = u n +1 γ n +1 ; (5.40) r n +1 = ∆ t v n +1 + r n +1 / , (5.41)where as usual quantities indicated in boldface are vectors and in normal fontscalars. Building the grid in the PM code
In our implementation of the particle mesh algorithm, essentially based on the superbox code, [54], the grid on which the charge density is displaced with anearest grid point method (Equations 5.9 and 5.10) is a cubic cartesian one. Since Chapter 5. Numerical methodswe are principally describing systems undergoing expansion we allow its spacing∆ s to increase with time. This is done as follows, first at t = 0 the centre of massof the system (c.o.m.) is identified and the coordinates of the particles changedso that it coincides with (0 , , r max istaken as the side of the part of the grid containing the real system and the firstmesh constructed along the three axes with the same ∆ S = 4 r max /N g so that for l = x , y or z r li = − r max + ( i − s, (5.42)the space (real plus ghost zones) is partitioned as sketched in Fig. 5.2 in N g nodes( r xi , r yj , r zk ).The Green function is then computed as in Eq. (5.18). The procedure isrepeated at every time step without identifying the c.o.m.. Since a normalizeddiscrete Fourier transform is used, and the grid spacing is uniform in all directions,one does not need to recompute the Green function, but just to multiply for theopportune function of ∆ s the normalized value of the grid based potential Φ i,j,k obtained by the Poisson solver routine before interpolating it at the particlespositions.Using a time dependent mesh has implications on the choice of the simulationstime step, we require that at every time∆ t < ∆ s ∆ σ v (5.43)where ∆ σ v is the variation of the particles’ velocity dispersion between two con-secutive time steps. In general this is obtained by taking ∆ t = 0 . s/ ∆ σ v . We have discussed until this point, only the numerical techniques used to treat thedynamics of systems of charged particles. However, in this thesis, we are mainlyinterested on the response of such systems when exposed to strong laser radiation,in particular short and intense x-ray impulses with typical durations of a few fem-toseconds, intensities in the range 10 − W / cm and photon energies of theorder of a few keV.When a cluster is irradiated with such strong laser fields (for monographicreviews see Refs. [136], [143] and [145], see also [164] and [159]), its constitut-ing atoms or molecules experience many photoionization events (i.e. electronsare stripped and released with kinetic energies K photo = ¯ hω − E bind where ω is the incoming photon frequency and E bind is the electron binding energy). Inaddition, when the inner electronic shells are ionized first, and the ion has stillother electrons, such unstable ionic configuration (say state a ) decays to a sta-ble configuration (state b ) via (in general chains of) the so called Auger process[140] that refills the inner shell and release a second electron with kinetic energy.3. Including quantum processes in classical dynamical simulations: a MonteCarlo approach K auger = E a − E b to the continuum. To carry out a full quantum description is Photoelectron Auger electron X-ray photon
Figure 5.5: Sketch of K − shell photoionization followed by an Auger decay in acarbon atom.impossible even on modern machines since it involves the solution a many electronSchr¨odinger equations in three spatial dimensions. On the other hand, using semi-classical methods or techniques based on mean field approaches such as DensityFunctional Theory (DFT, see [22]), implies a high degree of simplification thatprevent to study in detail the individual dynamics of the electrons, unless oneintroduces additional artifacts [9].In order to find an acceptable compromise, we follow the same line of Refs.[59], [60], [142] and [62], where the above mentioned quantum processes are treated“statistically”; roughly speaking, for the i − th atom or ion of the simulation firsta transition probability P is computed (for instance from the laser intensity, thephotoionization cross-section and the time step in the case of ionization), then arandom number p is extracted from a normal distribution and finally the ion’selectronic configuration is changed accordingly, whenever p ≤ P . Subsequentlythe electrons released by photoionization or Auger decay are initialized in thesimulation with their kinetic energies K photo or K auger and then treated as classi-cal particles.By doing this, each many-particle system treated in a single simulation is justa Monte Carlo realization of the problem considered, and therefore integrated aswell as time dependent quantities (such as the final charge states or the time evo-lution of the total kinetic energy), have to be averaged over many realizationsstarting with the same initial conditions, in order to have statistical meaning or Chapter 5. Numerical methodsto be compared for different parameter sets.In the next subsections we show how the photoionization and Auger decayprobabilities P photo and P auger are calculated. Photoionization
The probability used in the Monte Carlo sampling, that a photoionization eventhappens between times t and t + ∆ t for the j -th shell of the i -the atom or moleculein the cluster is given by P photo j,i = I ( t, r i ) σ j,i ∆ t ¯ hω ; j = 1 s, s, p, ... (5.44)where I ( t, r i ) is the laser intensity at time t at the ion’s position r i , and σ j,i thephotoionization cross section for which we refer to the values given in Ref. [12].If the j − th shell has one or more vacancies, P photo j,i is reduced by the factor χ = N e − N v N e , (5.45)where N e is the number of electrons in the full case and N v the number of vacancies.For the theoretical laser pulses we use Gaussian time envelopes, and constantspace profiles, so that I ( t, r ) = I e − ( t − t σ , (5.46)where I is the peak intensity and σ is the Gaussian’s standard deviation related tothe full width half maximum by FWHM = 2 √ σ , and t is the time at whichthe laser has its peak (i.e. I ( t ,
0) = I ). Alternatively, if the dimension of theirradiated sample is comparable to the radius of the the laser focus z a Lorentzianenvelope is used along the direction z in the focus plane and a Gaussian envelopewith standard deviation ω , along the pulse propagation direction r . The spaceenvelope therefore reads I ( r, z ) = I t
11 + ( z/z ) e − r ω [ z/z ] . (5.47)Note that, at photon energies of few keV, for the elements of the first row thatwe are considering in the present study, the photoionization happens mainly via K -shell photoionization (1 s orbitals), for ¯ hω = 1 keV in fact, the cross section σ s , given (see e.g. Ref. [50]) by σ s = 256 π α (cid:16) a Z (cid:17) (cid:18) E bind , s ¯ hω (cid:19) / , (5.48)where Z is the nuclear charge of the element, α = 1 /
137 and a = 5 . × − mare the fine structure constant and the Bohr radius respectively, is of the order of.3. Including quantum processes in classical dynamical simulations: a MonteCarlo approach × barn, while σ s and σ p are roughly 10 and 10 barn respectively. There-fore, In the numerical calculations we assume σ s = σ p = 0. Test simulations withreal values for all the cross sections, revealed a percentage of photo-electron com-ing from the outer shells smaller than the 3% at peak intensities ≥ W / cm ,making the latter assumption rather reasonable.Note also that in our simulations, when an electron is released by photoioniza-tion, the direction of its velocity is extracted randomly from a normal distributionwhile its mother ion recoils consistently with m ion v , ion hω − E bind = m e v , e m ion v , ion m ion v , ion = m e v , e + m ion v , ion , (5.49)where m ion and m e are the masses of the ion and the electron and v , ion and v , ion , v , e their velocities before and after the photon absorption. Auger decay
The K − shell photoionization produces molecular or atomic ion states prone todecay via Auger processes (see e.g. Ref. [4]) in which an external electron refillthe inner vacancy while a second is emitted with energy corresponding to theenergy difference between the initial (unstable) and the final configurations.For a given ion i , the transition probability from the state a to the state b viaan Auger decay between times t and t +∆ t is computed from the energy amplitudebetween the two configurations Γ abi and reads P auger ab,i = Γ abi ¯ h n e ( n e − n ( n −
1) ∆ t, (5.50)where n e and n are the total number of electrons in principle able to take part tothe transition, and the number of electrons present in the external shells in theequivalent neutral respectively. For the species of interest in this thesis, we usedthe values of Γ given in [93].It must be also mentioned that an analogous approach for the calculation ofAuger transition rates has been recently independently developed in [79] and [80],and used to investigate the dynamics and the fragmentation of water and methanemolecules under intense X-ray lasers. Recombination
The spatial finiteness of the plasmas studied here and the short time scales consid-ered ( < . Chapter 5. Numerical methodsFollowing the approach of Ref. [59], between every two time steps the revo-lution angle θ i,j ( t ) of the i th electron around the nearest ion j is determined. Ifsuch electron revolves two times around the same ion j (i.e. θ i,j ( t + ∆ t ) = 4 π ) itis considered localized to that ion and the electronic configuration of the latter isupdated with an electron in an excited state (i.e. a “classical recombination” hap-pened). A simpler procedure to assign an electron to the neighbouring ion wouldbe instead checking whether the total energy of the i th electron in the referencesystem of the nearest ion of charge q j E ij = eq j | r i − r j | + m e (cid:2) ( v x,i − v x,j ) + ( v y,i − v y,j ) + ( v z,i − v z,j ) (cid:3) + eU ( r j ) (5.51)is negative, together with the fact that the two particles are at a distance | r i − r j smaller than a certain fiducial scale length. In the definition above, U ( r j ) is thecluster electrostatic potential at the ion’s position. In both cases sketched here,the process is obviously strongly influenced by the way in which the Coulombinteraction is regularized for small separation. To conclude the chapter, we give a brief description of the tools used to analyzethe end products of our numerical simulations such as the routines to extract thefinal shapes of the charge distributions and the averaging procedures in the MonteCarlo simulations.
The shape of the final charge distribution
To extract the aspect ratio of a given non-spherical distribution of charges weuse a procedure based on tensor diagonalization. Following the algorithm used in[116] for the end products of gravitational N − body simulations, we compute thesecond order tensor I ij ≡ N (cid:88) k =1 r ( k ) i r ( k ) j (5.52)for the particles inside the sphere of radius r containing the 98% of the totalnumber of particles in the system, where r i are the Cartesian components of theposition vector in the reference frame with origin in the centre of mass. I ij isrelated to the inertia tensor by I ij = Tr( I ij ) δ ij − I ij . (5.53)The matrix I ij is diagonalized iteratively, requiring that the percentage differenceof the largest eigenvalue between two iterations to be smaller than 10 − . Thisprocedure requires on average 10 iterations, and we call I ≥ I ≥ I the threeeigenvalues. We finally apply a rotation to the system in order to have the three.4. The analysis of the end products eigenvectors oriented along the coordinate axes. For a heterogeneous ellipsoid ofsemi-axes a, b and c , we would obtain I = Aa , I = Ab and I = Ac , where A is a constant depending on the density profile. Consistently, for the end-productswe define b/a = (cid:112) I /I and c/a = (cid:112) I /I , so that the ellipticities in the principalplanes are (cid:15) = 1 − (cid:112) I /I and (cid:15) = 1 − (cid:112) I /I . The realization averages of distributions and global quantities
As we have previously mentioned, when using the Monte Carlo approach to treatquantum phenomena, each simulation represents only a single statistical realiza-tion of a given chain of events. On the other hand, in laboratory experiments oflaser-cluster interaction, one measures the yield of integrated quantities over theproducts of several clusters irradiated by the same laser pulse. To have a mean-ingful comparison between measured quantities and their simulated counterparts,or to have a reasonable picture of their scaling with the laser intensity I , thelatter have to be processed with an averaging procedure. When one performs aset of N r realizations of a given numerical experiment initialized with identicalinitial conditions (i.e. cluster size, laser intensity etc.), but different seeds for themachine’s random number generator, a time dependent quantity S ( t ) such as forinstance the number of absorbed photons n γ is averaged as (cid:104) S ( t ) (cid:105) = 1 N r N r (cid:88) i =1 S i ( t ) (5.54)while a quantity at fixed time t that depends itself on an average over the totalnumber of ions N ion such as their average charge state (cid:104) q (cid:105) ion is given by (cid:104) S (cid:105) = 1 N r N r (cid:88) i =1 N ion N ion (cid:88) j =1 S j,i . (5.55)We have observed that typically, for clusters containing ∼ (cid:104) S (cid:105) vs I curves.It must be also pointed out that in real experiments involving jets of clustersirradiated by a short and intense laser fields, it is unlikely that all the clusters ofthe jet are hit by the pulse at focus, if one needs for some reason to take thatinto account in the numerical calculations, a pulse average is required. Usually,in order to do so, following Ref. [124], one has first to chose a time and spatialprofile for the theoretical laser pulse and sample them to perform N r = N t × N s independent simulations with constant intensity I i,j in time and space from whichthe sampled quantity S ij is extracted. In this case, a time and space averagedquantity (cid:104) S (cid:105) s,t is finally obtained as (cid:104) S (cid:105) s,t = N s (cid:88) i =1 ln (cid:18) I i +1 I i (cid:19) i (cid:88) j =1 S i,j (cid:34)(cid:115) ln (cid:18) I i I j (cid:19) − (cid:115) ln (cid:18) I i I j +1 (cid:19)(cid:35) . (5.56) Chapter 5. Numerical methodsAs an alternative, one can instead assume an effective space profile for the laserpulse, for instance that given in Eq. (5.47), perform different runs where thetargets are displaced with respect to the laser focus according to a normal distri-bution, and finally average the wanted quantities with Eqs. (5.54) or (5.55). Inthe results presented in this work we have employed only the latter method. hapter 6
Summary and outlook
In this thesis we have presented an exploratory study of clusters exposed to shortand intense x-ray pulses. We have considered the implications of the massivecharging (i.e. one charge per atom/molecule) happening on a time scale of theorder of few femtoseconds on the dynamics of ions and electrons inside the cluster,and the process of expansion (i.e. Coulomb explosion). Such extreme conditionsare nowadays experimentally accessible due to high intensity femtosecond x-raypulses attainable in modern free-electron laser facilities. For the microscopic de-scription of such conditions a numerical code has been developed, integratingatomic processes such as photoionization, Auger decay and recombination with N − body molecular dynamics.Before studying specific systems, which are of particular interest in the viewof recent experiments, we have performed simulations of the explosion of chargedspherical clusters under different conditions of composition and density profile,and studied their end products by means of idealized analytic models based oncontinuum approximation.It appeared clear that, even though the expansion processes considered hereare mainly collisionless, the discrete nature of the systems still induces effects onthe energy distributions which are impossible to be accounted for in the contin-uum model. In particular, for systems starting with a homogeneous charge densityprofile, the energy distribution shows a peaked structure close to the cutoff en-ergy and therefore strongly departs from that which is predicted when assuminga continuum system.We have verified that this feature, discussed in Ref. [146], persist for differ-ent particle arrangements in the initial condition, with a general trend of gettingmilder for more randomized particle “displacement” and sharper with increasingnumber of particles.Remarkably, the differences between numerical and theoretical energy distri-butions of exploding clusters seem to be enhanced by the regularization of theCoulomb interaction. The latter, although leaving the long-range behavior basi- Chapter 6. Summary and outlookcally unaltered, induces for particles close to the surface of the system an unphys-ical “external field effect” (i.e. such particles suffer an extra force pushing themoutside). Such a fact obviously affects the analogous gravitational N − body codesand has not been documented in the literature.In view of molecular clusters or large molecules we have extended our studyto two-component (i.e. light and heavy atoms) and non-spherical systems. Wefind that the energy distribution for the light component of two species systemssignificantly diverges from that predicted by analytical (continuum) models, as itsdynamics is strongly influenced by the system’s discreteness at early stages of theexplosion.In the Coulomb explosion of homogeneous systems with ellipsoidal shapes, itturns out that the initial aspect ratio is not conserved during the expansion butinstead the initially prolate systems become oblate as they expand and vice versa.Expressions for the energy spectra are given for the case of axisymmetric systemsand limit aspect ratios are derived. Molecular dynamics simulations show remark-ably good agreement with our analytic predictions.Prompted by experiments with methane clusters at the x-ray free electron laserLCLS, we have studied molecular clusters. To shed some light on the surprisingexperimental findings we made a systematic study over a broad range of physicallyrelevant intensities (from 5 × to 10 W / cm ) and pulse lengths between 3 and100 fs. More importantly, a whole series of hydride clusters with iso-electronic con-figuration but different photoabsorption and Auger transition and different atomicmasses was analyzed.Spanning such ranges of intensities and pulse lengths, we observed a remark-ably non-monotonic trend of the ratio between the average kinetic energies ofprotons and heavy ions showing a maximum at a well determined laser intensityfor fixed photon energy and pulse length. This effect is in apparent contrast withthat which is found in the case of the pure Coulomb explosion of multi-species clus-ters where such a ratio is markedly monotonic, for every combination of chargeratios and percentages of the two species.This has been interpreted as the presence of three distinct electron-inducedregimes of energy redistribution between the two species of ions. In the first case(at low intensities), the cluster is weakly charged having suffered few photoion-izations and expands as a whole, as a quasi-neutral plasma. In the intermediateregime enough charge is created via photoionization. The electrons stripped bythe (time dependent) cluster electrostatic field do not leave the system but adjustto screen its charge in the core, leaving a charged shell at the surface from whichprotons escape with large kinetic energies due to their lower mass.Finally, at large intensities the whole cluster is strongly charged and expandsvia pure Coulomb explosion as in idealized models neglecting the contribution ofelectrons. It must be pointed out that, in contrast to what is observed in clusters irradiated by long wavelength lasers, for the conditions studied here, the charg-ing processes happen on a time scale that is of roughly two orders of magnitudesmaller than that of any appreciable (heavy) ion dynamics. The trapped electronsthat do not have enough time to thermalize are, however, inducing a previouslyunobserved species segregation channel.We stress the fact that the different scenarios described above, are obtainablesimply by tuning the laser intensity and, having noted that for the intermediateintensities the backbone of the cluster structure stays basically unaltered, may beof some relevance for x-ray based single molecule imaging.In conclusion, we point out that studying the dynamics of laser irradiatedclusters might be relevant for other areas of physics and their technological appli-cation. For instance, the possibility to control (by tuning the laser intensity) thefraction of energy absorbed by the cluster which is transferred to the light ions’kinetic energy, and therefore to a certain extent, control their energy distribution,opens to applications in particle acceleration or intra-cluster nuclear fusion, wherea narrow ion velocity spectrum is needed [98]. In this case, deuterated molecu-lar hydride clusters (i.e. hydrogen in CH or H O is substituted by its isotopedeuterium) are to be used. The nanoplasmas produced by the interaction of in-tense lasers with clusters are characterized by kinetic and potential energies K and U such that their Coulomb coupling parameter Γ, expressed as function ofits minimum inter-particle distance a , temperature T and the electron charge e by Γ = U/K = e /ak B T , is of the order of unity [141]. This places the lasergenerated cluster plasmas at the edge between ideal plasmas (Γ (cid:29)
1) [11], andstrongly correlated plasmas (Γ (cid:28)
1) [56]. The latter are thought to exist in naturein astrophysical environments such as the cores of giant planets or white dwarfsstars. Even if their intrinsic parameters such as densities and temperatures maydiffer by orders of magnitude, plasmas of equal Γ are expected to behave simi-larly and show analogous dynamical properties ( see e.g. [141]). Thus the studyof laser generated nanoplasmas is very likely to shed some light on the physicsof non-experimentally accessible plasmas. Nevertheless, due to their small sizesand relatively short life-times (being overall non-neutral, they expand under theirself consistent Coulomb potential doubling their initila radius in roughly 10 fs),the direct investigation of cluster plasmas dynamics is still prone to a number ofcomplications, and one relies in general on the indirect information obtained bythe fragmentation products and on numerical simulations.From the point of view of non-equilibrium thermodynamics, finite and globallynon-neutral Coulomb systems, where the electrons are confined by the space chargeof the ions, share some peculiar properties with gravitational systems [35], [163](or more generally with systems interacting via long-range forces, see [32], [20] and[44]), such as for example a negative specific heat in the microcanonical ensemble[109], [158], that is the more the system is heated, the more it cools. The laser Chapter 6. Summary and outlookinduced electron-ion plasmas treated in Chap. 4, likely fall in this regime andtherefore their study could be particularly relevant to test experimentally recentresults in non-equilibrium statistical mechanics and thermodynamics. ppendix A
Cluster production
Clusters can be classified in different ways according to their physical and struc-tural properties, in the first instance we distinguish between atomic an molecularcluster whereas they are composed by atoms or molecules. The first are held to-gether by either metallic, covalent or ionic bonds, or by London forces as in thecase of rare gas clusters, while the second prevalently by van der Waals forcesbetween the induced electric dipoles on molecules [156] or ionic bonds. We re-call that throughout this thesis we used the additional classification, speaking of homogeneous clusters when they are composed by a single species of atoms, and heterogeneous clusters when they are instead composed of more then one atomicspecies, or by hetero-nuclear molecules as in the case of molecular hydride clus-ters.The number N ∗ of constituting particles and the type of bonds supportingtheir structure are responsible for the shape of the cluster. Rare gas clusters forinstance (but also the molecular clusters of our interest) are characterized by theso called isocaederical structure (see Fig. A.1) where particles are arranged onconcentric shells so that with increasing N ∗ , the fraction of atoms sitting at thesurface n s decreases proportionally to N − / ∗ , which implies that for small sizesthe majority of the cluster’s constituents is located at its surface.Clusters are artificially produced by letting expand into void a sonic jet of gaswith pressure P and temperature T , coming out of a conical nozzle with aperture α and orifice diameter d of the order of 1 µ m (see [69], [84] and references therein).The average number of atoms or molecules in clusters formed in this way, is givenby (cid:104) N ∗ (cid:105) = A ( H ∗ / ) γ (A.1)where A and γ are constants and H ∗ is the so called (reduced) Hagena parameter[68] given by H ∗ = P d q eq T . q − . K ∗ ; d eq = 0 . d/ tan( α/ , (A.2)in which the constant K ∗ depends on the species and it is related to the minimuminter-particle distance in the solid phase r min and the sublimation enthalpy h s and Appendix A. Cluster productionreads K ∗ = 1 / ( r q − T . q − . ∗ ); T ∗ = h s /k B . (A.3)in all formulas above, the parameter q is found experimentally to fall in the interval0.5-1, see [37] and [69]. The distribution of sizes d P/ d N ∗ in a get of clusters isgiven by a lognormal distribution (see e.g. [133] and [112]) and readsd P d N ∗ = 1( N ∗ − N ) σ √ π exp (cid:20) − (ln( N ∗ − N ) − µ ) σ (cid:21) , (A.4)where N is the position of the “zero” and µ and σ the logarithms of the geometricmean and standard deviation. Detection of ultrafast plasma dynamics
We will discuss the scenario of a possible detection of the ultrafast plasma dynamicsin the context of laser pulses from unconventional light sources. One such sourceis a free electron laser (FEL) at soft X-ray and soon hard X-ray frequencies. It isexpected to deliver light pulses with frequencies of v
12 eV–12,000 eV over lessthan 50 fs and with unprecedented intensities of more than 10 W/cm . Light pulseswith such fantastic properties have their price: they are generated by first acceleratingelectrons to almost the speed of light over 2 km in a tunnel, then sending theelectrons through a series of magnets to force them into bent trajectories. Thereby,they radiate and organize themselves in bunches, leading overall to a dramaticcoherent amplification of the emitted light over a certain, narrow frequency range.This principle is called SASE (Self-Amplified Stimulated Emission) and machinesworking with it are built at DESY in Hamburg and at SLAC in Stanford.The second new development is attosecond pulses (1 as s) whoseduration comes close to the period of about 100 as an electron needs, in theground state of hydrogen, to travel around the proton and whose photon energy isbetween 20–150 eV. These pulses are generated from certain fractions of intenselonger laser pulses at longer wavelengths and, consequently, their intensity ismuch lower than the light from FELs. More details can be found in Ref. 11.Not only is the detection itself difficult, but the monitoring of the plasma mustbe synchronized with its creation; in other words, one needs to know when theclock started. These set-ups are called pump-probe experiments: the first ‘pump’pulse creates the dynamics one is interested in while the second ‘probe’ pulse,probes the dynamics at a well-defined time delay D t . Figure 5.
Closed shell configurations of rare gas clusters form so called isocaedersand are particularly stable. They contain the number of atoms as indicated
Complex Non-equilibrium Dynamics in Plasmas
Figure A.1: Rare gas clusters closed shell structure according to their (indicated)number of atoms N ∗ , so called “magic numbers”. The figure is reproduced fromRef. [141]. ppendix B Dynamical friction
Following the classical treatment of Chandrasekhar [29] and Spitzer [150], we de-rive here the formula of the dynamical friction, defined as the drag effect on a testparticle of charge q t and mass m t moving through a background of particles withcharge q f and mass m f due to long range collisions (i.e. Rutherford scattering[84]).Let us consider first as in the original derivation (see e.g. [15]) a binary col-lision in the so called impulsive approximation , for which the deflection of thetrajectories due to a long-range force is approximated as due to a finite-time in-teraction, (see [139]). Single orbital deflections, like that sketched in Fig. B.1,are computed independently, and their contributions are summed vectorially overall the encounters in the field distribution assuming no correlation between con-secutive collision (Markovian hypothesis). Over the relative orbit we define therelative velocity as v = v t − v f and the reduced mass as usual as µ = m f m t m f + m t . (B.1)According to Newton’s Third Law of Dynamics, the velocity change ∆ v t alongthe test particle’s (unbound) orbit is exactly given by∆ v t = µm t ∆ v , (B.2)where ∆ v is the vectorial change of the relative velocity, obtained from the solutionof the two body problem. We define in impulsive approximation the modulus ofeffective acceleration felt by the test particle as a = q t q f µb (B.3)where b is the so called impact parameter is defined as the smallest distance be-tween the unperturbed trajectory of the test particle and the target, or morerigorously as b = Lµv t (+ ∞ ) (B.4) Appendix B. Dynamical frictionwhere L is the modulus of the total angular momentum on the hyperbolic orbitFigure B.1: Unperturbed orbit (dashed dotted line) of a particle of charge Q andrelative orbit (solid line) when deflected by a field particle of charge q for the case Qq < v t (+ ∞ ) is the modulus of the asymptotic velocity of the test particle. If wenow take as the effective time on which the interaction takes place during thebinary collision t eff = 2 b | v | (B.5)we can express the velocity changing in the direction perpendicular to the initialrelative trajectory as ∆ v ⊥ ∼ q t q f b | v | . (B.6)Note that, in the latter the above expression becomes an exact equality in thelimit of large b or relative velocity. From impulsive approximation we can statethat || v || = || v + ∆ v ⊥ + v (cid:107) || ; now combining the latter with equation (B.2) forthe component v t (cid:107) of the test particle’s velocity parallel to the initial v t , we easilyobtain at first order ∆ v t (cid:107) = µ ∆ v (cid:107) m t ∼ − µ || ∆ v ⊥ || m t || v || v . (B.7)By substituting equation (B.6) on the r.h.s. of relation (B.7), after easy algebrawe finally get ∆ v t (cid:107) (cid:39) − q t q f ( m t + m f ) b m t m f || v || v . (B.8)We define now the number of binary collisions with impact parameter rangingfrom b to b + ∆ b in the time interval ∆ t as∆ N = 2 πbdb || v t − v f || ∆ tnf ( v f ) d v f (B.9) where n is the number density of field particles and f ( v f ) is their velocity distri-bution function that here we will assume to be a Maxwellian with dispersion σ ,that reads f ( v f ) = e − v f / (2 σ ) (2 π ) / σ . (B.10)It is now easy to write down as finite difference the velocity change on a timeinterval ∆ t for the test particle in the parallel direction after ∆ n encounters inthe form ∆ v t (cid:107) ∆ t = − πn ( m t + m f ) m t m f q t q f f ( v f ) v b || v t − v f || dbd v f . (B.11)The integration over b can be carried out separately since it does not involves thevelocity, obtaining (cid:90) b max b min dbb = ln( b max ) − ln( b min ) = ln b max b min = ln Λ (B.12)that is the well known Coulomb logarithm. It must be pointed out that, sincewe are working in impulsive approximation, a divergence appears for b min → b min → + ∞ (infrared divergence); we can get rid of them setting b min to its minimuminter-particle separation and b max equal to the Debye length λ D of the system whenthe test particle moves through particles with charge of the opposite sign, or to thesystem’s size in the opposite case. Note that, in unscreened “infinite systems” theinfrared divergence remains, due to the long range nature of Coulomb interaction.Note also that, due to this, the contribution to the deflection of the far particlesis in general larger than that due to short impact parameter encounters.The last step is now to integrate over the velocity. In order to do so accordingto the classical approach, we assume that the velocity distribution for the fieldparticles is spherically symmetric around the test particle and we take the velocityaveraged Coulomb logarithm ln ¯Λ. The integration in d v f gives finally d v t (cid:107) dt = − πn ( m t + m f ) m t m f q t q f ln ¯Λ Θ( v t ) v t v t . (B.13)where the so called fractional volume function Θ( v t ) is given assuming f ( v f ) fromEq. (B.10) byΘ( v t ) = 4 π (cid:90) v t f ( v f ) v f dv f = Erf( v t / √ σ ) − v t exp( − v t / σ ) σ √ π , (B.14)and the term ω coll = 4 πn ( m t + m f ) m t m f q t q f ln ¯Λ v t (B.15) Appendix B. Dynamical frictionis sometimes referred to as collision frequency.Equation (B.13) gives us the effective deceleration suffered by the test particlecrossing the sea of identical particles of charge q f and mass m f . It can be imme-diately noticed that, for different charge states of the same species, (same massdifferent charges), the particles with higher charge states suffer the largest decel-eration. This can be read for example considering an initially coherent beam ofdifferent ions of the same atomic species shot through a gas of identical particles,in the fact that the highest charged ions are more likely spread radially aroundthe beam’s initial direction of propagation.Recently, this formalism has been extended in the gravitational case to systemscharacterized by a mass spectrum Ψ( m f ) [34], so that their total number densityin the definition of ∆ N is given by n = (cid:90) ∞ Ψ( m f )d m f , (B.16)and formally each component has its fractional volume functionΘ( v t , m f ) = 4 π (cid:90) v t f ( v f , m f ) v f dv f . (B.17)In our case of charged particles, two case are physically relevant, namely thediscrete mass spectrum and discrete charge spectrum given byΨ( m f ) = n δ ( m f − m ) + n δ ( m f − m ); Ψ( q f ) = n δ ( q f − m ) + n δ ( q f − m ) , (B.18)corresponding respectively to systems where all particles have the same charge q f and different masses, or same mass m f and different charges. Here we haveconsidered for simplicity only two species with number densities n and n andmasses (charges) m and m ( q and q ). Obviously this can be generalized tosystems with an arbitrarily number of species, and to the more complicated casefor both mass and charge spectra. Following [34], by integrating over all thedifferent species m f , Eq. (B.13) becomes d v t (cid:107) dt = − πn ( m t + (cid:104) m f (cid:105) ) m t (cid:104) m f (cid:105) q t q f (cid:104) ln ¯Λ (cid:105) ˜Θ( v t ) v t v t . (B.19)where n (cid:104) m f (cid:105) = (cid:90) ∞ m f Ψ( m f )d m f = n m + n m , (B.20)the term ˜Θ( v t ) is the total velocity volume factor depending on the choice of f ( v f , m f ) and (cid:104) ln ¯Λ (cid:105) is the mass averaged Coulomb logarithm depending insteadon the choice of n and n . Analogous considerations and steps can be carried outfor the case of a charge spectrum Ψ( q f ). ppendix C Some remarks on kinetic theory
Another theoretical treatment of the Coulomb explosion of spherical nanoplasmasbased this time on Kinetic Theory [106], accounting for different initial densityprofiles and easily generalizable to multi-component systems, is that developed inRefs [94], [96] and [128], see also [119]. This approach makes use of the concept of one-particle phase space distribution function f ( r , v , t ), defined as the (differential)fraction of system at time t in given point ( r , v ) of the (six-dimensional) oneparticle phase space [11], in a way that once the time is fixed, the number density n ( r ) is obtained as n ( r ) = (cid:90) Ω f ( r , v )d v , (C.1)and the velocity field as v ( r ) = 1 n ( r ) (cid:90) Ω vf ( r , v )d r , (C.2)where the integrals are extended the domain Ω in phase space that is “accessible”for the system. Instead of following in time the trajectories of particles or fluid el-ements, here one studies the time evolution of f through the so-called CollisionlessBoltzmann Equation, also known among Plasma physicists as Vlasov Equation D f D t ≡ ∂f∂t + v · ∂f∂ r + ∇ Φ · ∂f∂ v = 0 , (C.3)where the self consistent electrostatic potential Φ is related to f through n ( r ) bythe Poisson equation now written as∆Φ( r ) = 4 πq (cid:90) Ω f ( r , v )d v . (C.4)When the system is in an equilibrium state or in an asymptotic state (i.e. ∂f /∂t =0), the phase space distribution f can be written as function of energy and angular Precisely speaking, one should refer to Eq. (C.3) coupled to the Maxwell equations as
Vlasovequations , see e.g. [153] se also [73] and references therein for an extended discussion. Appendix C. Some remarks on kinetic theorymomentum through the phase space coordinates r and v . If f is expressible as afunction of the energy E only, it is related to the differential energy distribution n ( E ) introduced in Chapter 2 by the relation f ( E ) = n ( E ) · g − ( E ) , (C.5)where g ( E ) = 4 π (cid:90) Ω | v (cid:112) E − Φ( r ))d r , (C.6)is the phase space volume with energy in the range ( E , ( E ) + d E ).In our case of spherically symmetric systems, where the density is assumed tobe angularly isotropic, if the distribution function at t = 0 can be factorized as f = f rad ( r, v rad , f tan ( v , , (C.7)it evolves under (C.3) so that its tangential part f tan is stationary. Therefore in thiscase it is sufficient to study only the evolution of the radial part of the distributionfunction f rad ( r, v rad , t ). For simplicity hereafter we will drop the suffix rad .Adding an initial condition, in our case f ( r, v, t = 0) = δ ( v ) n ( r ), to thesystem of Eqs. (C.3) and (C.4) one has a Cauchy problem whose explicit solutionin integral form for f is, see [94] for the derivation, f ( r, v, t ) = 1 r (cid:90) ∞ n ( η ) δ [ r − R ( η, t )] δ [ v − U ( η, t )] η d η, (C.8)where we the functions R and U are defined as solutions of ∂R ( η, t ) ∂t = U ( η, t ) ∂U ( η, t ) ∂t = qη ∇ Φ( η, t ) mR ( η, t ) R ( η, t = 0) = η ; U ( η, t = 0) = 0 . (C.9)Following the approach of [96], it is possible to obtain an (explicit) approximateexpression for f at any time as a function of the initial number density profile n . By substituting in (C.9) r ∇ Φ( r, t ) with its value for t = 0 and using thefactorization of the δ distribution δ [ g ( x )] = (cid:80) k δ ( x − x k ) / | ∂g/∂x k | , Equation(C.8) is rewritten as f ( r, v, t ) = (cid:88) k (2 s k − n ( η k ) s k | ∂R ( η, t ) /∂η | η = η k δ ( v − U ( η k , t )) , (C.10)where the s k are implicit solutions of t (cid:115) qη ∇ Φ( η ) mη = s ( s − s − s −
1) (C.11) for specific times t and coordinate η , R ( η, t ) = ηs ( η, t ) / (2 s ( η, t ) −
1) and U ( η, t ) =( s ( η, t ) − (cid:112) qη ∇ Φ( η ) / ( mη ) /s ( η, t ) and the summation in k runs over all thesolutions of s η s − − r = 0 . (C.12)Finally, the density and velocity profiles at time t are the obtained from f ( r, v, t )solving the integrals C.1 and C.2. Knowing v ( r ), from the usual expression n ( K ) = 4 πr n ( r, t )d K/ d r , (C.13)it is possible to compute the time dependent kinetic energy distribution n ( K )which is given in Ref. [94] as n ( K ) = 2 π √ K (cid:90) + ∞−∞ f ∗ ( v )d v × (cid:88) l η l n ( η l ) | ∂U ( η, v, t ) /∂η | η = η l ,U = ±√ K , (C.14)where f ∗ ( v ) is the part in velocity of the distribution function and the summationruns over the l solutions of U ( η, v, t ) = ±√ K. (C.15)In the limit of t → ∞ , (implying also that a ( η, t ) → ∞ ) when all the energy isconverted in kinetic energy so that n ( E ) ∞ = n ( K ) ∞ , Equation (C.14) is rewrittenaccording to [119] as n ( E ) = 2 π (cid:88) l η l n ( η l ) | q ∇ Φ( η l ) /m − η l ω ( η l ) | , (C.16)where this time the sum in l is over the solutions of2 q ∇ Φ( η l ) mη l − E = 0 (C.17)and the radial dependent plasma frequency is defined by ω ( r ) = (cid:114) πq n ( r ) m . (C.18)It must be pointed out that the model introduced here allows one to take intoaccount multi-stream motion and particle overtaking as well as to treat regimeswhere the hydrodynamical model is inapplicable, see [25].The approach discussed above can be easily extended to homogeneous twocomponent systems with densities n and n , and unit charges q and q , with“kinematic parameters” a and b defined by Eq. (2.60), and described by twodistribution functions f and f (see e.g. ref. [3], [96], [135] and [2]). Followingthe approach of [3], let us suppose that q n (cid:28) q n (corresponding to low values of Appendix C. Some remarks on kinetic theory b ). In this way, the contribution of the “more mobile” component 1 to the cluster’sradial electrostatic field E ( r, t ) is negligible, and thus the radial component of theelectric field can be written as E ( r, t ) = (cid:40) r (1 − b ) /
3; for r < (1 − b ) − r − /
3; for r ≥ (1 − b ) − , (C.19)where t = √ (cid:82) u d u (cid:48) / (1 − u (cid:48) ) . According to [95] and [96], the radial part of f has the same implicit form given by Eq. (C.8), where this time in the definitionof R and U as solutions of Eq. (C.9), enters the electric field given by (C.19).By integrating numerically Equation (C.14) in the limit of large times, one finds[135], that the asymptotic n ( E ) for the fast component has a power law tail atlow energies and an universal singular cutoff at maximal energy E max proportionalto ( E max − E ) − / . The cut off energy can be estimated as function of the clustersparameters as E max (cid:39) Q R ( a − / ζ, (C.20)where Q is the total charge of the fast component and ζ is the ratio of the unitmasses of the slow and fast component.As a final remark, it should be pointed out that the numerically obtained n ( E ) of Sect. 2.3 reproduce the power-law decay at low energies but do not showsharp cutoffs at E max , having instead a seemingly universal decay (when energiesare normalized with respect to their maximum value). The reason of that is thefact that for the values of a and b used in our simulations, the dynamics of theslower ions influences that of the fast component and, in addition, effects due tothe system’s particulate nature are not entirely negligible. ibliography [1] S.J. Aarseth, Direct methods for N -Body simulations , in Galactic Dynamicsand N -Body Simulations, Lecture Notes in Physics , 275, Springer Verlag,Berlin (1994).[2] A. A. Andreev, P. V. Nickles, & K. Y. Platonov, Quasi-Coulomb explosion ofmulticomponent laser cluster plasma
Phys. Plasmas , 023110 (2010).[3] I. A. Andriyash, V. Y. Bychenkov, & V. F Kovalev, Coulomb explosion of acluster with a complex ion composition , JETP Lett. , 623 (2008).[4] V. Averbukh, U. Saalmann, & J. M. Rost, Auger decay in the field of a positivecharge
Phys. Rev. A , 063405 (2012).[5] W. Bang, G. Dyer, H. J. Quevedo, A. C. Bernstein, E. Gaul, M. Donovan, &T. Ditmire, Optimization of the neutron yield in fusion plasmas produced byCoulomb explosions of deuterium clusters irradiated by a petawatt laser , Phys.Rev. E , 023106 (2013)[6] M. Barbui, W. Bang, A. Bonasera, K. Hagel, K. Schmidt, J. Natowitz, R.Burch, G. Giuliani, M. Barbarino, H. Zheng, G. Dyer, H. J. Quevedo, E. Gaul,A. C. Bernstein, M. Donovan, S. Kimura, M. Mazzocco, F. Consoli, R. DeAngelis, P. Andreoli, & T. Ditmire, Measurement of the Plasma AstrophysicalS Factor for the He3(d,p)He4 Reaction in Exploding Molecular Clusters , Phys.Rev. Lett. , 082502 (2013).[7] J. Barnes & P. Hut,
A hierarchical O(N log N) force-calculation algorithm ,Nature , 446 (1986).[8] Yu. K. Batygin,
Self-consistent analysis of three-dimensional uniformly chargedellipsoid with zero emittance , Phys. Plasmas , 3103 (2001).[9] D. Bauer & A. Macchi, Dynamical ionization ignition of clusters in intenseshort laser pulses , Phys. Rev. A , 033201 (2003).[10] A. Beck & F. Pantellini, Spherical expansion of a collisionless plasma intovacuum: self similar solution and ab initio simulations , Plasma Phys. Control.Fusion , 15004, (2009). Bibliography[11] P. M. Bellan,
Fundamentals of Plasma Physics , Cambridge University Press,(2009).[12] M. Berger,
XCOM: Photon Cross Sections Data-base , (1998).[13] M. Bergh, N. Tˆımneanu & D. van der Spoel,
Model for the dynamics of awater cluster in an x-ray free electron laser beam , Phys. Rev. E , 051904,(2004).[14] J.J. Binney, Anisotropic gravitational collapse , Astroph. J. , 492B (1977).[15] J.J. Binney & S. Tremaine
Galactic Dynamics , 2 nd ed., Princeton UniversityPress NJ, (2008).[16] J. Bogan, W. H. Benner, S. Boutet, U. Rohner, M. Frank, A. Barty, M. M.Seibert, F. Maia, S. Marchesini, S. Bajt, B. Woods, V. Riot, S. P. Hau-Riege,M. Svenda, E. Marklund, E. Spiller, J. Hajdu, & H. N. Chapman, SingleParticle X-ray Diffractive Imaging , Nano Lett. , 310 (2008).[17] C.M. Boily, C.J. Clarke & S.D. Murray, Collapse and evolution of flattenedstar clusters
Mon. Not. R. Astron. Soc. , 399 (1999).[18] R. Boll, D. Anielski, C. Bostedt, J. D. Bozek, L. Christensen, R. Coffee,S. De, P. Decleva, S. W. Epp, B. Erk, L. Foucar, F. Krasniqi, J. K¨upper,A. Rouz´ee, B. Rudek, A. Rudenko, S. Schorb, H. Stapelfeldt, M. Stener, S.Stern, S. Techert, S. Trippel, M. J. Vrakking, J. Ullrich & D. Rolles,
Femtosec-ond photoelectron diffraction on laser-aligned molecules: Towards time-resolvedimaging of molecular structure , Phys. Rev. A , 061402 (2013).[19] J. P. Boris, Relativistic Plasma Simulations - Optimization of a Hybrid Code ,in Proceedings of the 4th Conference on Numerical Simulations of Plasmas, Naval Research Lab. eds, Washington DC (1970).[20] F. Bouchet, S. Gupta& D. Mukamel,
Thermodynamics and dynamics of sys-tems with long-range interactions , Phys. A , 4389 (2010).[21] S. Boutet & G. V. Williams,
The Coherent X-ray Imaging (CXI) instrumentat the Linac Coherent Light Source (LCLS) , New J. Phys. , 035024 (2012).[22] B. H. Brandsen & C. J. Joachim, Physics of Atoms and Molecules , 2 nd edit.Pearson Education associated, Essex (2003).[23] P. H. Bucksbaum & J. M. Glownia, Ultrafast AMO physics at the LCLS x-rayFEL , EPJ Web of Conferences, , 04001 (2013).[24] F. Buersgens, K. W. Madison, D. R. Symes, R. Hartke, J. Osterhoff, W.Grigsby, G. Dyer, & T. Ditmire, Angular distribution of neutrons from deuter-ated cluster explosions driven by femtosecond laser pulses
Phys. Rev. E ,016403 (2006).ibliography [25] V. Y. Bychenkov & V. F. Kovalev, Coulomb explosion in a cluster plasma ,Plas. Phys. Rep. , 178 (2005).[26] V. Y. Bychenkov & V. F. Kovalev, Relativistic coulomb explosion of sphericalmicroplasma , JETP , 97 (2011).[27] A. M. Bykov, Shocks and particle acceleration in SNRs: theoretical aspects ,Advances in Space Research , 366 (2004).[28] A. Cavaliere & A. Messina, Propagation of blast waves , Astroph. J. , 424(1976).[29] S. Chandrasekhar,
Principles of stellar dynamics , Dover publications Inc.,New York (1943).[30] S. Chandrasekhar,
Ellipsoidal figures of equilibrium , Dover publications Inc.,New York (1987).[31] H. N. Chapman, A. Barty, M. J. Bogan, S. Boutet, M. Frank, S. P. Hau-Riege,S. Marchesini, B. W. Woods, S. Bajt, W. H. Benner, R. A. London, E. Pl¨onjes,M. Kuhlmann, R. Treusch, S. D¨usterer, T. Tschentscher, J. R. Schneider, E.Spiller, T. M¨oller, C. Bostedt, M. Hoener, D. A. Shapiro, K. O. Hodgson,D. van der Spoel, F. Burmeister, M. Bergh, C. Caleman, G. Huldt, M. M.Seibert, F. R. Maia, R. Lee, A. Sz¨oke, N. Tˆımneanu, & J. Hajdu,
Femtoseconddiffractive imaging with a soft-X-ray free-electron laser , Nature Physics , 839(2006).[32] P. H. Chavanis, Dynamics and thermodynamics of systems with long-range in-teractions: interpretation of the different functionals , in “Dynamics and Ther-modynamics of Systems with Long Range Interactions: Theory and Experi-ments”, Campa et al eds., AIoP Conf. Series , 39 (2008).[33] G. Chen, C. Wang, H. Lu, J. Liu, R. Li, G. Ni & Z. Xu,
Pure Coulomb explo-sions of highly charged methane clusters investigated by a simple electrostaticmodel , J. Phys. B , 105601, (2008).[34] L. Ciotti, Dynamical Friction from field particles with a mass spectrum , in“Plasmas in the laboratory and in the Universe”, Bertin et al. eds., AIoP Conf.Series , 117 (2010).[35] D. Comparat, T. Vogt, N. Zahzam, M. Mudrich & P. Pillet,
Star clusterdynamics in a laboratory: electrons in an ultracold plasma , Mon. Not. R. Astr.Soc. , 1227 (2005).[36] J. E. Crow, P. L. Auer & J. E. Allen,
The expansion of a plasma into avacuum , J. Plas. Phys. , 65 (1975). Bibliography[37] O. G. Danylchenko, S. I. Kovalenko & V. N. Samovarov,
Experimental verifi-cation of the Hagena relation for large clusters formed in conical nozzle , Tech.Phys. Lett. , 1037 (2008).[38] R.C. Davidson Physics of non-neutral plasmas , Imperial College Press, Lon-don (2001).[39] W. Dehnen,
A Family of Potential-Density Pairs for Spherical Galaxies andBulges , Mon. Not. R. Astron. Soc. , 250 (1993)[40] W. Dehnen,
Towards optimal softening in three-dimensional N-body codes -I. Minimizing the force error , Mon. Not. R. Astron. Soc. , 273 (2001).[41] W. Dehnen & J. I. Read,
N-body simulations of gravitational dynamics , EPJPlus , 55 (2011).[42] P.T. de Zeeuw,
A generalization of Kuzmin’s theorem , Mon. Not. R. Astron.Soc. , 599 (1985)[43]
P. F. Di Cintio & L. Ciotti,
Relaxation of spherical systems with long-rangeinteractions: a numerical investigation , Int. J. Bif. Cha. , 2279 (2011).[44] P. F. Di Cintio , L. Ciotti & C. Nipoti,
Relaxation of N-body systems withadditive r − α interparticle forces , Mon. Not. R. Astr. Soc. , 3177 (2013).[45] P. F. Di Cintio , U. Saalmann & J.M. Rost,
Proton ejection from molecularhydride clusters exposed to strong X-ray pulses , Phys. Rev. Lett. , 124401(2013).[46] T. Ditmire, T. Donnelly, A. M. Rubenchik, R. W. Falcone M. D. Perry,
Interaction of intense laser pulses with atomic clusters , Phys. Rev. A , 3379(1996).[47] T. Ditmire, M. Hohenberger, D. R. Symes, K. W. Madison, F. Buersgens,R. Hartke, J. Osterhoff, A. Henig, & A. Edens, Explosions of Methane andDeuterated Methane Clusters Irradiated by Intense Femtosecond Laser Pulses ,in Superstrong Fields in Plasmas, American Institute of Physics ConferenceSeries , 109, M. Lontano & D. Batani eds. NY (2006).[48] T. Ditmire, P. K. Patel, R. A. Smith, J. S. Wark, S. J. Rose, D Milathianaki,R.S. Marjoribanks & M. H. R. Hutchinson, keV x-ray spectroscopy of plasmasproduced by the intense picosecond irradiation of a gas of xenon clusters , J.Phys. B , 2825 (1998).[49] T. Ditmire, J. W. G. Tisch, E. Springate, M. B. Mason, N. Hay, R. A. Smith,J. Marangos & M. H. R. Hutchinson, High-energy ions produced in explosionsof superheated atomic clusters , Nature , 54 (1997).ibliography [50] G. W. F. Drake (editor),
The Springer Handbook of Atomic, Molecular andOptical Physics , Springer Verlag, Berlin (2006).[51] B. Erk, D. Rolles, L. Foucar, B. Rudek, S. W. Epp, M. Cryle, C. Bostedt, S.Schorb, J. Bozek, A. Rouzee, A. Hundertmark, T. Marchenko, M. Simon, F.Filsinger, L. Christensen, S. De, S. Trippel, J. K¨upper, H. Stapelfeldt, S. Wada,K. Ueda, M. Swiggers, M. Messerschmidt, C. D. Schr¨oter, R. Moshammer, I.Schlichting, J. Ullrich, & A. Rudenko,
Ultrafast Charge Rearrangement andNuclear Dynamics upon Inner-Shell Multiple Ionization of Small PolyatomicMolecules ”, Phys. Rev. Lett. , 53003 (2013).[52] T. Z. Esirkepov,
Exact charge conservation scheme for Particle-in-Cell simu-lation with an arbitrary form-factor , Computer Physics Communications ,144 (2001).[53] S.A.E.G. Falle,
Stability of a uniform, non-rotating spheroid collapsing underits own gravitation , Mon. Not. R. Astron. Soc. , 265E (1972).[54] M. Fellhauer,
Particle-Mesh Techniques and Superbox in The Cambridge N − body Lectures, Lecture Notes in Physics , 159, Springer Verlag, Berlin(2008).[55] T. Fennel, K. H. Meiwes-Broer, J. Tiggesb¨aumker, P. G. Reinhard, P. M.Dinh & E. Suraud, Laser-driven nonlinear cluster dynamics , Rev. Mod. Phys. , 1793 (2010).[56] V. S. Filinov, P. R. Levashov, A. V. Bot¸an, M. Bonitz & V. E. Fortov, Ther-modynamic properties and electrical conductivity of strongly correlated plasmamedia , J. Phys. A , 214002 (2009).[57] G. Fubiani, J. Qiang, E. Esarey, W.P. Leemans & G. Dugan, Space chargemodeling of dense electron beams with large energy spreads , Phys. Rev. Spec.Top. - Accelerators and Beams , 4402 (2006)[58] K. G. Gaffney & H. N. Chapman, Imaging Atomic Structure and Dynamicswith Ultrafast X-ray Scattering , Science , 1444 (2007).[59] I. Georgescu, U. Saalmann & J. M. Rost,
Clusters under strong vuv pulses: Aquantum-classical hybrid description incorporating plasma effects , Phys. Rev.A , 43203 (2007),[60] C. Gnodtke, U. Saalmann & J. M. Rost, Ionization and charge migrationthrough strong internal fields in clusters exposed to intense x-ray pulses , Phys.Rev. A , 41201, (2009).[61] C. Gnodtke, U. Saalmann & J. M. Rost, Dynamics of photo-activatedCoulomb complexes , New J. Phys. , 013028 (2011). Bibliography[62] C. Gnodtke, U. Saalmann & J. M. Rost,
Slow and fast multi-photon ionizationof clusters in strong XUV and X-ray pulses , Chem. Phys. 414, 65 (2013)[63] T. Gorkhover, M. Adolph, D. Rupp, S. Schorb, S. W. Epp, B. Erk, L. Foucar,R. Hartmann, N. Kimmel, K. U. K¨uhnel, D. Rolles, B. Rudek, A. Rudenko,R. Andritschke, A. Aquila, J. D. Bozek, N. Coppola, T. Erke, F. Filsinger, H.Gorke, H. Graafsma, L. Gumprecht, G. Hauser, S. Herrmann, H. Hirsemann,A. H¨omke, P. Holl, C. Kaiser, F. Krasniqi, J. H. Meyer, M. Matysek, M.Messerschmidt, D. Miessner, B Nilsson, D. Pietschner, G. Potdevin, C. Reich,G. Schaller, C. Schmidt, F. Schopper, C. D. Schr¨oter, J. Schulz, H. Soltau,G. Weidenspointner, I. Schlichting, L. Str¨uder, J. Ullrich, T. M¨oller & C.Bostedt,
Nanoplasma Dynamics of Single Large Xenon Clusters Irradiated withSuperintense X-Ray Pulses from the Linac Coherent Light Source Free-ElectronLaser , Phys. Rev. Lett. , 245005 (2012).[64] M. Grech, S. Skupin, R. Nuter, L. Gremillet & E. Lefebvre,
High-quality ionbeams from nanometric double-layer targets and their application to hadron-therapy , Nucl. Instr. and Meth. in Phys. Res. A , 63 (2010).[65] M. Grech, R. Nuter, A. Mikaberidze,
P. F. Di Cintio , L. Gremillet, E. Lefeb-vre, U. Saalmann, J. M. Rost & S. Skupin,
Coulomb explosion of uniformlycharged ellipsoids , Phys. Rev. E , 56404 (2011).[66] H. Grubm¨uller, H. Heller, A. Windemuth & K. Schulten, Generalized Verletalgorithm for efficient molecular dynamics simulations with long-range inter-actions , Mol. Sim. , 121 (1991).[67] F.J. Gr¨uner, C.B. Schroeder, A.R. Maier, S. Becker & J.M. Mikhailova, Space-charge effects in ultrahigh current electron bunches generated by laser-plasmaaccelerators , Phys. Rev. Spec. Top. - Accelerators and Beams , 701 (2009)[68] O. F. Hagena, Nucleation and growth of clusters in expanding nozzle flows ,Surf. Sci. , 101 (1981).[69] O. F. Hagena,
Cluster ion sources , Rev. Sci. Inst. , 2374 (1992).[70] S. P. Hau-Riege, R. A. London & A. Sz¨oke, Dynamics of biological moleculesirradiated by short x-ray pulses , Phys. Rev. E , 51906 (2004).[71] S. P. Hau-Riege, H. Sz¨oke, H. N. Chapman, A. Sz¨oke, S. Marchesini, A.Noy, H. He & M. Howells, SPEDEN: Reconstruction single particles from theirdiffraction patterns , Acta Crystall. A , 294 (2004).[72] S. P. Hau-Riege, R. A. London, H. N. Chapman, A. Sz¨oke & N. Tˆımneanu, Encapsulation and Diffraction-Pattern-Correction Methods to Reduce the Ef-fect of Damage in X-Ray Diffraction Imaging of Single Biological Molecules ,Phys. Rev. Lett. , 198302 (2007).ibliography [73] M. H´enon, Vlasov equation ? , Astron. & Astrophys. , 211 (1982).[74] L. Hernquist,
Performance Characteristics of Tree-Codes , Astrophys. J. Supp.Series , 751 (1987).[75] R. W. Hockney & J. W. Eastwood, Computer simulation using particles ,Taylor & Francis eds. (1998).[76] M. Hohenberger, D. R. Symes, K. W. Madison, A. Sumeruk, G. Dyer, A.Edens, W. Grigsby, G. Hays, M. Teichmann, & T. Ditmire,
Dynamic Accel-eration Effects in Explosions of Laser-Irradiated Heteronuclear Clusters , PhysRev. Lett. , 195003 (2005).[77] M. Hoener, C. Bostedt, H. Thomas, L. Landt, E. Eremina, H. Wabnitz,T. Laarmann, R. Treusch, A. R. B. de Castro & T. M¨oller, Charge recombina-tion in soft x-ray laser produced nanoplasmas , J. Phys. B , 181001 (2008).[78] Z. Huang, Brightness and coherence of Synchrotron Radiation and FELs ,SLAC Pub. (2013).[79] L. Inhester, C. F. Burmeister, G. Groenhof & H. Grubm¨uller,
Auger spectrumof a water molecule after single and double core ionization , J. Chem. Phys. ,144304 (2012).[80] L. Inhester, G. Groenhof & H. Grubm¨uller,
Core hole screening and decayrates of double core ionized first row hydrides , J. Chem. Phys. , 164304(2013).[81] M. R. Islam, U. Saalmann & J. M. Rost,
Kinetic energy of ions after Coulombexplosion of clusters induced by an intense laser pulse , Phys. Rev. A , 41201(2006).[82] B. Iwan, J. Andreasson, M. Bergh, S. Schorb, H. Thomas, D. Rupp, T.Gorkhover, M. Adolph, T. M¨oller, C. Bostedt, J. Hajdu, & N. Tˆımneanu, Explosion, ion acceleration, and molecular fragmentation of methane clustersin the pulsed beam of a free-electron laser , Phys. Rev. A , 33201 (2012).[83] J. D. Jackson, Classical Electrodynamics rd edition, Wiliey & Sons eds.(1999).[84] R. L. Johnston, Atomic and Molecular Clusters , Taylor & Francis eds. NY(2002).[85] Z. Jurek & G. Faigel,
The effect of inhomogenities on single-molecule imagingby hard XFEL pulses , EPL , 68003 (2009).[86] N. Kandadai, K. Hoffmann, A. Helal, H. Thomas, J. Keto, J. Andreasson, P.F. Di Cintio , U. Saalmann, J. M. Rost, B. Iwan, M. Seibert, N. Tˆımneanu, J. BibliographyHajdu, M. Adolph, T. Gorkhover, D. Rupp, S. Schrob, T. M¨oller, G. Doumy,L. F. Di Mauro, M. Hoener, B. Murphy, N. Berrah, J. Bozek, C. Bostedt& T. Ditmire,
Explosion of Methane Clusters in Intense X-ray FEL Pulses ,submitted to Phys. Rev. Lett. (2014).[87] H. E. Kandrup,
Chaos and Chaotic Phase Mixing in Galaxy Evolution andCharged Particle Beams in Galaxies and Chaos, Lecture Notes in Physics ,154 Contopoulos, G. & Voglis, N. eds., Springer Verlag Berlin (2003).[88] H. E. Kandrup, I.V. Sideris, C.L. Bohn,
Chaos and the continuum limit incharged particle beams , Phys. Rev. Spec. Top. - Accelerators and Beams ,4202 (2004).[89] H. E. Kandrup, C. L. Bohn, R. A. Kishek, P. G. O’Shea, M. Reiser, & I.V. Sideris, Chaotic Collisionless Evolution in Galaxies and Charged-ParticleBeams , Annals of the New York Academy of Sciences , 12 (2005).[90] A. E. Kaplan, B. Y. Dubetsky & P. Shkolnikov,
Shock Shells in CoulombExplosions of Nanoclusters , Phys. Rev Lett. , 143401 (2003).[91] O.D. Kellogg, Foundations of Potential Theory , Dover publications Inc., NewYork (1954).[92] Y. Kiwamoto, J. Aoki & Y. Soga,
Potential distribution of a nonuniformlycharged ellipsoid , Phys. Plasmas , 4868 (2004).[93] P. Kolorenˇc, & V. Averbukh, K-shell Auger lifetime variation in doubly ion-ized Ne and first row hydrides , J. Chem. Phys. , 134414 (2011).[94] V. F. Kovalev & V. Yu. Bychenkov,
Kinetic description of the Coulomb Ex-plosion of a Spherically Symmetric Cluster , JETP , 212 (2005).[95] V. F. Kovalev, V. Y. Bychenkov & K. Mima
Quasimonoenergetic ion bunchesfrom exploding microstructured targets , Phys. Plasmas , 103110 (2007).[96] V. F. Kovalev, K. I. Popov, V. I. Bychenkov & W. Rozmus, Laser triggeredCoulomb explosion of nanoscale symmetric targets ”, Phys. Plasmas , 53103(2007).[97] S.R. Krishnan, L. Fechner, M. Kremer, V. Sharma, B. Fischer, N. Camus, J.Jha, M. Krishnamurthy, T. Pfeifer, R. Moshammer, J. Ullrich, F. Stienkemeier,M. Mudrich, A. Mikaberidze, U. Saalmann & J.M. Rost, Dopant-Induced Ig-nition of Helium Nanodroplets in Intense Few-Cycle Laser Pulses , Phys. Rev.Lett. , 3402 (2011).[98] I. Last, & J. Jortner,
Nuclear Fusion induced by Coulomb Explosion of Het-eronuclear Clusters , Phys. rev. Lett , 033401 (2001).ibliography [99] I. Last, & J. Jortner, Electron and nuclear dynamics of molecular clustersin ultraintense laser fields. I. Extreme multielectron ionization
J. Chem. Phys. , 1336, (2004).[100] I. Last, & J. Jortner,
Electron and nuclear dynamics of molecular clustersin ultraintense laser fields. II. Electron dynamics and outer ionization of thenanoplasma
J. Chem. Phys. , 1348, (2004).[101] I. Last, & J. Jortner,
Electron and nuclear dynamics of molecular clusters inultraintense laser fields. III. Coulomb explosion of deuterium clusters
J. Chem.Phys. , 3030, (2004).[102] I. Last, & J. Jortner,
Electron and nuclear dynamics of molecular clustersin ultraintense laser fields. IV. Coulomb explosion of molecular heteroclusters
J. Chem. Phys. , 8329, (2004).[103] S. Lee, W. Roseker, C. Gutt, B. Fischer, H. Conrad, F. Lehmk¨uhler, I.Steinke, D. Zhu, H. Lemke, M. Cammarata, M. D. Fritz, P. Wochner, M.Castro-Colin, S. O. Hruszkewycz, P. H. Fuoss, G. B. Stephenson, G. Gr¨ubel,& A. Robert,
Single shot speckle and coherence analysis of the hard X-ray freeelectron laser LCLS , Optics Express , 24647 (2013).[104] E. Lefebvre, N. Cochet, S. Fritzler, V. Malka, M. N. Al´eonard, J. F. Chemin,S. Darbon, L. Disdier, J. Faure, A. Fedotoff, O. Landoas, G. Malka, V. M´eot,P. Morel, M. Rabec LeGloahec, A. Rouyer, C. Rubbelynck, V. Tikhonchuk, R.Wrobel, P. Audebert, & C. Rousseaux, Electron and photon production fromrelativistic laser plasma interactions , Nuclear Fusion , 629 (2003).[105] Y. Levin, R. Pakter & T. N. Teles, Collisionless Relaxation in Non-NeutralPlasmas , Phys. Rev. Lett. , 040604 (2008).[106] E. M. Lifshitz & L. P. Pitaevskii,
Course of Theoretical Physics, vol. 10,Physical Kinetics , Butterworth-Heinemann eds., Burlington Ma (1981).[107] C.C. Lin, L. Mestel & F.H. Shu,
The Gravitational Collapse of a UniformSpheroid , Astrophys. J. ,1431 (1965).[108] M. Lundwall, W. Pokapanich, H. Bergersen, A. Lindblad, T. Rander, G.¨Ohrwall, M. Tchaplyguine, S. Barth, U. Hergenhahn, S. Svensson & O. Bj¨orne-holm,
Self-assembled heterogeneous argon/neon core-shell clusters studied byphotoelectron spectroscopy , Jour. Chem. Phys. , 214706 (2007).[109] D. Lynden-Bell,
Negative Specific Heat in Astronomy, Physics and Chem-istry , Phys. A , 293 (1999).[110] A. Macchi, M. Borghesi, M.Passoni,
Ion acceleration by superintense laser-plasma interaction , Rev. Mod. Phys. , 751 (2013). Bibliography[111] A. McPherson, B. D. Thompson, A. B. Borlsov, K. Boyer & C. K. Rhodes,
Multiphoton-induced X-ray emission at 4-5 keV from Xe atoms with multiplecore vacancies , nature , 631 (1994).[112] K. J. Mendham, N. Hay, M. B. Mason, J. W. G. Tisch & J. P. Marangos
Cluster-size distribution effects in laser-cluster interaction experiments , Phys.Rev. A , 055201 (2001).[113] D. Merritt, & T. Fridman, Triaxial Galaxies with Cusps , Astrophys. J. ,136 (1996).[114] D. Merritt,
Optimal Smoothing for N-Body Codes
Astronom. Jour. , 2462(1996).[115] R. E. Meyer,
An introduction to mathematical Fluid Dynamics , John Wiley& Sons eds. NY (1971).[116] A. Meza & N. Zamorano,
Numerical Stability of a Family of Osipkov-MerrittModels
Astrophys. J. ,136 (1997).[117] A. Mikaberidze, U. Saalmann & J.M. Rost
Energy absorption of xenon clus-ters in helium nanodroplets under strong laser pulses , Phys. Rev. A , 1201(2008)[118] A. Mikaberidze, U. Saalmann & J.M. Rost
Laser-Driven Nanoplasmas inDoped Helium Droplets: Local Ignition and Anisotropic Growth , Phys. Rev.Lett. , 128102 (2009).[119] A. Mikaberidze, Atomic and molecular clusters in intense laser pulses , Ph.D.Thesis, Technische Universit¨at Dresden, (2011).[120] P. Mora & R. Pellat,
Self-similar expansion of a plasma into a vacuum ,Phys. Fluids , 2300 (1979).[121] M. Murakami & M. M. Basko, Self-similar expansion of finite-size non-quasi-neutral plasmas into vacuum: relation to the problem of ion acceleration , Phys.Plasmas , 12105 (2006).[122] M. Murakami, nanocluster explosions and generation of quasimonoenergeticions in “laser-driven Relativistic Plasmas applied to Science Industry andMedicine - the 2 nd International Symposium”, P.R. Bolton, S. V. Bulanov &H. Daido editors, American Institute of Physics conf. series , 113 (2009).[123] M. Murakami & K. Mima, Efficient generation of quasimonoenergetic ionsby Coulomb explosions of optimized nanostructured clusters, Phys. Plasmas , 103108 (2009).ibliography [124] M. J. Nandor, M. A. Walker, L. D. van Woerkom & H. G. Muller, Detailedcomparison of above-threshold-ionization spectra from accurate numerical in-tegrations and high-resolution measurements , Phys. Rev. A , 1771 (1999).[125] A. J. Nelson, S. Toleikis, H. Chapman, S. Bajt, J. Krzywinski, J. Chalupsky,L. Juha, J. Cihelka, V. Hajkova, L. Vysin, T. Burian, M. Kozlova, R. R.F¨austlin, B. Nagler, S. M. Vinko, T. Whitcher, T. Dzelzainis, O. Renner, K.Saksl, A. Khorsand, P. A. Heimann, R. Sobierajski, D. Klinger, M. Jurek,J. Pelka, B. Iwan, J. Andreasson, N. Tˆımneanu, M. Fajardo, J. S. Wark, D.Riley, T. Tschentscher, J. Hajdu & R. V. Lee, Soft x-ray free electron lasermicrofocus for exploring matter under extreme conditions , Optics Express ,18271 (2009).[126] R. Neutze, R. Wouts, D. van der Spoel, E. Weckert, & J. Hajdu, Potential forbiomolecular imaging with femtosecond X-ray pulses
Nature , 752 (2000).[127] K. Nishihara, H. Amitani, M. Murakami, S. V. Bulanov & T. Z. Esirke-pov,
High energy ions generated by laser driven Coulomb explosion of cluster ,Nuclear Instruments and Methods in Physics Research A , 98 (2001).[128] V. N. Novikov, A. V. Brantov,V. Y. Bychenkov, & V. F. Kovalev,
Coulombexplosion of a heated cluster , Plasma Physics Reports , 920 (2008).[129] F. Peano, R. A. Fonseca & L. O. Silva, Dynamics and Control of Shock Shellsin the Coulomb Explosion of Very Large Deuterium Clusters , Phys. Rev. Lett. , 33401 (2005)[130] F. Peano, R. A. Fonseca, J L. Martins & L. O. Silva, Controlled shock shellsand intracluster fusion reactions in the explosion of large clusters , Phys. RevA , 53202 (2006).[131] F. Peano, J. L. Martins, R. A. Fonseca, L. O. Silva, G. Coppa, F. Peinetti,& R. Mulas, Dynamics and control of the expansion of finite-size plasmasproduced in ultraintense laser-matter interactions , Phys. Plasmas , 56704(2007).[132] S. Pellegrini, L. Ciotti & J. P. Ostriker, X-Ray Properties Expected fromActive Galactic Nucleus Feedback in Elliptical Galaxies , Astrophys. J. , 21(2012).[133] I. P´ocsik,
Lognormal distribution as the natural statistics of cluster systems
Z. Phys. D , 395 (1991).[134] L. Pommier & E. Lefebvre, Simulations of energetic proton emission in laser-plasma interaction , Laser and Particle Beams , 573 (2003).[135] K. I. Popov, V. Y. Bychenkov, W. Rozmus, V. F. Kovalev, & R. D. Sydora, Mono-energetic ions from collisionless expansion of spherical multi-speciesclusters ”, Laser and Particle Beams , 321 (2009). Bibliography[136] J. Posthumus (editor),
Molecules and Clusters in intense Laser fields , Cam-bridge University Press, Cambridge (2001).[137] W. H. Press, S. A. Teukolsky, W. T. Vetterling & B. P. Flannery,
NumericalRecipes: The Art of Scientific Computing , 3 rd edit. Cambridge University Press(2007).[138] H. Rein, & S. Tremaine, Symplectic integrators in the shearing sheet model
Mon. Not. R. Astron. Soc. , 3168 (2011)[139] H. M. Robbins,
An analytical study of the impulsive approximation. , AIAAJournal , 1417 (1966).[140] N. Rohringer & R. Santra, Resonant Auger effect at high x-ray intensity ,Phys. Rev. A , 53404 (2008).[141] J. M. Rost, Complex non-equilibrium dynamics in plasmas , Euro. Rev. ,249 (2009).[142] U. Saalmann, I. Georgescu & J. M. Rost, Tracing non equilibrium plasmadynamics on the attosecond time scale in small clusters , New J. Phys. ,025014 (2008).[143] U. Saalmann & J. M. Rost, Ionization of Clusters in Strong X-Ray LaserPulses , Phys. Rev. Lett. , 143401 (2002).[144] U. Saalmann & J. M. Rost, Ionization Dynamics of Atomic Clusters inIntense Laser Pulses in Photonic, Electronic and Atomic Collisions, 160P. D. Fainstein, M. A. P. Lima, J. E. Miraglia, E. C. Montenegro, & R. D. Ri-varola eds. (2006).[145] U. Saalmann, C. Siedschlag, & J. M. Rost,
TOPICAL REVIEW: Mecha-nisms of cluster ionization in strong laser pulses , J. Phys. B , 39 (2006).[146] U. Saalmann, A, Mikaberidze & J.M. Rost, Spatial correlations in FiniteSamples revealed by Coulomb explosion , Phys. Rev. Lett. , 133401 (2013).[147] C. Sack & H. Schamel,
Plasma expansion into vacuum - a hydrodynamicalapproach , Phys. Rep. , 311 (1987).[148] Y. L. Shao, T. Ditmire, J. W. G. Tisch, E. Springate, J. P. Marangos & H.M. R. Hutchinson,
Multi-keV Electron Generation in the Interaction of IntenseLaser Pulses with Xe Clusters , Phys. Rev. Lett. , 3343 (1996).[149] E. Skopalov´a, I.C. El-Taha, A. Za¨ır, M. Hohenberger, E. Springate,J.W.G Tisch, R.A. Smith & J.P. Marangos, Pulse-Length Dependence of theAnisotropy of Laser-Driven Cluster Explosions: Transition to the ImpulsiveRegime for Pulses Approaching the Few-Cycle Limit , Phys. Rev. Lett. ,3401 (2010).ibliography [150] L. Spitzer Jr.
Physics of Fully Ionized Gases , 2 nd edition John Wiley & Sonsed. New York, (1961).[151] E. Springate, N. Hay, J. W. Tisch, M. B. Mason, T. Ditmire, H. M. Hutchin-son & J. P. Marangos, Explosion of atomic clusters irradiated by high-intensitylaser pulses: Scaling of ion energies with cluster and laser parameters , Phys.Rev. A , 063201 (2000).[152] P. St¨ackel, ¨Uber die Bewegung eines Punktes in einer n-fachen Mannig-faltigkeit , Mathematischen Annalen , 537 (1893).[153] T. H. Stix The theory of plasma waves , Mc-Graw Hill ed. New York, (1962).[154] A. Sugishima, H. Iwayama, S. Yase, H. Murakami, K. Nagaya, M. Yao, H.Fukuzawa, J. X. Liu, M. Motomura, K. Ueda, N. Saito, L. Foucar, A. Rudenko,M. Kurka, K. A. K¨uhnel, J. Ullrich, A. Czasch, R. D¨orner, R. Feifel, M.Nagasono, A. Higashiya, M. Yabashi, T. Ishikawa, T. Togashi, H. Kimura & A.Ohashi.
Charge and energy transfer in argon-core-neon-shell clusters irradiatedby free-electron-laser pulses at 62 nm , Phys. Rev. A , 33202 (2012).[155] E. Suraud, D. Cussol, C. Gr´egoire, D. Boilley, M. Pi, P. Schuck, B. Remaud& F. S´ebille, Explosions in Landau-Vlasov dynamics , Nuclear Phys. A , 73(1989).[156] H. Takeuchi,
The structural investigation on small methane clusters de-scribed by two different potentials , Comp. Theo. Chem. , 48 (2012).[157] S. Ter-Avetisyan, M. Schn¨urer, H. Stiel, U. Vogt, W. Radloff, M. Karpov, W.Sandner & P.W. Nickles,
Absolute extreme ultraviolet yield from femtosecond-laser-excited Xe clusters , Phys. Rev. E , 036404 (2001).[158] W. Thirring, Systems with negative specific heat , Z. Phys. , 339 (1970).[159] H. Thomas, A. Helal, K. Hoffmann, N. Kandadai, J. Keto, J. Andreasson, B.Iwan, M. Seibert, N. Tˆımneanu, J. Hajdu, M. Adolph, T. Gorkhover, D. Rupp,S. Schorb, T. M¨oller, G. Doumy, L. F. DiMauro, M. Hoener, B. Murphy, N.Berrah, M. Messerschmidt, J. Bozek, C. Bostedt, & T. Ditmire,
Explosionsof Xenon Clusters in Ultraintense Femtosecond X-Ray Pulses from the LCLSFree Electron Laser , Phys. Rev. Lett. , 133401 (2012)[160] N. Tˆımneanu, B. Iwan, J. Andreasson, M. Bergh, M. Seibert, C. Bostedt, S.Schorb, H. Thomas, D. Rupp, T. Gorkhover, M. Adolph, T. M¨oller, A. Helal,K. Hoffmann, N. Kandadai, J. Keto, & T. Ditmire,
Fragmentation of clus-ters and recombination induced by intense and ultrashort x-ray laser pulses , inSociety of Photo-Optical Instrumentation Engineers (SPIE) Conference Seriesvol. , (2013). Bibliography[161] N. G. van Kampen
Stochastic Processes in Physics and Chemistry , NorthHolland eds. Amsterdam (1981).[162] J. L. Vay,
Simulations of beams or plasmas crossing at relativistic velocity ,Phys. Plasmas , 56701 (2008).[163] D. Vrinceanu, G. S. Balaraman & L. A. Collins, The King model for electronsin a finite-size ultracold plasma , J. Phys. A , 5501P (2008).[164] H. Wabnitz L. Bittner, A. R. B. de Castro, R. D¨ohrmann, P. G¨urtler, T.Laarmann, W Laasch, J. Schulz, A. Swiderski, K. von Haeften, T. M¨oller, B.Faatz, A. Fateev, J. Feldhaus, C. Gerth, U. Hahn, E. Saldin, E. Schneidmiller,K. Sytchev, K. Tiedtke, R. Treusch, & M. Yurkov, Multiple ionization of atomclusters by intense soft X-rays from a free-electron laser , Nature , 482(2002).[165] J. Yang, C. J. Hensley, & M. Centurion,
Ultrafast 3D imaging of isolatedmolecules with electron diffraction , in Society of Photo-Optical Instrumenta-tion Engineers (SPIE) Conference Series , 884502 (2013).[166] Z. Zhou, J. Liu, H. Lu, Z. Wang, J. Ju, C. Wang, C. Xia, W. Wang, A.Deng, Y. Xu, Y. Leng, G. Ni, R. Li & Z. Xu,