Theoretical Description of Coulomb Balls - Fluid Phase J. Wrighton,1 J. W. Dufty,1 H. Kählert,2 and M. Bonitz2 arXiv:0909.0775v1 [cond-mat.stat-mech] 3 Sep 2009 1 Department of Physics, University of Florida, Gainesville, FL 32611 2 Institut für Theoretische Physik und Astrophysik, Christian-Albrechts Universität zu Kiel, 24098 Kiel, Germany (Dated: September 3, 2009) Abstract A theoretical description for the radial density profile of a finite number of identical charged particles confined in a harmonic trap is developed for application over a wide range of Coulomb coupling (or, equivalently, temperatures) and particle numbers. A simple mean field approximation neglecting correlations yields a density profile which is monotonically decreasing with radius for all temperatures, in contrast to molecular dynamics simulations and experiments showing shell structure at lower temperatures. A more complete theoretical description including charge correlations is developed here by an extension of the hypernetted chain approximation, developed for bulk fluids, to the confined charges. The results reproduce all of the qualitative features observed in molecular dynamics simulations and experiments. These predictions are then tested quantitatively by comparison with new benchmark Monte Carlo simulations. Quantitative accuracy of the theory is obtained for the selected conditions by correcting the hypernetted chain approximation with a representation for the associated bridge functions. PACS numbers: 52.27.Lw,52.27.Gr,52.65.Yy 1 I. INTRODUCTION Strong correlation effects in ensembles of charged particles are important in many fields of physics, including plasmas, the electron gas in solids, and electron-hole plasmas, e.g. [1, 2] and references therein. In recent years charged particles spatially confined in trapping potentials have attracted considerable interest. Examples are ultracold ions [3, 4], dusty plasmas [5, 6], and electrons in quantum dots [7]. Related studies of expanding neutral plasmas have been reported as well [8, 9]. Recently, the ground state structure of Coulomb crystals formed by classical charges in a harmonic trap has been studied in detail via molecular dynamics (MD) simulation for Coulomb interaction in [10, 11] and for Yukawa interaction in [6, 12]. The same system is studied here for its fluid phase at finite temperatures. For most experimental situations, e.g. in ultracold ions and dusty plasmas, finite temperature effects are expected. The objective here is to provide a theoretical description for the structure of this system over the full range of temperatures for the fluid phase. Among the primary properties of interest is the radial density profile of the particles in the trap. A theoretical description of the ground state has been obtained via an energy minimization for both Coulomb and Yukawa systems without correlations [13], and with correlations in a local density approximation [14]. Good agreement with simulations for the spatial average density profile was obtained in this way, but a more complete representation of correlations is required for the detailed spatial modulation (shells) expected at all except the highest temperatures. This is confirmed by the theoretical analysis developed here, following [15]. This theory is parameterized by the number of particles in the trap, N, and the plasma coupling constant, Γ, measuring the strength of the Coulomb energy for a pair of charges in the trap relative to the kinetic energy (defined below). Alternatively, the average kinetic energy relative to trap energy or to Coulomb energy can be used to define a dimensionless temperature T ∗ ≡ 1/Γ. The range of values explored is 1 < N ≤ 300 and 0 < Γ ≤ 40 (or 0.025 ≤ T ∗ < ∞). The primary results described here are: 1) an extension of the HNC approximation developed for bulk fluids to systems with localized densities, within the context of density functional theory; 2) a confirmation that mean field theory (no correlations) predicts a structureless, monotonically decreasing radial profile approaching a step function at T ∗ = 0; 3) a demonstration that HNC correlations lead to a shell structure (local radial peaks) for 2 strong coupling Γ ≥ 10 (T ∗ ≤ 0.1); 4) new Monte Carlo simulations to test this theory; and 5) a proposal for corrections to the HNC that leads to predictions in good quantitative agreement with the Monte Carlo simulations over the range of N, Γ tested. The qualitative features of the shell structure arising from correlations are: the number of shells increases with N, with new shells appearing at special values of N; a sharpening of the shells (more narrow, higher amplitude) with increasing Γ (decreasing T ∗ ); and shell populations that grow monotonically with N but which are almost independent of Γ (T ∗ ) (as observed in previous MD simulations [16, 17]). The HNC approximation provides agreement on all qualitative features of the density profile, e.g. the location of the shells and number of particles within the shells is well-described at all coupling values. However, discrepancies with the Monte Carlo simulations of the order of 20 − 40% for the widths and amplitudes of the shell structure are found. These discrepancies are removed by including a representation of the ”bridge functions” neglected in the HNC approximation, referred to in the following as the adjusted HNC (AHNC). The AHNC theory has the simplicity of the HNC, depending only on well-known correlations of the bulk one component plasma, and shows excellent agreement with simulations. This paper is organized as follows. First the Hamiltonian for the system is defined and the density profile is formulated in the Canonical ensemble. The relevant energy and length scales are introduced for a dimensionless representation in terms of N and Γ (or T ∗ ). This is the form most appropriate for Monte Carlo simulation. In Section III, the corresponding formulation is given in the Grand Canonical ensemble to allow exploitation of density functional theory. A formally exact representation for the density profile is obtained in Appendix A in terms of correlations in a bulk one component plasma (OCP), together with a corresponding exact equation for the OCP. The HNC approximation is defined in Section III as the neglect of the ”bridge functions” in these two sets of equations. A further neglect of all correlations, the mean field approximation, is explored first over the whole temperature range for reference in Section IV. The effects of correlations within the HNC approximation are illustrated in Section V. Selected results are compared with those from Monte Carlo simulations in Section VI, where the AHNC is defined and shown to give the necessary corrections to HNC required for quantitative agreement. Finally the results are summarized and discussed in the last Section. In closing this Introduction, it is worth emphasizing that the interesting radial structures 3 studied occur for strong plasma coupling, Γ & 10. Thus the system presents a good example of a strongly coupled Coulomb system. For very large Γ the shells become effectively systems of particles uniformly distributed on the surfaces of spheres of specific radii. At some point the rotational symmetry on these spherical surfaces is broken and these particles form spherical Wigner crystals [6]. An adequate description of correlations within density functional theory should provide both the fluid phase considered here and the freezing transition. The crystal phase is more subtle than in a bulk plasma since it entails packing on a non-Euclidean surface, with necessary disclinations depending on the value of N considered (related to Thomson’s problem [18]). The determination of the non-uniform density for charges confined in a trap is similar to the complementary problem of determining the enhanced density in an OCP with an impurity ion of the opposite charge. In the latter case the density profile of the bound states corresponds to that for the charges in the trap. The formalism here extends to this case in a natural way [15]. Both the structure and dynamics for this problem have been studied in some detail recently [19], and will not be considered further here. II. MONTE CARLO SIMULATION A system of N identical charges in a spherical container of radius R is considered. An attractive central force is applied at the origin. The Hamiltonian is H= N  X 1 i=1 2 mvi2 N 1 X V (rij ) + VW . + V0 (ri ) + 2 i6=j=1  (1) Here ri and vi are the position and velocity of charge i. The repulsive interaction potential between charges i and j is denoted by V (rij ) where rij ≡ |ri − rj |. The single particle potential V0 (ri ) represents the central ”confinement” potential (referred to as the ”trap” below) and VW is the wall potential (zero inside the container and infinite otherwise). The focus here is on the average density of charges n(r) for the equilibrium state, with spherical symmetry in the fluid phase. In the classical Canonical ensemble this is defined by R dr2 ..drN e−βV (r1 ,..,rN ) n(r1 ) = N R . dr1 ..drN e−βV (r1 ,..,rN ) 4 (2) Attention will be focused on a harmonic trap and Coulomb interactions among charges. The total potential energy in dimensionless form is given by # " N N 2 3 X X 1 1 mω r 0 r ∗2 + + VW∗ . V ∗ (r∗1 , .., r∗N ) ≡ βV (r1 , .., rN ) = Γ 2q 2 i=1 i 2 i6=j=1 r∗i − r∗j (3) Here, r∗i = ri /r0 , and Γ is the Coulomb coupling constant Γ = βq 2 /r0 . The usual choice for the length scale r0 is the ion sphere radius, or mean distance between charges, given by 4πr03 n/3 = 1, where n is the spatially averaged density for the confined particles. This uniform density can be characterized by the volume of a sphere inside of which the particles are localized. This localization can be due to either the walls or the trap alone, depending on conditions. Confinement by the trap is determined by the condition that the net force on the outermost charge at the surface of this sphere should be zero Nq 2 = mω 2 R0 , R02 ⇒ R03 = N q2 . mω 2 (4) This applies only if R0 < R; otherwise the particles are confined by the wall. In the case of interest here, R0 < R and the average density is therefore Z 3mω 2 q2 3 3 drn(r) = , r = . n≡ 0 4πR03 4πq 2 mω 2 (5) Note that the average density is independent of Γ and N, although the profile n(r) depends strongly on both. The dimensionless potential (3) is now " N # N X X 1 1 1 ri∗2 + V ∗ (r∗1 , .., r∗N ) = Γ , ∗ 2 i=1 2 i6=j=1 ri − r∗j (6) and the dimensionless number density is ∗ n (r1∗ ) ≡ n(r1 )r03 R ∗ ∗ ∗ ∗ dr2 ..dr∗N e−V (r1 ,..,rN ) = NR ∗ ∗ . ∗ ∗ dr1 ..dr∗N e−V (r1 ,..,rN ) (7) The parameters characterizing the profile in these units are Γ and N. The coupling constant Γ can be interpreted as the inverse temperature in units of the Coulomb energy q 2 /r0 T∗ ≡ r0 1 kB T = kB T 2 = . 2 2 mω r0 q Γ (8) The dimensionless density profile n∗ (r ∗ ) can be calculated numerically by direct Monte Carlo integration of (7) for given Γ and N. 5 III. THEORY AND APPROXIMATIONS For the theoretical analysis of the next sections, it is more convenient to work in the grand Canonical ensemble −1 n(r1 ) = Ξ  N Z ∞ X N zr03 ∗ ∗ ∗ dr∗2 ..dr∗N e−V (r1 ,..,rN ) , 3 N! λ N =0 with the normalization constant  N Z ∞ X 1 zr03 ∗ ∗ ∗ Ξ= dr∗1 ..dr∗N e−V (r1 ,..,rN ) . 3 N! λ N =0 1/2 Here z = eβµ , µ is the chemical potential, and λ = (h2 β/2πm) (9) (10) is the thermal wavelength. This choice for the ensemble allows application of methods from density functional theory [20, 21]. The relevant exact equations are summarized in Appendix A. In particular, the density profile n(r) is related to the excess free energy Fex (β | n) expressed as a functional of that density δβFex (β | n) n(r)λ3 = −βV0 (r) − . (11) ln z δn (r) The excess free energy is a universal functional of the density such that each choice for the density corresponds to a different applied external potential. The case of interest here is the non-uniform density resulting from the harmonic trap, βV0 (r) → Γr ∗2 /2. A different case is a uniform density, which results from an applied potential of a uniform neutralizing charge density. This latter case is the familiar one component plasma (OCP) [22]. Both are described by (11) with the same excess free energy functional; only the external potential is different. In Appendix A two exact representations for Fex (β | n) are given in terms of a uniform reference density. In the first case V0 (rα ) is chosen to be the Coulomb interaction for an ion placed at the origin, plus the uniform neutralizing background charge. The density in (11) ∗ then becomes the pair distribution function gOCP (r ∗ ) in the OCP. The appropriate reference density is its uniform value far from the ion, for which (11) becomes (see Appendix A) Z ∗ ∗ ∗−1 ∗ ∗ ln gOCP (r ) = −Γr + dr∗′ {gOCP (r ∗′ ) − 1} cOCP (|r∗ − r∗′ | ; Γ) − ΓBOCP (r ∗ | gOCP ), (12) together with the Ornstein-Zernicke equation for the direct correlation function cOCP Z ∗ ∗ ∗ ∗ (r ∗′ ) − 1} cOCP (|r∗ − r∗′ | ; Γ) . (13) (gOCP (r ) − 1) = cOCP (r ) + dr∗′ {gOCP 6 ∗ The function BOCP (r ∗ | gOCP ) is referred to as the OCP ”bridge function” [22]. A similar rearrangement of (11) is possible when V0 is chosen to be a harmonic trap using a different representation for Fex (β | n), to obtain an exact equation for the trap density n(r). In this case the relevant reference density far from the trap center is zero. The analysis is outlined in Appendix A and the resulting form for Eq. (11) is Z 1 ∗2 n(r)λ3 = −Γ r + dr∗′ n∗ (r ∗′ )cOCP (|r∗ − r∗′ | ; Γ) − ΓB (r ∗ | n∗ ) . ln z 2 (14) The second term is expressed in terms of the correlations for the uniform OCP, cOCP (r ∗ ; Γ), evaluated at the average trap density n of (5). The function B (r ∗ | n∗ ) differs from ∗ BOCP (r ∗ | gOCP ), but by analogy will be referred to as the trap bridge function. Equations (12) - (14) provide the formally exact description from which approximations are formulated for the trap density. Further details of their context are given in the Appendix. Correlations among charges occur in these expressions through the bridge functions and cOCP . The simplest approximation is the ”mean field theory” resulting from the neglect of all correlations, B → 0, BOCP → 0, and cOCP (r ∗ ) → −Γr ∗−1 , for which (14) and (13) become Z n(r)λ3 1 ∗2 −1 Γ ln = − r − dr∗′ n∗ (r ∗′ ) |r∗ − r∗′ | , z 2 Z −1 ∗ ∗ ∗−1 ∗ ∗′ ∗ ∗′ −1 Γ (gM F (r ) − 1) = −r − dr∗′ {gM . F (r ) − 1} |r − r | −1 (15) (16) Equation (16) gives the Debye-Hückel result for a weakly coupled OCP. The solution to (15) for the mean field description of the trap density is discussed below. This mean field equation was studied for Coulomb and Yukawa traps in references [13, 14] for the special case Γ−1 = T ∗ = 0. To account for correlations the HNC approximation is used, which results from neglecting the bridge functions B → 0, BOCP → 0 but retaining the effects of cOCP (which in this approximation is denoted cHN C ) Z 1 ∗2 n(r)λ3 = −Γ r + dr∗′ n∗ (r ∗′ )cHN C (|r∗ − r∗′ | ; Γ) , ln z 2 Z ∗ ∗ ∗−1 ∗ ∗′ ∗ ∗′ ln gHN C (r ) = −Γr + dr∗′ {gHN C (r ) − 1} cHN C (|r − r | ; Γ) , Z ∗ ∗ ∗ ∗ ∗′ ∗ ∗′ (gHN C (r ) − 1) = cHN C (r ) + dr∗′ {gHN C (r ) − 1} cHN C (|r − r | ; Γ) . 7 (17) (18) (19) ∗ Equations (18) and (19) are a closed set of equations to determine gHN C and cHN C for the OCP, a well-studied problem [22]. Equation (17) is an extension of the HNC to an approximation for the localized trap density. It is shown below that the inclusion of correlations leads to qualitatively different results for the density compared to results from the mean field theory. To go beyond the HNC approximation requires some estimate of the bridge functions. This is a difficult problem in general. For the OCP, effects of BOCP are important only at very strong coupling. However, for the trap density the bridge functions have a quantitative importance even at moderate coupling. This is discussed in Section VI where a phenomenological estimate is explored and shown to remove the observed discrepancies between HNC and Monte Carlo simulation results. Further work in this direction is planned for future studies. In the following it is convenient to write the above equations using a representation for the density in terms of a dimensionless effective potential U ∗ (r ∗ ), defined by ∗ e−ΓU ∗ n (r ) ≡ N 4π R R∗ 0 ∗ (r ∗ ) dr ∗ r ∗2 e−ΓU ∗ (r∗ ) . (20) Here the chemical potential µ has been eliminated in favor of the average number of particles R N = drn(r). Substitution of (20) into (14) gives the corresponding exact equation for U ∗ (r ∗ ) U ∗ 1 r , Γ, N = r ∗2 + N 2 ∗
R dr∗′ e−ΓU (r ,Γ,N ) c(|r∗ − r∗′ |, Γ) + B (r | n) , R ∗′ ∗ dr∗′ e−ΓU (r Γ,N ) ∗ ∗′ (21) with the notation c(r ∗ , Γ) ≡ −Γc(r ∗ , Γ). The solution is parameterized by only two dimensionless constants, Γ, N . IV. MEAN FIELD THEORY To develop a qualitative understanding of the dependence of the profile on the parameters Γ and N it is useful to explore this simple mean field approximation in the dimensionless variables of the last section. Taking the mean field limit of (21), and performing the angular 8 integration in the second term gives the mean field equation U ∗ r ∗ , Γ, N
= 1 ∗2 N r + R R∗ ∗ ∗ 2 dr ∗ r ∗2 e−ΓU (r ,Γ,N ) 0   ZR∗ Zr∗ ∗ ∗ 1 ∗ ∗ ×  ∗ dr ∗′ r ∗′2 e−ΓU (r ,Γ,N ) + dr ∗′ r ∗′ e−ΓU (r ,Γ,N )  . r (22) r∗ 0 For high temperatures T ∗ = Γ−1 >> 1 (22) gives  
3N 1 N ∗ ∗ 1 − ∗3 r ∗2 + , U r , Γ, N → 2 R 2R∗ (23) where an outer hard wall at R∗ has been restored. Consequently the trap density becomes ∗ n ∗
N r , Γ, N −→ 4π R R∗ 0 −Γ 12 dr ′ r ′ 2 e “ ” “ −Γ 12 1− N∗3 r ∗2 R ” . e 1− N∗3 r ′ 2 (24) R It is understood that r ∗ ≤ R∗ and the density is zero otherwise. For large R∗ (24) is just the expected Boltzmann factor for uncorrelated particles in an external harmonic potential. However, as the number of particles in the trap increases their Coulomb repulsion competes with the confining trap potential, enhancing the density toward the outer wall. For N < R∗3 the density increases toward the center of the trap, while for N > R∗3 it increases in the direction of the outer wall. The density is uniform at the special value R∗ = N 1/3 = R0∗ , the point where the Coulomb repulsive force and attractive confinement force exactly balance. At low temperatures T ∗ = Γ−1 << 1 and e−ΓU ∗ (r ∗ ) → 0 in (22) unless U ∗ (r ∗ ) is constant, in which case these factors cancel in the numerator and denominator. Therefore, a solution is sought of the form U ∗ ∗
r , Γ, N =   U0∗ , ∗ r ∗ < rmax  U ∗ (r ∗ ) > 0, ∗ r ∗ > rmax , ∗ ∗ for some constant U0∗ and rmax < R∗ . Then, for Γ >> 1 (22) gives for r ∗ < rmax  
N 3N 1 ∗ ∗ 1 − ∗3 r ∗2 + ∗ U r , Γ, N → 2 rmax 2rmax (25) (26) ∗ The first term must vanish for consistency with (25), which determines rmax ∗ rmax = R0∗ = N 1/3 . (27) Next, for Γ >> 1 and R0∗ < r ∗ ≤ R∗ (22) gives
1 N U ∗ r ∗ , Γ, N → r ∗2 + ∗ . 2 r 9 (28) 0.3 Γ=0.1 Γ=1 Γ=10 Γ=20 Γ→∞ n*(r*) 0.2 0.1 0 0 2 4 r* 6 8 FIG. 1: (color online) Mean field results for N = 100 for various values of Γ. The Coulomb limit (Γ → ∞) is also shown. Consequently, the trap density in this limit is  1, r ∗ < R0∗ 3N  ” “ . n∗ (r ∗ ) → 4πR0∗3  e−Γ 21 r∗2 + rN∗ , r ∗ > R0∗ (29) The density is uniform up to R0∗ which depends only on N , and is exponentially small for larger r ∗ ; the density effectively has a step profile. More generally, (22) can be solved numerically for the full range of Γ and N . The dependence on Γ is illustrated in Figure 1 for N = 100 and Γ = 0.1, 1, 10, and 20. Also shown is the limiting step function for Γ → ∞. The edge of the profile is the radius at which the confinement force equals the total Coulomb force, R0∗ = N 1/3 . Although instructive as a reference, the low temperature behavior within the mean field approximation is relevant only for the average density, since charge correlation effects become dominant for Γ & 10, as shown in the next section. 10 3 Γ=0.1 Γ=1 Γ=10 Γ=100 Coulomb limit 2.5 c(r*,Γ) 2 – 1.5 1 0.5 0 0 0.5 1 r* 1.5 2 FIG. 2: (color online) Dependence of the scaled direct correlation function on Γ. Also shown is the small Γ Coulomb limit (mean-field approximation ). V. HNC APPROXIMATION The density in the HNC approximation is given by (17) - (19). The effective potential is determined from the equation U ∗ 1 r , Γ, N = r ∗2 + N 2 ∗
R dr∗′ e−ΓU (r ,Γ,N ) cHN C (|r∗ − r∗′ |, Γ) , R ∗′ ∗ dr∗′ e−ΓU (r Γ,N ) ∗ ∗′ (30) The only difference from the mean field form is the replacement of the Coulomb interaction by the direct correlation function cHN C . This is first calculated from the corresponding HNC approximation for the OCP, (18) and (19), and the results are illustrated in Figure 2 for Γ = 0.1, 1, 10, and 100. It is seen that the Coulomb limit is recovered for all r ∗ in the limit Γ → 0, and for any Γ at sufficiently large r ∗ . With cHN C known, (30) can be solved to determine the trap density. Figure 3 shows the results for the same conditions as in Figure 1, illustrating the effects of correlations relative to the mean field results. The monotonic decrease with r ∗ is now replaced by the appearance of structure (”shells”) for Γ = 10 and 20. The latter two curves also suggest that the location of the shells is insensitive to the value of Γ, as is confirmed below. Figure 4, for Γ = 40, 11 0.4 Γ=0.1 Γ=1 Γ=10 Γ=20 0.35 0.3 n*(r*) 0.25 0.2 0.15 0.1 0.05 0 0 1 2 3 4 5 r* FIG. 3: (color online) Radial density profile with correlations included, for the same conditions as Figure 1. 0.6 N=50 N=100 N=300 0.5 n*(r*) 0.4 0.3 0.2 0.1 0 0 1 2 3 4 5 6 7 r* FIG. 4: (color online) Dependence of radial density profile on N̄ for Γ = 40. 12 80 Γ=10 Γ=20 Γ=40 MC (Γ=12.6) MC (Γ=126) Ni (number in each shell) 70 60 50 40 30 20 10 0 20 30 40 50 60 70 80 N (total number in system) 90 100 FIG. 5: (color online) Shell population for two-shell conditions. The higher trend curve is for the the outer shell. Data is shown for HNC and MC calculations. shows that the number of ”shells” increases with N . Figure 5 illustrates the formation and filling of the first three shells, giving the number of particles in each shell (defined as the number between two spheres at successive radial minima) as a function of N for Γ = 10, 20, and 40. The initial dependence is linear as only a single outer shell is populated. However, after formation of the second shell the dependence is more complex as additional particles can go into either shell. Remarkably, this process of filling is almost independent of Γ (in agreement with experiment and simulations [16, 17]). VI. COMPARISON TO MONTE CARLO SIMULATION To assess the quality of predictions based on the HNC approximation some benchmark Monte Carlo simulations have been performed for comparison. For moderate coupling, Γ = 1, correlations are important (see Figure 2, for example) but shell structure does not appear. Figure 6 shows the comparison with Monte Carlo results for Γ = 1, and N = 10, 100, and 1000. Not shown are results for N < 10 where the agreement is poor even at weak coupling. The first conclusion is that the HNC approximation is accurate at weak 13 0.3 HNC N=10 MC N=10 HNC N=100 MC N=100 HNC N=1000 MC N=1000 n*(r*) 0.2 0.1 0 0 2 4 6 8 r* 10 12 14 16 FIG. 6: (color online) Comparison with Monte Carlo results for Γ = 1. to moderate coupling, Γ ≤ 1, and N ≥ 10. The formation of shell structure at larger values of Γ is tested in Figure 7 for Γ = 20, N = 25 and 100, and in Figure 8 for Γ = 40, N = 300. It is seen that the formation and location of the peaks is very well described, but systematically their widths are too large and amplitudes too small by as much as 20 − 40%. Nevertheless, the particle number in the shells is quite accurate, as shown in Figure 5. VII. ADJUSTED HNC The HNC approximation provides a surprisingly good representation of correlation effects on the density profile of strongly coupled charges in a trap: formation, population, location, and number of shells. Its primary limitation is the significant discrepancy in the amplitudes and widths of the shells. This is due to neglect of important contributions from the bridge functions B (r | n) and BOCP (r | n) in (12) and (14) at strong coupling. This problem was addressed some time ago for the OCP by Ng [23], who observed that the peaks of the ∗ HNC pair distribution function gOCP (r ∗ ) appear at the correct position but with incorrect amplitudes and widths. This is clearly analogous to the current errors of the trap density. 14 HNC N=25 MC N=25 HNC N=100 MC N=100 0.6 n*(r*) 0.4 0.2 0 0 1 2 3 r* 4 5 6 FIG. 7: (color online) Comparison with Monte Carlo results for Γ = 20. Ng corrected the HNC by an estimate of the OCP bridge function in the form BOCP (r | n) → λ (Γ) V0 (r ∗ ) , V0 (r ∗ ) = 1 , r∗ (31) where λ (Γ) is a chosen function to interpolate between λ (0) = 0 and some constant value λ (∞). The advantage of the chosen form is that when inserted in (12) the form of the HNC approximation is recovered, but with a renormalized coupling constant Γ′ = (1 + λ (Γ)) Γ Z ∗ ∗ ′ ∗−1 ∗ (r ∗′ ) − 1) cOCP (|r∗ − r∗′ | ; Γ′ ) . (32) ln gOCP (r ) = −Γ r + dr∗′ (gOCP Note that cOCP (|r∗ − r∗′ | ; Γ′ ) also becomes a function of the renormalized coupling constant ∗ since it inherits that dependence only from gOCP (r ∗ ) via the Ornstein-Zernicke equation (13). ∗ Ng found that the Monte Carlo data for gOCP (r ∗ ) could be fit within two percent with the choice λ (∞) = 0.6. Alternatively, (31) can be viewed as a representation of the OCP bridge ∗ function determined from (13) using gOCP (r ∗ ) from Monte Carlo as input. It is plausible to use this same scheme for the trap bridge function B (r | n) as well B (r | n) → λ (Γ) V0 (r ∗ ) , 1 V0 (r ∗ ) = r ∗2 , 2 (33) with the same choice for λ (Γ). Accordingly, (14) gives the HNC form for the trap density Z n(r)λ3 ′ 1 ∗2 = −Γ r + dr∗′ n∗ (r ∗′ )cOCP (|r∗ − r∗′ | ; Γ′ ) . (34) ln z 2 15 0.7 HNC N=300 MC N=300 0.6 n*(r*) 0.5 0.4 0.3 0.2 0.1 0 0 1 2 3 4 r* 5 6 7 8 FIG. 8: (color online) Comparison with Monte Carlo results for N̄ = 300 and Γ = 40. This is the same as the HNC approximation of Section V, except with the bridge function corrected value Γ′ . Based on the results of Ng, the same value of λ (∞) = 0.6 for the strong coupling conditions is used here. Equations (32) and (34) will be referred to as the adjusted HNC (AHNC) approximation. Figures 9 and 10 show a revised comparison for the conditions of figures 7 and 8 but now with the AHNC calculated with Γ′ = 1.6Γ. Most of the discrepancies are seen to be removed by this AHNC approximation. It describes accurately the complexities of the trap density over a range of Γ, N including strong coupling, with the structural simplicity of the HNC approximation and only OCP correlations as input. Its domain of applicability and limitations will be explored in more detail elsewhere. It is somewhat remarkable that the simple forms (31) and (33) for the bridge functions work so well. This success should motivate future work to address the theoretical basis for such practical approximations to the bridge functions. VIII. DISCUSSION The objective here has been the development and exploration of a theoretical framework for analysis of the density profile of charges in a harmonic trap. Previous extensive work 16 AHNC N=25 MC N=25 AHNC N=100 MC N=100 0.6 n*(r*) 0.4 0.2 0 0 1 2 3 r* 4 5 6 FIG. 9: (color online) Comparison of adjusted HNC with Monte Carlo results for Γ = 20. 0.7 AHNC N=300 MC N=300 0.6 n*(r*) 0.5 0.4 0.3 0.2 0.1 0 0 1 2 3 4 r* 5 6 7 8 FIG. 10: (color online) Comparison of adjusted HNC with Monte Carlo results for N̄ = 300 and Γ = 40. 17 on this problem has focused on experimental and simulation studies of the low temperature crystal structure (and its melting) for Coulomb, Yukawa, and Lennard-Jones interactions. Here the complementary theoretical study starting from the high temperature fluid phase has been described, showing the evolution of shell structure at lower temperatures due to correlations that are precursors of the observed structure in the crystal phase. A first order representation of correlations in the HNC approximation is shown to capture the radial structural features of the density profile. Particularly remarkable is the formation of ”blurred” shells for moderate coupling at the magic numbers for tiling the sharp crystal shells. This anticipation in the fluid phase of the low temperature crystal structure will be discussed elsewhere. The analysis above provides justification for a simple and practical description within the fluid phase for charged particle confinement, the AHNC approximation. Although the discussion here has been limited to the case of Coulomb interactions there is no complication involved in extending it to other more practical forms such as Yukawa particles. The formulation of this approximation within the context of density functional theory also provides a complete thermodynamics for this system. In this way it is hoped that the AHNC could be improved to include a description of freezing and crystal structure. IX. ACKNOWLEDGEMENTS This work is supported by the Deutsche Forschungsgemeinschaft via SFB-TR 24, and by the NSF/DOE Partnership in Basic Plasma Science and Engineering under the Department of Energy award DE-FG02-07ER54946. APPENDIX A: REPRESENTATIONS FROM DENSITY FUNCTIONAL THEORY Consider a given system in an external potential. For simplicity only the case of spherically symmetric central potentials are considered. Two classes of external potential are identified. In the first class, the potential goes to zero or a constant at large r, resulting in an average density that becomes uniform at large distances. This is the usual case of localized perturbations in a bulk fluid. The second class is potentials that become unbounded 18 for large r, so that the resulting average density is localized in some region about the origin. This is the case of interest here for charges in a trap. In the following, different exact representations of the inhomogeneous density appropriate for each of these cases are given. The approach is based on the formal structure of density functional theory [20, 21]. To review briefly the relevant equations, the grand potential is defined first by βΩ (β | µ − V0 ) = − ln Ξ ∞ X T re−β(H−µN ) , (A1) N =0 where Ξ is the grand partition function of (10). The notation makes explicit the functional dependence of Ω on the external potential V0 . The density is obtained from the functional derivative of Ω n(r) = δΩ . δV0 (r) The roles of n(r) and V0 (r) can be exchanged by the Legendre transformation Z F (β | n) ≡ Ω (β | µ − V0 ) + dr (µ − V0 (r)) n (r) , δF (β | n) = µ − V0 (r) . δn (r) (A2) (A3) (A4) The new functional F (β | n) is the free energy expressed as a universal functional of the density. Equations (A2) and (A4) make explicit the fact that the density is a functional of the applied potential and vice versa. The functional F (β | n) is universal in the sense that each choice for V0 corresponds to evaluation of this functional at a corresponding specific density field. Equation (A4) will be applied here as an exact equation for the density in terms of the chosen V0 . The free energy can be divided into an ”ideal gas” form F0 (β | n), (its form in the absence of interactions among the particles), plus the ”excess free energy” due to interactions, Fex (β | n) F (β | n) = F0 (β | n) + Fex (β | n). (A5) The first term can be calculated directly so that (A4) takes the more practical form of (11) ln δβFex (β | n) n(r)λ3 = −βV0 (r) − . z δn (r) (A6) A further decomposition of Fex is useful for the introduction of approximations. The central idea is to exploit knowledge of correlations in uniform systems to characterize the 19 non-uniform density profile due to V0 . First, note the identity for any functional X(g) of the function g (r) X(g) − X(gr ) = Z 1 dλ∂λ X [λg + (1 − λ) gr ] Z 1 Z δX (g) = dλ dr |λg+(1−λ)gr (g (r) − gr (r)) δg (r) 0 0 (A7) where gr (r) is an arbitrary reference function. Application of this identity twice to the excess free energy Fex gives Z Fex (β | n) = Fex (β | nr ) − drβ −1c(1) (r |nr ) (n (r) − nr ) Z 1 Z λ Z ′ − dλ drdr′β −1 c(2) [r, r′ |λ′ n + (1 − λ′ ) nr ] (n (r) − nr ) (n (r′ ) − nr ) . dλ 0 0 (A8) Here nr is an arbitrary reference density, and the direct correlation functions have been introduced by the definitions c(m) (r1 , r2 ..rm | n) ≡ −β δ m Fex (β | n) . δn (r1 ) δn (r2 ) ..δn (rm ) (A9) For external potentials V0 that vanish for large r there is a natural reference density n corresponding to the uniform system at distances far from the perturbation. Then (A2) becomes Z 1 dλ (1 − λ) Fex (β | n) = Fex (β | n) − 0 Z × drdr′β −1 c(2) [r, r′ |λ′ n + (1 − λ′ ) n] (n (r) − n) (n (r′ ) − n) . (A10) Here the second term of (A8) vanishes by requiring that nV = N , the average total particle number. Again applying (A7) to express c(2) (r, r′ |λ′ n + (1 − λ′ ) n) in terms of the uniform density gives Z 1 Fex (β | n) = Fex (β | n) − drdr′β −1 c(2) (|r − r′ | |n) (n (r) − n) (n (r′ ) − n) 2 Z 1 Z λ Z ′ − dλ (1 − λ) drdr′ dr′′ (n (r) − n) (n (r′ ) − n) (n (r′′ ) − n) dλ 0 0 × c(3) [r, r′ , r′′ | λ′ n + (1 − λ′ )n] 20 (A11) Next, consider external potentials that become unbounded for large r so that the density vanishes at large distances. In this case the appropriate choice for the reference density in (A8) is nr = 0, giving Fex (β | n) = Z 1 dλ (1 − λ) 0 Z drdr′n (r) n (r′ ) β −1 c(2) (r, r′ | λn) , (A12) The first two terms on the right side of (A8) vanish at zero density, and the first λ integral has been performed. The integration domains are now controlled by the density, which restricts them to the bounded domain. In this domain, it is reasonable to identify an appropriate uniform reference density n and to express c(2) (r, r′ | λn) in terms of this uniform density by a further application of (A7) c (2) ′ (r, r | λn) = c (2) ′ (r, r | n) + Z 1 dx∂x c(2) [r, r′ | xλn + (1 − x)n] 0 = c(2) (r, r′ | n) Z 1 Z + dx dr′′ (λn (r′′ ) − n) c(3) [r, r′, r′′ | xλn + (1 − x)n] . (A13) 0 Then the excess free energy becomes Z 1 drdr′n (r) n (r′ ) β −1 c(2) (|r − r′ | | n) Fex (β | n) = − 2 Z 1 Z 1 Z   − dλ (1 − λ) dx dr′′ β −1 c (3) r, r′ , r′′ | xλn + (1 − x)n ((λn (r′′ ) − n)) , 0 0 (A14) Equations (A11) and (A14) are both formally exact and equivalent, but each is formulated so that the leading terms are reasonable approximations for different classes of external potentials. Use of (A11) in (A6) leads to the form (12) of the text; use of (A14) in (A6) leads to the form (14). The Ornstein-Zernicke equation (13) is an identity that follows from the definition of c(2) in (A9). The last terms of (A11) and (A14) give rise to what is referred to as bridge functions in the text, and provide the formal definitions for each. The above is quite general. For the case of interest here, a system of charges, the uniform density case corresponds to the one component plasma. [1] D.H.E. Dubin, and T.M. O’Neill, Rev. Mod. Phys. 71, 87 (1999). 21 [2] 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, Physics of Plasmas 15, 055704 (2008). [3] D. J. Wineland, J. C. Bergquist, W. M. Itano, J. J. Bollinger, and C. H. Manney, Phys. Rev. Lett. 59, 2935 (1987). [4] M. Drewsen, C. Brodersen, L. Hornekær, J. S. Hangst, and J. P. Schiffer, Phys. Rev. Lett. 81, 2878 (1998). [5] O. Arp, D. Block, A. Piel, and A. Melzer, Phys. Rev. Lett. 93, 165004 (2004). [6] M. Bonitz, D. Block, O. Arp, V. Golubnychiy, H. Baumgartner, P. Ludwig, A. Piel, and A. Filinov, Phys. Rev. Lett. 96, 075001 (2006). [7] A. V. Filinov, M. Bonitz, and Y. E. Lozovik, Phys. Rev. Lett. 86, 3851 (2001); A.V. Filinov, Yu.E. Lozovik, and M. Bonitz, Phys. Stat. Sol. (b) 221, 231 (2000). [8] T. Pohl, T. Pattard, and J.M. Rost, Phys. Rev. Lett. 92, 155003 (2004). [9] T. C. Killian, S. Kulin, S. D. Bergeson, L. A. Orozco, C. Orzel, and S. L. Rolston, Phys. Rev. Lett. 83, 4776 (1999); S. Kulin, T.C. Killian, S. D. Bergeson, and S. L. Rolston, Phys. Rev. Lett. 85, 318 (2000); T. C. Killian, M. J. Lim, S. Kulin, R. Dumke, S. D. Bergeson, and S. L. Rolston, Phys. Rev. Lett. 86, 3759 (2001). [10] P. Ludwig, S. Kosse, and M. Bonitz, Phys. Rev. E 71, 046403 (2005). [11] O. Arp, D. Block, M. Bonitz, H. Fehske, V. Golubnychiy, S. Kosse, P. Ludwig, A. Melzer, and A. Piel, J. Phys. Conf. Series 11, 234 (2005). [12] H. Baumgartner, D. Asmus, V. Golubnychiy, P. Ludwig, H. Kählert, and M. Bonitz, New Journal of Physics 10, 093019 (2008). [13] C. Henning, H. Baumgartner, A. Piel, P. Ludwig, V. Golubnichiy, M. Bonitz, and D. Block, Phys. Rev. E 74, 056403 (2006). [14] C. Henning, P. Ludwig, A. Filinov, A. Piel, and M. Bonitz, Phys. Rev. E 76, 036404 (2007). [15] J. Wrighton, J.W. Dufty, C. Henning, and M. Bonitz, J. Phys. A: Math. Gen. 42, 214052 (2009). [16] V. Golubnychiy, H. Baumgartner, M. Bonitz, A. Filinov, and H. Fehske, J. Phys. A: Math. Gen. 39, 4527 (2006). [17] H. Baumgartner, H. Kählert, V. Golubnychiy, C. Henning, S. Käding, A. Melzer, and M. Bonitz, Contrib. Plasma Phys. 47, 281 (2007). 22 [18] M. Bowick, A. Cacciuto, D. Nelson, and A. Travesset, Phys. Rev. Lett. 89 185502 (2002); M. Bowick and L. Giomi, arXiv:0812.3064 (submitted to Advances in Physics 2009). [19] B. Talin, A. Calisti, J. Dufty, and I. Pogorelov, Phys. Rev. E 77, 036410 (2008); J. Dufty, I. Pogorelov, B. Talin, and A. Calisti, J. Phys. A 36, 6057 (2003); J. Dufty, B. Talin, and A. Calisti , in Theory of Energy Deposition, Adv. Quant. Chem. 46, 293 (2004); J. Wrighton and J. Dufty, J. Stat. Mech. P10021 (2008). [20] R. Evans, in Fundamentals of Inhomogeneous Fluids, edited by D. Henderson (Marcel Dekker, NY, 1992). [21] J. Lutsko, Recent Developments in Classical Density Functional Theory, preprint (2000). [22] J-P Hansen and I. MacDonald, Theory of Simple Liquids, (Academic Press, San Diego, CA, 1990). [23] K-C Ng, J. Chem. Phys. 61, 2680 (1974). 23 View publication stats