O Bernard

J.-P Simonin

Association of counterions on polyelectrolytes: Thermodynamic properties in the binding mean spherical approximation

We present a model for polyelectrolyte solutions within the binding mean spherical approximation. An approach developed previously, to describe polyelectrolytic chain solutions and on the other hand to describe the association of counterions on spherical polyions is generalized, considering both the polyelectrolytic chain formation and the association of counterions on the chains. Thermodynamic properties deduced from this model are presented. The associative part of the Helmholtz energy is deduced from the thermodynamic perturbation theory. Analytic expressions for the electrostatic contributions to the internal and Helmholtz energies are established.

Introduction

Understanding the properties of polyelectrolyte solutions is essential to biological processes and in many industrial applications. Extensive efforts have been made to develop theories and models for this type of solution. Linear polyelectrolytes are generally seen as a sequence of charged subunits linked together to form a chain. Counterions are also present in solution to electrically neutralize the chains. Because of the electrostatic repulsion between the charges on the chain, polyelectrolytes are often described as charged rods at low concentrations. On the other hand, when their concentration is high, or in the presence of added salts, the polyelectrolytes are deemed to adopt the conformation of random coils. Moreover, in the same way as in colloidal solutions, in some cases, a part of the counterions is expected to be bound to the charged sites on the chains. This association can be of purely electrostatic nature when the density of charge along the chains is sufficiently high. In this case, association is usually called condensation of the counterions. In order to account for this electrostatic association at low concentrations, a polyelectrolyte is often represented as a uniformly charged infinite cylinder. The distribution of counterions in the vicinity of a polyelectrolyte chain can be described using the Poisson-Boltzmann equation [1,[START_REF] Fuoss | Proc. Natl. Acad. Sci. U. S[END_REF] in the framework of the cell model [START_REF] Wennerström | [END_REF]. Additionally, the degree of condensation of counterions can be assessed analytically [4][5][START_REF] Oosawa | Polyelectrolytes[END_REF]. However, specific interactions exist also between charged subunits and counterions. This is particularly the case when charged monomers along the chains have acid-base properties or with certain counterions that can form specific complexes. These interactions are in particular necessary to account for the differences in the degree of association observed for different counterions of similar charges and sizes. Thus, for a long time, the polyelectrolytes have been described as chains of charged monomers bearing sites, in particular to describe the titration of their solutions [START_REF] Harris | [END_REF][8][START_REF] Rice | Polyelectrolyte Solutions[END_REF][START_REF] Zhou | [END_REF][11][12]. On the other hand, when the solutions are sufficiently concentrated, the screening length between the ions can become much smaller than the length of the stretched polyion. As a result, the repulsion between the charged monomers along the chains, becomes weaker and the polyions may adopt a more collapsed conformation. Moreover, for an infinitely long polyion, the association of a counterion on a given monomer is independent of the position of this monomer along the polyion. On the other hand, for small polyions, the electrostatic interaction between a counterion and a given monomer depends on the distance of this monomer to the center of the polyion [13][14][15][16]. Under these conditions, it is desirable to explicitly take into account the various subunits constituting a polyion.

Moreover, an integral equations theory has been developed to describe polymers consisting of hard spheres linked together by sticky points [17,18]. These equations have been solved and used to represent the thermodynamic properties for uncharged chain solutions [19][20][21][22][23][24]. Charged chains have also been described within this framework [25][26][27][28][29][30], using the mean spherical approximation [31,32]. However, the strong interactions between counterions and charged monomers are not properly taken into account in this approximation. Besides, a thermodynamic perturbation theory (TPT) has been developed to describe the thermodynamic properties of polymers [18,21] and polyelectrolytes [33]. Expressions of thermodynamic quantities have been established and used to describe the polyelectrolyte solutions [33][34][35]. Interestingly, the association between sites carried by counterions and monomers was included in these models [35]. However, only interactions between nearest neighbors were taken into account in order to evaluate the associative contribution to the thermodynamic properties. Besides, a previous study of polyelectrolyte solutions has been carried out using integral equation [28]. In this case, when calculating the thermodynamic properties, the interactions between all the subunits of the chain are taken into account. Therefore, in order to evaluate the associative contribution to thermodynamic properties, other interactions than those between the nearest neighbors could be taken into account within the framework of integral equations. We therefore decided to continue this study of polyelectrolyte solutions, within this approach, taking into account the association between the polyion subunits and the counterions. Integral equations are expected to take more systematic account of all interactions.

In the next section, the polyelectrolyte model is presented. In the following section the thermodynamic properties deduced from this model are established. The various contributions to the Helmholtz energy are presented. The contribution of association between neutral species is introduced first. Then, the calculation of the screening parameters is carried out. This allows one to derive the electrostatic contributions to the internal and Helmholtz energies in the last subsection. Finally, a discussion and a conclusion are presented in the last section.

The model of polyelectrolyte

We consider an ionic system with an arbitrary number of components m in a solvent. The solvent is regarded as a continuum of relative permittivity 4 r . The temperature of the system is T, Boltzmann's constant is denoted as k B and we use b = 1/k B T throughout. The ionic components have number densities q i , charge ez i (e is the elementary charge) and hard core diameter s i . As previously [28], we consider that the particles 1, 2,

• • • n with z 1 = z 2 = • • • = z n and s 1 = s 2 = • • • = s n
are the constituents of the polyelectrolyte chain. The first element of the chain (taken arbitrarily as head) will be denoted with index 1. It possesses a bonding site of type S 1 . The other end of the chain (the tail), of index n, will have a bonding site of type S 2 . All elements with indices 2, • • • , n -1 have two different sites S 1 and S 2 and we allow only bonds between sites S 1 and S 2 of components that have neighboring indices. Initially, the bond between two consecutive atoms is described using an association equilibrium between the sites S 1 and S 2 of these two atoms. The set of particles with indices 1, 2, • • • n participate in the bonding to form the polyelectrolyte chains.

In order to consider only the totally formed chain, it is considered that the association constant between these sites is large enough so that the balance between the free sites and the linked sites is completely displaced in favor of the linked sites. The chain thus formed of constituents {1, • • • , n}, will be called a polyion, and it will be denoted by p. There is also counterions c of charge z c , diameter s c and number density q c such that the system is electroneutral. The system may also contain coions of arbitrary sizes and charges. In order to allow association of counterions on the elements of the chains, we also consider that each constituent of the chain bears an additional site of type S 4 which can be bound solely to counterions. Each counterion bears only one site of type S 4 . The system is depicted in Fig. 1.

The pair potential for this model is given by

u ij (r) = u HS ij (r) + u el ij (r) + S k ,S l U (S k S l ) ij (r 12 ) ( 1 )
where the indices i and j denote the ionic species, u HS ij (r) is the hard sphere potential, u el ij (r) is the electrostatic Coulomb potential, and U (S k S l ) ij (r 12 ) is a short-ranged site-site potential responsible for association. r is the vector joining the centers of ions i and j, r 12 is the vector joining the sites S k and S l . The Coulomb potential, u el ij (r), is given by

u el ij (r) = e 2 z i z j er ( 2
)
where e = 4p4 0 4 r with 4 0 the permittivity of vacuum. By averaging the correlation functions over the orientations of each site, the orientation dependence of the site-site potentials is eliminated to obtain the total orientation-averaged pair correlation function, h ij (r).

The short-ranged site-site potential U

(S k S l ) ij
(r) between sites of type S k and S l on ions i and j are defined in terms of the Mayer function,

f (S k S l ) ij (r), as f (S k S l ) ij (r) = exp -bU (S k S l ) ij (r) -1 = B (S k S l ) ij d r -s - ij (3)
where

B (S k S l ) ij
is an element of the matrix of the sticky interactions. The orientation-averaged pair correlation function h ij (r) can be split into orientation-averaged partial contributions h

(ab) ij (r) as, h ij (r) = h (00) ij (r)+ a a (a) i h (a0) ij (r)+ b a (b) j h (0b) ij (r)+ a,b a (a) i a (b) j h (ab) ij (r). (4)
In this relation, the superscripts denote states of bonding for the sets of sites a and b on ions i and j, respectively. The case a = 0 means a non-bonded particle and a (a) i represents the fraction of ion of type i for which the sites of the set a are not bonded. As previously [45], we label the states (a) by using consecutive numbers in the following way: (1) ≡ S 1 , (2) ≡ S 2 , (3) ≡ S 1 S 2 , (4) ≡ S 4 , (5) ≡ S 1 S 4 , (6) ≡ S 2 S 4 , (7) ≡ S 1 S 2 S 4 . The case a = S k means that only the site S k is bonded, a = S k S l means that only the sites S k and S l are bonded, etc. The partial correlation functions h (r) through the following Wertheim-Ornstein-Zernike (WOZ) equation [17,18]

h ij (r) = c ij (r) + k q k dr 1 c ik (|r 1 -r|) a k h kj (r 1 ) ( 5 )
where h ij (r) and c ij (r) are matrices containing the partial correlations h (r). In order to solve this set of equations, different closure relations can be used. In order to obtain analytical results, the associative mean spherical approximation (AMSA), also called binding mean spherical approximation (BiMSA), closure relation has been used. For the pair correlation function matrix, the closure condition is,

h (ab) ij (r) = -d a0 d b0 +(1 -d a0) (1 -d b0 ) t (ab) ij 2ps ij d r -s - ij for r ≤ s ij (6)
where the constants t (ab) ij are proportional to the sticky parameters

B (S k S l ) ij
(defined by Eq. ( 3)) involved in the association of the ions i A site S 1 of a subunit can associate with the site S 2 of the following subunit. The chain is formed when all these bonds are achieved. Each of the sites S 4 on these subunits can associate with a site S 4 carried by a counterion.

and j and to the contact values of the cavity function [25,36]. For the direct correlation function, the closure condition is [37,38],

c (ab) ij (r) = -
be 2 e z i z j r d a0 d b0 for r > s ij [START_REF] Harris | [END_REF] with s ij = (s i + s j )/2 the i-j distance of closest approach and d a0 the Kronecker's delta function, d np = 1 if n = p and d np = 0 otherwise. Eq. ( 7) means that only the c (ab) ij (r) with a = b = 0 has a contribution for r > s ij . This approximation is a generalization of the MSA, used previously to describe simple electrolytes [31,32,39,40].

A formal solution of the OZ equations with these closure conditions has been obtained previously for ionic mixtures in which an arbitrary finite number of bonds is allowed [36,41]. The solution of this model had led to relations between various parameters and the excess internal energy of the studied system. Previously, this approach was used, on the one hand, to describe polyelectrolytic chain solutions without association [28], and on the other hand to describe the association of counterions on spherical polyions [42][43][44][45]. In the present work, we establish a generalization of these studies considering both the polyelectrolytic chain formation and the association of counterions on the subunits of the chains.

Thermodynamic properties of charged chains

Background

Similar to our previous work [28,34,46], the excess Helmholtz energy is expressed as being constituted of four contributions, namely,

D A = D A HS + D A chain + D A ass + D A el (8)
where the first, D A HS is a hard sphere term, the second, D A chain , is a term related to the connectivity in the uncharged chain, the term D A ass is a contribution arising explicitly from the association between the counterions and the subunits of the polyion and D A el is the electrostatic contribution arising explicitly from the charges on the associated chains and counterions. The hard sphere contribution can be computed from the expressions of Boublik or Mansoori et al. [47]. From the solution of Eq. ( 5) for uncharged chains, the chain contribution D A chain can be determined using the compressibility or virial equations [18,[21][22][23]. As an alternative the chain contribution can be deduced from the thermodynamic perturbation theory [17,18,21]. From the thermodynamic relation bD A = b 0 DE db , the electrostatic contribution to the Helmholtz energy DA el , can be deduced by thermodynamic integration of the excess electrostatic energy DE el . The electrostatic energy D E el can be deduced from the solution of the WOZ Eq. ( 5), for the closures defined by Eqs. ( 6)-( 7), using the following expression:

DE el = 2p ij q i q j ∞ 0 dr r 2 u el ij (r) h ij (r). (9)
However, in view of the specific form of the functions h ij (r) given by Eq. ( 4), it is obvious that the electrostatic contributions to these properties are dependent of the model used to express the association between counterions and polyions. In our previous study of association on a spherical polyion [45], all the sites were equivalent. In this case, the description of association was expressed either in terms of species in chemical equilibrium, or in a mean field approximation in terms of the association on the different sites of the polyion. The description of association in terms of different species in chemical equilibrium is more general. However, when the number of sites on each polyion becomes large, the solution of the set of equations describing all the possible chemical equilibrium becomes increasingly tedious. As an alternative, in the mean field approximation, it has been assumed that equations defining the fractions of independent sites and counterions are those defined in the thermodynamic perturbation theory (TPT) [17,18] or in the statistical associating fluid theory (SAFT) [48][49][50]. Now, in the present study, the subunits of the polyions are clearly not equivalent. It is therefore much more complicated to describe all the possible species obtained with many counterions at different position along the polyelectrolyte chain. This is why, in this work, we have again used the mean field approximation deduced previously in order to get a tractable description of this system.

Associative part of the Helmholtz energy

The thermodynamic perturbation theory was developed in order to calculate the free energy of linear polymers from the reference free energy of its dissociated constituents. The statistical associating fluid theory is a generalization of this approach for associating fluids with an arbitrary number of bonding sites. In these theories, the associative part of the Helmholtz energy D A ass is a function of the set of fractions of free sites a i of the counterions and subunits of the polyion. Initially, the subunits are involved both in polymerization equilibria and in association with the counterion. Then, we consider that the chains are completely formed. The subunits are then involved only in equilibria of association with the counterions. Subsequently, we shall consider only the association between the subunits and the counterions:

bD A ass = bD A TPT1 = q c ln(a c ) + 1 2 (1 -a c) + q p n k=1 ln(a k ) + 1 2 (1 -a k ) . ( 10
)
The fractions a k of free sites on ions k is related to the association constant k ck by the following relation

1 = a k + a k q c a c k ck with: 1 ≤ k ≤ n. ( 11
)
In the same way, the fractions a c of free counterions is given by

1 = a c + a c q p n k=1 a k k ck . ( 12
)
Now, in our mean field approximation, the association constant k kc are related to the various interactions by the following exponential relation [45]

k ck = k 0 ck g ck = k 0 ck g HS kc (s ck ) exp(-bW ck ) ( 13
)
where ĝck expresses the departure of the equilibrium constant k ck from that of the corresponding uncharged system, k 0 ck , at infinite dilution. This term is approximated by the product of the pair distribution g HS ck (s ck ) for uncharged hard spheres and of an electrostatic contribution exp( -bW ck ). By differentiation of Eq. ( 10) and utilizing Eqs. ( 11)-( 12), one gets,

∂bD A ass ∂b = -q p n k=1 (1 -a k ) ∂ ln k ck ∂b = -q p n k=1 (1 -a k ) ∂ ln k 0 ck ∂b - ∂bW ck ∂b . ( 14
)
In view of the Eq. ( 13), the partial derivative ∂lnk ck /∂b depends on the derivative of bW ck . As in previous work, this relation will be used in a step of the thermodynamic integration used to obtain the Helmholtz energy from the energy. Then the suitable function W ck allowing this integration will be determined.

Screening parameters used for the electrostatic energy

A formal solution of the WOZ equation, Eq. ( 5), for the class of closures defined by Eqs. ( 7)-( 6) was obtained earlier [41]. The thermodynamic properties, such as the internal energy, the pressure, and the compressibility, can be deduced from this solution. The electrostatic contributions to these quantities are all expressed in terms of the screening parameter, C B , given by the following equation:

C B 2 = p be 2 e 0 k q k X T k a k X k ( 15
)
where X T k stands for the transpose of matrix X k which are obtained by solving the following linear system

k M ik X k = Z i ( 16
)
with

M ik = 1 + C B s k d ik I + p 2D s k s 2 i 1 0 0 0 q k a k -s i t ik q k a k ( 17
)
where the zeros in bold, in a matrix of the right-hand side, indicate that the set of the other elements of the rows and columns are null.

X i = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ X (0) i X (1) i X (2) i . . ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ Z i = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ z i 0 0 . . ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ ( 18
)
where I is the identity matrix and,

D = 1 -pf 3 /6 ( 19
)
f n = k q k s n k . ( 20
)
Moreover, in order to solve the system of Eqs. ( 15)-( 20), one needs to define the matrices a i and t ij . The elements of these two matrices are related to the fractions a S k i of free site S k on element i and to the sticky parameters B (S k S l ) ij respectively. When only two sites S 1 and S 2 are allowed on each particle there are simple expressions for the matrices a i and t ij .

When there are more than two sites per particle, an extension of this approach called ideal network approximation (INA) has been introduced [24,[51][52][53][54]. This approximation considers that association on each site is independent of the association on other sites. We will use this approximation in this study to define the matrices a i and t ij . First, we recall the forms taken by these matrices when only two sites per atom can be present, for polyelectrolytes without association.

We have previously considered a model of a fixed length polymer chain where the individual beads have two different sites S 1 and S 2 , and we allow only bonds between S 1 and S 2 , and only between particles that have contiguous indices.

In the case of only two different sites on each particle, the matrix a i appearing in Eqs. ( 5), ( 15) and ( 17) had the form [20,26],

a i = ⎡ ⎢ ⎢ ⎢ ⎣ 1 a (1) i a (2) i a (3) i a (1) i 0 a (3) i 0 a (2) i a (3) i 0 0 a (3) i 0 0 0 ⎤ ⎥ ⎥ ⎥ ⎦ (21)
in which the matrix element a

(1) i = a S 1
i is the fraction of free site S 1 , a

(2) i = a S 2
i , and a

(3) i is the fraction of ions i which have the sites S 1 and S 2 free at the same time. Upon the simplifying assumption that the sites are independent, the INA leads to use the relation, a

(3) k = a (1) k a (2)
k , which permitted us to represent the matrices a i as a tensor product of smaller matrices [45]

a i = a (1) i ⊗ a (2) i (22)
where the matrix a (a) i are defined by

a (a) i = 1 a Sa i a Sa i 0 ( 23
)
with (a) = ( 1) or ( 2). In the INA, the mathematical form of the matrices a i expressed by Eq. ( 22) can be extended to the case where three sites are present, which leads to the remarkably simple relation,

a i = a (1) i ⊗ a (2) i ⊗ a (4) i ( 24
)
a i = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 1 a (1) i a (2) i a (3) i a (4) i a (5) i a (6) i a (7) i a (1) i 0 a (3) i 0 a (5) i 0 a (7) i 0 a (2) i a (3) i 0 0 a (6) i a (7) i 0 0 a (3) i 0 0 0 a (7) i 0 0 0 a (4) i a (5) i a (6) i a (7) i 0 0 0 0 a (5) i 0 a (7) i 0 0 0 0 0 a (6) i a (7) i 0 0 0 0 0 0 a (7) i 0 0 0 0 0 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (25)
where a

(1)

i , a

(2) i and a

(3) i have the same meaning as before and:

a (4) i = a S 4 i ; a (5) i = a S 1 i a S 4 i ; a (6) i = a S 2 i a S 4 i ; a (7) i = a S 1 i a S 2 i a S 4 i . ( 26
)
In the case of two different sites on each particle, the matrices t ij appearing in Eq. ( 17) had the form [28],

t ik = ⎡ ⎢ ⎢ ⎣ 0 0 0 0 0 0 t 12 ik 0 0 0 0 0 0 0 0 0 ⎤ ⎥ ⎥ ⎦ when k = i + 1 and t ik = ⎡ ⎢ ⎢ ⎣ 0 0 0 0 0 0 0 0 0 t 21 ik 0 0 0 0 0 0 ⎤ ⎥ ⎥ ⎦ when k = i -1. ( 27
)
In other words, our t ij ( or the sticky interaction B ij matrix) had nonzero elements only when j = i ± 1. In the INA, only the elements t 12 i,i+1 = t 21 i+1,i of the matrix t ij have non-zero values. The elements of the matrix t ij with upper label (3) ≡ S 1 ∩S 2 have been neglected. Now, when a third site S 4 is added on the subunits i of the polyion, we need to consider also a matrix t ic which contain only one element t 44 ic . With only one site per particle, the matrix t ik was a 2 × 2 matrix [37,38,42].

In the case of two different sites on each particle, we had to use 4 × 4 matrices. With a third site, we have 8 × 8 matrices. The matrices a i are given by Eq. ( 24) or Eq. ( 25). As before, all the elements of the matrices t ij are null except the terms t 12 i,i+1 = t 21 i+1,i for the constituents i of the chain.In the matrices t ic = t ci between the constituents i of the chain and the counterions c, all the elements are null except the terms t 44 ic at the intersection of the fifth line and of the fifth column.

Using these results, we can deduce general expressions from Eqs. ( 16)- (17). The results are given in Eqs. ( 30)-( 32) further. We note that only the functions X (a) i with the exponents (a) = 0,1,2 and 4 are different from zero. This is due to the choice of the functions t ij used. To solve the expressions given by Eqs. ( 16)-( 17), we have used 8 × 8 matrices a i and t ik , which is tedious. However, we note that the matrix t ik has a great number of rows and columns filled with zeros. This particularity can be used to simplify the problem. First, we define a compact form of the matrix t ik given by:

t ik = ⎡ ⎢ ⎢ ⎣ 0 0 0 0 0 0 t 12 ik 0 0 t 21 ik 0 0 0 0 0 t 44 ik ⎤ ⎥ ⎥ ⎦ ( 28
)
where we have kept the first row and column of the initial matrix t ik , which are related to the unassociated part of the MSA in our system of equations. We recall that: the term t 12 ik = 0 only when k = i + 1, the term t 21 ik = 0 only when k = i -1 and the term t 44 ik = 0 only when i and k are a subunit of the polyion and a counterion. By keeping only the rows and columns which contain elements which can interact with the kept elements of the compacted matrix t ik , the compact form of the matrix a i is obtained:

a i = ⎡ ⎢ ⎢ ⎢ ⎣ 1 a (1) i a (2) i a (4) i a (1) i 0 a (3) i a (5) i a (2) i a (3) i 0 a (6) i a (4) i a (5) i a (6) i 0 ⎤ ⎥ ⎥ ⎥ ⎦ . ( 29
)
We have now a set of 4 × 4 matrices in place of 8 × 8 matrices.

Extension can be made easily to ions with more sites. We have verified by substitution that we obtain the same results with the big and the smaller matrices. So, from Eqs. ( 16)-( 17), we find:

X 0 i = z i -g T s 2 i 1 + C B s i ( 30
)
where now g T is given by

g T = p 2D k q k s k 4 a=0 a (a) k X (a) k (31)
and for (a) > 0

X (a) i = s i 1 + C B s i k q k 4 b=1 t ab ik a (b) k ⎛ ⎝ X 0 k + 4 c =b a (c) k X (c) k ⎞ ⎠ . ( 32
)
Using Eqs. ( 30) and (32) in Eq. ( 15), we find:

C B 2 = p be 2 e ⎡ ⎣ m j=1 q j X 0 j 2 + 2 n k=1 q * c (k) X 0 k X 0 c + n-1 k=1 n-k l=1 q * (k, k + l) X 0 k X 0 k+l + S * C (k, k + l) + n k=2 k-1 l=1 q * (k, k -l) X 0 k X 0 k-l + S * C (k, k -l) ⎤ ⎦ . ( 33
)
The functions q * c (k) are related to the condensation of counterions on the unit k of the polyelectrolyte chain

q * c (k) = q c a c t 44 c,k a k q k s * c + s * k ( 34
)
where the following notation has been used

s * i = s i 1 + C B s i ( 35
)
with i = c or k. The functions q * (k, k ± l) are equal to q * k,k±l defined in the previous article without condensation (Eq. ( 36) of Ref. [28]). The explicit expressions of these functions for l = 1.2 and 3 were given in Eq. ( 36) of this previous article. The terms S C (k, k ± l) used in Eq. ( 33) are given by

S * C (k, k±l) = X 0 c q * c (k, k ± l) q * (k, k ± l) X 0 k±l + q * c (k ± l, k) q * (k, k ± l) X 0 k + q * cc (k, k ± l) q * (k, k ± l) X 0 c ( 36
)
where the ratios of functions q * c (k, k ± l)/q * k,k±l , and q * c (k ± l, k)/q * k,k±l are related to the association of one counterion on subunits k and k±l respectively. The ratio of functions q * cc (k, k ± l)/q * k,k±l , is related to the association of two counterions on subunits k and k ± l. The functions q * c (k, k±l), q * c (k±l, k) and q cc (k, k±l) have mathematical forms similar to those of the functions q * (k, k ± l). The explicit expressions of these functions for l = 1, 2 and 3 are recalled in Appendix A. From these expressions, we deduce by induction that the ratios appearing in Eq. ( 36) are given by

q * c (k, k + l) q * (k, k + l) = q c a c t 44 c,k s * k s * c + S * k+l S * k+l q * c (k + l, k) q * (k, k + l) = q c a c t 44 c,k+l s * k+l s * c + S * k+l S * k+l q * cc (k, k + l) q * (k, k + l) = q c a c t 44 c,k q c a c t 44 k+l,c s * k s * k+l 2 s * c + S * k+l S * k+l ( 37
)
with

S * k+l = l j=0 s * k+j . ( 38
)
In previous work without association, we established links between the q * (k, k ± l) functions and the densities of the various aggregates formed during polymerization of the chain. We considered the limit of large dilution of various electrostatic quantities such as the screening parameter. In this limit, we studied the q(k, k ± l) functions -noted without the star exponent -which are deduced from the q * (k, k ± l) functions considering that the C B parameter is zero in this limit. Subsequently, for finite concentrations, the C B dependence is reintroduced when the link between the q(k, k ± l) functions and the various aggregates has been established. We have imposed that electrostatic quantities satisfy the limiting laws of Debye and Hückel. Moreover, when the constants of association become very large, we have also studied the limit of complete polymerization obtained for chains formed with various numbers of subunits. When the polymerization is incomplete, chains of various lengths are present. The q(k, k ± l) functions must account for this mixture of chains of various lengths. Therefore, without association to recover this limiting law, we have found that q(k, k + 1) must be the sum of the densities of all the polymers containing at least subunits k and k + 1. In the same way, q(k, k + 2) must be related to the set constituted by the trimer formed with k, k + 1 and k + 2, and all the higher polymers (tetramers, pentamers and so on) containing k, k + 1 and k + 2. The same deduction must be made for higher polymers appearing in the equation of [C B ] 2 . Now, when the association between subunits and counterions occurs, in the limit of high dilution, we studied the functions q c (k, k ± l) and q cc (k, k + l)noted without the star exponent -which are also deduced from the functions q * c (k, k ± l) and q * cc (k, k + l), considering again that the C B parameter is zero in this limit. We noted that the functions q c (k, k + l) and q cc (k, k + l) have the same mathematical form as the functions q(k -1, k + l) and q(k -1, k + l + 1), respectively. They must therefore represent the same set of aggregates as these functions, replacing the subunits k -1 and k + l + 1 by counterions. Then the same deduction must be made for the functions q c (k, k + l) and q cc (k, k + l), containing counterions at their ends, as for the functions q(k, k + l). Therefore, when association occurs, q c (k) must be the sum of densities of all the polymers containing c bonded on subunit k; q c (k, k ± l) must be the set of polymers containing at least the subunits comprised between k and k ± l and with a counterion c bonded on subunit k. Finally, q cc (k, k ± l) must be the set of polymers containing at least the subunits comprised between k and k ± l and with counterions c bonded on subunits k and k ± l. Now, to obtain more explicit expressions, we consider only the case where the subunits of the chains and the other ions all have the same diameter s. When the chain is completely formed q(k, k ± l) = q p = q 1 = q 2 = • • • = q n . Then in this case from Eq. (66), we have:

q * (k, k ± l) = q(k, k ± l) (1 + C B s) l = q p (1 + C B s) l = q p y l ( 39
)
where

y = 1 (1 + C B s) . ( 40
)
Now, we consider the functions related to the association between counterions and the subunits of the chain. When the chain is completely formed, the functions q c (k), must be the density of the chains containing the counterion c bonded on subunit k

q c (k) = q p (1 -a k ) = q p āk (41)
where

āk = 1 -a k ( 42
)
is the fraction of bonded site k. In order to obtain similar expressions for the other functions related to the association of counterions on the chains, we used the INA and assumed the simplifying hypothesis of apparent independence of the sites. Consequently, all the functions q c (k, k±l) must be equal to q c (k). Furthermore, all the functions q * cc (k, k ± l), related the interactions between two associated counterions, bound to subunits k and k±l, are proportional to (1a k ) and (1a k±l ). Consider first, the expressions of q * c (k, k ± l), q c (k ± l, k) and q * cc (k, k ± l) when C B tends to zero (i.e., without the star in exponent).

q c (k, k ± l) = q(k, k ± l)(1 -a k ) = q p āk q c (k ± l, k) = q(k, k ± l)(1 -a k±l ) = q p āk±l q cc (k, k ± l) = q(k, k ± l)(1 -a k ) (1 -a k±l ) = q p āk āk±l . (43)
In the above equations, the first equality is a consequence of the hypothesis of independence of the sites. The second equality results from the complete polymerization limit, when all the chains have reached their maximum size.

Let us consider now the expressions of q * c (k), q * c (k, k±l), q c (k ± l,k) and q * cc (k, k ± l), when C B is different from zero (that is to say with the star in exponent)

q * c (k) = q p (1 -a k ) 1 + C B s = q p āk y ( 44
)
q * c (k, k ± l) = q c (k, k ± l) (1 + C B s) l+1 = q p āk y l+1 (45) q * c (k ± l, k) = q c (k ± l, k) (1 + C B s) l+1 = q p āk±l y l+1 (46) q * cc (k, k ± l) = q cc (k, k ± l) (1 + C B s) l+2 = q p āk āk±l y l+2 . ( 47
)
Then from the Eq. ( 33) for (C B ) 2 we deduce:

C B 2 =p be 2 e ⎡ ⎣ m j=1 q j X 0 j 2 + 2 n k=1 q c (k) X 0 k X 0 c y + q p n-1 k=1 n-k l=1 X 0 k X 0 k+l + S C (k, k + l) y l +q p n k=2 k-1 l=1 X 0 k X 0 k-l + S C (k, k -l) y l ⎤ ⎦ ( 48
)
where X 0 k is given by Eq. ( 30) and with

S C (k, k±l) = X 0 c q c (k) q p X 0 k±l y + q c (k ± l) q p X 0 k y + q cc (k, k ± l) q p X 0 c y 2 . ( 49
)
Eq. ( 48) for [C B ] 2 can be simplified by expressing the q c (k) and q cc (k, k ± l) functions in terms of the āk (assuming that q cc (k, k ± l) are simply proportional to the product āk āk±l ) :

C B 2 = p be 2 e ⎡ ⎣ m j=1 q j X 0 j 2 + 2q p n k=1 āk X 0 k X 0 c y + q p n-1 k=1 n-k l=1 X 0 k + āk X 0 c y X 0 k+l + āk+l X 0 c y y l +q p n k=2 k-1 l=1 X 0 k + āk X 0 c y X 0 k-l + āk-l X 0 c y y l ⎤ ⎦ . ( 50
)
In the same way, using again Eqs. ( 30) and (32) in Eq. ( 31), the parameter g T can be explicitly expressed as a function of the charges, diameters and of the C B parameter. The expressions obtained are given in Appendix B.

Internal and Helmholtz energies

The electrostatic contributions to the thermodynamic properties of the present model can be expressed with the parameters defined in the previous section. An expression for the electrostatic contribution to the internal energy DE el has been obtained as a function of the parameters involved in the solution of Eqs. ( 5)-( 6) [36,41,45]:

DE el = e 2 e k q k 1 s k Z T k a k (X k -Z k ) ( 51
)
where Z T stands for the transpose of matrix Z given by Eq. ( 18) and in which the elements of the vectors X k are given by Eqs. ( 30) and (32). We find

DE el = e 2 e ⎡ ⎣ m j=1 q j z j M 0 j + q p 2s n k=1 āk z c X 0 k + z k X 0 c y + q p s n-1 k=1 n-k l=1 z k X 0 k+l (l + 1) + āk (z k+l X 0 c + z c X 0 k+l )y (l + 2) + āk āk+l z c X 0 c y 2 (l + 3) y l + q p s n k=2 k-1 l=1 z k X 0 k-l (l + 1) + āk (z k-l X 0 c + z c X 0 k-l )y (l + 2) + āk āk-l z c X 0 c y 2 (l + 3) y l ⎤ ⎦ (52)
with

M 0 k = X 0 k -z k s k . ( 53
)
The energy route is known to be an accurate way to determine the thermodynamic properties of charged systems for MSA-type closures [39]. In order to determine the excess Helmholtz energy, one starts from the thermodynamic relation ∂(bD A)/∂b = DE. By integrating the electrostatic contribution by parts, one gets [41,45] bD

A = bDE el - b 0 db b ∂ ∂b DE el . ( 54
)
First, we remark that bDE el depends on C B and on the fractions of bonded sites āk . Then we deduce

bD A = bDE el - b 0 db b ∂DE el ∂C [ ā] ∂C ∂b - b 0 db b n k=1 ∂DE el ∂q p āk C B ∂q p āk ∂b (55)
where [ ā], in the first integral after the equality, is the set of all the āk 's. By differentiation of Eq. ( 52), we find as usual the relation

pb ∂DE el ∂C B [ ā] = -[C B ] 2 . ( 56
)
Using the simplifying Eq. ( 56) in Eq. ( 55) and integrating by parts, we get:

bD A = bDE el + [C B ] 3 3p -b n k=1 q p āk ∂DE el ∂q p āk C B + b 0 db n k=1 q p āk ∂ ∂b b ∂DE el ∂q p āk C B . ( 57
)
It is found that the derivatives ∂DE el /∂q p āk can be written as the following products

∂DE el ∂q p āk C B = X 0 c V k . ( 58
)
Starting from Eq. ( 52), by differentiation of DE el , one gets for the subunits k for which 1 < k < n,

V k = e 2 es ⎡ ⎣ X 0 k + 2 n-k l=1 X 0 k+l y l (l + 2) + āk+l X 0 c y l+1 (l + 3) + 2 k-1 l=1 X 0 k-l y l (l + 2) + āk-l X 0 c y l+1 (l + 3) ⎤ ⎦ . (59)
In the same way, for the subunits 1 and n, one gets

V 1 = e 2 es X 0 1 + 2 n-1 l=1 X 0 1+l y l (l + 2) + ā1+l X 0 c y l+1 (l + 3) (60)
V n = e 2 es X 0 n + 2 n-1 l=1 X 0 n-l y l (l + 2) + ān-l X 0 c y l+1 (l + 3) . ( 61
)
In order to integrate the last term of the right-hand side of Eq. ( 57), one can relate, as before, the expression to be integrated with the derivative of the Helmholtz energy, bD A ass . As in a previous work [45], we introduce the relation

∂DE el ∂q p āk C B = W ck (62)
between the above derivatives of DE el and the functions W ck , introduced in Eq. ( 13), which makes a link with our mean field approximation. It can be shown that the term in the right-hand side of Eq. ( 14) is equal to the term inside the integral in Eq. (57). Therefore, by utilizing Eq. ( 14) in Eq. ( 57), we deduce

bD A = bDE el + [C B ] 3 3p -bq p k āk V k X 0 c + bD A TPT1 . ( 63
)
It is noticed that this D A, obtained by the latter thermodynamic integration, comprises both an electrostatic part and an associative part. The other contributions (ideal, hard spheres and chains) must be added to obtain the complete expression of A.

Discussion and conclusion

Previously, to describe linear polyelectrolytes in solution, we considered the polyions as being formed of charged subunits bonded together to form chains. The electrostatic properties of these ions in solution were calculated in the framework of the mean spherical approximation. The results obtained with this approximation take into account the fact that the polyions consist of subunits linked together. In order to take into account the association of counterions on polyelectrolytes, we considered that each of the subunits constituting the polyions has an additional attracting site allowing the formation of a bond with a counterion. As before, we considered that the intensity of the attraction between the consecutive subunits is very high and that the chains constituting each of the polyions are totally formed. On the other hand, we left unchanged the equilibria of association between the subunits and the counterions.

Analytical expressions of internal energy and screening parameters have been established. Due to the electrostatic interactions, the association between a counterion and a given site on a polyion depends on the possible association of counterions with the other sites of the polyion. This is particularly obvious in the expression of internal energy Eq. ( 52). This expression includes three types of terms:

• terms that do not directly depend on the quantities āk , • terms proportional to āk , • terms proportional to āk āk±l .

The first terms are those obtained previously for the polyions without association of the counterions. The terms proportional only to āk account for the electrostatic interactions between a counterion c linked to k and all the elements of the chain. Finally, the terms proportional to the products āk āk±l account for the electrostatic interaction between a counterion c linked to k and a counterion c linked to k ± l.

In order to highlight the influence of the association between the counterions and the monomers of the polyion, we calculated the electrostatic energy DE el and its dependence on concentration in two cases:

• without counterion-polyion association ( āk = 0 for 1 ≤ k ≤ n) • with complete association of the counterions on the charged monomers of the polyion ( āk

= 1 for k 1 ≤ k ≤ n).
Experimentally, these two states can be encountered during the titration of a polyelectrolyte whose charged groups have acid-base properties. Depending on the pH, the protons will be free or bound to the acid-base groups of the polyion. Eqs. ( 50) and (75), for C B and g T respectively, were solved iteratively as previously. Knowing these two quantities and all the āk 's, the internal energy can be deduced from Eq. (52). As the energy increases sharply with the amount of polyion, we chose to consider bDE el /q where q is the total number density (monomers + counterions). Practically, we have plotted the difference bDE exc /q defined as bDE exc q = bDE el q -bDE 0 q (64) in which

bDE 0 q = lim q→0 bDE el ( āk = 0) q (65)
is the electrostatic energy due to the bonds between the charged monomers on the chains at infinite dilution. We present in Fig. 2 the result for a polyelectrolyte chain with 32 monomers of equal valency z 1 = 1 and diameters s = 0.857 nm. The counterions have a valency z c = 1 and the same diameter. The Bjerrum length L B = e 2 b/4 is related to the diameter by: L B = 0.833s.

Without association, a strong variation of bDE exc /q with concentration is observed. This variation is due to two contributions: one that is related to the screening of electrostatic interactions of all charged subunits, monomers and counterions (the terms proportional to M j in Eq. ( 52), that already are present for dissociated electrolytes); and another one which is due to the weakening with concentration of interactions between the monomers within each chain (terms proportional to z k X 0 k±l ). When association between the counterions and the monomers is complete, only a very small Fig. 2. Variation of bDE exc /q as a function of log(p q s 3 /6). Continuous line: without association of counterions. Dashed and dotted line: with complete association of counterions on the charged monomers along the chains. variation of bDE exc /q is observed at the highest values of volume fraction (p/6)q s 3 . The value obtained at low concentrations results from the electrostatic energy of interaction between the counterions bound to the monomers and, besides, to the energy between counterions bonded to different monomers along each polyelectrolyte chain. When all the counterions are linked to the monomers, the polyelectrolyte chains are overall neutral. Then, the energy variation is much lower than when the polyelectrolytes are charged.

In order to evaluate the Helmholtz energy by thermodynamic integration, we have chosen to describe association using Wertheim's thermodynamic perturbation theory. Consequently, the association energy W ck between a counterion and a subunit k also depends on the potential association of counterions on the other elements of the chain (see Eqs. (58) and ( 62)). This has the consequence of complicating the evaluation of the fractions of associated counterions, and of the electrostatic energy of these associated counterions. Even if the constants of association between the subunits and the counterions are taken equal initially, due to the formation of links between the subunits along the chain and due to the charges carried by these subunits, the sites are clearly not equivalent and the association is not the same at the center and at the ends of a chain. This complication accounts for electrostatic interactions existing within an associated polyelectrolyte. It is therefore perfectly justified. In short, our model takes into account the interactions between all the constituents of an associated polyelectrolyte and not only the interactions between first neighbors. It was one of the objectives that motivated the development of this study. As an alternative, the Helmholtz free energy can be obtained using exponential approximation and numerical integration of the expression for the internal energy. However, assumptions must be made to define the association energy W ck between a counterion and a subunit k, used in the exponential approximation. The use of the association energies obtained in this study to deduce the Helmholtz free energy using numerical integration should lead to the same result.

Other thermodynamic quantities, such as osmotic and activity coefficients, could be deduced from this model in the future. Moreover, in order to evaluate the Helmholtz energy, more elaborate expressions of the free energy of association D A ass , taking into account the cooperativity when associating different sites [18,56,57], could be used. The functions q * (k, k ± l), intervening in the expression for C B in Eq. (33), are related to the electrostatic interactions between the subunits k and k ± l. For l = 1, 2 and 3, they are given by

q * (k, k + 1) = q k T 12 k,k+1 q k+1 S * k+1 q * (k, k + 2) = q k T 12 k,k+1 q k+1 T 12 k+1,k+2 q k+2 s * k+1 S * k+2 q * (k, k + 3) = q k T 12 k,k+1 q k+1 T 12 k+1,k+2 q k+2 T 12 k+2,k+3 q k+3 s * k+1 s * k+2 S * k+3 . . . ( 66
)
where we have used the following notation to represent terms relate to the fraction of bond segments

T ab ij = a (a) i t ab ij a (b) j (67)
and

S * k+l = l j=0 s * k+j . ( 68
) We recall that s * k = s k /(1 + C B s k ).
Accordingly, when all the subunits have the same diameter s, the functions q

* (k, k + l) are proportional s l /(1 + C B s) l .
The functions q * c (k, k ± l) and q * c (k ± l, k) characteristic of the interactions between the condensed counterions in Eq. ( 36), are given by

q * c (k, k + 1) = q c T 44 c,k q k T 12 k,k+1 q k+1 s * k s * c + S * k+1 q * c (k, k + 2) = q c T 44 c,k q k T 12 k,k+1 q k+1 T 12 k+1,k+2 q k+2 s * k s * k+1 s * c + S * k+2 q * c (k + 1, k) = q k T 12 k,k+1 q k+1 T 44 k+1,c q c s * k+1 s * c + S * k+1 . . . (69)
Accordingly, when the counterions and all the subunits have the same diameter s, the functions q * c (k, k + l) and q

* c (k + l, k) are proportional s l+ 1 /(1 + C B s) l+ 1 .
The functions q * cc (k, k ± l) characteristic of the interactions between two condensed counterions in Eq. ( 36), are given by Accordingly, when the counterions and all the subunits have the same diameter s, the functions q * cc (k, k + l) are proportional s l+ 2 /(1 + C B s) l+ 2 . The functions q * c (k, k±l), q * c (k±l, k) and q cc (k, k±l) have mathematical forms similar to those of the functions q * (k, k ± l). We note that the mathematical expression of the function q * c (k, k + l) can be obtained from that of the function q * (k -1, k + l) by replacing the subunit k -1 by a counterion c. Similarly, the mathematical expression of the function q * c (k + l, k) can be obtained from that of the function q * (k -1, k + l + 1) by replacing the subunit k + l + 1 by a counterion c. Finally, the mathematical expression of the function q * cc (k, k + l) can be obtained from that of the function q * (k -1, k + l + 1) by replacing the subunits k -1 and k + l + 1 by two counterions c. The ratio of functions (q * c (k, k ± l)/q * k,k±l ), and (q * c (k ± l, k)/q * k,k±l ) are related to the association of one counterion on subunits k and k±l respectively. For l = 1 or 2, they are given by q * c (k, k + 1) q * (k, k + 1) = q c T The ratio of functions (q * cc (k, k ± l)/q * k,k±l ), are related to the condensation of two counterions on subunits k and k ± l . For l = 1 or 2, they are given by q * cc (k, k + 1) q * (k, k + 1) = q c T 44 c,k T 30) and (32) in Eq. ( 31), the parameter g T can be expressed as a function of the charges, diameters and of the C B parameter. We make the same substitutions given by Eqs. ( 39)-( 47) for all the q(k • • •). Then we obtain: (75)

q * cc (k, k + 1) = qcT 44 c,k q k T 12 k,k+1 q k+1 T 44 k+1,c qc s * k s * k+1 2 s * c + S * k+1 q * cc (k, k + 2) = qcT 44 c,k q k T 12 k,k+1 q k+1 T
g T = p 2D s ⎡ ⎣ m j=1 q j X 0 j + q p 2
Using again Eq. ( 30), the term X 0 k can be split in a term which depends on z k and C B , and in a term proportional to g T . By grouping together the terms depending on g T , we deduce: (77)

  (a,b) ij (r) can be coupled with the partial direct correlation functions c (a,b) ij

Fig. 1 .

 1 Fig.1. Sketch of chain, consisting of five subunits. The subunit 1 has a bonding site of type S 1 and one of type S 4 . The subunit 5 (the last in this chain), has a site S 2 and a site S 4 . The other subunits have the three types of sites: S 1 , S 2 and S 4 . A site S 1 of a subunit can associate with the site S 2 of the following subunit. The chain is formed when all these bonds are achieved. Each of the sites S 4 on these subunits can associate with a site S 4 carried by a counterion.

  Appendix A. The first clusters intervening in the expression for C B in Eq.(33)

  of the interactions between two condensed counterions. In view of the relation between the function for k±l and those for k ± l + 1, we deduceq * c (k, k + l) q * (k, k + l) = q c T 44 c,k s * k s * c + S * k+l S * k+l q * c (k + l, k) q * (k, k + l) = q c T 44 c,k+l s * k+l s * c + S * k+l S * k+l q * cc (k, k + l) q * (k, k + l) = q c T 44 c,k T 44 k+l,c q c s * k s *k+l of the parameter g T By using Eqs. (