PHYSICAL REVIEW E 91, 062140 (2015)

Counterion density profile around a charged disk: From the weak to the strong association regime Juan Pablo Mallarino* Facultad de Ciencias, Laboratorio Computacional HPC, Universidad de los Andes, Bogot´a, Colombia

Gabriel T´ellez† Departamento de F´ısica, Universidad de los Andes, Bogot´a, Colombia (Received 18 February 2015; published 29 June 2015) We present a comprehensive study of the two-dimensional one-component plasma in the cell model with charged boundaries. Starting from weak couplings through a convenient approximation of the interacting potential we were able to obtain an analytic formulation to the problem deriving the partition function, density profile, contact densities, and integrated profiles that compared well with the numerical data from Monte Carlo simulations. Additionally, we derived the exact solution for the special cases of  = 1,2,3, . . . , finding a correspondence between those from weak couplings and the latter. Furthermore, we investigated the strong-coupling regime taking into consideration the Wigner formulation. Elaborating on this, we obtained the profile to leading order, computed the contact density values as compared to those derived in an earlier work on the contact theorem. We formulated adequately the strong-coupling regime for this system that differed from previous formulations. Ultimately, we calculated the first-order corrections and compared them against numerical results from our simulations with very good agreement; these results compared equally well in the planar limit, whose results are well known. DOI: 10.1103/PhysRevE.91.062140

PACS number(s): 05.20.Jj, 05.10.−a, 02.70.−c

I. INTRODUCTION

In this work we present a thorough analysis of the condensation phenomenon of counter-ions around a charged disk. The model under consideration is a two-dimensional (2D) system formed by an impenetrable disk of charge Q1 surrounded by ions dispersed freely in a larger disk with external charged boundary Q2 and no dielectric discontinuities between the regions delimited by the geometry. The problem resembles an annulus with particles moving freely between the inner and outer radii as in Fig. 1. The N ions have, respectively, a charge −q in such a way that neutrality yields Q1 + Q2 = N q.

(1.1)

We will assume a pointlike geometry for the free charges, which is not a problem due to electrostatic repulsion alone. This model is seemingly the one-component plasma (2D-OCP) with a small variation. First, the neutralizing charge is not distributed homogeneously in the background and, second, the inner core is impenetrable. This is the two-dimensional analog of the Manning counterion condensation phenomenon around charged cylinders [1–3]. Unlike the 3D situation where the Coulomb potential shapes as 1/|r|, the partition function for two-dimensional Coulomb systems is written as a product of contributions which, in some cases, may be computed exactly. Moreover, the logarithmic nature of the potential motivated the theoretical computation of the abundant static and dynamic properties of electrolytes for two-dimensional systems.

* †

[email protected] [email protected]

1539-3755/2015/91(6)/062140(28)

The interaction between two unit charges separated by a distance r is given by the two-dimensional Coulomb potential − ln(r/L), where L is an irrelevant arbitrary length scale. We are interested in the equilibrium thermal properties of the system at a temperature T . As usual, we define β = 1/(kB T ), where kB is the Boltzmann constant. There are three important dimensionless parameters which characterize the system. The 2D equivalent of the Manning parameter is ξ = βqQ1 /2, which characterizes the strength of the Coulomb coupling between the inner disk and the counterions. The Coulomb coupling between counterions is  = βq 2 /2. Here the familiarized reader may have noticed that for the two-dimensional systems it is typical to use  for the notation of the coupling constant [4–8]; the previously defined coupling equates the standardized one as  = /2, but we have kept  motivated by a discussion with the three-dimensional case. Finally, the third parameter, defined as  = ln(D/R), measures in log scale the size of the system. For our problem, we acknowledge the published works from Refs. [9–11], which have covered thoroughly the basics of the 2D construction of the cell model. They observed that condensation obeys the well-known threshold at ξ = 1 where below this value there is none. In the mean field, it is known that a disk or cylinder with a dimensionless charge inferior to unity is unable to bind counterions. When the charge is above unity, it attracts an ion cloud in such a way that it neutralizes partially the disk or cylinder so the effective dimensionless charge of the disk or cylinder and the cloud is unity. Hence, the fraction of condensed ions in the cloud corresponds to the excess above unity by the total, also known as the Manning fraction fM = 1 − 1/ξ . They also showed that the system with two dimensions differs from the three-dimensional one in several ways. First, when the outer boundary is uncharged, Q2 = 0, the relationship between the coupling parameter and the Manning parameter is dictated by neutrality, which, unlike

062140-1

©2015 American Physical Society

´ JUAN PABLO MALLARINO AND GABRIEL TELLEZ

PHYSICAL REVIEW E 91, 062140 (2015)

Finally, we study, in Sec. VI, the strong-coupling situation when   1 (or   2) with special attention to condensation and evaluate the profile and other statistical properties. Ultimately, we will retake the contact theorem to obtain the value of the densities at contact comparing to the values derived by the models and numerical Monte Carlo simulations. II. THE MODEL

The Hamiltonian for the system considering the logarithmic Coulomb potential interaction between charges as the solution to the Poisson equation ∇r2 v(r) = 2π δ(r) reads as H : = βH

in 3D where they both enjoyed independence,1 it fixes one or the other in such a way that  = ξ/N,

(1.2)

as will be discussed in Sec. II. Additionally, if we were to look at a system with a fixed number of particles going from a temperature above the critical temperature for condensation, i.e., ξ < 1, and allow it to undergo a cooling process, the energy, heat capacity, and other quantities will present successive transitions due to iterative localization phenomena occurring when a counterion is condensed [9–11]. The interesting regime for the two-dimensional construction is that which corresponds to a small number of counterions. Due to Eq. (1.2), a large number N of ions is equivalent to  → 0, which is unquestionably the mean-field regime. Therefore, as mentioned in Ref. [9], our interest is in the range for weak and strong couplings controlled by the ratio between ξ and the number of counterions. The outline for the following work begins with the presentation of the two-dimensional model (Sec. II), followed in Sec. III by an analysis of the regime when  < 1 (or  < 2), where we extend previous works [9–12]. Then, in Secs. IV and V, the study of the special cases  = 1,2,3, . . . (or  = 2,4,6, . . . ), which are exactly solvable using expansions of powers of Vandermonde determinants in basis of symmetric or antisymmetric monomials [4,5,13,14]. Through these analytic models we will be able to obtain the profiles and the energy along with all the quantities associated with them. Furthermore, we will present a model for condensation in each case.

1

The dimensionless parameters in 3D are as follows:

N 

|rj | D + NβQ2 q ln L L j =1     rj − rk   − βq 2 ln   L 1j |, {2 ln |r1 − r2 |}eff := 2π 0 (3.2)

Burak et al. realized that the Boltzmann factor could be written conveniently as a product of functions and the form of the integral involved is a successive convolution of functions. Using the Laplace transformation we can evaluate the partition function analytically. Arranging the set of {yj }’s by the [BO] H reads, in simplified form, H =

(3.3)

N 

aj (yj +1 − yj ) − aN  =

j =0

N 

(aj −1 − aj )yj , (3.7)

j =1

where y0 = 0 and yN+1 = , and the {aj } are conveniently defined as  ξ −1 aj = j − + 

positive constants 1 2

 2

ξ −1 − 2(j − 1). (3.8)  The choice for the set of {aj } is not unique since we have a total of N + 1 variables and N conditions; then an arbitrary aj −1 − aj = 2

2

where |r> | ≡ max(|r1 |,|r2 |); consequently, 2 ln |r1 − r2 |  2 ln |r> |.

(3.4)

With the partition function defined as  1 Z= d 2N r e−H(r1 ,r2 ,...,rN ) . N!

Notice that the excess free energy is then

This approximation is equivalent to computing the angular correlations using a cumulant expansion of the Boltzmann factor 062140-3

exc = − ln Z + ln(V N ). F

´ JUAN PABLO MALLARINO AND GABRIEL TELLEZ

PHYSICAL REVIEW E 91, 062140 (2015)

choice for any term will define the rest. Anticipating the following steps, choosing each of the terms to be positive will be convenient in the calculation of the Laplace transformation of the partition function. Note that aj has a minimum sitting on √

jmin =

a0 =

ξ −1 1 + ,  2

j* j*sup j*inf

(3.9) 1

thus indicating that the smallest term of the set is that in which j is the integer closest to jmin . Through this construction we are now able to evaluate the partition function writing the integral as a convolution of functions {fj }. From Eq. (3.6),   2π R 2 N aN −E0  e dyN Z=  0  yN  y2 N  × dyN−1 . . . dy1 fj (yj +1 − yj ) 

0

 =

0

2π R 

FIG. 2. (Color online) The value for j  as a function of fM N /(1 + ξB /ξ ). Notice that the solution is the composition of the red (dark) and green (light gray) parts which tells us that j  =

fM N/(1 + ξB /ξ ) + 1 and signals the discontinuities precisely at the points where fM N /(1 + ξB /ξ ) is a whole number. The dashed line represents j  = fM N /(1 + ξB /ξ ).

j =0

2 N

eaN −E0

× [fN ⊗ fN−1 ⊗ · · · ⊗ f1 ⊗ f0 ](),

(3.10)

with fj (x) = e−aj x , whereas the Laplace transform of the convolution part is then  T{f[] (s) = N ⊗···⊗f0 }

∞

[fN ⊗ fN−1 ⊗ · · · ⊗ f1 ⊗ f0 ]

0

× () e =

−s()

d()

N 

1 . s + aj j =0

(3.11)

The milestone of this analysis is the approximation of the interacting potential term, valid for a large box size (  1). A small box size surfaces other effects which cannot be neglected. The route of the simplification comes from transforming the two-dimensional gas into a one-dimensional mean interacting strip of ions [11]. As we proceed to invert the Laplace transform to obtain the partition function, factors of the form e−aj  will emerge in the function where the dominating contributions will come from the smallest of the {aj }’s; this fact enforces choosing a positive value for aj . Hence, the prevailing√term is the smallest of the set, indicated by j  , and closest to a0 , or3  jinf

=

fM N 1+

ξB ξ

1 + 2

 j  

fM N 1+

ξB ξ

1 + 2

However, choosing the lower or upper bounds depends on whether

fM N 1  ajinf  ajsup ⇒ j = + 2 1 + ξξB   fM N 1    . ajinf > ajsup ⇒j = + 2 1 + ξξB   Since jsup = jinf + 1 we can rewrite the previous condition  using Eq. (3.8) in such a way that the sign of [(ξ − 1)/ − jinf ] will determine the minimal parameter; a representation of this situation is given in Fig. 2. We can summarize,

fM N  + 1. (3.13) j = 1 + ξξB

The mathematical relationship for j  hides even so a more important fact about this particular system in relation to condensation. Inspecting the above-mentioned equation, N/(1 + ξB /ξ ) corresponds to the number of counterions that neutralize the center disk charge (namely Nn ); or, by using Eq. (2.5), Nn = ξ/ = N − ξB / (i.e., Nn = N − 2ξB / ). Furthermore, we know from previous works in the mean field [9–11,15] that the ratio of condensed ions is the well-known Manning relationship fM . In other words, from Eq. (3.13), j Nn

 =

fMN 1+ξB ξ

1/ 2

 jsup .

(3.12)

N→∞

= fM ,

(3.14)

thus recovering the celebrated Manning condensation fraction in the thermodynamic limit. As a result, j  will be intrinsically related to the number of condensed ions; this will be further clarified in the following sections. A. The partition function

The  and   notations are reserved for the floor and ceiling functions. 3

In order to invert the transform for the partition function Z we need to revise three possible cases according to the set of {aj }. These cases depend on the values of  and ξ . In fact,

062140-4

COUNTERION DENSITY PROFILE AROUND A CHARGED . . .

PHYSICAL REVIEW E 91, 062140 (2015)

aj

aj

j*

j*-1

j

j*

−1 given that the aj = (j − c)2 , with c = ξ  + 12 , the set will be degenerate if given two values k and j , aj = ak . That is the case of either j = k or that 2c ∈ Z. In other words, c could be either a semi or a whole number. In terms of our variables ξ , ξB , , and N it implies that

2 (1 + ξB ) ∈ N  or

(3.15)

2 + 2(N − Nn ) ∈ N.  There are three cases that correspond to different kinds of solutions as follows: (i) The nondegenerate case For this case, the values of aj are nondegenerate, as shown in Fig. 3. Given that we expect a simple separation of the product terms from Eq. (3.11) into sums as (s) = T{f[] N ⊗···⊗f0 }

N 

j

j*

FIG. 3. (Color online) Representation of aj as a function of j with j  the location of the minimum of the {aj } for the nondegenerate case.

 1 1 = Ck , (3.16) s + a s + ak k k=0 k=0 N

j*

FIG. 4. (Color online) Representation of aj as a function of j with j  the location of the minimum of the {aj } for the degenerate case with 2(1 + ξB )/ an even number.

(ii) The degenerate case: 2(1 + ξB )/ even The even degeneracy case is represented graphically in Fig. 4 where the value of j  is degenerate and the number of degenerate aj ’s equates the minimum between (j  − 1) and (N − j  + 1). If {j † } denotes the set of degenerate values, then the partition function reads from Eq. (3.11),5 T{f[] (s) N ⊗···⊗f0 } ⎡ ⎤⎡ ⎤   1 1 ⎦⎣ ⎦ =⎣ s + ak (s + ak )2 k ∈{j / †} k∈{j † } ⎧ ⎫ N  ⎨  1 ⎬ 1 = ⎩ al − ak ⎭ s + ak k ∈{j / † } l=0,l=k ⎧ ⎫ N  ⎨  1 ⎬ + ⎩ al − ak ⎭ † † k∈{j }

⎡

l=0,l=k,k

⎛

N 

⎞

⎤

1 1 ⎠ 1 ⎦ −⎝ a − ak s + ak (s + ak )2 l l=0,l=k,k †

  1 1 1 † = , Ck + Ck,k† − S k,k 2 s + ak s + ak (s ) + a k † † k ∈{j / } k∈{j } ×⎣

 1 4 with coefficients Ck = N l=0,l=k al −ak simplified to  −1 1    −1 1   + 2  k − 2 ξ + k − ξ 2(−1)k Ck =   ξ −1 1  2 , k!(N − k)!  N +1+k−2  + 2

(3.17) which inverted gives for the partition function,   N 2π R 2 N aN −E0  e Ck e−ak  . Z=  k=0

(3.18)

(3.20)  with Ck,k† = N l=0,l=k,k † coefficients simplified to

If the size of the box is large such that  → ∞, then the partition function scales as Z

4

∼

→∞

Cj  e−aj   .

Ck,k†

1 al −ak

and Sk,k† =

N

1 l=0,l=k,k † al −ak

2 † (−1)k+k k − k † =− k!(k † )!(N − k)!(N − k † )!

(3.21)

(3.19)

A key remark on the notation used to avoid confusion between the coupling constant  and the Gamma function (x).

5 The  symbol represents a product or sum where only one of the two degenerate indeces (k or k † such that ak = ak† ) is taken into account.

062140-5

´ JUAN PABLO MALLARINO AND GABRIEL TELLEZ

PHYSICAL REVIEW E 91, 062140 (2015)

argument is to use the nondegenerate partition function from Eq. (3.18) for all future calculations.

aj

B. Density profile

The density profile is given by  N ρ(r) = d N r δ(r − r1 ) e−H N! Z  2 N  R N δ(R ey/ − R ey1 / ) = d N yd N θ N! Z  2π R ey/ × e−H+  2

j

j*

=

FIG. 5. (Color online) Representation of aj as a function of j with j  the location of the minimum of the {aj } for the degenerate case with 2(1 + ξB )/ an odd number.

and6 Sk,k† =

=

#

2 1 + 0 (N − k + 1) − 0 (N − k † + 1) † k − k k − k† $ + 0 (k † + 1) − 0 (k + 1) . (3.22)

Ultimately, the partition function  Z=

2π R 2  

+

N

⎛

eaN −E0 ⎝



k ∈{j / †}

Ck,k† [ − Sk,k† ] e−ak  ⎠.

(3.23)

k∈{j † }

(iii) The degenerate case: 2(1 + ξB )/ odd Different from the previous case, j  is not degenerate as shown in Fig. 5 and so we would expect that the partition function gives ⎛  2 N  2π R Z= eaN −E0 ⎝Dj  + Ck e−ak   k ∈{j / † }∧j  ⎞  + Ck,k† [ − Sk,k† ] e−ak  ⎠, (3.24)

 Dk =

1 k!(N − k)!

 ρy , 2π R 2 e2y/

(3.26)

(3.27)

N δ(y − y1 ){yj } = δ(y − y1 )T{yj } + δ(y − y2 )T{yj } + · · · + δ(y − yN )T{yj } ,

(3.28)

where the T stands for a truncation of the integration to the subregion of the phase space delimited by the [BO] leading to

(3.25)

The partition function is given by three different expressions, (3.18), (3.23), and (3.24), depending on whether 4(1 + ξB )/ is a whole number. Since the free energy has to be continuous, the partition function requires the same property. Despite that fact, it is interesting to verify that we can recover (3.23) and (3.24) from an appropriate limit of (3.18). This is done in Appendix A. The rationale proceeding to this

6

 N 2π R 2 N  e−E0 2 2y/ 2π R e  N! Z 

× d N y δ(y − y1 ) e−H

which speaks equally for ρ. The average for ρy is easily calculated using the procedure for the partition function considering that in the partitioning of the phase space δ(y − y1 ) takes (N − 1)! permutations with y1 at a given position in the group of {yj }’s. Since each arrangement is obtained through a series of permutations of the [BO] (y1 < y2 < · · · < yN denoted by [BO]), then the average yields that

2 .

yj

ρy : = N δ(y − y1 ){yj } with   2π R 2 N −E0 1 δ(y − y1 ){yj } : = e  N! Z 

× d N y δ(y − y1 ) e−H ,

k∈{j † }

with

j

where

Ck e−ak  ⎞



δ(y − yk )T{yj } N −E0   e 2π R 2

= d N y δ(y − yk )e−H  Z [BO] N aN −E0  e 2π R 2 =  Z  N × d N y δ(y − yk )e− j =0 aj (yj +1 −yj ) .  =

[BO]

2π R 2  

×

eaN −E0 Z

d N y δ(y − yk ) %

We have used the 0 ( ) notation for the Digamma function.

N

[BO]

&'

N  j =0

recalled as δkT

062140-6

fj (yj +1 − yj ) . (3.29) (

COUNTERION DENSITY PROFILE AROUND A CHARGED . . .

The N ! steps in to account for all permutations of any given arrangement of the {yj }’s deriving from the (N − 1)! intrinsic permutations times N from the average of the density. Let us evaluate δkT term using the Laplace transformation. From the steps followed in Eq. (3.10), #   yN  ∞ [y] T{δT } (s) = dye−s y dyN dyN−1 . . . k

0

⎛



y2

×

dy1 ⎝

0



N 

×

⎛ y2

dy1 ⎝

0

N 

[y]

T{δT } (s) = k

⎭

e−(s+aj ) s + aj − am

k−1 N  

Ck,N;j C0,k−1;m

e−aj  , s − aj + am

Ck,N;j C0,k−1;l e−aj −(al −aj )y N −aj  j =0 Cj e .- k−1 . - N −aj (−y) −aj y j =k Ck,N;j e j =0 C0,k−1;j e = . N −aj  j =0 Cj e =

[y] k

⊗ g1 ⊗ g0 ](),

dj,k (s) =

C0,k−1;j Ck,N;m

δ(y − yk )T{yj } N k−1

T{δT } (s) = [fN ⊗ fN−1 ⊗ · · · ⊗ fk ⊗ gk−1 ⊗ gk−2 ⊗ · · ·

 N  1 1 l=0,l=j −aj +al l=k −aj −s+al ,  k−1  N  1 1 l=0 −aj +s+al l=k,l=j −aj +al ,

k−1  N 

which is invertible noting that the first term is proportional to (y − ), the Heaviside step function,7 hence trivial since y  . Finally, the truncated density reads

l=0

which, again, admits a solution via the Laplace transformation on ; therefore, ⎞⎛ ⎞ ⎛ k−1 N   1 1 ⎠⎝ ⎠ ⎝ T [] [y] (t) = {T T } t + s + aj t + aj {δk } j =k j =0 ⎛ ⎞ k−1 N   d (s) d (s) j,k j,k ⎠, =⎝ + t + s + a t + a j j j =k j =0

j =k

l=0

(3.32) The previous formulation for the density teaches us many things about the behavior of the density close to R and D. First, inspecting the above relationship, we find that δ(y − yk )T{yj } ∝

Z[,ξ,k − 1, − y/] × Z[,ξ,N − k,y/] ; Z[,ξ,N,] (3.33)

the truncated profile at a given position equates the product of partition functions corresponding to k − 1 particles before that position and N − k particles beyond the latter. This trait is characteristic of decorrelated fluids as presumed in the meanfield regime. Summarizing, if we define ρ (r) = 2π R 2 ρ(r)/(N ξ ), the density profile is given by

if j < k if j  k.

The presentation of the previous constants is not convenient for the final inversion. Expanding the second and first products in, respectively, the lower and upper bounds of j we obtain +  , ⎧ 1 ⎨C0,k−1;j − N , if j < k C k,N;m m=k s+aj −am + , dj,k (s) = k−1 1 ⎩Ck,N;j if j  k, m=0 C0,k−1;m s−aj +am ,

ρ (r)

 2 R r  D −aj  .-k−1  r −aj  . N N  j =k Ck,N;j r j =0 C0,k−1;j R × .  D −aj  N k=1 j =0 Cj R

 = Nξ

(3.30)

(3.34)

n

1 with the constant Cm,n;k = l=m,l=k al −a (m,n ∈ N  m < k n ∧ k ∈ [m,n]), evaluated equally as Ck in Eq. (3.17), given by

Cm,n;k =

dj,k;s e−aj 

j =k

j =k m=0

⎞) * k−1   −s(y −y ) fj yj +1 − yj ⎠ e l+1 l .

) k−1

N 

j =0 m=k

Defining gj (x; s) = e−s x fj (x),

with

dj,k;s e−(s+aj ) +

j =0

=−



j =0

k−1 

+

dyN−1 . . . 0



fj (yj +1 − yj )⎠δ(y − yk )

The first inversion yields

⎫ ⎬

yN

dyN 0

⎞

j =0





=

0

0

PHYSICAL REVIEW E 91, 062140 (2015)

2(−1)k−m (k − m)!(n − k)!   −1 1    −1 1  k − ξ + 2  m + k − 2 ξ + ×   ξ −1 1  2 .  n+1+k−2  + 2

From the previous relationship we can extract the leading behavior of the weakly coupled regime, which turns out to be precisely that of mean field; ρ ∝ 1 − 2x/μ with x the perpendicular distance from the disk and μ is the Gouy length.

7

#

(3.31) 062140-7

(x) =

0 1

x0 . x>0

´ JUAN PABLO MALLARINO AND GABRIEL TELLEZ

PHYSICAL REVIEW E 91, 062140 (2015)

On the other hand, as  becomes large, the functional form of the profile simplifies close to both boundaries. Knowing that the first term (that which sums from k to N ) will be nonzero, or at least relevant, if k  j  compared to the partition function, which scales as e−aj   from Eq. (3.19). Conversely, for the outer shell, or close to D, the situation is such that the second term is nonzero if k > j  . This provides an interesting parallel rationale to condensation.

δ(y − yk )T{yj }

Those particles which contribute to the profile near the surface, i.e., the condensed population, are j  . On the other hand, the remaining N − j  contribute to the profile at the outer shell. Approximating from Eq. (3.32) and looking at the profile near R we are examining the truncated profiles for the case of k  j  since all profiles of index beyond j  will correspond to evaporated counterions. Hence,

⎡ k−1 k−1  C0,k−1;j −(aj −aj  )y ⎣ C0,k−1;j Ck,N;j  ±1 −(aj −aj  ±1 )y  e + e  C C0,N;j  j =0 0,k−1;j j =0 ⎤ k−1 Cj  ±1  C0,k−1;j −(aj −aj  )y ⎦ −(aj  ±1 −aj  ) − e + O(e−2 ). e Cj  j =0 C0,k−1;j 

(3.35)

Conversely, the contribution to the exterior shell reads, for k > j  , δ(y − yk )T{yj }

⎡ N N  Ck,N;j −(aj −aj  )(−y) ⎣ C0,k−1;j  ±1 Ck,N;j −(aj −aj  ±1 )(−y)  e + e Ck,N;j  C0,N;j  j =k j =k ⎤ N Cj  ±1  Ck,N;j −(aj −aj  )(−y) ⎦ −(aj  ±1 −aj  ) − e + O(e−y− ), e Cj  j =k Ck,N;j 

(3.36)

sharing the same functional form at both edges. Notice that the density, according to the procedure we have followed, is written as an expansion of powers of the radial distance. Now we turn to Fig. 6, where the comparison between the analytic prediction and the simulation results is displayed. From the results, the deviations from the data, as expected, increase with higher coupling (see ξ = 9 in the plot). Notice that despite the differences the contact density as the profile approaches to r = R matches that reached by the analytic profile. C. Integrated charge

The integrated charge amounts to the number of condensed counterions from the density profile Eq. (3.32) . The integrated charge reads, in terms of the centrifugal variables {yj =  ln(rj /R)}, as follows:  y N  (3.37) Q(r) = ρy (y ) dy = Qk (y), 0

k=1

which, split into contributions of each truncated density, reads  k−1 C0,k−1;l . N −aj  −(al −aj )y ) j =k Ck,N;j e l=0 al −aj (1 − e , Qk (r) = N −aj  j =0 Cj e with particular emphasis on a functional form suited for distances close to R, or  e−(al −a−j )(−y) Ck,N;l . k−1 - −aj   Cj − C0,k−1;j N j =0 e l=k al −aj Qk (r) = 1 − , N −aj  C e j =0 j

corresponding to that closest to D. Notice how the two functional forms consistently show that when y → 0 (close to R) Qk → 0 and y →  (close to D) Qk → 1 as expected since the truncated densities correspond to the contribution of a single particle in the [BO]. In order to understand how many ions condense, let us look at Eq. (3.38) and notice that when  is large compared to y, Qk → 0 if k > j  and so the condensed counter-ions

(3.38)

(3.39)

will equate j  , an argument consistent with what we have assumed constructing the density profile for infinite box sizes. With regards to condensation, it means that the fraction of  condensed ions is f = jN , with special attention to a minor, yet relevant, detail of the density profile. For simplicity we have omitted the degeneracy issue here considering that the pressure is a continuous function. Therefore, the density profiles should not exhibit any particular behavior at the troublesome values.

062140-8

COUNTERION DENSITY PROFILE AROUND A CHARGED . . .

ξ=3, Δ=20 ξ=6, Δ=20 ξ=9, Δ=20 2 ξ=3, Δ=10 2 ξ=6, Δ=10 2 ξ=9, Δ=10

0.8

ρ

0.6 0.4 0.2 0 1

1.1

1.2

1.3

1.4

1.5

r/R

FIG. 6. (Color online) The density profile ρ near the charged disk for different values of the Manning parameter for N = 10. The data display the results for different box sizes. The dotted and dashed curves represent the analytic prediction Eq. (3.34). The values for the Manning parameter are chosen for small  and very close to unity where the theory is no longer valid.

As a matter of fact, there is one problem when aj  = aj  −1 , as discussed earlier, where the profile behaves as ρ ∼ 1/r 2 and the corresponding integrated quantity Qj  (r)  y/(). This form tells us that the condensing counterion is evenly distributed in the box (using y coordinates). An important rationale then for condensation is that this particle is neither condensed nor evaporated and, thus, the integrated charge will exhibit a constant slope form when this occurs. This situation does not take place in the three-dimensional case [9,10,15]. With this matter clarified, the fraction of condensed counterions corresponds then to 0 /   1 fM N 1 ξ −1 = f = . (3.40) N  N 1 + ξξB The numerical results from Monte Carlo simulations are presented in Fig. 7, highlighting the evaporated counterions as indicated by the plateau. We observe that, at the transitions, the integrated charge displays a constant slope form associated to the counterion which lies at the borderline of condensation

N − Q(r)

3

PHYSICAL REVIEW E 91, 062140 (2015)

(in Fig. 7,  = 1/2, 1). Away from the transitions the plateau is flat, similarly to the three-dimensional case. The onset for condensation in two dimensions coincides with the well-known result ξ = 1. This signals a transition between the regime of full evaporated to partially condensed counterions. However, the details of the behavior beyond this point are unique to the two-dimensional problem. Unlike the situation in three dimensions, the thermodynamic limit (N → ∞), at a fixed Manning parameter and no charge in the exterior boundary (ξB = 0), takes the coupling to zero, which is the mean-field regime investigated thoroughly in previous works. In such a case, f → fM = 1 − ξ1 . The situation maintaining a fixed, nonvanishing, coupling constant  and N → ∞ will imply an infinite Manning parameter ξ → ∞ and, consequently, a strongly bound set of counterions. However, despite the situation, the number of condensed counterions increases with N while the remaining evaporated remains all the same. This is easy to see from Eq. (3.13) where the number of evaporated ions is

1 + ξB Nevap = . (3.41)  In other words, evaporation, looking at ξ > 1, is dominated by the coupling parameter. This is the reason for the division of the two-dimensional case among  < 1,  = 1, and  > 1; assuming that ξB = 0 if  < 1, there is evaporation,  > 1 full condensation, leaving  = 1 the critical value for the coupling parameter where there is only one unbound counterion. From the earlier considerations, the problem with  = 1 will have, exactly, one free ion while the problem for greater couplings is equivalent to full condensation. This helps to find an analytic route towards the profile for large couplings and the contact theorem. One question that arises with the previous result is in regard to the effective charge of the center disk and the ion cloud. For simplicity we assume that Q1 > 0. The condensed ions will screen the charge of the center disk, thus reducing the mean potential interaction between the center region and the ions in the outer region of the disk. In other words, what is the value of ξeff := ξ − Nc ?

1 + ξB ξeff =  − ξB , (3.42)  which tells us that ξeff  1 for any coupling and exterior charge Q2 ; also, if   1, then ξeff  0. Last, if ξB = 0, then ξeff  0. However, if  > 1, then the effective charge could be negative only if Q2 = 0. This striking effect, similar to charge inversion, is possible only at strong couplings.

2

1

D. Contact density

0 0

10

20

30

40

50

ln(r/R)

FIG. 7. (Color online) The integrated charge N − Q(r) as a function of the logarithmic distance for  = 102 , ξB = 0, and N = 10 for various . The plots read for the coupling parameter, from top to 10 1 10 2 10 , 2 , 19 , 3 , 11 ,1, 109 , and 2. bottom,  = 25 , 21

For the contact density, we turn to Eq. (3.32) at y = 0 and y = . In both cases, we encounter that the truncated densities have vanishing values except for two particular cases. For y = 0, ) N −aj  j =1 C1,N;j e  N k=1 −aj  T  C e δ(y − yk ){yj } y=0 = , (3.43) j =0 j 0 k>0

062140-9

´ JUAN PABLO MALLARINO AND GABRIEL TELLEZ

(a)

PHYSICAL REVIEW E 91, 062140 (2015)

(b)

1

1.5

Δ=1/2 Δ=2 Δ=5 Δ=10

0.6 0.4

Δ=1/2 Δ=2 Δ=5 Δ=10

0.2 0

0.1

1

ξ 2 e2Δ ρ|r=D

ρ|r=R

0.8 1

0.5

0

2

0.1

Ξ

1

2 3 45

Ξ

FIG. 8. (Color online) The density ρ  at contact in r = R (a) and r = D (b) as a function of the coupling parameter  with ξB = 0; here  varies and N = 10. The range of  extends below the onset for condensation up to high couplings. The dashed curves represent the exact contact density taken from evaluating the density from Eq. (3.26). We considered the analytic exact result since  is small.

likewise for y = , δ(y − yk )T{yj }

 

y=

=

and

⎧ ⎨0 ⎩

N−1

−aj  j =0 C0,N−1;j e N −aj  j =0 Cj e

k=N

1 (aN − aj  ) Nξ



  1 1 1 j j − , = − fM − + fM − ξ N N ξ N

(D) = e2 ρ

k 1, coinciding with the onset for condensation discussed earlier. As expected, when all ions condense, or j  /N → 1, the density at the exterior shell vanishes. The values compare very well with the simulation data from Figs. 8(a) and 8(b) [Figs. 9(a) and 9(b) in logarithmic scale] considering a small . Notice that the contact densities are continuous functions; for r = R the function does not present any appreciable or qualitative changes while its exterior counterpart shows a succession of bumps coming from the evaporated counterions. Seen that increasing  bolsters the transitions we plot the contacts at  = 102 in Fig. 10. The prediction gives a very accurate estimate of the contact densities seen in the figures at all ranges. Observe that the data

-7

0.1

1

2

ξ 2 e2Δ ρ|r=D

m 

10 10

Ξ

-1

-2

10 10

0

Δ=1/2 Δ=2 Δ=5 Δ=10

-3

-4

0.1

1 Ξ

FIG. 9. (Color online) Same as in Fig. 8 but with logarithmic scaling in the y axis.

062140-10

2 3 45

COUNTERION DENSITY PROFILE AROUND A CHARGED . . .

PHYSICAL REVIEW E 91, 062140 (2015)

1 0.8

N=10 0.6 0.5

0.4

E/Δ

ρ|r=R

N=10

1

ξ 2 e2Δ ρ|r=D

1

0

0.2 0

0.1

1

Δ=20 2 Δ=10

0.5

-0.5 0.1

0

1

Ξ FIG. 10. (Color online) The density ρ  at contact in r = R (squares, left vertical axis scale) and r = D (up triangles, right vertical axis scale) as a function of the Manning parameter; here  = 100 and N = 10. The range of ξ extends below the onset for condensation up to high couplings. The dashed curves represent the  → ∞ formulation from Eqs. (3.45) and (3.46). The arrows indicate the location of  = 1, the borderline to strong coupling.

FIG. 11. (Color online) The energy as a function of the Manning parameter; here N = 10 and  = 20,102 . The dashed curve corresponds to the  → ∞ prediction from Eq. (3.47). The dotted curve corresponds to MF’s prediction, or N → ∞, likely for  → ∞.

Then μ(2D) evap

displayed consider values of the Manning parameter which fall out of the scope of the present course; at  = 1 ( = 2) we have the drastic change of behavior to strong coupling, as indicated in Figs. 8–10.

   N − j N + 1 − j = N ξj   N N   

 N +1−j −1 N − j − 1 − ξ N N = 2,

The energy for the system can be found with the derivative of the partition function from Eq. (3.19) with respect to ξ . Then  = − ξ ∂ (ln Z) , E N ∂ξ

(3.47)

This equation for the energy has a saw-like shape, as was acknowledged in Refs. [9,10] and explored via the aborementioned procedure for  < 1 in Ref. [11], due to the transitions at temperatures below the critical temperature (ξ > 1) when N/ξ becomes a whole number. Figure 11 presents the numerical results for two values of  displaying, as anticipated, the effect of the box size to the transitions. Notice that for  > 1 ( > 2) the energy goes below zero due to the interaction among the counterions. Related to the problem of condensation and minimal free energy, the energy shift at any transition is given for the change of energy when an unbound ion is condensed.8

IV. THE  = 1 CASE

The case when  = 1, or  = 2, represents the borderline before full condensation and also could be understood as the limiting situation before the strong-coupling regime. This situation admits an exact formulation that is interesting to explore. To begin with, let us look at the Boltzmann factor of the Hamiltonian Eq. (2.3) as it reads,   N      rj − rk 2  rj −2ξ −H −E0   e =e   (4.1)  R  . R j =1 1j 1, then γ (,ξ,μk ) > 0 and we notice that each term in the sum (5.25) is positive. Therefore ∂ 2 Q(r) > 0 and never vanishes: The curve Q(r) versus ln r (∂ ln(r/R))2 does not have an inflexion point.

The prefactor is ρ(N−1) =

1 π R 2 Z(,ξ,N,∞) N 

×

 μκ (N) with μ1 =(N−1)

cμ(N) ()2  i mi !

(ξ − μk − 1)−1 .

(5.21) 2. Case Q 2 > 0: Unbinding of one counterion

k=2

The above sum involves only partitions μ = (μ1 ,μ2 , . . . ,μN ) such that μ1 = (N − 1). Consider the partition μ˜ = (μ2 , . . . ,μN ), it is a partition with N − 1 elements with |μ| ˜ = (N − 1)(N − 2)/2. One can think of μ as being composed by two partitions, one with one element and the other one with N − 1 elements: μ = ((N − 1),μ). ˜ Using a factorization property of the coefficients of the partitions (). Therefore, shown in Ref. [23] one has cμ(N) () = cμ(N−1) ˜ the numerator of (5.21) is the partition function of a system with N − 1 particles. Then ρ(r) ∼

r→∞

1 Z(,ξ,N − 1,∞) −2 r . π R 2 Z(,ξ,N,∞)

(5.22)

Since we are considering only the case  > 1, this means that the density decays faster that r −2 when r → ∞. This behavior differs from the prediction of the Poisson-Boltzmann equation: the mean-field regime does not apply even at large distances from the inner disk. This can be contrasted to the situation for a three-dimensional system of a charged cylinder, where at large distances the density profile behaves as the mean-field prediction [15], a fact that can be understood as a sign that the coupling is independent of the density in two dimensions. The decay of the density faster than r −2 is also an indication that all counterions condense into the inner charge disk. Indeed, the integrated charge (5.12) has a limit for  → ∞, Q(N) c (r) =

(N−1) 

Zn∞ (1 − (r/R)−2(ξ −n−1) )

with N  cμ(N) ()2  1  (ξ − μk − 1)−1 , Z(,ξ,N,∞) (N) m ! i i k=1 μκ with n∈μ

(5.24) which has the limit Q(N) c (r) → N , when r → ∞: All counterions are condensed. The notation Q(N) c with the subscript c (condensed) and the superscript (N ) has been chosen to recall that this is the integrated charge of a system with N condensed ions. Another indication of the complete condensation of counterions when  > 1 and Q2 = 0 can be noticed in the absence of an inflexion point in the curve of Q(r) as a function of ln r. Indeed, returning to the case D < ∞, we have N  cμ(N) ()2  ∂ 2 Q(r) 1  = γ (,ξ,μk ) [∂ ln(r/R)]2 Z(,ξ,N,) i mi ! k=1 (N) μκ

×

ξ −1 ξ −1  Nc < + 1 that is Nc =  



 ξ −1 , (5.26) 

(5.23)

n=0

Zn∞ =

If the outer shell is charged Q2 > 0, then one can choose independently the charge of the inner disk  and the number of counterions N , these are no longer restricted by the relation ξ = N . The global electroneutrality now reads N  − ξ = Q2 /q. Starting from the situation of the previous section ξ = N, we will see that if ξ decreases keeping fixed N , some counterions will progressively start to unbind from the inner disk. The existence of unbound ions can be traced back to a divergence in the partition function when D → ∞. In the expression (5.6) for the partition function we notice that Z will diverge if at least one of the functions γ (,ξ,μk ) diverges when  → ∞. This occurs if it exists at least one μk such that μk + 1  ξ . By construction, the largest possible value for μk is (N − 1). Therefore unbinding of ions will occur as soon −1 as (N − 1) + 1  ξ , i.e., N  ξ  + 1. We can use this argument as a basis for the definition of the number of condensed ions Nc . Let us define Nc as the number of ions such that the partition function of Nc ions converges, lim→∞ Z(,ξ,Nc ,) < ∞, but the partition function with Nc + 1 ions diverges, lim→∞ Z(,ξ,Nc + 1,) = ∞. Then, from the previous analysis, Nc is the integer such that

N  4(ξ − μ − 1)2 (r/R)2(μ+1−ξ ) . (5.25) 1 − e−2(ξ −μ +1) =1

where x is the ceiling function. Notice that, in general, the number of condensed ions Nc and unbound ones N − Nc differ from the number of charges at the inner and outer disk,   Q1 ξ ξ −1 and = , (5.27) Nc =  q    Q2 ξ −1 ξ N − Nc = N − and = N − . (5.28)  q  For a fixed number of counterions N , as ξ is decreased from the value N  (situation where Q2 = 0), the number of condensed counterions will decrease, in a piecewise fashion, due to the ceiling function. While N  + 1 −  < ξ  N  all counterions remain condensed Nc = N , the partition function of the system has the same form as in the previous section given by Eq. (5.19) in the limit  → ∞. The first evaporation of one ion will occur when ξ = N  + 1 − , and then Nc = N − 1. Let us consider this case in some detail. If ξ = (N − 1) + 1, then for some partitions such that its largest member is μ1 = (N − 1), we have μ1 = ξ − 1. Then in the partition function (5.6), the corresponding function γ (,ξ,μ1 ) = 2 diverges when  → ∞. Thus the partition function will

062140-16

COUNTERION DENSITY PROFILE AROUND A CHARGED . . .

diverge when  → ∞. Its leading order is

PHYSICAL REVIEW E 91, 062140 (2015)

ions plus an additional contribution from the evaporated ion,

Z(,ξ,N,) ∼ Z∞ (,ξ,N )

(r) + Q(r) = Q(N−1) c

→∞



= 2

μκ (N) with μ1 =(N−1)

N cμ(N) ()2 



i

mi !

→∞

(ξ − μk − 1)−1 .

(5.29)

k=2

Using the same argument that leads to (5.22), we recognize in the above expression the partition function of a system with N − 1 particles, Z∞ (,ξ,N) = 2Z(,ξ,N − 1,∞).

− ln[Z∞ (,ξ,N)] = − ln(2) − ln[Z(,ξ,N − 1,∞)]. (5.31) From this, we see that the free energy of the system has a dominant contribution, −kB T ln(2), due to the free energy from the unbound ion, plus a subdominant contribution {− ln[Z(,ξ,N − 1,∞)], finite as  → ∞} from the N − 1 condensed ions. More precisely, adding the contributions from the charged boundaries, the total free energy (in units of kB T , F = − ln Z) of the system is F (,ξ,N,) = − ln(2) + ( − 2) + F (,ξ,N − 1,∞) R +  ln − ln(π R 2 ) + O(−1 ), (5.32) L where F (,ξ,N − 1,∞) is the free energy of a system with N − 1 = Nc condensed particles when  → ∞. A similar analysis can be done for the density profile by separating the contributions of partitions with μ1 = (N − 1) ∞ which are dominant. We have ρ(r) ∼ ρ(N) (r) with

Q(r) = N − 1 +

×

N 

 μκ (N) with μ1 =(N−1)

ln(r/R) + O[e−2 ln(r/R) ], ln(D/R)

(5.36)

where we recalled the fact that  = ln(D/R). From (5.36), we see that when Q(r) is plotted as a function of ln(r/R) it should be a function which varies fast (in the log scale) from 0 to N − 1 and then linearly to the value N . This behavior can be observed in Fig. 13 for the case  = 1. In this limiting situation when (ξ − 1)/ is an integer (equal to N − 1), the unbound ion is in fact “floating” between the inner and outer disks. The curve Q(r) versus ln(r) does not exhibit yet an inflexion point but rather a linear tendency. Suppose now that we decrease ξ below the previous value [(N − 1) + 1] when one ion has unbound from the inner disk, (N − 2) + 1 < ξ < (N − 1) + 1. We can repeat the previous analysis, separating the partitions for which μ1 = (N − 1) from the rest. For those partitions, we have 2((N−1)+1−ξ ) γ (,ξ,μ1 ) ∼ e(N−1)+1−ξ → ∞ when  → ∞, while for any other element of the partition μk (with k = 1), γ (,ξ,μk ) → (ξ − μk − 1)−1 < ∞. Using this and following similar steps to the one that lead to (5.30) we obtain − ln[Z(,ξ,N,)]

→∞

2 1 ∞ (r) = ρ(N) π R 2 Z∞ (,ξ,N )

(5.35)

where Q(N−1) (r) is the integrated charge of a system with N − c 1 ions which are all condensed (with  = ∞), and it is given by an expression similar to (5.23) but with the replacement of N by N − 1 particles. In particular, we notice from (5.23) that if r  R, Q(N−1) (r) = N − 1 + O((r/R)−2 ). Therefore, when c r  R,

(5.30)

Then

ln(r/R) , 

∼ − ln[Z(,ξ,N − 1,∞)] − 2[(N − 1) + 1 − ξ ]

→∞

cμ(N) ()2  i mi !

+ ln[(N − 1) + 1 − ξ ].

(5.37)

The free energy is

(ξ − μk − 1)−1

F∞ (,ξ,N ) = ( − 2) + F (,ξ,N − 1,∞) +  ln

k=2

3

4 N  (r/R)2(μ −ξ ) (r/R)−2 × + . (ξ − μ − 1)−1 2 =2

R L

− ln(π R 2 ) + ln[(N − 1) + 1 − ξ ]. (5.33)

(5.38)

Introducing the density ρc(N−1) (r) for a system with N − 1 particles (all condensed), but with the same values of  and ξ and  → ∞, we notice that

From this expression we see once again that the leading contribution to the free energy is given by the unbound ion, which here contributes with a term ( − 2), plus some subleading terms from the N − 1 condensed ions. Similarly as above, the density profile appears as a sum of a contribution from the N − 1 condensed ions and the unbound ion, (N − 1) + 1 − ξ + r ,2[(N−1)−ξ ] , ρ(r) = ρc(N−1) (r) + π D2 D (5.39)

∞ (r) = ρc(N−1) (r) + ρ(N)

1 . 2π r 2 

(5.34)

If r  D, then the last term vanishes as  → ∞, and the density of the system with N is the same as the one with N − 1 particles. This is an explicit indication that one counterion has evaporated. The contribution to the density from this evaporated ion is the term 1/(2π r 2 ), which vanishes in the limit  → ∞. From the above expression we deduce that the integrated charge has also one contribution from the N − 1 condensed

when  → ∞. Close to the inner disk, r  D, the second term is negligible and ρ(r)  ρc(N−1) (r). The second term becomes important only close to the outer disk when r → D.

062140-17

´ JUAN PABLO MALLARINO AND GABRIEL TELLEZ

PHYSICAL REVIEW E 91, 062140 (2015)

In the limit  → ∞, the integrated charge is

because μk + 1  (N − Nu ) + 1 = Nc  + 1 > ξ . Furthermore, by the nature of the squeezing operations one has that the sum

Q(r) = Q(N−1) (r) + e−2[(N−1)+1−ξ )(−ln(r/R)] c + O(e−2[(N−1)+1−ξ ] ).

(5.40)

Nu 

When R  r  D, we have

k=1

Q(r) = N − 1 + O{e−2[(N−2)+1−ξ ] ln } r R

+ O{e−2[(N−1)+1−ξ ](−ln R ) }, r

μk =

(5.41)

and as r → D [i.e., ln(r/R) → ], the integrated charges approaches the value N , exponentially fast in ln scale, i.e., as e−2[(N−1)+1−ξ ][−ln(r/R)] .

Nu 

(N − ) =

=1

(N − 1 + Nc )Nu  2

is fixed and equal to the the same sum for the root partition. Then Nu 

γ (,ξ,μk ) ∼ eNu [(N−1)+Nc +2(1−ξ )] →∞

k=1

Nu 

1 . μ +1−ξ k=1 k

3. Case Q 2 > 0: Unbinding of many counterions

(5.45)

In this subsection we study the general case where many counterions unbind from the inner disk. Suppose that the charge of the inner disk ξ is such that (N − Nu − 1) + 1 < ξ < (N − Nu ) + 1, with Nu an integer, Nu  2, which is the number of unbound counterions. Indeed, if we recall that the −1 number of condensed ions is Nc =  ξ   < N − 1, then we have Nu = N − Nc . To put in evidence the counterion condensation from the analytical expressions for the partition function, the density, and integrated charge profiles, it is convenient to recall some properties of the coeficients cμ(N) of the expansion (5.2) of the power of the Vandermonde determinant. The partitions μ present in the expansion are dominated by the root partition κ (N) , that is, μ can be obtained from κ (N) by “squeezing” operations: μj → μj + n and μk → μk − n with j < k and n > 0. Consider that we divide the root partition into two parts κ (N) = (κ˜ (Nu ) ,κ (Nc ) ), where κ (Nc ) = ((Nc − 1),(Nc − 2), . . . ,0) and κ˜ (Nu ) = ((N − 1),(N − 2), . . . ,(N − Nu )). The latter can be a partition of Nu elements κ (Nu ) = ((Nu − 1),(Nu − 2), . . . ,0), with all its parts shifted by Nc : κ˜ k(Nu ) = Nc  + κk(Nu ) . Suppose that squeezing operations are performed on κ (Nc ) to obtain a partition μ(Nc ) and, separately, squeezings operations are performed on κ (Nu ) to obtain μ(Nu ) . We now define the shifted u) u) = Nc + μ(N and consider the composite partition μ˜ (N k k (Nu ) (Nc ) partition μ = (μ˜ ,μ ) of N = Nu + Nc particles. Then it is shown in Refs. [23,26] that the corresponding coefficients of the expansion of the power  of the Vandermonde determinant of these partitions satisfy the factorization relation c((N) () = cμ(N(Nuu)) () cμ(N(Ncc)) () . μ˜ (Nu ) ,μ(Nc ) )

(5.42)

In the analysis of the previous sections we used a special case (N) of this factorization property where Nu = 1: c((N−1),μ (N−1) ) = (1) (N−1) cμ(N−1) c(0)

(1) c(0)

with = 1. With the aid of the factorization property (5.42) we will be able to factorize the leading order, when  → ∞, of the partition function of N particles into partitions functions of Nu and Nc particles. Indeed, consider the contribution to the partition function (5.6) from partitions μ = (μ˜ (Nu ) ,μ(Nc ) ) constructed as explained above. For the parts μk of μ with k = 1, . . . ,Nu that belong to μ˜ (Nu ) , we have γ (,ξ,μk ) ∼

2(μk +1−ξ )

e →∞ μk + 1 − ξ

when  → ∞

(5.44)

(5.43)

Notice that the coefficient of the exponential is independent of the partition considered. On the other hand, the contribution from the other parts of the partition are finite when  → ∞ N 

γ (,ξ,μk ) →

→∞

k=Nu +1

N  k=Nu

1 . ξ − μk − 1 +1

(5.46)

Consider the contribution from a partition ν that is not constructed from the mechanism below. This means that at some point a squeezing operation was performed with a part from the first Nu terms and another part from the remaining Nc terms. For example, consider the following partition ν constructed from the following squeezing of the root partition κ (N) : κj(N) → κj(N) − n = νj and κ(N) → κ(N) + n = ν , j  Nu and   Nu + 1.

with (5.47)

Then, if ν < (N − Nu ), the leading-order contribution of this partition, when  → ∞, is Nu 

γ (,ξ,νk ) ∼ eNu [(N−1)+Nc +2(1−ξ −n)] →∞

k=1

Nu 

1 , ν + 1−ξ k=1 k (5.48)

which is a subdominant contribution compared to (5.46). In the case where ν  (N − Nu ) the contribution is γ (,ξ,ν )

Nu 

γ (,ξ,νk )

k=1 Nu    1 1 exp Nu (N − 1) + Nc  →∞ ν + 1 − ξ ν +1−ξ k=1 k .   (5.49) + 2 1 − ξ − ξ − κ(N) + 1  ,

∼

but since ξ − κ(N) + 1 > 0, this contribution is again subdominant compared to (5.46). Thus, to leading order in  only partitions of the form μ = (μ˜ (Nu ) ,μ(Nc ) ) contribute to the

062140-18

COUNTERION DENSITY PROFILE AROUND A CHARGED . . .

partition function. Then, using the factorization property of the coefficients (5.42), the partition function can also be factorized,

PHYSICAL REVIEW E 91, 062140 (2015)

with the contribution from the condensed ions, ρc(Nc ) (r)

Z(,ξ,N,) ∼ Zu (,ξ,Nu ,)Z(,ξ,Nc ,∞) →∞

1 1 = π R 2 Z(,ξ,Nc ,∞)

(5.50)



cμ(N(Ncc)) ()2  (Nc ) ! i mi

μ(Nc ) κ (Nc )

Nc   −1 ξ − μk(Nc ) − 1 ×

with

k=1

×

Zu (,ξ,N,) =e

=1



Nu [(N−1)+Nc +2(1−ξ )]

μ(Nu ) κ (Nu )

×

Nc 

cμ(N(Nuu)) ()2  i∈μ(Nu ) mi !

Nu   (Nu ) −1 μ˜ k + 1 − ξ

(5.56)

and the contribution from the unbound ions, ρu(Nu ) (r) =

1 1 2 ∗ π D Zu (,ξ − Nc ,Nu )

(5.51)

k=1

×

Zu (,ξ,N,) = eNu [(N−1)+Nc +2(1−ξ )] Zu∗ (,ξ − Nc ,Nu )

×

 μ(Nu ) κ (Nu )

cμ(N(Nuu)) ()2  i mi !

Nu   (Nu ) −1 μk + 1 + Nc  − ξ k=1

and Z(,ξ,Nc ,∞) given by (5.19) but for a system of Nc (condensed) particles. The partition function (5.51) gives the u) = contribution from the unbound ions. Recalling that μ˜ (N k (Nu ) Nc + μk , we notice that this contribution (5.51) can be rewritten as

Nu 

 (N )  u (r/D)2 μ +Nc −ξ

=1

u) [μ(N + 1 + Nc  − ξ ]−1 

(5.57)

.

From this, it follows that the integrated charge is (Nu ) c) Q(r) = Q(N c (r) + Qu (r)

(5.58)

with

(5.52)

c) Q(N c (r) = Nc −

with

×

n  cμ(n) ()2   (μk + 1 − z)−1 . Zu∗ (,z,n) = m ! i i (n) k=1

(Nc )

(r/R)2(μ −ξ )  −1 , c) ξ − μ(N −1 

(5.53)

×

μκ



1 Z(,ξ,Nc ,∞)

Nc 

1

k=1

c) μ(N k

Nc 

ξ−

−1

(Nc )

e−2(ξ −μ

μ(Nc ) κ (Nc )

cμ(N(Ncc)) ()2  (Nc ) ! i mi

−1) ln(r/R)

(5.59)

=1

Notice that, formally, Zu∗ (,z,n) = (−1)n Z(,z,n,∞) if Z(,z,n,∞) is analytically continued using the left-hand side of (5.19) when z < (n − 1) + 1. Thus, in a loose sense, Zu∗ (,ξ − Nc ,Nu ) is the partition function of a system of Nu particles and the inner disk with its charge reduced by the Nc charges of the condensed ions. The free energy is

and u) Q(N u (r)

×

F∞ (,ξ,N) = F (,ξ,Nc ,∞) − ln[Zu∗ (,ξ − Nc ,Nu )]

R 2 + Nu  ln + ( − 2) − ln(π R ) . L (5.54)

1 = ∗ Zu (,ξ − Nc ,Nu )

×

Nu 

 μ(Nu ) κ (Nu )

cμ(N(Nuu)) ()2  (Nu ) ! i mi

1

(Nu ) k=1 μk Nu  -

+ 1 + Nc − ξ (Nu )

e−2[μ

+1+Nc −ξ ](−ln

=1 (Nu )

− e−2[μ

+1+Nc −ξ ]

.

r R)

(5.60)

.

We notice that if R  r  D, then −2[ξ −(Nc −1)−1] ln R c) Q(N } c (r) = Nc + O{e

(5.61)

−2[(N−1)+1−ξ ](−ln R ) u) }. Q(N u (r) = O{e

(5.62)

r

Using the factorization of the partition function, one can find that, at leading order in , the density profile is the sum of two contributions, ρ(r) = ρc(Nc ) (r) + ρu(Nu ) (r)

(5.55)

and r

Thus the integrated charge density increases from 0 when r = R to Q(r)  Nc in that intermediate region R  r  D. Then, when r is close to D, u) Q(N u (r) → Nu ,

r→D

062140-19

(5.63)

´ JUAN PABLO MALLARINO AND GABRIEL TELLEZ

PHYSICAL REVIEW E 91, 062140 (2015)

13

exponentially fast in log scale, to finally recover the total number of particles, Q(r) → Nc + Nu = N.

Q(r)

(5.64)

r→D

The case when (ξ − 1)/ is an integer is a special limiting case. The number of condensed counterions is Nc = (ξ − 1)/. The dominant terms in the partition function, when  → ∞, are due to partitions of the form μ = (μ˜ (Nu −1) ,Nc ,μ(Nc ) ), where μ˜ (Nu −1) is a partition of Nu − 1 parts squeezed from the corresponding root partition κ˜ (Nu −1) . The Nu -th part of the partition is fixed μNu = Nc , and μ(Nc ) is, as above, a partition of Nc parts squeezed from the root partition κ (Nc ) . The results (5.50) and (5.52) for the partition function become Z(,ξ,N,) ∼ 2eNu (Nu −1) Zu∗ (,1,Nu − 1)Z(,ξ,Nc ,∞).

→∞

(5.65)

− ln(2).

(5.66)

This ln  correction, which is a contribution coming from the Nu th part of each partition (μNu = Nc ), is the fingerprint of the existence of a “floating” counterion. Indeed, the density profile is now the sum of three contributions, one from Nc condensed ions ρc(Nc ) , one from Nu − 1 unbind counterions ρu(Nu −1) , and an additional contribution from one floating ion, proportional to r −2 , ρ(r) = ρc(Nc ) (r) + ρu(Nu −1) (r) + →∞

1 . 2π r 2 

(5.67)

10 0

ln(r/R) 

50

100

ln(r/R) FIG. 15. (Color online) The integrated charge Q(r) as a function of the radial distance r/R for  = 2,  = 100, N = 13, and different ξ as indicated. The solid curves represent the exact result from Eq. (5.10) as compared to numerical results from Monte Carlo simulations (symbols). VI. STRONG COUPLINGS: THE CASE   1

The strong-coupling regime, indicated for   1 (  2), is characterized by full condensation for ξB = 0, as we have seen throughout the previous sections. We have insisted that the behavior at strong couplings is attributed to the small fluctuations of the counterions around the ground state [27], a perspective of the minimal energy configuration at zero temperature (as shown in Fig. 16) that matches a Wigner crystal like a pebbled necklace√ with a minimum interparticle distance a; notice that a = R 2 − 2 cos θ2 . This has been evidenced in the investigations on strong coupling in Ref. [28] on the two-dimensional case with a charged plate. Then, the ground-state positions are r(0) j = R[cos θj ıˆ + sin θj jˆ],

(6.1)

with 2π (j − 1) . (6.2) N In the T → 0 image, the profile is dominated by the leading order, or the contribution to the energy due to the θj =

Similarly, the integrated charge can be cast as (Nu −1) c) Q(r) = Q(N (r) + c (r) + Qu

12

11

Thus, the free energy acquires an additional ln  correction F∞ (,ξ,N) = F (,ξ,Nc ,∞) − ln[Zu∗ (,1,Nu − 1)]

R + Nu  ln + ( − 2) − ln(π R 2 ) L

ξ=22.8 ξ=23 ξ=23.2

(5.68)

with the integrated charge corresponding to the condensed ions c) Q(N c (r) given by (5.59) and the one corresponding to Nu − 1 u −1) (r) given by (5.60) with the replacement unbind ions Q(N u Nu → Nu − 1. The charge of the condensed ions converges exponentially fast (in log scale) to Nc when r  R [Eq. (5.61)]. The term ln(r/R)/, in log scale, varies linearly from 0 to 1 when r varies from R to D, thus linearly increasing the total integrated charge from Nc (close to r = R) to Nc + 1 (r = D). Close to r = D, the charge corresponding to the Nu − 1 unbind ions varies exponentially fast (in log scale) from 0 to Nu − 1, thus completing at r = D the total charge Q(D) = N . This is illustrated in Fig. 15, for the case  = 2. When ξ = 22.8, the number of condensed ions is Nc = 11, and for ξ = 23.2, Nc = 12. The results from Fig. 15 for these two cases are compatible with the approximation Eq. (5.58). In the case where ξ = 23, there is a “floating” ion because (ξ − 1)/ = 11 is an integer. In that case the integrated charge shown in Fig. 15 follows Eq. (5.68).

FIG. 16. (Color online) The ground state for the cylindrical cell model in two dimensions. The counterions sit on a position with θk = 2π (k − 1)/N .

062140-20

COUNTERION DENSITY PROFILE AROUND A CHARGED . . .

PHYSICAL REVIEW E 91, 062140 (2015)

charged disk alone. Corrections to the leading order arise from the interaction between charges which are attributed to their arrangement at the crystalline positions. In this sense, looking at the ground state in Fig. 16, we can write the shift of energy when an ion moves to a new position from its ground state as follows; from Eq. (2.3), using Eq. (3.1),

This statement provides the starting point for the analysis of the strong coupling in two dimensions that is unique. Furthermore, it redefines the characteristic length of the diffusive layer, which corresponded to the Gouy-Chapman length like R/ξ . The redefined scale, according to the potential energy, is

N 

|rj | − 2 R j =1 8 7   rj − rk  1  × − ln[2 − 2 cos θj k ] , ln R 2 1j