text:
# Three Sampling Formulas

Abstract. Sampling formulas describe probability laws of exchangeable combinatorial structures like partitions and compositions. We give a brief account of two known parametric families of sampling formulas and add a new family to the list.

1 Introduction. By an integer composition of weight $`n`$ and length $`\mathrm{}`$ we shall mean an ordered collection of positive integer parts $`\lambda =(\lambda _1,\mathrm{},\lambda _{\mathrm{}})`$; we write $`\lambda n`$ for $`\lambda _j=n`$. It will be convenient to also use variables $`\mathrm{\Lambda }_k=\lambda _k+\mathrm{}+\lambda _{\mathrm{}},`$ $`k\mathrm{}`$, so that $`\lambda _j=\mathrm{\Lambda }_j\mathrm{\Lambda }_{j1}.`$

A composition structure is a nonnegative function $`q`$ on compositions such that for each $`n`$ the values $`\{q(\lambda ):\lambda n\}`$ comprise a probability distribution, say $`q_n`$, and the $`q_n`$’s satisfy the following sampling consistency condition. Imagine an ordered series of randomly many nonempty boxes filled in randomly with balls, so that the distribution of occupancy numbers from the left to the right is $`q_n`$. The condition requires that if some $`k<n`$ balls are sampled out uniformly at random then the distribution of the reduced occupancy numbers in nonempty boxes (in same order) must be exactly $`q_{nk}`$ (without loss of generality we can take $`k=1`$).

Ignoring the order of boxes yields Kingman’s partition structure  (see  and  for systematic development of the theory of partition structures, their relation to exchangeability and many references). But the relation cannot be uniquely inverted, because for a given partition structure there are many ways to introduce the order in a consistent fashion.

Gnedin  showed that all composition structures can be uniquely represented by a version of the Kingman’s paintbox construction . Let $`U`$ be a paintbox – a random open subset of $`[0,1]`$. With a paintbox we associate an ordered partition of $`[0,1]`$ comprised of the intervals of $`U`$ and of individual elements of the complement $`U^c`$, with the order of blocks induced by the order on reals. Suppose $`n`$ independent uniform random points are sampled from $`[0,1]`$ independently of $`U`$. The sample points group somehow within the partition blocks and we obtain a random composition by writing the nonzero occupation numbers from the left to the right. With probability one there is no tie among the sample points and the consistency for various sample sizes follows from exchangeability in the sample.

From a topological viewpoint, the representation establishes a homeomorphism between the space of extreme composition structures and the compact space of open subsets of $`[0,1]`$ (endowed with a weak topology), and also identifies the generic composition structure with a unique mixture of extremes. Thus already the set of extremes is intrinsically infinite-dimensional, not to say about the mixtures. It is therefore a question of interest to find smaller parametric families which admit a reasonably simple description.

In this note we discuss briefly three such families: one is an ordered modification, due to Donnelly and Joyce , of the ubiquitous Ewens sampling formula (corresponding to $`(0,\theta )`$-partition structure from the Ewens-Pitman two-parametric family ); another one, due to Pitman , is an ordered (symmetric) version of the $`(\alpha ,0)`$-partition structure, and the third composition structure is new. Despite the fact that the new composition structure is, in a sense, constructed from the beta distributions, like the first two, the corresponding partition structure does not fit in the Ewens-Pitman family. All three belong to a large infinite-dimensional family of regenerative compositions introduced and characterised in  and all three are mixed, i.e. generated by genuinely random paintboxes.

2 Ordered ESF. In their encyclopaedical exposition of the multivariate Ewens distribution Ewens and Tavaré presented the ordered version of ESF (see , Eqn. 41.6),

(1)
$$e(\lambda )=\frac{\theta ^{\mathrm{}}n!}{[\theta ]_n}\underset{j=1}{\overset{\mathrm{}}{}}\frac{1}{\mathrm{\Lambda }_j}\theta >0,$$

in connection with a size-biased permutation of the Ewens partition structure (here and forth, $`[]`$ is the Pochhammer factorial). The special case $`\theta =1`$ is well known to combinatorialists as the distribution of cycle lengths in a uniform random permutation, provided the cycles are ordered by increasing of their least elements.

Donnelly and Joyce  observed that the formula also defines a composition structure, i.e. that (1) determines a consistent sequence of random partitions taken together with an intrinsic ordering of classes, based neither on the sizes of classes nor on labeling of ‘balls in boxes’. They argued that the ordered structure is of some significance for biological applications, and proved the following paintbox representation of $`e`$.

Let $`Z_j`$ be independent random variables with beta density

$$\mathrm{d}\omega =\theta (1z)^{\theta 1}\mathrm{d}z.$$

Let $`U_e`$ be the open set complementary to the stick-breaking sequence

$$1\underset{j=1}{\overset{k}{}}(1Z_j)k=1,2,\mathrm{}$$

taken together with the endpoints of $`[0,1]`$. Then rephrasing Theorem 10 from  we have

###### Theorem 1

The composition structure $`e`$ can be derived from the paintbox $`U_e`$.

The proof of this result given in  relied on the twin fact about weak convergence of the size-biased permutation of ESF. Next is a direct argument which offers some more insight and exemplifies the approach taken in this paper.

Proof. Introduce the binomial moments

$`w(n:m)`$ $`=`$ $`\left({\displaystyle \genfrac{}{}{0pt}{}{n}{m}}\right){\displaystyle _0^1}z^m(1z)^{nm}d\omega (z)`$
$`=`$ $`\theta \left({\displaystyle \genfrac{}{}{0pt}{}{n}{m}}\right)\mathrm{B}(m+1,nm+\theta )`$
For $`I`$ the leftmost interval of $`U`$ (adjacent to $`0`$) the size of $`I`$ equals $`Z_1`$. Denoting $`\widehat{e}`$ the composition structure derived from $`U_e`$ we aim to show that $`\widehat{e}=e.`$

The argument is based on two facts. Firstly, suppose $`n`$ uniform points have been sampled from $`[0,1]`$ and $`I`$ occured to contain $`m`$ sample points, then conditionally given $`I`$ the configuration of other $`nm`$ points is as if it were a uniform sample from $`I^c`$. The second fact is that given $`I`$ the set $`U_eI`$ is a scaled distributional copy of $`U_e`$, as it is clear from the definition of the paintbox via stick-breaking.

Composition $`(\lambda _1,\mathrm{},\lambda _{\mathrm{}})`$ can only appear if the interval $`J`$, defined to be the leftmost of the intervals of $`U_e`$ discovered by the sample, contains exactly $`\lambda _1`$ sample points. The chance that $`J`$ coincides with $`I`$ is $`w(n:m)`$ and in this case the composition derived from the piece of $`U_e`$ to the right from $`J`$ must be $`(\lambda _2,\mathrm{},\lambda _{\mathrm{}}).`$ Otherwise $`\lambda `$ can appear only if $`I`$ contains no sample points and all $`n`$ group within $`U_eJ`$ in accord with $`\lambda `$.

Combining these facts we get equation

$$e(\lambda _1,\mathrm{},\lambda _{\mathrm{}})=w(n:\lambda _1)e(\lambda _2,\mathrm{},\lambda _{\mathrm{}})+w(n:0)e(\lambda _1,\mathrm{},\lambda _{\mathrm{}})$$

leading to the recursion

$`e(\lambda _1,\mathrm{},\lambda _{\mathrm{}})={\displaystyle \frac{w(n:\lambda _1)}{1w(n:0)}}e(\lambda _2\mathrm{},\lambda _{\mathrm{}}).`$
which is solved as

(3) $`\widehat{e}(\lambda )={\displaystyle \underset{j=1}{\overset{\mathrm{}}{}}}q(\mathrm{\Lambda }_j:\lambda _j)`$
where

$`q(n:m)`$ $`:=`$ $`{\displaystyle \frac{w(n:m)}{1w(n:0)}}`$
$`=`$ $`{\displaystyle \frac{\theta }{n}}{\displaystyle \frac{n!}{(nm)!}}{\displaystyle \frac{[\theta ]_{nm}}{[\theta ]_n}}.`$
Cancelling common factors we arrive at (1), thus $`\widehat{e}=e.`$ $`\mathrm{}`$

There is a canonical correspondence between composition structures and probability distributions of exchangeable compositions of an infinite set $`\{\underset{¯}{1},\underset{¯}{2},\mathrm{}\}`$ (see ). In terms of the paintbox representation the composition derived from $`U`$ is obtained by sampling infinitely many uniform points and then assigning objects $`\underset{¯}{i}`$ and $`\underset{¯}{j}`$ to distinct classes if the closed interval spanned on the $`i`$th and the $`j`$th sample points has a nonempty intersection with $`U^c`$.

The infinite composition associated with $`e`$, call it $``$, has a simply ordered collection of blocks, and the law of large numbers says that the asymptotic frequencies of the blocks (in a growing sample) coincide with the sizes of stick-breaking residuals, from the left to the right. When we view $``$ from the perspective of restrictions $`_n`$ on finite sets $`\{\underset{¯}{1},\mathrm{},\underset{¯}{n}\}`$’s, the collection of blocks stabilises (with probability one) in the sense that for any $`k`$ no new block appearing in $`_n^{}`$, for $`n^{}>n`$, will interlace with the collection of the first $`k`$ blocks represented in $`_n`$, provided $`n`$ is sufficiently large (a zero-one law). Compositions with this property were called ‘representable’ in  and the class of such compositions generated by a general stick-breaking paintbox was characterised in .

3 PSF. Pitman’s composition structure is given by Eqn. (30) in :

(5)
$$p(\lambda )=\frac{n!\alpha ^{\mathrm{}}}{[\alpha ]_n}\underset{j=1}{\overset{\mathrm{}}{}}\frac{[1\alpha ]_{\lambda _j}}{\lambda _j!}0<\alpha <1.$$

This sampling formula was derived from the following paintbox representation.

###### Theorem 2

The paintbox $`U_p`$ for $`p`$ is the union of excursion intervals of the Bessel bridge of dimension $`22\alpha `$.

Equivalently, the complement $`U_p^c`$ is the set of zeroes of the Bessel bridge on $`[0,1]`$. The case $`\alpha =1/2`$ corresponds to the Brownian bridge.

In fact, $`p`$ is a conditional version of another Pitman’s composition structure $`p^{}`$ derived from the set of zeroes of a Bessel process (which has final meander interval adjacent to the rightpoint of $`[0,1]`$). Pitman obtained a formula for $`p^{}`$ akin to (5) (see , Eqn. (28)) using selfsimilarity of the Bessel process and distribution of the length of meander interval. The relation between the structures is that

$$p(\lambda )=\mathrm{const}(n)p^{}(\lambda ,1)\lambda n.$$

Gnedin and Pitman  give a characterisation of $`p`$ related to the observation that this composition structure is also of the product form (similar to (3)) with

$$q(n:m)=\frac{\left(\genfrac{}{}{0pt}{}{\alpha }{m}\right)\left(\genfrac{}{}{0pt}{}{\alpha }{nm}\right)}{\left(\genfrac{}{}{0pt}{}{\alpha }{n}\right)}$$

For $`\mathrm{}`$ fixed, $`p`$ is a symmetric function of the parts. This reflects in that $`U_p`$ is symmetric, that is has component intervals ‘in random order’ (in  the open sets with this kind of invariance are called ‘exchangeable interval partitions’). Summing $`p(\lambda )`$ over distinct permutations of parts yields a function on integer partitions which is the $`(\alpha ,0)`$-partition structure from the Ewens-Pitman family. It follows that $`p`$ could be obtained from the partition structure by permuting the parts in uniform random order (this is the general device allowing to derive symmetric composition structures and symmetric open sets from their unordered relatives ).

Blocks of the Pitman’s composition $`𝒫`$ on $`\{\underset{¯}{1},\underset{¯}{2},\mathrm{},\}`$ are ordered like the set of rational numbers and a such cannot be labeled by integers consistently with their intrinsic order. This happens each time a composition has infinitely many blocks (almost surely) and is symmetric. A consequence is that the infinite composition $`𝒫`$ has no definite first, second, etc or the last block, in particular the first (hence $`k`$th) block in $`𝒫_n`$ does not stabilise as $`n`$ grows.

4 A new sampling formula. Here is a new composition structure

(6)
$$g(\lambda )=\frac{n!}{[\theta ]_n}\underset{j=1}{\overset{\mathrm{}}{}}\frac{1}{\lambda _jh_\theta (\mathrm{\Lambda }_j)}\theta >0$$

where

$$h_\theta (n)=\underset{k=1}{\overset{n}{}}\frac{1}{\theta +k1}$$

are the generalised harmonic numbers which coincide with the partial sums of the harmonic series when $`\theta =1`$.

To explain the genesis of the formula consider stick-breaking with the general beta density

(7)
$$\mathrm{d}\omega (z)=\mathrm{const}z^{\alpha 1}(1z)^{\theta 1}\alpha ,\theta >0.$$

The resulting paintbox generates a composition structure given by the RHS of (3) with

(8)
$$q(n:m)=\left(\genfrac{}{}{0pt}{}{n}{m}\right)\frac{[\alpha ]_m[\theta ]_{nm}}{[\alpha +\theta ]_n[\theta ]_n}$$

where (8) is obtained like (Three Sampling Formulas ) from the binomial moments of the beta density (7) (to see this just follow the lines in the proof of Theorem 1).

For general $`\alpha `$ and $`\beta `$ the induced composition structure cannot be expressed by a simple product formula, because the denominator has no good factorisation. One notable exception is the ESF appearing when $`\alpha =1`$. Another exception is the case $`\alpha =0`$ giving rise to $`g`$; but this should be interpreted properly because measure $`\omega `$ becomes infinite.

###### Theorem 3

When $`\alpha 0`$ the stick-breaking composition structure directed by the beta density (7) converges to $`g`$.

Proof. Expansions in powers of $`\alpha `$ start with

$$[\alpha ]_m=\alpha (m1)!+\mathrm{},[\alpha +\theta ]_n[\theta ]_n=\alpha h_\theta (n)+\mathrm{}$$

therefore when $`\alpha `$ approaches $`0`$ we get

(9)
$$q(n:m)=\frac{n!}{(nm)!}\frac{[\theta ]_{nm}}{[\theta ]_n}\frac{1}{mh_\theta (n)}.$$

which yields $`g`$ as in (3). $`\mathrm{}`$

Distribution (8) underlying $`g`$ is especially simple for $`\theta =1`$ when it gives a weight proportional to $`m^1`$ to each $`m=1,\mathrm{},n`$.

To determine the paintbox representation for $`g`$ we will extend the classical stick-breaking procedure by embedding the process into continuous time and allowing infinitely many breaks within any time interval. Note that defining a composition structure via the RHS of (3), through the binomial moments of some measure $`\omega `$ and

$$q(n:m)=\frac{w(n:m)}{w(n:1)+\mathrm{}+w(n:n)}$$

we need not require that the measure $`\omega `$ be finite and do need to only impose the condition

$$_0^1zd\omega (z)<\mathrm{}$$

to have all binomial moments finite for $`1mn<\mathrm{}`$.

In particular, our $`g`$ appears when we take improper density

(10)
$$\mathrm{d}\omega (z)=z^1(1z)^{\theta 1}\mathrm{d}z$$

(see  for more examples). For this $`\omega `$ consider a planar Poisson process (PPP) in the infinite strip $`[0,\mathrm{}]\times [0,1]`$ with $`\mathrm{Lebesgue}\times \omega `$ as intensity measure. The PPP has countably many atoms $`(\tau _j,\xi _j)`$ (we adopt the conventional fake labeling of atoms which is not intended to say that $`\tau _j`$ or $`\xi _j`$ is a definite random variable for particular $`j`$), and each location on the abscissa is a concentration point for the set of atoms. Define a pure-jump process with increasing cadlag paths

$$S_t=1\underset{(\tau _j,\xi _j):\tau _jt}{}(1\xi _j)$$

where the product is over all PPP atoms to the left from $`t`$. For any $`t`$ the product converges because $`z\omega (\mathrm{d}z)`$ is a finite measure. The process $`(S_t)`$ is a geometric subordinator: for $`t^{}>t`$ the ratio $`(1S_t^{})/(1S_t)`$ is independent of the partial path on $`[0,t]`$ and has same distribution as $`1S_{t^{}t}`$.

(The reader feeling more comfort with breaking sticks from the right to the left should translate paintbox formulas using involution $`z1z`$ and also mirror the sampling formulas.)

###### Theorem 4

The paintbox $`U_g`$ representing $`g`$ is the complement to the closure of the random set $`\{S_t:t>0\}`$, which is the range of the geometric subordinator.

Proof. Fix $`\lambda n`$ and consider a uniform sample of size $`n`$. The composition $`\lambda `$ appears when for some $`\tau _j`$ the interval $`[0,S_{\tau _j}]`$ contains $`m`$ sample points grouped in one component interval of $`U_g[0,S_{\tau _j}]`$ and the composition on the remaining $`(nm)`$ sample points is $`(\lambda _2,\mathrm{},\lambda _{\mathrm{}})`$. From the properties of uniform distribution and because PPP is ruled by a product measure follows that the composition structure induced by $`U_g`$ is of the product form as in (3) and we only need to justify the formula (9) for $`q`$ which is the distribution of the first part of composition of $`n`$.

To that end, let $`\pi (t)`$ be the probability that some $`m`$ sample points group in one interval of $`U_g[0,S_t]`$ and denote $`ϵ_1,\mathrm{},ϵ_n`$ the increasing order statistics of uniform sample. Considering a small time interval $`[0,\mathrm{d}t]`$ it is not hard to see that $`\pi `$ satisfies the differential equation

$$\pi ^{}=a\pi +b,\pi (0)=0$$

with constant coefficients

$$a=E\omega [ϵ_1,1]=w(n:1)+\mathrm{}+w(n:n)\mathrm{and}b=E\omega [ϵ_m,ϵ_{m+1}]=w(n:m)$$

(with 1 in place of $`ϵ_{m+1}`$ in case $`m=n`$) where $`w(n:m)`$’s are the binomial moments of (10). Solving the equation we obtain $`\varphi (t)=(b/a)(1e^{at})b/a=q(n:m)`$, as $`t\mathrm{}`$ whence $`q(n:m)=b/a`$ and this is (9). $`\mathrm{}`$

The infinite composition $`𝒢`$ associated with $`g`$ has infinitely many blocks, and the set of blocks is order isomorhic to the set of rational numbers. Unlike Pitman’s $`𝒫`$ it is not symmetric, i.e. $`g`$ is sensible to permutation of parts $`\lambda _j`$ when $`\mathrm{}>1`$, and the representing paintbox $`U_g`$ is not an ‘exchangeable interval partition’. A combinatorialist might find natural to view $`g`$ as a function on Young diagrams $`(\mathrm{\Lambda }_1,\mathrm{},\mathrm{\Lambda }_{\mathrm{}})`$ with strictly decreasing parts.

Ignoring the order in $`𝒢`$ yields a novel partition structure. For no $`\theta `$ belongs this partition structure to the Ewens-Pitman two-parameter family, which had covered practically all explicit sampling formulas known to date. The distinction can be seen immediately by comparing the probability of one-class partition, our $`g(n)=q(n:n)`$ given by (8) versus the analogous quantity computed via Eqn. (16) in  (the formulas do not match for $`n>4`$ whatever the values of parameters).

Taking other integer values of $`\alpha `$ in (7) leads to formulas involving products of stereotypic polynomial factors, e.g. for $`\alpha =2`$ we have

$$g_2\left(\lambda \right)=\frac{n!\theta ^{\mathrm{}}\left(1+\theta \right)^{\mathrm{}}}{\left[\theta \right]_n}\underset{j=1}{\overset{\mathrm{}}{}}\frac{\lambda _j+1}{\mathrm{\Lambda }_j+2\theta +1}.$$

The resulting infinite compositions have simply ordered blocks and thus are more in line with $``$.

Department of Mathematics, Utrecht University
Postbus 80010, 3508 TA Utrecht, The Netherlands
gnedin@math.uu.nl


format: text

text:
# Monte-Carlo simulations of star-branched polyelectrolyte micelles

## 1 Introduction

Charged diblock copolymers with a long and fully charged hydrophilic part and a short hydrophobic head are attracting a great attention. In water, their hydrophobic cores tend to aggregate and they form colloidal solutions of star-branched spherical micelles. These physical systems have many applications in surface adhesion, stabilization of colloidal solutions and are also used in drug delivery processes. Due to the presence of long-range Coulomb interactions, simple scaling arguments are generally unsufficient to describe their behavior. Quantitative predictions concerning the concentration profiles of monomers and counterions have been obtained through Self-Consistent-Field (SCF) models Borisov1 ; Borisov2 ; Borisov3 . Numerical simulations can also provide useful information. The Molecular Dynamics (MD) approach has been widely used in the study of neutral star-branched micelles Kremer1  and charged linear polyelectrolytes Kremer2 ; Kremer3 . Recently Jusufi, Likos and Löwen (JLL) Likos1 ; Likos2  applied it to star-branched polyelectrolyte micelles. They studied, at relatively low linear charge densities of the chains, the monomer and counterion profiles of one micelle and the effective interaction between two micelles. They also proposed analytical approximations to be compared with the results of MD simulations.

We focus here our attention on the problem of one isolated micelle, in a different strongly-charged regime, which is relevant for some dilute solutions of (sodium-polystyrene-sulfonate)(NaPSS)– (polyethylene– propylene) (PEP) that have been investigated through neutron scattering in our laboratory Muller1 .

We use the Monte-Carlo method which, through a stochastic approach, should provide the same information as the Molecular Dynamics. Our final goal is to reach real-system sizes like stars with $`f=54`$ arms of $`N=251`$ monomers (27108 particles), corresponding to experimental conditions and compare with our neutron-scattering data, without any adjustable parameter. As a prelude, we first present a systematic investigation of smaller micelles within a more general theoretical framework and compare our numerical Monte-Carlo results with the predictions of scaling theories and analytical approximations.

After the description of our model and Monte-Carlo method (section 2), we study in section 3 the evolution of the structure of a micelle with $`f=27`$ arms of $`N=130`$ monomers as a function of the following dimensionless parameter measuring the strength of the Coulomb interaction:

$$\zeta =|z_i|\xi =|z_i||z_m|l_B/a$$
(1)
where

$$l_B=(e^2/4\pi ϵϵ_0)(1/k_BT)$$
(2)
is the Bjerrum length, $`\xi =|z_m|l_B/a`$ represents the “Manning parameter” Manning  which characterizes the linear charge density of polymers; $`a`$ is the distance between consecutive charged monomers, $`T`$ is the temperature, $`ϵ`$ the dielectric constant of the solvent and $`z_m`$, $`z_i`$ represent the valence of monomers and counterions respectively. Throughout our study $`z=z_m=z_i>0`$ and all monomers are charged.

The Bjerrum length is $`l_B0.71nm`$ in water at room temperature and it can be varied through the temperature and dielectric constant. The monomer length is usually of the order of a fraction of a nanometer. $`\zeta `$ is approximately 3 for sodium polystyrene sulfonate in water. It can be increased by choosing other solvents with lower dielectric constants and$`/`$or by taking polyvalent monomers and counterions ($`z>1`$). Higher values of $`\zeta `$ are thus easily accessible. Although, in our fully charged model small values of $`\zeta `$, much less than 1, are at present experimentally unrealistic, it is nevertheless interesting to extend our theoretical investigations down to that range.

When the thermal energy is much larger than the Cou-lomb interaction $`(\zeta <<1)`$, most counterions are outside the micelle corona. This corresponds to the “unscreened regime” Borisov1 . Above some threshold value $`\zeta _{sph}1/f`$, a spherical condensation of counterions inside and in the vicinity of the corona is expected Borisov1 ; Borisov2 ; Borisov3 ; Luc1 . In this so-called “osmotic regime” Borisov1 ; Borisov2 ; Borisov3 , the spherically condensed counterions lower the interaction between branches. Both unscreened and osmotic regime have been extensively studied through mean-field analytical models Borisov1 ; Borisov2 ; Borisov3 ; Pincus1 . Finally, when the “Manning criterion” Manning  $`\xi _M1/z_i`$, i.e $`\zeta _M1`$, is reached, a “cylindrical condensation” of part of the counterions around the branches is superimposed to the spherical condensation. Scaling and SCF theories Borisov1 ; Borisov2 ; Borisov3 ; PGDG  have predicted a linear extension of the branches with corona radius $`R_c`$ proportional to $`N`$ and monomer density varying roughly as $`1/r^2`$, where $`r`$ represents the radial distance to the center of the micelle. This has justified a “rod-like” picture of the branches. This picture has been confirmed by clear experimental evidence of the stretching of the chains Guenoun . For highly charged polymers, a Poisson-Boltzmann analysis PB  has been proposed to interpret the distribution of counterions around the rod-like branches in terms of a ionic condensation Manning ; Luc2 .

In the three latter regimes, we have compared the results of our simulations to the expectations from approximate analytical theories. The simulation of a real system of (sodium-polystyrene-sulfonate) (NaPSS)– (polyethylene– propylene)(PEP) with $`f=54`$ branches of $`N=251`$ monomers, is described in section 4, with a direct comparison to our neutron-scattering data Muller1 .

## 2 The Model and Monte-Carlo method

### 2.1 The model

We use the simple “freely-jointed chain” (or “pearl necklace”) model. Each monomer is represented by a hard core of diameter $`\sigma _m=a`$. The angle between two successive monomer-monomer segments is free. We consider that the hydrophilic chain is fully ionized in water and each monomer carries a charge $`(z_me)`$. This is relevant for sodium- polystyrene- sulfonate, in our experimental conditions, where the degree of ionization is generally of order 0.95. The counterions are modeled by hard spheres of diameter $`\sigma _i=a/2`$, they carry a charge $`(z_ie)`$. For sodium- polystyrene- sulfonate $`z_i=z_m=1`$.

One isolated micelle and the counterions are treated in the “primitive cell model” i.e. they are confined inside a hard-wall cube of edge $`A`$, surrounding the center of the micelle. The dielectric constant is assumed to be the same inside and outside the cubic cell. This approach is relevant for low concentrations with respect to the overlap concentration $`c^{}`$ of micelles. The weak dependence of our physical results on the size $`A`$ of this cubic box (provided it is larger than the diameter of the micelle) means that a more sophisticated treatment of Coulomb interactions (e.g. periodic boundary conditions and Ewald sums) is not required here.

The hydrophobic core is modeled by a hard sphere of radius $`R_{core}`$ of order 5 to $`20a`$. The hydrophilic parts are attached at fixed points chosen on this sphere. Two distributions of these fixed points have been considered:

* a random distribution
* a regularized distribution with approximate equidistance between first-neighbor pairs.

The results obtained with these different distributions are practically identical and it is thus not useful to perform mean values of Monte-Carlo results over different realizations of random distributions of points on the sphere.

All distances are expressed in units of the monomer length $`a`$ and the strength of the Coulomb interaction compared to $`k_BT`$ is measured by the dimensionless parameter $`\zeta `$ defined in equation (1).

### 2.2 Monte Carlo moves

For a counterion, an elementary move is chosen at random in a cube of edge $`\delta x10a`$ centered around it.

For polymers we use the “pivot algorithm”, which has been shown to satisfy ergodicity pivot . A monomer $`P`$ is chosen at random, and a unit vector $`\widehat{u}`$ is taken randomly in some solid angle $`\delta \mathrm{\Omega }`$ centered about the axis defined by $`P`$ and the next monomer $`(P+1)`$. A rigid rotation of random angle around $`\widehat{u}`$ is applied to the part of the polymer going from $`(P+1)`$ to its end. However, beyond Manning condensation $`(\zeta >1)`$, such a pivot algorithm involving only the monomers becomes inefficient since many counterions are localized around each polymer and the cost in energy to move a group of monomers away from these counterions prohibits any Monte-Carlo move.

To circumvent this problem, we have used the “partially clothed” pivot algorithm proposed by Gordon and Valleau Valleau . Having chosen a pivot $`P`$ and the corresponding end part of the polymer, we look for all counterions at a distance less than $`d`$ of this polymer part (the distance from the polymer is defined as the distance to the closest monomer). Each of these counterions is rotated in block with the part of the polymer, with a probability $`p1`$. The detailed-balance equation is satisfied if, in the Metropolis scheme, we take as the probability of acceptance of the Monte-Carlo move Valleau

$$P_{acc}=min\{1,(1p)^{n_{new}n_{old}}\mathrm{exp}[(U_{new}U_{old})/k_BT]\},$$
(3)
where $`n_{new}`$ and $`n_{old}`$ are the number of counterions at a distance less than $`d`$ after and before the temptative rotation. $`U_{new}`$ and $`U_{old}`$ are the Coulomb energies after and before the trial move.

In practice, we have observed that “clothed pivots” have a reasonable probability to be accepted in the Metro-polis scheme if $`p>0.9`$. Since we did not observe very important differences in the convergence to equilibrium for $`0.9<p1`$, we chose the value $`p=1`$ which makes the algorithm simpler and the calculation shorter. (The choice $`p=1`$ has also been made by Lobaskin and Linse Lobaskin  in similar “cluster moves” of macro-ions surrounded by a shell of small counterions). The maximum ion-polymer distance for clothed pivots has been optimized to $`d=4a`$.

To facilitate stretching and contraction of the polymer, we have also applied some small “clothed block translations” of parts of polymer in the following way. A monomer $`P`$ is chosen at random. The next monomer $`(P+1)`$ is moved from $`\stackrel{}{r}_{(P+1)}`$ to $`\stackrel{}{r}_{P+1}^{}`$ within a solid angle $`\delta \mathrm{\Omega }`$ centered around $`(\stackrel{}{r}_{(P+1)}\stackrel{}{r}_P)`$. The rest of the polymer from P to its end is translated in block by a vector $`(\stackrel{}{r}_{(P+1)}^{}\stackrel{}{r}_{(P+1)})`$. Such moves are also combined with “cluster moves” of counterions at distance less than $`d`$. Compared to pivots they have a higher probability to be accepted and a suitable combination of these both types of cluster-moves accelerates the convergence.

Our main criterion for equilibrium is the convergence of the Coulomb energy to the same stationary value starting from two drastically different initial configurations:

* Fully extended polymers with counterions distributed uniformly in the volume of the cell.
* Contracted star with polymers in the equilibrium configuration of uncharged polymer and counterions distributed uniformly in the volume of the cell.

As a typical example, Fig. 1 illustrates the convergence of the Coulomb energy, starting from these two extreme configurations for a micelle with 27 branches of 130 monomers ($`2\times 3510`$ particles) at $`\zeta =2`$. The equilibrium was reached after $`10^3`$ steps per particle, representing 25 hours of monoprocessor time on a Compaq RS232 parallel computer. The relevant statistical properties (distribution of monomers and counterions, form factors etc…) were measured over $`10^3`$ further steps. A “cluster move” of a polymer part was attempted every 3 iterations (the two others being devoted to single counterion moves). The cluster moves were either pivots (3/4 of them) or translations (1/4 of them) as previously described. During the measure of statistical averages, each counterion was moved approximately $`1300`$ times, about $`80`$ pivots and $`200`$ translations per monomer were accepted.

The number of accepted pivots as a function of the position of the pivoting point $`P`$ is represented in Fig. 2. Close to the core of the micelle, the number of accepted moves is low (of the order of 10 per monomer). Nevertheless, it grows roughly exponentially as a function of the position of monomer $`P`$ and statistics becomes acceptable beyond a few tens of monomers.

The size $`\delta x`$ of the cubic box for elementary counterion displacements was $`10a`$, and the elementary solid angle $`\delta \mathrm{\Omega }`$ for pivot rotations was fixed to $`2\pi 10^2`$ steradians. We have also tried some simulation with $`\delta \mathrm{\Omega }`$ increasing with the distance to the core of the micelle, without any substantial improvement of the rate of convergence to equilibrium.

## 3 General results

### 3.1 The theoretical model studied

We present here some theoretical results concerning the variations of the monomer and counterion profiles as a function of the dimensionless parameter $`\zeta `$ \[cf. equation (1)\] measuring the strength of the Coulomb interaction with respect to $`k_BT`$.

We consider two micelles with numbers of monomers and branches smaller than these corresponding to the experimental system studied in the next section:

$`N=130`$ and $`f=6`$

$`N=130`$ and $`f=27`$

The radius of the core $`R_{core}`$ is chosen to be small compared to the radius of the micelle corona and has thus a little influence on the corona. We take in all cases $`R_{core}=5a`$. The influence of the boundary conditions was investigated by comparing the results obtained with two different sizes of the cubic cell: $`A=320a`$ and $`A=640a`$.

We vary the dimensionless parameter $`\zeta `$ keeping the ratios $`\stackrel{~}{A}=A/a`$, $`\stackrel{~}{R}_{core}=R_{core}/a`$ and $`\stackrel{~}{\sigma }_i=\sigma _i/a`$ constant. We start from a low value $`\zeta _0=0.005`$ and increase it by steps. At each new step $`\zeta _{n+1}`$ we take as initial state of the Monte-Carlo calculation the equilibrium state obtained at the preceding one $`\zeta _n`$. The system is thus progressively “annealed”, which provides a faster convergence to equilibrium at each step. As a typical test for equilibrium we compare in Fig. 1 the energy obtained in that way at $`\zeta =2`$ starting from the equilibrium configuration at $`\zeta =1`$ to the results obtained with far-from-equilibrium arbitrary initial configurations. After 1000 iterations the differences between all curves are of the order of the fluctuations.

We calculate the mean values of the radial distributions of monomers and counterions. The deviations of these distributions from the monomer distribution corresponding to fully extended arms are better emphasized by considering, in terms of $`\stackrel{~}{r}=r/a`$, the dimensionless quantities:

$$\rho _{\{m,i\}}(\stackrel{~}{r})=\frac{4\pi }{f}\stackrel{~}{r}^2a^3c_{\{m,i\}}(\stackrel{~}{r})$$
(4)
where $`c_{\{m,i\}}(\stackrel{~}{r})`$ represent the local volume densities.

For a micelle with fully extended arms, $`\rho _m`$ should be 1 for $`\stackrel{~}{R}_{core}<\stackrel{~}{r}<\stackrel{~}{R}_{core}+N`$ and 0 elsewhere.

### 3.2 Monomer profile as a function of $`\zeta `$

For $`N=130`$, $`f=27`$ and $`\stackrel{~}{A}=320`$, typical monomer and counterion profiles at represented in Figs. 3a, 3b and 3c at weak coupling ($`\zeta =0.02`$), intermediate coupling ($`\zeta =0.2`$) and strong coupling ($`\zeta =3`$) respectively. Corresponding snapshots of the micelle are shown in Fig. 4.

The variations of the radius of gyration $`R_g`$ and of the corona radius $`R_c`$ as a function of $`\zeta `$ are shown in Fig. 5 (circles). The corona radius $`R_c`$ is defined here as the abscissa of the maximum of the free-end distribution, which is approximately gaussian (see Fig. 6). It also roughly coincides with the rightmost inflexion point in the concentration profile $`\rho _m(\stackrel{~}{r})`$.

We have studied the influence of boundary conditions by doubling the size $`A`$ of the cubic cell (crosses). $`R_c`$ and $`R_g`$ are very weakly dependent on boundary conditions. Relative difference of the order of a few percent at most are observed at maximum stretching ($`\zeta `$ of order 1). We conclude that it is not useful to investigate more sophisticated treatments of Coulomb interaction like these corresponding to periodic boundary conditions with Ewald sums.

The influence of the number of arms is seen in Fig. 5 by comparing our results with a small number of arms $`f=6`$ (diamonds) to the results with $`f=27`$ (circles). It is only important at weak coupling.

The Figs. 3, 4 and 5 show clearly, as a function of $`\zeta `$, three different regimes. These three regimes are governed by changes in the behavior of the counterions, which are expected from mean-field theories Borisov1 ; Borisov2 ; Borisov3 ; Manning ; Luc1 .

### 3.3 The “unscreened regime”.

At weak coupling $`\zeta <<1`$, the contribution of the counterions to the entropy dominates the electrostatic energy. The local volume density $`c_i`$ of counterions is roughly constant in the cell: $`c_ifN/A^3`$ and from equation (4) $`\rho _i4\pi N\stackrel{~}{r}^2/\stackrel{~}{A}^3`$. If $`A`$ is larger than $`R_c`$, the fraction of counterions in the micelle corona is small and the electrostatic interactions between branches are weakly screened.

In the unscreened limit, the self-consistent field arguments proposed a long time ago for linear polyelectrolytes by de Gennes et al. PGDG  have been straightforwardly generalized to star-branched polyelectrolytes Borisov1 ; Borisov2 . Let us consider a chain whose first unit is placed at the origin. The $`n`$-th unit is at an average distance $`r(n)`$ and feels two forces:

* In a simple gaussian model, there is an elastic force which tends to contract the chain:

$$F_e=\frac{3k_BT}{a^2}\frac{d^2r}{dn^2}.$$
(5)
* If we assume equal stretching of the $`f`$ chains of a star, i.e. all chain ends are at the same distance $`R_c`$ from the origin, the spherical region at distance smaller than $`r(n)`$ contains a mean charge $`Q=nfz_me`$. Applying Gauss theorem, there is an electrostatic force:

$$F_c=\frac{f(z_me)^2}{4\pi ϵϵ_0}\frac{1}{r(n)^2}.$$
(6)

Writing the equilibrium between the two forces, one obtains the following dimensionless equation:

$$\frac{d^2u}{dx^2}+\frac{1}{3}\frac{x}{u^2}=0$$
(7)
with $`x=n/N`$ and $`u=r/(Nd)`$ where the length scale $`d`$ is:

$$d=a\left(\frac{fz_m^2l_B}{a}\right)^{1/3}$$
(8)
The boundary conditions are $`u(0)=0`$ and $`du/dx=0`$ at $`x=1`$. Taking $`f=1`$ in the previous algebra is nothing else than De Gennes et al. equations PGDG .

The center to end distance of the chains is thus proportional to $`Nd`$:

$$R_cf^{1/3}\zeta ^{1/3}Na$$
(9)
(we take $`z_i=z_m`$).

These simple scaling arguments lead to a stretching of the chains which is proportional to the number $`N`$ of monomers. They have justified a so called “rod-like” picture of the branches.

In Fig. 5, the slope $`d\mathrm{log}R_c/d\mathrm{log}\zeta `$ is smaller than $`1/3`$, the value expected through the previous scaling arguments. The dependence of $`R_c`$ on $`f`$ is also weaker than predicted by equation (9). At very weak coupling, we obtain (see Fig. 5) $`R_c(f=27)/R_c(f=6)1.35`$ compared to the expected scaling: $`(27/6)^{1/3}=1.65`$. The explanation is that when $`\zeta `$ is low enough to ensure that very few counterions are inside the corona, the excluded volume effects (which are not taken into account in the previous scaling arguments) become significant. The value of $`R_c`$ has a lower limit which is given by the scaling law for neutral stars Daoud :

$$R_c^{neutral}af^{1/5}N^{3/5}$$
(10)
This limiting value for $`N=130`$ and $`f=27`$, within our pearl-necklace model is indicated in Fig. 5. Our lowest value for $`R_c`$ shows an inflexion towards this asymptote.

Close to the core, our monomer profiles can be roughly scaled on a universal curve by plotting $`\rho _m(R_cR_{core})/Na`$ as a function of $`(rR_{core})/(R_cR_{core})`$ (i.e. the integral of each curve is normalized to 1) and compared to the solution of Eq. 7 given in Fig. 2 of Ref. Borisov1  (see Fig. 7). Close to the core, our Monte-Carlo results agree qualitatively with the solution of equation (7) which gives a roughly constant profile of $`\rho _m`$ (i.e. $`a^3c_m(\stackrel{~}{r})\stackrel{~}{r}^2)`$. At $`R_c`$, equation (7) leads to some unphysical divergence, which is smoothed in more refined SCF models taking into account the fluctuations of the chain ends Borisov1 ; Borisov3 . Nevertheless a substantial increase of $`\rho _m`$ near $`R_c`$ remains in these SCF models, as well as in our simulation.

### 3.4 The “osmotic regime” with spherical condensation of counterions

Increasing $`\zeta `$ we expect a “spherical condensation” of a part of the counterions inside and in the vicinity of the corona. This occurs when the Coulomb energy of a counterion in the field created by the charged monomers:

$$(z_mNf)z_ie^2/(4\pi ϵϵ_0R_c)$$

at the limit of the corona radius is of the order of $`kT`$, i.e. for

$$\zeta _{sph}\frac{1}{f}$$

(at this level of approximation we take $`R_cNa`$). The condensation of ions near the surface of an impenetrable charged sphere Luc1  has been formulated in an analogous way to the cylindrical condensation of ions around a charged rod studied by Manning Manning . The more complex problem of counterion condensation in a penetrable micelle has been studied through Self-Consistent Field methods Borisov1 ; Borisov2 ; Borisov3 . It can also be formulated in a simpler way through a Poisson-Boltzmann analysis. The Poisson-Boltzmann equation can be solved numerically for a penetrable sphere of radius $`R_c`$ with charge density

$$c_m(r)=fN/(4\pi R_cr^2)$$

inside a spherical cell of volume $`A^3`$. The Fig. 8 compares the Monte-Carlo to the Poisson-Boltzmann results for the osmotic coefficient $`\varphi `$. The osmotic coefficient is defined as the ratio of the osmotic pressure to its value corresponding to a uniform repartition of the ions in the cell. The osmotic pressure is equal to $`k_BT`$ times the local ionic density at the edge of the cell.

There is a clear change of slope of $`\varphi (\zeta )`$ at a point which can be defined as the crossover to the condensed regime. For different numbers of arms and cell sizes, we deduce the crossover to the spherically condensed regime at the following values $`\zeta _{sph}`$:

* $`\zeta _{sph}0.15`$ for $`f=6`$, $`\stackrel{~}{A}=320`$
* $`\zeta _{sph}0.03`$ for $`f=27`$, $`\stackrel{~}{A}=320`$
* $`\zeta _{sph}0.08`$ for $`f=27`$, $`\stackrel{~}{A}=640`$

($`N=130`$ in the three cases). At fixed $`\stackrel{~}{A}`$, $`\zeta _{sph}`$ varies in $`1/f`$.

In Fig. 9, we show that at the thresholds defined above, there is a drastic change of the counterion profile. Below $`\zeta _{sph}`$, the local density $`c_i`$ of counterions is roughly constant $`c_ifN/A^3`$ leading to a quadratic behavior of $`\rho _i`$. Above the threshold, the curves $`\rho _i`$ show an inflexion point at $`\stackrel{~}{r}R_c/a`$. Since in that $`\zeta `$ range $`R_c`$ is roughly constant, the abscissa of this inflexion point remains constant.

For $`\zeta >>1/f`$, most of the counterions are inside the micelle corona (see Fig. 9). This regime has often been called “osmotic” by reference to the osmotic pressure of counterions inside the corona Borisov2 ; Borisov3 .

The range $`1/f<<\zeta <1`$ corresponds to the largest stretching of the star with a flat maximum in the curve $`R_c(\zeta )`$ (see Fig. 5).

### 3.5 The “Manning condensed” regime

For larger $`\zeta `$ values, beyond the “radial condensation” of counterions in the micelle corona, we expect the superposition of a “cylindrical condensation” of counterions around each branch of the micelle. This phenomenon is clearly seen in Fig. 3c and Fig. 4c. The Manning criterion for this “cylindrical condensation” is Manning ; Luc1 :

$$\xi _M=1/|z_i|$$

which corresponds to the threshold value $`\zeta _M=1`$.

The Fig. 10 represents an enlargement of the last decade in Fig. 5. The corona radius $`R_c`$ presents a flat maximum at $`0.5<\zeta <1`$ and decreases at $`\zeta >1`$. At this stage, it is interesting to compare our approach with the Molecular Dynamics study by JLL Likos1 ; Likos2  where $`\zeta `$ has been fixed to 3, but where only one monomer out of three, four or six is charged. The “charge fraction” $`\alpha `$ is thus varied from $`\alpha =1/3`$ to $`\alpha =1/6`$ whereas it is $`\alpha =1`$ in our model. The relevant Manning parameter which characterizes the linear charge density of the polymer is $`\xi =z_ml_B/b`$, where $`b=a/\alpha `$ is the distance between consecutive charged monomers Footnote1 . In Fig. 10 we have plotted some of JLL results, as a function of $`\zeta ^{}=\alpha \zeta `$. The fact that from $`\alpha =1/6`$ to $`\alpha =1/3`$ their corona radius increases indicates that they have not reached the Manning condensed regime at $`\alpha =1/6`$ and $`\alpha =1/4`$. They are just reaching it near $`\alpha =1/3`$.

At $`\zeta >>1/f`$, (near the broad maximum of $`R_c`$ and in the Manning condensed regime), the monomer profiles are only very weakly dependent on $`f`$. In Fig. 11, we compare the monomer profiles for different values of the number of branches $`f=6`$, $`27`$, and $`54`$. The number of monomer per branch is $`N=130`$, the radius of the hydrophobic core is $`\stackrel{~}{R}_{core}=5`$ and the size of the cell is $`\stackrel{~}{A}=320`$. The Fig. 12 represents the variation of the gyration radius as a function of the number of branches $`f`$. It increases only by a few percent when the number of branches vary by a factor of 10. The MD results of JLL (see Tables I and II in Ref. Likos2 ), also show a weak dependence of $`R_c`$ on $`f`$.

The Fig. 13 represents a blow up over the last 80 monomers of one branch in Fig. 4c. At a length scale of the order of 10 monomers, a local “rod-like” picture is relevant. To illustrate the phase space explored by a given branch, we have registered 100 snapshots of the same branch (Fig. 14), at equal time intervals in thermodynamic equilibrium. The phase space explored corresponds roughly to a cone of solid angle $`4\pi /f`$ steradians which shows that each branch occupies a fraction $`1/f`$ of the sphere and rarely “mixes” with neighboring branches. An “urchin-like” picture Muller2  is thus relevant.

The last decade in Fig. 5 and its enlargement in Fig. 10 show that, in our fully charged model ($`\alpha =1`$), beyond the Manning threshold, the corona radius $`R_c`$ decreases in $`\zeta ^{1/3}`$. This is quite a new result. Both this exponent and the independence of $`R_c`$ on $`f`$ can be interpreted with the following simple arguments.

For $`\zeta >1`$, there is a cylindrical condensation of a fraction of counterions along the branches. The model proposed by Manning Manning , considers that a fraction $`(11/\zeta )`$ of the counterions condenses around the branches in such a way that the monomers get an effective charge $`\alpha _{eff}z_m`$ with an effective charge fraction

$$\alpha _{eff}=1/\zeta .$$

The remaining fraction $`1/\zeta `$ of counterions are only spherically condensed and distributed between the branches. We assume that $`f`$ is sufficiently high ($`\zeta >>1/f`$) to neglect the fraction of “free” counterions outside the corona. Taking a spherical distribution in $`1/r^2`$, we obtain the following effective ion density:

$$c_i^{eff}(\stackrel{~}{r})=\frac{1}{\zeta }\frac{fN}{4\pi a^2R_c\stackrel{~}{r}^2}.$$

The corresponding “effective” screening length is

$$\kappa _{eff}^1(\stackrel{~}{r})=[4\pi l_Bz_i^2c_i^{eff}]^{1/2}=a\stackrel{~}{r}\left[f\frac{|z_i|}{|z_m|}\frac{Na}{R_c}\right]^{1/2}.$$

Since we restrict our study to $`|z_i|=|z_m|`$, it is independent of $`\zeta `$. The “effective” screening length roughly keeps the same value it reaches at the Manning condensation threshold $`\zeta =1`$. At a radial position $`r`$, the mean transverse distance between two neighboring branches is roughly $`4a\stackrel{~}{r}f^{1/2}`$. It exceeds the screening length by a factor larger than 4. Hence we can neglect the electrostatic interactions between branches. Since the polymers are relatively extended, the excluded volume interactions between branches are also weak. This is the reason why the monomer profiles are so weakly dependent on $`f`$.

To justify the $`(1/3)`$ exponent, we consider one branch in the familiar electrostatic “blob” picturePGDG . Let us represent a branch as a succession of $`N/m`$ blobs of size $`l`$, each containing $`m`$ monomers. As far as the size of a blob is smaller than the screening length $`\kappa _{eff}^1`$, the electrostatic energy of a blob is

$$E_c=(\alpha _{eff}z_mem)^2/(4\pi ϵϵ_0l).$$

Its elastic energy from entropic origin is

$$E_e=(3/2)k_BTl^2/(ma^2).$$

Minimizing the sum $`E_c+E_e`$ with respect to $`m`$, we obtain:

$$(l/ma)=\left[\frac{4}{3}\frac{(\alpha ^{eff})^2z_m^2l_B}{a}\right]^{1/3}=\left[\frac{4}{3}\frac{|z_m|}{|z_i|}\zeta \right]^{1/3}$$

and the length of the branch is

$$R_c=(N/m)l=Na\left[\frac{4}{3}\frac{|z_m|}{|z_i|}\zeta \right]^{1/3}.$$

Since in the present study $`|z_m|=|z_i|`$, $`R_c`$ is proportional to $`\zeta ^{1/3}`$ Footnote2 .

At high $`\zeta `$, Schiessel and Pincus Pincus2  predicted a collapse of the polymers with the following arguments. When the Coulomb energy dominates over the entropic contribution, each counterion is practically fixed to a monomer and forms an electric dipole. The attraction of dipoles should lead to a transition towards a collapse of the star. Up to $`\zeta =10`$ our results follow a $`\zeta ^{1/3}`$ power law, and we do not observe any precursor of this effect.

The Manning picture allows us to make some approximate mapping of our model at $`\alpha =1`$, $`\zeta =3`$ onto the model of JLL at $`\alpha =1/3`$ and the same $`\zeta =3`$. In our model, at $`\zeta =3`$, due to Manning condensation, $`(11/\zeta )=2/3`$ of the counterions condense on the monomers. The remaining effective charge fraction of monomers is thus $`\alpha _{eff}=1/3`$ and the other one third of the counterions is only spherically condensed inside the corona. This effective model is thus close to the real situation studied by JLL at $`\alpha =1/3`$ and $`\zeta =3`$. The arrow in Fig. 10 shows a very nice agreement of the corona radii corresponding to this mapping. This is a strong evidence of the relevance of Manning’s picture.

Up to now, radial distributions of counterions were only considered. In the Manning regime, it is of prime interest to study the cylindrical distribution of counterions around one branch. Defining the distance $`d`$ of a counterion to the closest branch as the distance to the closest monomer we have measured the mean value of the radial density $`n(r,d)`$, $`n(r,d)\delta r\delta d`$ representing the number of counterions at distance between $`d`$ and $`(d+\delta d)`$ from a branch and at distance between $`r`$ and $`(r+\delta r)`$ from the center of the micelle.

The Fig. 15 represents the density $`n(\stackrel{~}{r},\stackrel{~}{d})`$ of counterions at distance $`d=\stackrel{~}{d}a`$ from the closest branch and at a distance $`\stackrel{~}{r}a`$ from the center of the micelle for three values of $`\stackrel{~}{r}`$: close to the core $`(\stackrel{~}{r}=20)`$, in the middle $`(\stackrel{~}{r}=40)`$ and close to the edge of the corona $`(\stackrel{~}{r}=70)`$. We chose $`\zeta =3`$, well above the Manning threshold, $`N=130`$, $`f=27`$ and $`\stackrel{~}{A}=320`$.

It is interesting to compare these results to those expected from a “rod-like” picture of the arms and a Poisson-Boltzmann approach for the distribution of counterions.

Using the “FIDISOL/CADSOL” package developed at Rechenzentrum, Universität Karlsruhe fidisol  for the numerical solution of non-linear partial differential equations, we have solved the Poisson-Boltzmann equation for a rod in a cone of solid angle $`(4\pi /f)`$. The radius of the rod is chosen to $`0.75a`$, which corresponds to the sum of the hard-core radius of monomer and counter ion. The height of the cone is equal to the radius of a sphere of volume $`A^3`$. The distribution of charges in the rod corresponds to the monomer profile obtained in Fig. 3c through Monte-Carlo simulation. (In practice, we fit this distribution to an approximation by a rational function). Using a two-dimensional $`80\times 120`$ grid, we obtained an accuracy better than $`10^4`$ on almost every point. The calculation took only one minute of monoprocessor time. There is a nice agreement except in the vicinity of the corona radius.

Fig. 16 compares the overall density profile $`\rho _i(\stackrel{~}{r})`$ obtained through the Poisson-Boltzmann approach to the Monte-Carlo results of Fig. 3c.

An interesting quantitative way to characterize the Manning threshold is to plot in log-linear coordinates the number $`\nu _i(\stackrel{~}{d})`$ of counterions at distance less than $`\stackrel{~}{d}`$ from the closest monomer (i.e the integral over $`\stackrel{~}{d^{}}`$ of $`n(\stackrel{~}{r},\stackrel{~}{d^{}})`$ between $`3/2`$ and $`\stackrel{~}{d}`$) Luc2 . Within the Poisson-Boltzmann approximation, $`\nu _i(\stackrel{~}{d})`$ is expected to follow a logarithmic law at the Manning threshold PB . As shown in Fig. 17, there is a change of curvature of the lines between $`\zeta =0.5`$ and $`\zeta =1.5`$. A straight line (logarithmic behavior) is obtained at $`\zeta =1`$.

The Fig. 18 shows that, in the Manning condensed regime, $`\nu _i(\stackrel{~}{d})`$ is practically independent of the boundary conditions: the results obtained after doubling the size of the cubic cell ($`\stackrel{~}{A}=640`$) are practically identical to those obtained with $`\stackrel{~}{A}=320`$. $`\nu _i(\stackrel{~}{d})`$ depends on the number of arms but the threshold value $`\zeta _M1`$ corresponding to Manning condensation is unchanged.

## 4 Simulation of an experimental system

### 4.1 Description of the experimental system

We used asymmetric diblock copolymers made of fully deuterated sodium poly-(styrene-sulfonate) $`(\mathrm{NaPSS}_\mathrm{d})`$ and of (polyethylene-propylene) (PEP) Muller1 . These materials were synthetized by anionic polymerization. The degrees of polymerization were measured and found to be 251 for $`(\mathrm{NaPSS}_\mathrm{d})`$ moiety and 52 for PEP moiety. The aggregation number, i.e the number of branches per micelle, was found to be 54 and the core radius is $`R_{core}4.5nm18a`$. The monomer size for $`(\mathrm{NaPSS}_\mathrm{d})`$ moieties is $`a=0.25nm`$. Taking $`l_B=0.71nm`$ as the Bjerrum length of pure water at room temperature, the dimensionless parameter $`\zeta `$ is evaluated to $`\zeta =l_B/a=2.8`$. The concentration is $`0.25c^{}`$ where $`c^{}`$ represents the overlap concentration. Our cubic cell was chosen such that its volume represents the mean volume occupied by one micelle in the solution i.e. $`\stackrel{~}{A}=560`$.

### 4.2 Monomer and counterion profiles

The density profiles of monomers and counterions are plotted in Fig. 19. The corona radius is $`R_c=174a`$ leading to $`(R_cR_{core})/Na0.62`$ this value has to be compared to that obtained in section 3 with $`N=130`$ and $`f=54`$ at the same $`\zeta `$: $`(R_cR_{core})/Na0.57`$. We agree with the Molecular Dynamics results of JLL Likos1 ; Likos2  on smaller branch lengths showing that the ratio $`R_c/N`$ is not strictly constant but slightly decreases when $`N`$ increases. The monomer and ion profiles are quite similar to those observed in section 3 for $`N=130`$ and $`f=27`$ (compare with Fig. 3c).

### 4.3 Form factor, comparison to experimental results

The monomer-monomer form factor obtained through the Monte-Carlo simulation, with the parameters described in subsection 4.1 is shown in Fig. 20 (dashed line). It is close to the monomer-monomer form factor obtained though neutron scattering Muller1  except in the range $`q\pi /R_{core}`$ where the size of the core has a large influence on the spectrum. Another Monte-Carlo calculation with a slightly smaller core $`R_{core}=12a`$, (full line) was performed. It is in a better agreement with the neutron scattering results.

It is interesting to compare the simulated and experimental results to straightforward analytical approximations at low and high $`q`$.

1. At low $`q`$, ($`q<\pi /R_{core}`$), the form factor can be reasonably approximated by that of a $`1/r^2`$ spherical distribution of monomers from $`r=R_{core}`$ to $`r=R_c`$ and zero otherwise:

$$F^{sph}(q)=\left[\frac{Si(qR_c)Si(qR_{core})}{qR_c}\right]^2$$
(11)
where $`Si(z)`$ represents the SineIntegral function Stegun  The thin dashed and dotted lines represent this function for $`R_{core}=12a`$ and $`R_{core}=18a`$ respectively. Up to $`qa0.1`$ they closely follow the corresponding Monte-Carlo simulations (heavy full and dashed lines). At $`q\pi /R_{core}`$ large oscillations appear in the function $`Si(qR_{core})`$, the first of them being responsible for a dip observed in the simulation. Close to this dip, the form factor is very sensitive to the size of the core. We find that our simulation with $`R_{core}=12`$ is closer to the experimental results than that with $`R_{core}=18`$. The interface between the hydrophobic core and the hydrophilic branches might happen to be not so smooth and as well defined as assumed in our model. For example, some hydrophilic monomers may be trapped inside the outer part of the core, leading to a lower effective core radius.
2. At high $`q`$, the branches can be approximated by rods of length $`L=R_cR_{core}`$ with negligible thickness. The corresponding form factor is:

$$F^{rod}(q)=\frac{1}{f}\left\{\frac{2}{qL}Si(qL)\left[\frac{\mathrm{sin}(qL/2)}{qL/2}\right]^2\right\}$$
(12)
The thin full line represents this function with $`L=R_cR_{core}=156a`$ corresponding to our Monte-Carlo estimate. At large $`q`$, $`F^{rod}(q)[1/fq(RcR_{core})]`$ and the position of this line of slope (-1) in logarithmic coordinates gives a nice estimate of the corona radius. In conclusion the corona radius obtained through our Monte-Carlo simulation is in agreement with the experimental value.

The peak observed in the experimental data at $`qa0.02`$ arises from long-range correlations between micelles, mostly due to the remaining uncondensed charges. It is a collective effect which, of course, is absent in our model. The concentration is $`c=5\times 10^3g/cm^3`$. The corresponding $`q`$ wave vector characterizing the inter-micelle correlations is $`q2\pi [c𝒩/fM]^{1/3}`$ where $`M=5.7\times 10^4`$ is the molar mass of one diblock copolymer and $`𝒩`$ the Avogadro number. We obtain $`qa0.016`$ in reasonable agreement with the position of the peak.

## 5 Conclusion

Monte-Carlo simulations appear to be a useful tool for understanding the complex structure of highly charged copolymer urchin-like micelles. As the Coulomb coupling is increased, three different regimes can be distinguished. i) At low coupling, the freely-jointed polyelectrolyte branches become stretched due to the monomer-monomer repulsion and the micelle swells. The monovalent counterions are not sensitive to the low micellar charge and are uniformly distributed inside the cell. ii) At intermediate coupling, the arms are almost fully stretched and most of the counterions are trapped inside or in close vicinity of the micelle. This spherical condensation is due to the presence of many, although weakly charged, branches in the same micelle. iii) At large coupling. a classical ionic condensation appears locally around each highly charged branch which adds to the previous one. The effective post-condensation monomer charge decreases, the monomer-monomer repulsion weakens and the micelle deswells. These Monte-Carlo results agree with Poisson-Boltzmann analyses in different geometries and nicely reproduce small-angle neutron scattering data. The influence of intrinsic chain rigidity, and the role of ion-ion correlation in presence of multivalent counterions will be investigated in future papers.