Structure and Phase Transitions of Yukawa Balls

Contributions to Plasma Physics, 26/10/2018 Structure and Phase transitions of Yukawa balls Henning Baumgartner1 ∗ , Dietmar Block2 , and Michael Bonitz1 1 arXiv:0901.4309v1 [physics.plasm-ph] 27 Jan 2009 2 Institute for Theoretical Physics and Astrophysics, Christian Albrechts University Kiel, Leibnizstrasse 15, D-24098 Kiel, Germany Institute for Experimental and Applied Physics, Christian Albrechts University Kiel, Leibnizstrasse 19, D-24098 Kiel, Germany Received XXXX, revised XXXX, accepted XXXX Published online XXXX Key words dusty plasmas, dust crystals, Yukawa balls, strong Coulomb correlations PACS 52.27.Gr, 05.30.-d,52.27.Lw In this review, an overview of structural properties and phase transitions in finite spherical dusty (complex) plasma crystals – so-called Yukawa balls – is given. These novel kinds of Wigner crystals can be directly analyzed experimentally with video cameras. The experiments clearly reveal a shell structure and allow to determine the shell populations, to observe metastable states and transitions between configurations as well as phase transitions. The experimental observations of the static properties are well explained by a rather simple theoretical model which treats the dust particles as being confined by a parabolic potential and interacting via an isotropic Yukawa pair potential. The excitation properties of the Yukawa balls such as normal modes and the dynamic behavior, including the time-dependent formation of the crystal requires, in addition, to include the effect of friction between the dust particles and the neutral gas. Aside from first-principle molecular dynamics and Monte Carlo simulations several analytical approaches are reviewed which include shell models and a continuum theory. A summary of recent results and theory-experiment comparisons is given and questions for future research activities are outlined. Copyright line will be provided by the publisher 1 Introduction Complex plasmas are gas plasmas consisting of electrons, ions, and neutral atoms that additionally contain microscopic particles with sizes ranging from 10nm to some 10µm. This state of matter is ubiquitous in space, e.g. in the interplanetary medium, in interstellar clouds, in comet tails, and in the ring systems of the giant planets as well as in mesospheric noctilucent clouds [1, 2]. At the same time, in microchip manufacturing, avoiding particle contamination during the many production steps that involve plasmas is a technological challenge [3, 4]. On the other hand, the growth, transport, and deposition of nanoparticles is the central goal of many plasma deposition techniques, e.g. in the manufacturing of amorphous solar cells [5–8]. Since the early 1990s, research in complex plasmas is rapidly evolving. The original name dusty plasmas is nowadays often replaced by “complex plasmas” – in analogy to complex fluids, in order to emphasize many fundamentally different properties of this medium that distinguish it from ordinary gas plasmas. The field of complex plasmas is now maturing and the interested reader can find timely monographs, e.g. on dusty plasmas in space [9], technological applications [10], or on waves in dusty space plasmas [11]. The microparticles in complex plasmas are electrically charged by collection of plasma electrons and ions as well as by photoemission or secondary electron emission [1, 2]. In laboratory plasmas usually the collection processes dominate and the particles attain a high negative charge of a few thousand elementary charges. To some extent, complex plasmas behave like negative ion plasmas. In both cases a heavy, negatively charged species substitutes part of the electrons and thus affects the mobility and the contribution to shielding by negative charges. Most ordinary plasmas in space and laboratory are weakly coupled, which means that the interaction energy of nearest neighbours is much smaller than their thermal energy. This situation is completely different in complex plasmas with micrometer sized particles which then become strongly coupled [12]: the negative charge can ∗ Corresponding author E-mail: [email protected], Phone: + 49 431 880-4063, Fax: + 49 431 880-4094 Copyright line will be provided by the publisher 2 H. Baumgartner, D. Block, and M. Bonitz: Structure of Yukawa balls be so large that the interaction energy of nearest neighbours exceeds their thermal energy by several orders of magnitude. A consequence of strong coupling is the formation of liquid or solid phases. In 1934, a liquid to solid phase transition had been predicted by Eugene W IGNER for the three-dimensional (3D) electron gas in metals at low temperature and low density due to strong Coulomb repulsion [13]. This phase transition became known as Wigner crystallization, and the solid phase as Wigner or Coulomb crystals. This kind of crystallization has been studied for decades in a varity of systems ranging from electrons trapped on top of liquid helium [14, 15], electrons trapped in semiconductor quantum wells [16] and quantum dots [17], strongly coupled radio-frequency discharges of dusty plasmas [18], hole crystallization in semiconductors [19] to laser-cooled trapped ions in Pauland Penning traps [20,21]. It is remarkable that these physically completely different systems exhibit very similar collective behavior which is due to the governing role of the long-range Coulomb interaction. What makes the dust systems particularly attractive is that here collective phenomena, such as crystallization, occur at a comparatively large length scale of several millimeters and at room temperature which makes detailed experimental investigations rather simple. The formation of ordered solid phases in dusty plasmas was predicted by I KEZI [22] and first observed as plasma crystals in the early 1990s [23–25]. Since then crystalline structures and dynamical processes have been observed in a variety of confinement geometries ranging from 1D to 3D, see e.g. Ref. [26] for an overview and additional references. The present review will focus on the 3-dimensional spherical dust crystals recently observed by A RP et al. [27]. Their crystalline structure was found to be very similar to those observed before in ion crystals in Paul traps and the crystals which have been predicted to exist in expanding neutral plasmas [28, 29]. Compared to the latter, the advantage of these so called Yukawa balls is that the individual particle positions and trajectories can be directly visually observed and traced and recorded with video cameras. As a result, very detailed comparisons to theory and computer simulations can be performed which go far beyond comparisons in more traditional plasmas. In particular, the comparison is not limited to average quantities such as distribution functions on long time scales, but also covers correlation effects, fluctuations, and even short-time effects become accessible. Thus, the aim of this review is to show, on the example of Yukawa balls, that complex plasmas are an excellent object to study in great detail static and dynamic phenomena in finite systems such as ground and metastable states, transitions between configurations and phase transitions. These properties have, in the past few years, been studied in a very close connection between experiment, theory and computer simulation. For this purpose the review is organized as follows: In Section 2 we discuss the experimental realization of the Yukawa balls. Then, in Section 3 we give an introduction to the theory and computer simulations and discuss the low-temperature structural properties of Yukawa balls. Section 6 is devoted to finite temperature excitations of the crystals, to structural transitions and the melting behavior. We conclude this paper with a discussion of the results and an outline of further research directions in Sec. 7. 2 Experimental Realization of Yukawa balls To realize and investigate spherical 3D dust clouds two ingredients are indispensable: an isotropic confinement to prevent Coulomb explosion of the alike charged particles and powerful 3D diagnostics. For a detailed report on the progress in the field of 3D diagnostic tools we refer the reader to [30]. Here, we will limit the discussion to dust confinement, because the shape of the confinement potential is a central input parameter for theory and simulation. To confine dust particles inside a laboratory plasma reactor, the gravitational force has to be balanced. Thus, dust confinement is typically achieved in the plasma sheath region above an electrode where strong electric fields are present [31]. The trapping of the dust particles in horizontal direction is established by depressions of the electrode surface [32] or flat metal rings on the electrode [33]. However, this results in a very anisotropic confinement potential. The confinement in vertical direction is much stronger and thus these dust clouds are mostly 2D systems which nevertheless can form highly ordered crystals with a hexagonal lattice structure [23–25, 34]. Further, the supersonic ion flow towards the electrode is focused below each particle [35–39]. In multilayer systems, the resulting positive space-charge attracts particles in a lower layer. This process is responsible for chain formation observed in all dust clouds which are confined in regions of strong electric fields, i.e. regions with strong ion flows. This alignment vanishes if small particles and high gas pressure are used [40–42]. However, Copyright line will be provided by the publisher cpp header will be provided by the publisher 3 Fig. 1 (a) Side view of the discharge arrangement for Yukawa balls. The lower electrode is heated (T < 90◦ C). The vacuum vessel is grounded and kept at room temperature. The dust cloud is confined inside a glass cube where the upper and lower side are left open. The inset shows an image of a large dust cloud with 1 cm in diameter. (b) A thin slice at the front side of the cloud is illuminated. The particles basically arrange in a hexagonal lattice. After Ref. [27]. extended homogeneous 3D plasma crystals cannot be generated this way and investigations of dust dynamics are not feasible due to strong damping. To produce extended 3D dust clouds, various approaches were followed. M ERLINO and coworkers confined dust in a magnetized anodic plasma [43, 44] and investigated dust acoustic waves. Later, similar experiments on dust dynamics were performed in other discharges [45–51]. To produce 3D plasma crystals, a number of experiments have been performed under microgravity conditions. These experiments have provided many interesting observations, e.g. of localized crystalline structures [52], of complex plasma boundaries [53, 54], of coalescence of complex plasma fluids [55], of transport properties [56,57], and of low-frequency waves and instabilities [58–61]. However, the most striking observation was the formation of a dust-free zone (void) in the center of the discharge [52,62]. It was proposed that the ion drag force is responsible for the formation of voids [63–66], and a combination of simulations [67–69], experiments [70–73], and recent ion drag models [74–76] were able to verify this. Although it was shown very recently that a void closure can be achieved [77], the formation of void-free crystalline dust clouds is still an important issue. Such dust clouds in a well defined confinement would allow to extend a large fraction of the successful research on 2D plasma crystals to 3D. Interesting observations were reported by A NNARATONE and coworkers [78, 79] who observed spherical dust clouds with less than 50 particles in a secondary discharge in front of an adaptive electrode. Unfortunately, these clouds are rather in a liquid state and their confinement is not yet understood. A different approach to confine void-free dust clouds is motivated by experiments, simulations and theories on void formation [62, 63, 73, 80]. Those results showed that a void-free dust confinement at moderate plasma densities is hampered by the inherent ion flows towards the discharge boundaries. As a consequence, A RP and coworkers [27] have modified the usual setup of an asymmetric capacitively coupled rf-discharge [Fig. 2(a)]. To establish a vertical temperature gradient the rf-electrode was heated to T = (60 − 80)◦ C. The grounded vacuum vessel was kept at room temperature. The resulting thermophoretic force should balance gravity for particles with a radius rp < 3 µm, i.e. allow for dust levitation in the bulk plasma. In addition to the experimental setup used by ROTHERMEL [73], Arp et al. proposed to place a glass box with square cross section on the electrode that was left open at the top and at the bottom. For low discharge power (P < 10 W), they discovered that dust, which was dropped into the glass box, formed spherical void-free clouds [27]. A picture of a laser illuminated dust cloud (Yukawa ball) consisting of roughly 10, 000 dust particles is shown in the inset of Fig. 2(a). Already if only the front of such a dust cloud was illuminated with a laser sheet, interesting structural properties were observed [Fig. 2(b)]. The particle arrangement was static and particles on the surface of Yukawa balls seemed to have preferably 5 or 6 nearest neighbors. Moreover, the particles did not form vertical chains, i.e. first evidence was found that the Yukawa balls were in a solid or even crystalline state. To achieve a quantitative description of the dust confinement, fluid simulations with the S IGLO -2D code [81, 82] in combination with additional experiments were performed [83]. The results are summarized in Fig. 2. Copyright line will be provided by the publisher 4 H. Baumgartner, D. Block, and M. Bonitz: Structure of Yukawa balls Fig. 2 Experiments on the confinement of Yukawa balls: (a) Superposition of multiple video frames. As soon as the discharge is switched off, the trapped particles fall down. Their motion is only affected by gravity and the thermophoretic force. (b) Vertical component of the experimentally determined thermophoretic force field. In ’red’ regions the thermophoretic force exceeds gravity, in ’blue’ regions gravity is dominant. The solid line shows where both forces balance. (c) Trap potential obtained from PIV measurements and fluid simulations. (d) A horizontal section through the trap center reveals an almost parabolic confinement. From Ref. [83]. It was found that thermophoresis alone was not sufficient to explain the confinement [Fig. 2(a)]. By means of particle imaging velocimetry (PIV) experiments the thermophoretic force field was measured [Fig. 2(b)]. It was found that thermophoresis contributed about 70 percent to the vertical confinement and that the radial confinement and the remaining 30 percent in vertical direction were provided by plasma induced forces, namely electrostatic fields due to surface charges on the glass [83]. Further, the fluid simulations showed that the Yukawa balls were confined in a region where the plasma production was negligible. Due to the low plasma density (n ≈ 1013 m−3 ) and the weak electric fields the ion drifts were significantly below the ion sound speed. The resulting ion drag force was about two orders of magnitude smaller than the electrostatic confinement force and could, therefore, be neglected and chain formation is not expected. The confinement provided by the combination of gravity, thermophoresis and electric fields is plotted in Fig. 2(c) and is in very good agreement with the experimental observations [Fig. 2(a)]. It yields the correct levitation height, its depth is sufficient to trap particles with a few thousand elementary charges and even its asymmetry was observable for huge clouds in the experiment [83]. For further structural analysis and especially for comparison with simulations as well as with other strongly coupled systems, it is important to note that the trapping potential is essentially isotropic and parabolic for dust clouds with a radius of less than 2 mm, i.e. N < 1000 particles. Furthermore, for the trapping process an interaction of dust and plasma, as proposed by T OTSUJI [84, 85], is not required. The particle trap is provided by external forces only and thus this trap geometry is closely related to those of Penning and Paul traps [86]. This allows for a direct comparison of the structural properties of Yukawa balls and e.g. trapped laser-cooled ions. Copyright line will be provided by the publisher cpp header will be provided by the publisher 3 5 Theoretical Description of Yukawa balls A theoretical description of complex plasmas is very complicated, due to the large number of different particle species, which include electrons, ions, neutrals and dust particles, and their huge difference in mass. Thus the length and time scales of the motion of the different particles differ by many orders of magnitude, making a selfconsistent modeling unfeasible. The most advanced computer models use particle in cell (PIC) simulations with a full account of the time-dependence of the discharge. [87]. However, this currently allows to treat only a very small number of dust particles in the plasma environment. A simplified model which treats the streaming electrons and ions within linear response and computes the dynamically screened potential of the dust particles was developed by J OYCE and L AMPE, see [88] and references therein. Both simulations yield an anisotropic and non-monotonic (wake) potential around a dust grain if the ion streaming velocity exceeds the sound speed which is normally the case in the plasma sheath. At the same time, simulations for the conditions of Yukawa balls which are produced in the plasma bulk revealed that there the streaming velocities are much lower and the dust-dust interaction is nearly isotropic and statically screened. It is, therefore, well approximated by a Yukawa potential. These findings justify the use of a much simpler theoretical model which is discussed in the following sections. 3.1 Theoretical Model A standard model to describe the Yukawa balls and their dynamics is based on the Hamiltonian H(r1 . . . , p1 , . . . rN , . . . pN , t) = N N N N X X X p2i q2 α 2 X ext U (ri , t) + + ri + · e−κrij , (1) 2m 2 4πǫr ij i=1 i=1 i=1 i6=j containing kinetic energy, a confinement and the interaction energy. The isotropic and parabolic confinement is based on the experimental results, cf. Sec. 2. q is the particle charge and α the confinement strength. For both values the experiment sets rather narrow boundaries [83]. Finally, ǫ is the background dielectric constant and rij is the distance between particles i and j. The effect of the surrounding plasma on the dust is condensed in the static screening parameter κ which will be treated as a parameter below. A second important effect of the medium on the dust particles can be described by a friction force. It is beyond the hamiltonian model but can be included in the equations of motion, see Sec. 3.2. In writing the formula (1) we assumed that all particles have identical mass and charge which is a good approximation of the experimental conditions. In fact, commercially available dust particles have a diameter which fluctuates by less than one percent, and the remaining charge variations were shown to have a negligable effect on the structure of Yukawa balls [89]. Thus, the present model is expected to be well appropriate to reproduce the experimental conditions. In Eq. (1) we have also included a possibly time-dependent potential U ext which describes external manipulations of the dust particles. The Hamiltonian (1) is brought to a simpler dimensionless form suitable for numerical analysis by introducing units of length and energy, r0 = (q 2 /16πǫα)1/3 and E0 = (αq 4 /32π 2 ǫ2 )1/3 , where r0 corresponds to the ground state distance of two particles in the trap in the absence of screening and E0 is the corresponding pair interaction energy. As p the third scale parameter we introduce the unit of time T0 which is the inverse of the trap frequency T0−1 = ω = α/m. 3.2 Computer simulations The model (1) is well suited for accurate computer simulations. Since the dust particles consitute a classical system, exact simulation methods are available which include molecular dynamics (MD), Langevin dynamics (LMD) and Metropolis Monte Carlo (MC) which have been successfully applied in recent years. Since most of the theoretical results in this paper are based on these methods we give a brief discussion. Molecular dynamics solves Newton’s equations corresponding to the hamiltonian (1), ṗi (t) ri (0) ~ ext (ri , t) + = −αri − ∇U = r0i , pi (0) = p0i , X j6=i q2 −κrij , 3 (1 + κrij )rij · e 4πǫrij i = 1, . . . N, (2) Copyright line will be provided by the publisher 6 H. Baumgartner, D. Block, and M. Bonitz: Structure of Yukawa balls where the last line contains the initial conditions. The coupled system (2) of ordinary differential equations is efficiently solved for up to several thousand particles by standard techniques such as Runge-Kutta methods. The present version gives a solution in the microcanonical ensemble where the total energy H is conserved. The initial conditions are chosen as close as possible to the experimental situation. This scheme has been used to obtain the ground and metastable states of the Yukawa balls, see Sec. 4.1 To reproduce the equilibrium properties of Yukawa balls which are characterized by a given temperature it is more appropriate to perform simulations in the canonical ensemble. This can be realized by Langevin MD where, in addition to the total force Fi on particle i [r.h.s. of Eq. (2)] a friction force describing the momentum loss of the particles to the neutral gas and a random fluctuating force modeling the collisions of the particles with the plasma medium (Brownian motion), ṗi (t) = Fi (t) − νmṙi + yi (t) i = 1, . . . N, hyi i = 0, hyiα (t)yjβ (t′ )i = 2Dδi,j δα,β δ(t − t′ ), (3) α, β = x, y, z, to be supplemented by the same initial conditions as in Eq. (2). Choosing the noise intensity according to the Einstein relation, D = mkB T ν, the system relaxes towards a Maxwellian velocity distribution with temperature T – tpyically the temperature of the neutral gas. This method has been used to obtain the probability of metastable states, cf. Sec. 6.2 and for nonequilibrium dynamical simulations of Yukawa balls [90]. Finally, a first principle method to study the thermodynamic equilibrium properties of the dust particles is Monte Carlo. In the canonical ensemble where the dust temperature Td and the particle number are assumed to be fixed, the probability of a certain configuration s of the Yukawa ball with the coordinates Rs = {rs1 , . . . rsN } is defined by Z s 1 P (Rs ; Td , N, α) = e−H0 (R )/kB Td , d3 r1 . . . d3 rN P (Rs ; Td , N ; α) = 1 (4) Z where Z is the normalization constant (partition function). Note that in contrast to thermodynamics of macroscopic systems, here the volume is not fixed but its role is taken over by the confinement strength. The exponent H0 of the probability distribution P is the hamiltonian (1) with the momenta put equal to zero [the momentum and coordinate distributions factorize, and the former is given by a Maxwellian and does not influence the particle configurations]. Also, the effect of dissipation is not included in this description. The equilibrium state of a Yukawa ball, is obtained by maximizing the probability P which is efficiently done using the Metropolis algorithm, e.g. [91]. MC methods have been applied to find the ground and metastable states, see Sec. 4.1, 6.2, and finite temperature effects and melting points, Sec. 6.4. 4 Structure of Yukawa balls Let us start the analysis of Yukawa balls by considering their ground state properties. They are obtained in the zero-temperature limit where all particles settle in the minima of the total potential energy H0 given by the sum of confinement and pair interactions. While zero temperature seems to contradict the experimental situation – the experiments are performed at room temperature – this is in fact a very good approximation. The reason is the very high charge collected on the grain surface leading to a potential energy per particle which exceeds the kinetic energy by a factor of 500 or more. Therefore, kinetic energy has a negligible effect on the stationary properties. In most cases the potential energy has several minima. In this case, the deepest minimum corresponds to the classical ground state whereas the others correspond to metastable states. Due to the spherical symmetry of the hamiltonian H0 the low temperature stationary states have the same symmetry: they consist of spherically symmetric concentric shells with a very small width which is in excellent agreement with the experiment. Different states differ by the radii and populations Nk of the shells. The particle arrangement within the shells is highly symmetric with a combination of five-fold and six-fold symmetry (particles dominantly have five or six nearest neighbors) [27, 92]. 4.1 Ground and metastable states Much work has been devoted to first-principle simulations of spherical Coulomb and Yukawa crystals. This is efficiently done by MD simulations by solving Eqs. (2), combined with “simulated annealing”. There one starts Copyright line will be provided by the publisher cpp header will be provided by the publisher 7 with a random initial configuration and reduces the velocity of all particles proportionally in small steps, until they have zero kinetic energy, [92, 93]. Alternatively, one can use LMD simulations, Eq. (3), with a very small damping coefficient ν, e.g. [90] or MC. The ground state for the case of pure Coulomb interaction was reported many years ago by Hasse et al. [94] where also earlier references are given, for larger crystals, see also [95]. These data correctly describe the general trend previously observed in spherical ion crystals. However, with the first experimental results on Yukawa balls in dusty plasmas very detailed comparisons of the precise shell structure became possible. Then it turned out that most data of Ref. [94] refered to metastable states rather than to the ground state. The first accurate data for spherical Coulomb crystals were obtained by L UDWIG et al. [92] who performed MD with simulated annealing. Extensive tables for N ≤ 160 and some larger clusters were presented in Ref. [96] and the results were later confirmed by other groups, e.g. [97, 98]. Fig. 3 (color online) Stationary states observed in the LMD simulations for N = 31, κ = 1.4 for cooling to a minimum ˙ min ¸ kinetic energy Ekin = 10−8 . The runs are sorted by the energy or the final state, see also table 1. For slow cooling (black bars, ν = 0.05) one can clearly distinguish different states, cf. the horizontal lines. The length of the bold lines is proportional to the occurrence frequency of a state in a total of 5000 runs. In the case of strong friction (red, dashed line, ν = 5.3) the particles often lose their kinetic energy before they can settle into the equilibrium positions and the fine structure cannot be resolved. From Ref. [90] ∆ E/N 3.030266 0.000006 0.000009 0.000291 0.000372 config. (27, 4) (27, 4) (27, 4) (26, 5) (26, 5) ∆ E/N 0.000479 0.000499 0.000530 0.000656 0.000669 config. (26, 5) (26, 5) (26, 5) (25, 6) (25, 6) Table 1 Energy difference between metastable states and the ground state (the ground state and its energy is given by bold numbers) as seen in Fig. 3. States with the same shell configuration but different energy differ only by the arrangement of the particles on the same shell (fine structure) [92]. First structural data for Yukawa balls have only recently been obtained, e.g. [99, 100]. The first systematic analysis was performed in Ref. [101] for small crystals with N ≤ 60 for a range of screening parameters of 0 ≤ κ ≤ 5.0. An example of the energetically lowest stationary states is given in Fig. 3 and in table 1 for a cluster with N = 31 and κ = 1.4. The table also contains the shell configuration of the states. As one can see in some cases there exist several states with the same configuration which differ only in the arrangement of the particles on the shells. This “fine structure” was first observed in Ref. [92] and revealed by a Voronoi analysis [102]. It is remarkable how small the energetic spacing of these states is which underlines the high numerical accuracy needed to find the ground states. Since no method is able to reveal the ground state with one Copyright line will be provided by the publisher 8 H. Baumgartner, D. Block, and M. Bonitz: Structure of Yukawa balls hundred percent reliability it is necessary to perform a large series of repeated calculations, as shown in Fig. 3, where 5000 runs were performed. Further, the results sensitively depend on the cooling rate in the annealing or, equivalently, on the magnitude of the friction in an LMD simulation. The ground state and the metastable states can only be resolved for sufficiently small damping, i.e. for slow cooling [90]. Further note that the number of metastable states increases exponentially with N which makes it very difficult to find the ground state for clusters with more than about 1000 particles. 160 κ = 0.0 κ = 0.3 κ = 0.6 κ = 1.0 Ns 120 80 40 0 20 30 40 2/3 N Fig. 4 (color online) Experimental (symbols) and theoretical (lines) shell population NS of the ground states versus system size. The MD results are obtained for different screening parameters κ and show that the particles redistribute towards inner shells with increased κ. From Ref. [99] 4.2 Effect of screening A comparison of experimental and theoretical cluster configurations is shown in Fig. 4 for particle numbers in the range of N ≤ 500. In addition to the measured configurations of 43 clusters MD simulations of the ground states of the same particle numbers and with different values of κ were performed [99], using simulated annealing, as described in Sec. 4.1 above. Comparison of the results shows that the model of isotropic static screening in fact reproduces the measured shell occupations very well. Furthermore, it is obvious that the measurements cannot be explained by pure Coulomb interaction. Since κ is the only free parameter in the simulations within the model (1), comparison to the measurements allows for reliable predictions of the screening strengths in the experiments which is not directly accessible. The result of κ ≈ 0.6, corresponding to λD r̄ = 1.5, gives a value which agrees well with other estimates and the accuracy was estimated to be within 50 percent [99]. Sources of errors such as charge fluctuations and finite temperature were found to be small [89] and also the occurence of metastable state in the experiments, see Sec. 6.2, has a negligible influence on this result. For completeness, we note that Ref. [99] also provided detailed comparisons of the experimental and theoretical shell radii Rs . It turned out that the radii are almost equidistant and increase approximately as N 1/3 , as one would expect in bulk matter. Most importantly, while the absolute value of the radii decreases with screening (due to the reduced interparticle repulsion), no effect of screening is observed for the ratio of Rs to the mean interparticle distance within a shell, for fixed N . Since at the same time, the particle number on the outer shells decreases with κ, this has important consequences for the mean density distribution within the Yukawa balls which will be considered more in detail in Sec. 5.2. There are two general effects of screening on the structure of Yukawa balls: first, increase of screening leads to cluster compression and, second, to increased occupation of inner shells. However, this is only the average trend, and noteble deviations have been found. A detailed numerical analysis of small Yukawa balls in a broad range of screening parameters has revealed a number of interesting anomalies [101], the shell populations are shown in Fig. 5. In general the particle number on the inner shell increases one by one, when κ is increased. Copyright line will be provided by the publisher cpp header will be provided by the publisher 9 Fig. 5 (color online) Ground states of small Yukawa balls vs. screening parameter. The white numbers denote the number of particles on the inner shell(s). The black circles indicate anomalies of the 1st kind. The white circles indicate the end of the screening range, where anomalies of the 2nd kind appear. The ground states for a screening parameter κ = 20.0 are plotted above the diagram. The cyan bar for N = 44 at κ = 20.0 refers to a ground state of (11, 1) in the center region; it is the only time this configuration is part of a ground state. The small squares just below κ = 20.0 indicate anomalies of the 3rd kind, where a ground state configuration reappears with increased screening. From Ref. [101] However, in several cases, two particles move inward simultaneously, see for example N = 30 at κ ≈ 1.5 where the configuration changes according to (26, 4) → (24, 6). These correlated transitions have been called “anomaly of the 1st kind” and are observed 18 times. A second anomaly is observed when N is increased by one at constant κ. The normal trend is an increase of the inner or outer shell population by one. However, in some cases the inner shell population decreases by one whereas the outer shell population increases by two. This has been observed for Coulomb balls in a single case N = 59 → 60 [103] when the configuration changes according to (46, 12, 1) → (48, 12). In the case of Yukawa balls this “2nd anomaly” is observed frequently, see e.g. N = 11 → 12 for 2 ≤ κ ≤ 4 where the configuration changes according to (10, 1) → (12, 0). Finally, there are five cases (N = 35, 36, 37, 39, 54) of reentrant shell configuration changes where for a fixed N and increasing κ a particle moves from an inner to an outer shell restoring a configuration which existed before, at a smaller value of κ. The reason for these non-trivial configuration changes are special symmetries of the particle arrangement which allow to reduce the total energy of the Yukawa ball. A table containing a complete list of the anomalies is given in Ref. [101]. 5 Analytical Models for the structure of Yukawa balls While computer simulations provide virtually exact structural data for a certain hamiltonian a lot of additional information and physical insight can be gained from analytical models. Appropriate simplification of the system allows to observe the dominant trends and effects. Here we summarize the results of two such models – shell models and statistical models based on a fluid-type continuum description. We start by considering the Coulomb Copyright line will be provided by the publisher 10 H. Baumgartner, D. Block, and M. Bonitz: Structure of Yukawa balls limit which is simpler – the ground state configurations of small clusters with N ≤ 8 can be found analytically [104] – and then proceed to Yukawa interaction. 5.1 Shell models The first idea of a simplified picture of the stationary particle configuration is to use a shell model, where the particles are confined to concentric shells with zero width. The first shell model was derived by AVILOV and H ASSE [94] and was improved by T SURUTA and I CHIMARU [103] who included correlation effects. K RAEFT and B ONITZ further improved this model and presented detailed comparisons with MD simulations [105]. They gave the following form of the total energy per particle for a system of L homogeneously charged concentric spherical shells L X Emodel (N ) Nν 1/3 = 2 N (Ze)2 /r0 N xν µ=1 " ! # √ Nν − ǫ Nν X 9 1 3 + Nµ + ζ + xν − N 2/3 , 2 2 10 µ<ν (5) where ζ = 0, 1 accounting for a particle in the trap center and xµ is the radius of the shell µ in units of r0 . Here, the first term (proportional to Nν2 ) is the surface energy of a spherical capacitor of radius a, containing Nν charges, Esurf (N ) = (N (N − 1)e2 )/(2a). The contribution from the sum over µ accounts for the electrostatic interaction of the shell ν with all inner shells. The x3ν contribution describes the confinement energy and the last term takes into account a (possible) compensating background and is not relevant for the structure of the Coulomb ball. This model, without the term proportional to ǫ, can be rigorously derived from a mean field theory [106] which is discussed in Sec. 5.2 and was given by H ASSE et al. [94] although they missed the −1 correction to the surface energy Esurf . The cluster configuration can now be derived simply by an optimization procedure searching for those shell populations and radii which minimize the total energy (5) which is much simpler than to solve the exact problem. While this yields the correct qualitative trend, however, with ǫ = 0 there are large quantitative deviations. T SU RUTA and I CHIMARU corrected the Hasse model by introducing a correlation correction (they used ǫ = 1) which reflects the discreteness of the particles and takes into account that the area around a particle an a shell cannot be occupied by others. Using ǫ = 1 allows to reproduce the shell populations within 5 percent [105]. However, the exact expression for the correlation corrections is unknown, and there is no obvious reason why this coefficient should be equal to one. In fact, a further improvement of the quality of the model is achieved by leaving ǫ as an additional free parameter that is optimized to minimze the energy. This has been done in Ref. [105] and allows to reduce the deviations from the exact ground states to (1 . . . 2)%. The resulting values for the correlation parameters were slightly above one and converged to ǫ = 1.104 for large N , a result which could be recently derived from the Thomson model by C IOSLOWSKI [108, 109]. For Yukawa balls the situation is more complex. A mean field-type model (ǫ = 0) was recently obtained by T OTSUJI et al. [84] and derived from continuum theory by H ENNING et al. [106], L X α 2 e−κRν Rν + Q2 × 2 Rν ν=1 √  X sinh(κRµ )  sinh(κRν ) Nν − ǫν (N, κ) Nν Nµ +ζ + , κRν 2 κRµ µ<ν Emodel (N ; κ) = Nν  (6) where, compared to Eq. (5), the background term (last term) has been dropped and we multiplied with the energy unit (denominator on the l.h.s.). One readily confirms by performing the limit κ → 0 that this result includes the Coulomb case. The model (6) differs from Refs [84, 106] by the additional correlation corrections which generalize the Coulomb expression (5). This model was used in Ref. [107] and optimized to minimize the total energy. The detailed comparison to exact simulation data showed that, for Yukawa balls, the mean field model (neglecting the ǫ terms) performs rather poorly, but satisfactory agreement with an accuracy of the shell populations within 5% can be achieved only by including correlation corrections and using different ǫ-values for different shells. Yet a systematic derivation of the correlation corrections remains open. Copyright line will be provided by the publisher cpp header will be provided by the publisher 5.2 11 Continuum theory approach. Radial density profile of Yukawa balls An entirely different approach to the many-particle behavior of the Yukawa balls has been developed by K RAEFT and B ONITZ [110] and by H ENNING et al. in Refs. [106, 111]. It is based on a fluid-like statistical theory where the particles are treated as a continuum. This may seem far from reality but one can hope to describe several aspects of the Yukawa ball structure such as spatially averaged properties and achieve deeper physical insight. One such question is the mean radial density profile of the Yukawa balls. It is well known that, in a spherically symmetric parabolic potential, particles interacting via the Coulomb potential establish a radially constant-density profile (this follows from the fact that the potential profile inside a homgeneously charged sphere is a parabolic function of the distance from the center). However, in the case of Yukawa balls the pair interaction is screened, which leads to enhanced populations of inner shells with increaed κ while the shell radii (in units of the mean interparticle distance within the shells) does not change with κ, cf. Sec. 4.2. From this we may expect that radial distribution of matter in the Yukawa balls will be inhomogeneous and the density has to decay towards the edge. This has indeed been observed experimentally [112, 113] although the trend is weak, see Fig. 9. To analzye this question quantitatively we consider the model of a finite one-component plasma containing N identical particles described by the Hamiltonian (1). The classical ground state energy follows from neglecting in Eq. (1) all particle momenta and can be written as a functional of the density [106] Z H0 [n] = E[n] = d3 r u(r), (7) with the potential energy density   Z q2 N −1 d3 r2 n(r2 ) e−κ|r−r2 | + ucorr , u(r) = n(r) Φ(r) + 2N 4πǫ|r − r2 | (8) where the first term on the right is the mean-field contribution with the confinement potential Φ, and ucorr denotes the density of the correlation energy. The ground state density profile can be obtained from minimizing ! the R 3total energy E, i.e. from the solution of the variational problem 0 = δE[n]/δn(r) under the constraint d r n(r) = N =const.This problem can be solved for a general anisotropic confinement potential Φ(r) in a mean-field approximation, i.e., ucorr ≡ 0, with the result [106]: 4πq 2 n(r)(N − 1)/N = (∆ − κ2 )Φ(r) + κ2 µ, where µ is a Lagrange multiplier (chemical potential) assuring the normalization. 15 9 MF LDA MD (κdc = 0.0) MD (κdc = 1.0) MD (κdc = 2.0) 6 N = 1000 n(r)/nc 12 κdc = 2.0 κdc = 1.0 3 0 a) 0 1 κdc = 0.0 2 3 4 5 6 r/dc 7 8 9 10 Fig. 6 Radial density profiles of a 3D plasma of N = 1000 particles calculated with the mean-field model (solid lines) and with the LDA including correlation contributions (dashed lines) for Coulomb and Yukawa interaction with four different screening parameters: κdc = 1, 2, 3, and 5. Averaged shell densities of molecular dynamics simulations of a plasma crystal for the same parameters are shown by the symbols. dc is the two-particle equilibrium distance r0 , cf. Sec. 3 For a parabolic confinement, Φ(r) = α/2r2 , the ground state density profile is isotropic and parabolically decaying away from the trap center,   κ2 r 2 αN R2 κ2 3 + κR c − n(r) = Θ(R − r), c = 3 + , (9) 4π(N − 1)q 2 2 2 1 + κR Copyright line will be provided by the publisher 12 H. Baumgartner, D. Block, and M. Bonitz: Structure of Yukawa balls and the density drops to zero at a finite radius R(N, κ) which follows from the normalization Q2 Q2 (N − 1) − 15 κ(N − 1)R + 15R3 + 15κR4 + 6κ2 R5 + κ3 R6 = 0. (10) α α This equation has four complex and two real solutions, only one of which is non-negative, and thus constitutes the unique proper result entering equation (9). This density is equivalent to a local force balance [106, 114]. Figure 6 shows the density profiles derived from the mean field (MF) result (9) (solid lines) and compares them to the averaged shell densities of MD simulations (symbols). In the case of weak screening (a) the MF result agrees perfectly with the exact MD results, while for increased screening discrepancies show up, especially in the trap center. The origin of these deviations is the neglect of correlations. An extension of the previous analysis to include correlation effects is difficult as the form of the pair correlation function is unknown. One way around this problem is to apply a local density approximation (LDA) [111, 115]. There one replaces the nonlocal terms within the energy density at point r by local expressions using the known energy density of the homogeneous system [85]   4  −1/3 −1/3 corr 2 4/3 −5 u (n0 , κ) = −1.444Q n0 exp −0.375κn0 + 7.4 ∗ 10 , (11) κn0 −15 where n0 is the corresponding density of the homogeneous system. Using this expression in equation (8), leads to an improved equation for the density profile which has to be solved numerically. The resulting density profiles for a cluster with 1000 particles are included in Fig. 6. In the limit of strong screening, where the MF result failed to describe the simulation results, the LDA approximation with correlations is very accurate, see Fig. 6b. In contrast, for weak and moderate screening, κr0 ≥ 2, the LDA behaves very poorly. This is due to the long range of the interaction which cannot be accounted for within a local approximation. Thus, the mean-field model together with the presented local density approximation complement one another in the description of strongly correlated spatially confined one-component plasmas providing a complete theoretical description of the average ground state properties of Yukawa balls. An extension of these results to finite temperatures has been outlined in Ref. [116] 6 Structural transitions Up to this point we concentrated on the ground state properties of Yukawa balls, but as was shown, in Fig. 3 and table 1 also other stationary states – metastable states are frequently resolved in the simulations. Their energy is higher than that of the ground state but often only marginally. It may, therefore, be expected that these metastable states might be observable in the experiment and also transitions between stationary states might occur. We first analyze these questions theoretically and then compare with recent experimental results. 6.1 Coexistence of stationary states. Potential barriers As we have seen in Fig. 5, upon increase of κ the ground state configurations frequently change. This requires that, around some critical value κcr , both, the old and new configuration are stationary states, corresponding to minima of the total potential energy H0 . When κ departs from κcr one of the configurations (the metastable state) will eventually disappear. This is illustrated for N = 31 in the left part of Fig. 7. The right (left) arrows indicate the range of κ values where various metastable states appear (vanish). If stationary states co-exist one could expect to observe transitions between them which may occur spontaneously, driven by thermal fluctuations. However, this does not only depend on the local minima but also on the potential landscape between them, in particular on the potential barrier. The radial potential barriers can be “measured” in a computer simulation by moving, e.g. an outer particle to the inner shell, while allowing all other particles to relax. The results of such a MC simulation are shown in the right part of Fig. 7. They indicate that, at the chosen screening, transitions (6, 25) → (5, 26) and (5, 26) → (4, 27) occur practically spontaneously whereas the opposite transitions are hampered by a significantly higher barrier. This asymmetry his even more pronounced for the transition between the configurations (4, 27) and (3, 28). This is nicely confirmed by the experiments [113] where a slow transition (4, 27) → (5, 26) could be observed followed by a rapid return to the ground state (4, 27), see Fig. 9. Copyright line will be provided by the publisher cpp header will be provided by the publisher 12 (27,4)-(26,5) (26,5)-(25,6) ∆ E/N in 10-4 [E0] 4 10 2 0 E/N in [E0] 8 13 κcr = 1.61 2 κcr = 1.56 1 -2 6 (27,4) (26,5) (25,6) -4 1.4 1.5 1.6 κ 1.7 1.8 4 2 (27,4) (26,5) (25,6) 0 0 1 2 κ 3 4 5 Fig. 7 (colour online) Left: Energy per particle of a Yukawa ball with 31 particles vs. screenings strength. The ground state configuration is (27, 4) for κ ≤ 1.56, which is replaced by (26, 5) in the inverval 1.56 ≤ κ ≤ 1.61, followed by (25, 6) for larger κ. Inset shows the energy differences of the configurations. The vertical dashed lines with the arrows denote the beginning [blue for (26,5) and green for (25,6)] and the end [red for (27,4)] of occurance of these configurations in the simulations. Right: Potential energy landscape measured by a MC simulation where one particle is moved, while all other particles are allowed to relaxate into a new minimum energy configuration (N = 31, κ = 0.8). Curves show the potential barriers for three transitions (27, 4) → (26, 5) (solid line), (26, 5) → (25, 6) (dashed) and (27, 4) → (28, 3). The peak in the solid curve around 0.75 < r0 < 0.9 occurs when the number of MC steps is reduced simulating incomplete relaxation. Interestingly this peak matches the return point in the experiment, cf. Fig. 9 (a). 6.2 Probability of occurance of metastable states In the case of coexisting stationary states the question arises with what probability each of them will be observed when the Yukawa balls are formed. This subtle question could in fact be answered experimentally by performing a sequence of repeated measurements where the confinement was rapidly turned off and on again such that the plasma parameters were not altered [112, 113]. The statistics was sufficient to extract reliable occurrence frequencies of different states which are shown in the left part of Fig. 8 by the horizontal solid straight lines (the dashed lines indicate the statistical error). The measurements proved the existence of three stationary states. Most suprisingly, the metastable configuration (26, 5) was found to occure with an almost two times higher probability than the ground state (27, 4) [(62 ± 13)% compared to (35 ± 10)%]. We will now analyze these observations theoretically. The occurence probability of stationary states can be described analytically by a harmonic approximation of the partition function [90]. There, the potential energy of a given state s (ground or metastable state) is expanded around the local minimum Es0 . With the assumption that the vibrational and the three rotational modes of the whole system are decoupled, the full energy of the state s is, to second order in the displacements, ξs,i = ri − ris , Es = Es0 + f X i=1 ( p2ξs,i m 2 2 + ωs,i ξs,i 2m 2 ) + 3 X L2s,i . 2Is,i i=3 (12) Copyright line will be provided by the publisher 14 H. Baumgartner, D. Block, and M. Bonitz: Structure of Yukawa balls Fig. 8 (color online) Experimental and theoretical probability of stationary states for a Yukawa ball with N = 31. Left: Experimental results including mean probability (full lines) and standard deviations (dashes) corresponding to the states (from top to bottom) with 5, 4 and 6 particles on the inner shell compared to LMD results (lines with symbols) versus friction at κ = 0.8. Right: Analytical theory (lines) compared to MC results (symbols) as a function of temperature. The solid lines are obtained from Eq. (15) and the dashed lines neglect the statistical weight factor caused by the eigenfrequencies, (Ωs = 1) for all states. Agreement with the experimental data is obtained for T ∼ 0.025. The term in parentheses denotes the vibrational kinetic and potential energy of the f = 3N − 3 normal modes with the normal coordinates ξi . Even though the harmonic approximation is only applicable for low temperatures (or strong coupling) when the particles oscillate around the equilibrium positions with a small amplitude, it can be used to estimate phase transitions and compare the results with simulations and experimental data. With the factorized partition function Zs = ns Zsint Zsvib Zsrot , N! , with the degeneracy factor ns = QL s i=1 Ni ! (13) PL where L is the number of shells and Nis the occupation number of shell i with i=1 Ni = N . The degeneracy factor ns denotes the number of possibilities to form a configuration with shell occupation (N1 , N2 , ..., NL ) from 0 N distinguishable particles. The internal partition function is given by Zsint = e−βEs , and the vibrational and rotational partition functions are Zsvib (T ) =  kB T h̄Ωs f , Zsrot (T ) =  2πkB T Is h̄2 3/2 , (14) Qf where we introduced the geometric mean eigenfrequency Ωs = ( i=1 ωs,i )1/f and the mean moment of inertia Is = (Is,1 Is,2 Is,3 )1/3 . Then the occurence probability of the stationary state s is obtained as Ps = ns Zs ns Zs . = PM Z ′ ′ s′ =1 ns Zs (15) The moments of inertia cancel to a good approximation for clusters with 27-40 particles, while for lower particle numbers they should be included. Since the total partition function Z is not known, it is advantageous to compute instead the ratio of probabilites of two states s and s′ , using Eqs. (14) and (15) ns Ps = Ps ′ n s′  Ω s′ Ωs f  Is Is′ 3/2 0 0 e−β(Es −Es′ ) ≈ ns n s′  Ωs′ Ωs f 0 0 e−β(Es −Es′ ) . (16) Thus the probability ratio of two states depends on three factors: their energy difference Es0 − Es0′ , the ratio of 0 0 degeneracy factors ns /ns′ and the ratio of mean eigenfrequencies Ωs′ /Ωs . The Boltzmann factor e−β(Es −Es′ ) gives preference to states with a low energy. For low temperatures it will be the dominant factor but it becomes Copyright line will be provided by the publisher cpp header will be provided by the publisher 15 0 0 less important for higher temperatures when kB T >> Es0′ − Es0 and e−β(Es −Es′ ) ≈ 1. The degeneracy factor assigns large statistical weight to states with more particles on inner shells. For example, for N = 31 the ratio of the factors for two configurations with, respectively, 4 and 5 particles on the inner shell is n(27,4) /n(26,5) = (26!5!)/(27!4!) = 5/27. The energy differences are computed by MD simulations and, as an exmmple for N = 31 and κ = 1.4, are given in table 1. The effect of the mean eigenfrequency, i.e. the effect of the local curvature of the potential surface cannot be that easily predicted but the ratio of the eigenfrequencies are in the range of 0.8 − 30 and can therefore strongly influence the ratio of the probabilty of states. The analytical results can be compared with MC simulation results [90] which are practically exact. In particular, they include all anharmonic corrections. For a fixed temperature the occurrence probability of each stationary state is recorded. Both results are shown in Fig. 8 as a function of temperature, for the example N = 31. Evidently, the analytical results are accurate at temperatures below kB T < 0.02. Above this value anharmonic corrections become important. For other particle numbers and screening parameters simliar results are obtained. Let us now compare to the experimental results, cf. left part of Fig. 8. These data were obtained at room temperature, corresponding to approximately T = 0.0015. At this temperature the MC and analytical results both are in agreement with each other (right figure) predict practically 100 percent occurence of the ground state (27, 4) which is in striking contrast to the experiment. To resolve this contradiction, independent Langevin dynamics simulations have been performed [90, 113] which model the experimental situation of turning on and off the confinement by including a time-dependent potential U ext in Eq. (1) and which are included in the left figure part. Using a screening parameter which corresponds to the experimental parameters (the same as in the MC simulations) one finds a strong dependence on the chosen values of the friction coefficient ν. For small ν the simulation always favor the ground state which is easily understood: with small friction the system is cooled very slowly and has enough time to find the lowest potential minimum (this is nothing but simulated annealing which was applied in Sec. 4.1 to obtain the ground states). In contrast, with large damping the particles rapidly loose their kinetic energy even before settling in a potential minimum. The experimental damping parameters are well above ν = 2 where Fig. 8 indicates very good agreement between LMD and experiment. On the other hand, Monte Carlo and analytical theory are based on a hamiltonian model, cf. Sec. 3.2, and do not include damping effects. Thus, a direct comparison with the experimental probabilities is not possible. However, as was demonstrated in Ref. [90] due to the large friction in the experiments the plasma particles rapidly relax towards a Maxwellian distribution already at temperatures of about T = 0.05. During further cooling the system remains in equilibrium, only its temperature is being gradually reduced until it reaches T ∼ 0.002. The “decision” about what stationary state will be reached is “made” at an early stage, before the Yukawa ball freezes which occurs around T ∼ 0.02. In fact, Fig. 8 shows that the MC simulations at T ≈ 0.025 agree very well with the experiment and the LMD [90]. 6.3 Structural transitions in the experiment Transitions between two stationary configurations which were discussed above could, in fact, be observed experimentally [113]. Figure 9(a) shows the motion of a single particle of a Yukawa ball with 31 particles over a period of ten minutes. It can be seen that the particle moves from the outer shell inwards and returns to the outer shell within about 200 seconds. This transition is temperature induced (it occurs spontaneously) and is observed in other Yukawa balls as well [113]. This behavior is well explained by the MC simulation of the total potential profile [see the discussion in Sec. 6.1] which are shown in the right part of Fig. 7. The potential barrier for moving a particle inwards from (27, 4) to (26, 5) [lowest curve] shows where the particle returns; interestingly the sharp increase at r0 ≈ 0.85 is matching the return point of r = 0.39mm in the experiment, cf. Fig. 9. This peak corresponds to a MC run with fewer steps where the other particles have not enough time to fully relax in response to the moving paricle – this nonequilibrium situation is apparantly observed in the experiment. 6.4 Finite temperature effects. Solid-liquid phase transition So far we have considered only the low-temperature behavior of the Yukawa balls where the system was in the ground state or low-lying metastable state. Finite thermal energy caused only occasional transitions between stationary states. This behavior is typical for matter at very strong coupling characterized by very low kinetic Copyright line will be provided by the publisher 16 H. Baumgartner, D. Block, and M. Bonitz: Structure of Yukawa balls Fig. 9 (a) Radial component of the trajectory of one particle leaving the outer shell of a Yukawa ball with N = 31, the inset shows the trajectory in the ρ-z plane (black symbols). (b) Average Wigner-Seitz cell volume per particle, VW S , for all particles belonging to the inner shell as a function of the radial position of the traveling particle. For R = 0.55mm, the Wigner-Seitz cell of the traveling particle is included in the average. The solid lines are the average VW S of the configurations (4, 27) and (5, 26) obtained by repeating the experiment and averaging. energy, typical for a Wigner crystal. This physical situation can be quantified by introducing a coupling parameter, which is defined as the ratio of the mean interaction energy to kinetic energy. For a classical macroscopic Coulomb system in thermal equilibrium the definition is ΓC = q2 , 4πǫakB T (17) where a is the mean interparticle distance. The transition from strong coupling where the system is in a crystalline state to a liquid occurs at critical value Γcr which is of the order 175 (137) in 3D (2D). The transition to a gas-like state occurs around ΓC = 1, (for generalizations to quantum systems, see Ref. [117]). Here we are concerned with systems with Yukawa interaction and have to generalize the coupling parameter. The straighforwad way to do this would be to replace the Coulomb energy of a pair of particles by the mean Yukawa energy and ΓC → e−κa · ΓC , but this makes it impossible to systematically define a value for the melting transition in a Yukawa system, so for the phase diagram [118] one has to use other quantities. On the other hand, it is not the total potential energy which governs the melting transition, but the process is driven by the excitation energy of particles above the potential minimum which leads to increasing fluctuations of the particles around their stationary positions. These fluctuations are defined by the local curvature of the potential [99, 117], i.e. its second derivative in the minimum. This allows to define a proper effective Yukawa coupling parameter [119], Γ(κ) = ΓC eκa , 1 + κa + (κa)2 /2 (18) which allows to qualitatively correctly predict the melting point Γcr (κ) for any κ, based on the known critical value for a Coulomb system, i.e. by using, in Eq. (18), ΓC = Γcr . However, for the analysis of melting in Yukawa balls one has to take into account the peculiarity of small systems. Finite systems behave fundamentally different to macrososcopic systems: there exist no strict phase Copyright line will be provided by the publisher cpp header will be provided by the publisher 17 transitions. Nevertheless, a similar transition or crossover from solid-like to liquid-like behavior is well known and has been studied as a function of particle number e.g. by S CHIFFER [95]. For small N the melting temperatures are typically higher than in the macroscopic limit and the crossover proceeds over a finite range of coupling parameters (temperatures) making it often difficult to identify a melting point, see below. Furthermore, finite systems are strongly influenced by the symmetry of the confinement potential, and the shell structure of the Coulomb balls affects the melting behavior. The effect of screening on the stability and melting of Yukawa balls was analyzed in Ref. [99], where a simple analytical estimate can be observed for the example N = 2. Expanding the total potential energy to second order around the ground state, as in the macroscopic system, a formula analogous to (18) is obtained Γ(κ, N = 2) = ΓC (N = 2) eκr0Y 2 /3 . 1 + κr0Y + κ2 r0Y (19) Here r0Y (κ) is the two-particle distance for Yukawa interaction which is obtained by inverting the equation 3 2ω 2 q 2 eκr0Y r0Y = = r03 (r0Y ), 1 + κr0Y 4πǫm (20) where r0 is the equilibrium distance of two particles with Coulomb interaction, see Sec. 3. As in the macroscopic case, screening always destabilizes the crystal. Another quantity to locate the melting point are the mean position fluctuations of particles which are known to be small (large) in the solid (liquid) phase, independently of the interaction and system size [117] which is commonly known as Lindemann criterion. However, this parameter diverges in 2D systems, and an improved quantity are the mean relative interparticle distance fluctuations, i.e. the fluctuations normalized to the mean distance s N 2i X hrij 2 − 1. (21) urel = N (N − 1) hrij i2 1≤i≤j This quantity has been successully applied to many finite systems and also to Yukawa balls [89, 98, 100, 107], an example is shown in the left part of Fig. 10. One clearly sees the increase of the fluctuations typical for a melting process. It is also evident that this increase extends over more than one order of magnitude in the temperature. A rather detailed analysis of melting in Coulomb balls was performed by A POLINARIO [100] who observed that the melting often proceeds in several steps. Finally, in a recent study [120] it was observed that the quantity urel has a number of difficulties in small systems – the result strongly depends on the way of calculation. Instead it was suggested to use another quantity – the variance of the block averaged interparticle distance fluctuations (VIDF) σurel = K q 1 X hu2rel (s)i − hurel (s)i2 . K s=1 (22) This parameter measures the width of the probability distribution P (urel ) which has a maximum in the vicinity of the phase transition. In the solid phase the particles show only small distance fluctuations, and in the liquid phase all particles strongly fluctuate and P (urel ) has a peak centered around a large value, yet the width of the peak is again small. However, in the vicinity of the phase transition, the particles may be in one state for some time and frequently make transition over the potential barrier to another stable state and so on. This causes a broad peak of P , as can be seen in the rightmost part of Fig. 10. An experimental analysis of the melting behavior is difficult because it require a systematic scan of the plasma parameters in a sufficiently broad range around the expected melting point. To do this without changing the screening is difficult and requires future advances in the handling of the Yukawa balls. 6.5 Further results on Yukawa balls In this review we had to concentrate on a few properties of Yukawa balls leaving out many other interesting results some of which are mentioned in the following. Copyright line will be provided by the publisher 18 H. Baumgartner, D. Block, and M. Bonitz: Structure of Yukawa balls Fig. 10 Mean distance fluctuations (21), crosses and their variance (22), circles, for a Yukawa ball with κ = 0.4 in the vicinity of the melting point Tcr . Left: Temperature dependence of urel and σurel for a Yukawa ball with N = 40. Right: Block averaged interparticle distance fluctuations urel (s) vs block number s during an MC simulation and their accumulated probability P (rightmost figure) for N = 40 and four temperatures (from bottom to top): below, close to, at and above Tcr . Block length M = 1000. The normal mode spectrum of Yukawa balls which was mentioned in Sec. 6.2 has been analyzed in great detail. Among the most interesting collective excitations of these finite systems is the breathing mode (uniform radial expansion and contraction). Recently it could be shown [121] that Yukawa balls do not possess such a mode. Instead they can exhibit simultaneously several monopole oscillations (periodic increase and reduction of the mean radius). Also, first experimental observations of the normal modes of small Yukawa balls have been performed [122]. The confinement potential in the experiments is certainly never perfectly isotropic, so the analysis of anisotropy effects is of interest. Ground and metastable states in anisotropic parabolic confinements have been analyzed as a function of the anisotropy by A POLINARIO et al. [98, 100]. They also analyzed the normal modes of these systems. The results are of interest also for experiments where the confinement is made anisotropic on purpose. Another interesting theoretical investigation has been performed by P SAKHIE et al. They performed MD simulations for Yukawa balls containing two sorts of particles having different masses [97]. In their model the confinement forces, the termophoretic force, gravity, electric field force and ion drag have been explicitely included. Increasing the mass difference they found that beyond a critical mass asymmetry the Yukawa ball may split in two balls containing different particles which are vertically separated. Recently, first simulations of the short-time dynamics of Yukawa balls have been performed by K ÄHLERT et al. [90]. Starting from a random initial state and rapidly turning on the confinement showed a very interesting nonequilibrium behavior with a non-monotonic kinetic energy increase. This behavior is very different from the ultrafast dynamics of weakly coupled plasmas and is a very interesting subject of further research. 7 Discussion and Outlook In this review, and overview on the recently discovered spherical crystals in dusty plasmas – Yukawa balls – has been given. The experimental setup and the diagnostic methods have been briefly discussed. The main attention was then devoted to the theoretical modeling and to computer simulatinos of Yukawa balls, when possible, aiming at a detailed comparison to the experimental data. This close connection between experiment and theory has lead to substantial advances in our understanding of finite Yukawa crystals. It was shown that the Yukawa balls produced in the Kiel and Greifswald labs are well described by an isotropic statically screened pair potential and that wake effects are apparantly not of importance. This situation may change when larger balls are studied which is an interesting question for further research. Copyright line will be provided by the publisher cpp header will be provided by the publisher 19 It was shown that screening has a crucial effect on the internal structure of Yukawa balls, including the shell populations, the ground and metastable states. Comparisons between measurement and simulations allowed to accurately determine the value of the screening parameter. Besides the shell structure, also the mean distribution of charge throughout the Yukawa balls could be explained. This was achieved by an analytical theory based on a fluid-like description. It could be shown that, with increased screening, the average density is more and more concentrated in the central region of the cluster and rapidly decays towards the edge. This is an important difference to Coulomb systems, including spherical ion crystals, which exhibit a constant radial density profile. Not only the ground state of the Yukawa balls could be observed. The experiments were able to reliably detect metastable states and to reveal the probability of their occurence. These measurements could be well explained by theory and simulations where it became clear that the (in some cases) dominant role of metastable states – despite relatively low temperature – is due to the large friction experienced by the dust particles in the plasma. Besides the static properties also several dynamic aspects are now well understood. The experiments could demonstrate that spontaneous transitions between various stationary states occur which is well explained by a theoretical analysis of the potential energy landscape, in particular from the energy barriers between local potential minima. The theory further predicted a rather interesting melting behavior of the Yukawa balls upon further temperature increase. It will be of high interest to verify these predictions in future experiments. At the same time many fundamental theoretical questions remain to be answered. In the shell strucutre of the Yukawa balls correlation effects play a crucial role. This became clear by studying greatly simplified shell models. These effect were also important to explain the density profile of Yukawa balls at large screening using a local density approximation. Yet no analytical results for the pair correlations are known and no derivation of these effects from statistical mechanics has been given. The excellent experimental accuracy gives a unique opportunity to advance and verify modern theoretical approaches in strongly coupled plasmas. Finally, questions of further interest in the field of Yukawa balls include their normal mode spectrum, and their interaction with external electric and magnetic fields. This will require to take into account the modification of the pair interaction. Here also wake effects could play a more important role. Last but not least it is becoming possible to advance to an analysis of the short-time behavior Yukawa balls. Excerting a fast excitation, e.g. by laser manipulation, the time-dependent response of these strongly correlated systems can be systematically studied. Acknowledgements We acknowledge stimulating discussions with J.W. Dufty, Ch. Henning, H. Kählert, P. Ludwig, A. Melzer and A. Piel. This work is supported by the Deutsche Forschungsgemeinschaft via SFB-TR 24 Projects A3, A5, and A7. M.B. acknowledges hospitality of the Physics Department of the University of Florida. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] C K Goertz, 1989, Rev. Geophys. 27, 271–292 T G Northrop, 1992, Physica Scripta 45, 475–490 G S Selwyn, J Singh, and R S Bennett, 1989, J. Vac. Sci. Technol. A 7, 275827–65 G S Selwyn, J S McKillop, K Haller, and J Wu, 1990, J. Vac. Sci. Technol. A 8, 1726–1731 A Bouchoule, A Plain, L Boufendi, J P Blondeau, and C Laure, 1991, J. Appl. Phys. 70, 1991–2000 P Belenguer, J P Blondeau, L Boufendi, M Toogood, A Plain, A Bouchoule, C Laure, and J P Boeuf, 1992, Phys. Rev. A 46, 7923–7933 C Hollenstein, 2000, Plasma Phys. Control. Fusion 42, R93–R104 J Perrin, J Schmitt, C Hollenstein, and A Howling, 2000, Plasma Phys. Control. Fusion 42, B353–B363 P Bliokh, V Sinitsin, and V Yaroshenko, 1995, Dusty and Self-Gravitational Plasma in Space (Dordrecht: Kluwer) A Bouchoule, 1999, Dusty Plasmas: Physics, Chemistry, and Technological Impacts in Plasma Processing (New York: Wiley) F Verheest, 2000, Waves in Dusty Space Plasmas (Dordrecht: Kluwer) G E Morfill, B M Annaratone, P Bryant, A V Ivlev, H M Thomas, M Zuzic and V E Fortov, 2002, Plasma Phys. Control. Fusion 44, B263-B277 E P Wigner, 1934, Phys. Rev. 46, 1002 C C Grimes, and G Adams, 1979, Phys. Rev. Lett. 42, 795 Copyright line will be provided by the publisher 20 H. Baumgartner, D. Block, and M. Bonitz: Structure of Yukawa balls [15] For recent experimental results on electron crystals on superfluid helium, see the paper E Rousseau, D Ponarin, L Hristakos, O Avenel, E Varoquaux, and Y Mukharsky, 2009, Phys. Rev. B 79, 045406 and the discussion by F Peeters, 2009, Physics 2, 4 [16] E Y Andrei, G Deville, D C Glattli, F I B Williams, E Paris, and B Etienne, 1988, Phys. Rev. Lett. 60, 2765 [17] A Filinov, M Bonitz and Yu E Lozovik, 2001, Phys. Rev. Lett. 86, 3851 and 2000, phys. stat. sol. (b) 221, 231 [18] J H Chu and Lin I, 1994, Phys. Rev. Lett. 72, 4009 [19] M Bonitz, V S Filinov, V E Fortov, P R Levashov, and H Fehske, 2005, Phys. Rev. Lett. 95, 235006 and 2006, J. Phys. A: Math. Gen. 39, 4717 [20] D J Wineland et al, 1987, Phys. Rev. Lett. 59, 2935 [21] M Drewsen et al, 1998, Phys. Rev. Lett. 81, 2878 [22] H Ikezi, 1986, Phys. Fluids 29, 1764–1766 [23] J H Chu, J B Du, and I Li, 1994, J. Phys. D: Appl. Phys. 27, 296–300 [24] H Thomas, G E Morfill, V Demmel, J Goree, B Feuerbacher, and D Möhlmann, 1994, Phys. Rev. Lett. 73, 652 [25] Y Hayashi, and K Tachibana, 1994, Jpn. J. Appl. Phys. 33, L804–L806 [26] A Piel and A Melzer, 2002, Plasma Phys. Control. Fusion 44 [27] O Arp, D Block, A Piel and A Melzer, 2004, Phys. Rev. Lett. 93, 165004 [28] T Pohl, T Pattard, and J M Rost, 2004, Phys. Rev. Lett. 92, 155003 [29] T Killian, 2004, Nature 429, 815 [30] D Block and A Melzer, Imaging diagnostics in dusty plasmas, in: “Introduction to Complex Plasmas” edt. by M Bonitz, N Horing, J Meichsner and. P Ludwig, Springer, Berlin, Heidelberg 2009. [31] M A Liebermann and A J Lichtenberg, Principles of plasma discharges and material processing. John Wiley & Sons Inc., New York, 1994. [32] J B Pieper and J Goree, 1996, Phys. Rev. Lett. 77, 3137 [33] T Trottenberg, A Melzer, and A Piel, 1995, Plasma Sources Sci. Technol. 4, 450 [34] A Melzer, T Trottenberg, and A Piel, 1994, Phys. Lett. A 191, 301 [35] F Melandsø and J Goree, 1995, Phys. Rev. E 52, 5312 [36] F Melandsø and J Goree, 1996, J. Vac. Sci. Technol. A 14, 511 [37] A Melzer, V A Schweigert, I V Schweigert, A Homann, S Peters, and A Piel, 1996, Phys. Rev. E 54, R46 [38] M Nambu, S V Vladimirov, and P K Shukla, 1995, Phys. Lett. A 203, 40 [39] V A Schweigert, I V Schweigert, A Melzer, A Homann, and A Piel, 1996, Phys. Rev. E, 54, 4155 [40] Y Hayashi, 1999, Phys. Rev. Lett. 83, 4764 [41] J B Pieper, J Goree, and R A Quinn, 1996, Phys. Rev. E 54, 5636 [42] M Zuzic, A V Ivlev, J Goree, G E Morfill, H M Thomas, H Rothermel, U Konopka, R Sütterlin, and D D Goldbeck, 2000, Phys. Rev. Lett. 85, 4064 [43] A Barkan and R L Merlino, 1995, Phys. Plasmas 2, 3261. [44] C Thompson, A Barkan, N D’Angelo, and R L Merlino, 1997, Phys. Plasmas 4, 2331 [45] V E Fortov, A G Khrapak, S A Khrapak, V I Molotkov, A P Nefedov, O F Petrov, and V M Torchinsky, 2000, Phys. Plasmas 7, 1374 [46] S V Ratynskaia, S A Khrapak, A V Zobnin, M H Thoma, M Kretschmer, A D Usachev, V Yaroshenko, R A Quinn, G E Morfill, O Petrov, and V E Fortov, 2004, Phys. Rev. Lett. 93, 085001. [47] S V Ratynskaia, M Kretschmer, S A Khrapak, R A Quinn, G E Morfill, M H Thoma, A V Zobnin, A D Usachev, O Petrov, and V E Fortov, 2004, IEEE Trans. Plasma Sci. 32, 613. [48] E Thomas Jr. and R L Merlino, 2001, IEEE Trans. Plasma Sci. 29, 152. [49] T. Trottenberg, D. Block, and A. Piel 2006, Phys. Plasmas 13, 042105. [50] I. Pilch, A. Piel, T. Trottenberg, and M.E. Koepke 2007, Phys. Plasmas 14, 123704. [51] I. Pilch, T. Reichstein, and A. Piel 2008, Phys. Plasmas 15, 103706. [52] A P Nefedov, G E Morfill, V E Fortov, H M Thomas, H Rothermel, T Hagl, A V Ivlev, M Zuzic, B A Klumov, A M Lipaev, V I Molotkov, O F Petrov, Y P Gidzenko, S K Krikalev, W Shepherd, A I Ivanov, M Roth, H Binnenbruck, J A Goree, and Y P Semenov, 2003, New J. Phys. 5, 33. [53] B M Annaratone, M Glier, T Stuffler, M Raif, H M Thomas, and G E Morfill, 2003, New J. Phys. 5, 92. [54] P Bryant, 2004, New J. Phys. 6, 1 [55] A V Ivlev, H M Thomas, G E Morfill, V Molotkov, A M Lipaev, V E Fortov, T Hagl, H Rothermel, and S Krialev, 2006, New J. Phys. 8, 25. [56] V E Fortov, O S Vaulina, O F Petrov, V I Molotkov, A M Lipaev, G E Morfill, H Thomas, S A Khrapak, Y P Semenov, and A I Ivanov, 2004, Plasma Phys. Control. Fusion 46, 359 [57] V E Fortov, O S Vaulina, O F Petrov, V I Molotkov, A M Lipaev, V M Torchinsky, H M Thomas, G E Morfill, S A Khrapak, Y P Semenov, A I Ivanov, S K Krikalev, A Y Kalery, S V Zaletin, and Y P Gidzenko, 2003, Phys. Rev. Lett. 90, 245005. [58] S A Khrapak, D Samsonov, G E Morfill, H Thomas, V Yaroshenko, H Rothermel, V Fortov, A Nefedov, V Molotkov, O Petrov, A Lipaev, A Ivanov, and Yu M Baturin, 2003, Phys. Plasmas 10 1 [59] A Piel, M Klindworth, O Arp, A Melzer, and M Wolter, 2006, Phys. Rev. Lett. 97, 205009. Copyright line will be provided by the publisher cpp header will be provided by the publisher 21 [60] D Samsonov, G E Morfill, H Thomas, T Hagl, and H Rothermel, 2003, Phys. Rev. E 67, 036404. [61] V V Yaroshenko, B M Annaratone, S A Khrapak, H M Thomas, G E Morfill, V E Fortov, A M Lipaev, V I Molotkov, O F Petrov, A I Ivanov, and M V Turin, 2004, Phys. Rev. E 69, 066401. [62] G E Morfill, H M Thomas, U Konopka, H Rothermel, M Zuzic, A Ivlev, and J Goree, 1999, Phys. Rev. Lett. 83, 1598 [63] J Goree, G E Morfill, V N Tsytovich, and S V Vladimirov, 1999, Phys. Rev. E 59, 7055 [64] V Tsytovich, 2001, Physica Scripta 89, 89 [65] V N Tsytovich, S V Vladimirov, G E Morfill, and J Goree, 2001, Phys. Rev. E 64, 029902 [66] V N Tsytovich, S V Vladimirov, G E Morfill, and J Goree, 2001, Phys. Rev. E, 63, 056609. [67] M R Akdim and W J Goedheer, 2002, Phys. Rev. E 65, 015401. [68] V Land and W J Goedheer, 2006, New J. Phys. 8, 8. [69] V Land, W J Goedheer, and M R Akdim, 2005, Phys. Rev. E 72, 046403. [70] M Klindworth, A Piel, and A Melzer, 2004, Phys. Rev. Lett. 93, 195002. [71] G Praburam and J Goree, 1996, Astrophys. J. 441, 830 [72] M Wolter, A Melzer, O Arp, M Klindworth, and A Piel, 2007, Phys. Plasmas. [73] H Rothermel, T Hagl, G E Morfill, M H Thoma, and H M Thomas, 2002, Phys. Rev. Lett. 89, 175001 [74] I H Hutchinson, 2005, Plasma Phys. Control. Fusion 47, 71 [75] I H Hutchinson, 2006, Plasma Phys. Control. Fusion 48, 185. [76] S A Khrapak, A V Ivlev, G E Morfill, and H M Thomas, 2002, Phys. Rev. E 66, 046414. [77] A M. Lipaev, S A. Khrapak, V I Molotkov, G E Morfill, V E Fortov, A V Ivlev, H M Thomas, A G Khrapak, V N Naumkin, A I Ivanov, S E Tretschev, and G I Padalka, 2007, Phys. Rev. Lett. 98, 265006. [78] B M Annaratone, T Antonova, D D Goldbeck, H. M Thomas, and G E Morfill, 2004, Plasma Phys. Control. Fusion 46, 495 [79] T Antonova, B M Annaratone, D D Goldbeck, V Yaroshenko, H M Thomas, and G E Morfill, 2006, Phys. Rev. Lett., 96, 115001. [80] M R Akdim and W J Goedheer, 2004, IEEE Trans. Plasma Sci. 32, 680 [81] J P Boeuf and L C Pitchford, 1995, Phys. Rev. E 51, 1376 [82] J P Boeuf, Ph Belenguer, and T Hbid, 1994, Plasma Sources Sci. Technol. 3 407 [83] O Arp, D Block, M Klindworth, and A Piel, 2005, Phys. Plasmas 12, 122102. [84] H Totsuji, T Ogawa, C Totsuji, and K Tsuruta, 2005, Phys. Rev. E 72, 036406. [85] H Totsuji, C Totsuji, T Ogawa, and K Tsuruta, 2005, Phys. Rev. E 71, 045401 [86] W Paul, 1990, Rev. Mod. Phys. 62, 531 [87] K. Matyash et al., Poster P5.043 at EPS Conference Warsaw 2007, http : //epsppd.epf l.ch/W arsaw/pdf /P 5043.pdf [88] M Lampe, G Joyce, G Ganguli and V Gavrishchaka, 2000, Phys. Plasmas 7, 10 [89] V Golubnychiy, H Baumgartner, M Bonitz, A Filinov and H Fehske, 2006, J. Phys. A: Math. Gen. 39, 4527 [90] H Kählert, P Ludwig, H Baumgartner, M Bonitz, D Block, S Käding, A Melzer and A Piel, Phys. Rev. E 78, 036408 (2008) [91] Introduction to Computational Methods for Many-body physics, M. Bonitz and D. Semkat (eds.), Rinton Press, Princeton 2006 [92] P Ludwig, S Kosse and M Bonitz, 2005, Phys. Rev. E 71, 046403 [93] S Kirkpatrick and C D Gelatt, 1983, Science, Vol. 220, pp 671 [94] R W Hasse and V V Avilov, 1991, Phys. Rev. A 44, 4506 [95] J P Schiffer, 2002, Phys. Rev. Lett. 88, 205003 [96] O Arp, D Block, M Bonitz, H Fehske, V Golubnychiy, S Kosse, P Ludwig, A Melzer, and A Piel, 2005, J. Phys.: Conf. Ser. 11, 234 [97] S G Psakhie, K P Zolnikov, L F Skorentsev, D S Kryzhevich, and A V Abdrashitov, 2008, Phys. Plasmas 15, 053701 [98] S W S Apolinario, B Partoens and F M Peeters, 2007, New J. Phys. 9, 283 [99] M Bonitz, D Block, O Arp, V Golubnychiy, H Baumgartner, P Ludwig, A Piel and A Filinov, 2006, Phys. Rev. Lett. 96, 075001 [100] S W S Apolinario, PhD thesis, Universtiy Antwerp 2008 [101] H Baumgartner, D Asmus, V Golubnychiy, P Ludwig, H Kählert, and M Bonitz, 2008, New J. Phys. 10, 093019 [102] F Aurenhammer, 1991, ACM Computing Surveys 23, 3 [103] K Tsuruta and S Ichimaru, 1993, Phys. Rev. A 48, 1339 [104] T Kamimuraa, Y Suga and O Ishihara, 2007, Phys. Plasmas 14, 123706 [105] W D Kraeft and M Bonitz, 2006, J. Phys.: Conf. Ser. 35, 94 [106] C Henning, H Baumgartner, A Piel, P Ludwig, V Golubnychiy, M Bonitz and D Block, 2006, Phys. Rev. E 74, 056403 [107] H Baumgartner, H Kählert, V Golubnychiy, C Henning, S Käding, A Melzer and M Bonitz, 2007, Contrib. Plasma Phys. 47, 281-290 [108] J Cioslowski, 2008, Phys. Rev. E 78, 026416 [109] J Cioslowski, 2008, J. Chem. Phys. 128, 164713 [110] W D Kraeft and M Bonitz, 2006, J. Phys.: Conf. Ser. 35, 78 Copyright line will be provided by the publisher 22 H. Baumgartner, D. Block, and M. Bonitz: Structure of Yukawa balls C Henning, P Ludwig, A Filinov, A Piel and M Bonitz, 2007, Phys. Rev. E 76, 036404 S Käding, D Block, A Melzer, A Piel, H Kählert, P Ludwig, and M Bonitz, 2008, Phys. Plasmas 15, 073710 D Block, S Käding, A Melzer, A Piel, H Baumgartner and M Bonitz, 2008, Phys. Plasmas 15, 040701 A Piel, O Arp, D Block, I Pilch, T Trottenberg, S Käding, A Melzer, H Baumgartner, C Henning and M Bonitz, 2008, Plasma Phys. Control. Fusion 50, 124003 [115] H Totsuji, C Totsuji and K Tsuruta, 2001, Phys. Rev. E 64, 066402 [116] J Wrighton, J W Dufty, C Henning, and M Bonitz, 2009, accepted for publication in J. Phys. A, arXiv:0809.3071 [117] M Bonitz, P Ludwig, H Baumgartner, C Henning, A Filinov, D Block, O Arp, A Piel, S Käding, Y Ivanov, A Melzer, H Fehske and V Filinov, 2008, Phys. Plasmas 15, 055704 [118] S Hamaguchi, R T Farouki and D H E Dubin, 1997, Phys. Rev. E 56, 4671 [119] O S Vaulina and S A Khrapak, 2000, JETP 90, 287 [120] J Böning, A Filinov, P Ludwig, H Baumgartner, M Bonitz and Yu E Lozovik, 2008, Phys. Rev. Lett. 100, 113401 [121] C Henning, K Fujioka, P Ludwig, A Piel, A Melzer, and M Bonitz, 2008, Phys. Rev. Lett. 101, 045002 [122] C Henning, H Kählert, P Ludwig, A Melzer, and M Bonitz, 2009, accepted for publication in J. Phys. A, arXiv:0810.2592 [111] [112] [113] [114] Copyright line will be provided by the publisher