Solvation force between tethered polyelectrolyte layers. A density functional approach

We use a version of the density functional theory to study the solvation force between two plates modified with a tethered layer of chains. The chains are built of tangentially jointed charged spherical segments. The plates are immersed in an electrolyte solution that involves cations, anions and solvent molecules. The latter molecules are modelled as hard spheres. We study the dependence of the solvation force and the structure of chains and of solute molecules on the grafting density, length of chains, architecture of the chains and on concentration of the solute.


Introduction
Layers of charged polymers tethered onto solid surfaces have attracted great attention for the recent years due to their biological importance and numerous practical applications [1]. The behavior of charged polymers at interfaces is not only more complex as compared to charged polymers in bulk solution [2][3][4][5] but is also more complex than the behavior of systems involving uncharged tethered polymers [6][7][8]. Generally, these systems contain polymers carrying positive and/or negative charges, surface charges, mobile ions in solution and solvent molecules. All these components make the properties of a system intricately dependent on electrostatic interactions [9][10][11][12][13].
A special type of charged polymers are polyampholytes, i.e., polymers consisting of positively and negatively charged segments and, in some cases, also neutral segments. Differently charged, as well as neutral segments, can be distributed randomly or in a certain sequence [14][15][16][17]. A number of synthetic and natural (e.g., proteins) polymers can be classified as polyampholytes [18][19][20] A special example of such systems are the so-called polyzwitterionic polymers [21].
Due to their unique features, polyampholytes have a wealth of important and practical applications, e.g., as dispersing additives, gelling agents, rheology modifiers, etc. [18]. Nevertheless, polyampholytes have received much less attention, particularly in theoretical studies, as compared to neutral polymers. Consequently, our understanding of polyampholytes, especially tethered polyampholytes, is far from being complete.
In the case of uncharged solid surfaces and tethered chains whose segments bear the same charges, the electrostatic repulsion forces can be strong even at low grafting densities. Consequently, under such conditions, the chains are stretched and the system enters the brush regime [17]. In the case of tethered polyampholytes, the pinned layer behaves like a polyelectrolyte when the chains are globally charged and the resultant single chain charge is significant. However, for zero or for low net charge, the polyampholyte can coil and the tethered layer collapses. Low density tethered polyampholytes exhibit the highest deformation amplitude among the stretched and the collapsed states [22,23]. The structure and thermodynamic properties of tethered polyelectrolytes have been investigated using both theoretical [16,[24][25][26][27][28][29] and computer simulation [15,17,[30][31][32][33][34][35] methods. Density functional methods are among the most successful theoretical approaches [36][37][38][39][40][41].
Interaction forces between surfaces and colloidal stability are closely related. A stable colloidal suspension is characterized by repulsive forces between the colloidal particles, while attractive forces lead to an unstable suspension. One might expect that pinned polyelectrolyte layers induce repulsive forces and thus it leads to suspension stabilization. However, it has been also observed that polyelectrolytes may destabilize colloidal suspensions [42]. The problem of evaluation of the solvation forces between the surfaces modified with tethered charged chains and their dependence on the properties of a solvent and on the properties of the chains themselves was discussed in numerous papers [34,[43][44][45][46][47][48][49].
The purpose of this work is to use the density functional theory to calculate the force between two plates covered with tethered charged brushes. The approach used by us is an extension of the theory outlined in reference [39]. We will consider both the chains with the segments having identical charges (i.e., typical polyelectrolytes) and polyamphylic molecules with no net charge per chain. In the latter case, we wonder how the distribution of the charges along the chain (i.e., the chain architecture) effects the solvation force and the structure of a system.

Theory
We consider two identical surfaces, lower and upper, at a distance of H apart. Each surface is covered with charged tethered chains (P ). The chains are built of M tangentially jointed charged hard spheres of the same diameter, σ. The charge of the segment j , j = 1, . . . , M is Z (P ) j e, where e is the magnitude of an elementary charge. The chain connectivity is ensured by imposing a bonding potential [36][37][38][39] where R = (r 1 , . . . , r M ) is the vector describing the positions of all the segments. For tangentially bonded chains, the potential V b satisfies the relation The first segment of each chain, j = 1, is pinned at one of the surfaces. The surface bonding potential has the form where C is a constant and z ′ = z − σ/2 for lower and z ′ = z − (H − σ/2) for the upper surface. The interaction of all the remaining segments of a chain with the walls is described by the potential where v (P ) and

33801-2
is the Coulomb potential. In the above, Qe is the surface charge density of each wall and ǫ is the dielectric constant. According to equation (6), the surface charge density on each wall is assumed to be the same.
The total chain-surface potential is then V (P ) (z 1 , . . . , . The modified surfaces are in contact with a three-component fluid. The indices α = S, C and A abbreviate solvent, cation and anion species, respectively. The particles of the second and of the third species (C, A) are charged hard spheres as well, for the sake of simplicity, of a diameter σ, and carry the charges e Z C 1 and e Z A 1 . In other words, we consider the so-called solvent primitive model of an electrolyte. This is the simplest model of electrolyte solutions that takes into account the presence of non-zero size solvent molecules. The polar nature of a solvent is taken into account by retaining the dielectric constant in the Coulomb interaction between the ions. Note that the solvent primitive model was used by Henderson et al. to describe the forces between macroscopic particles in an electrolyte and to model electrolytemembrane systems [50,51]. Other examples of applications of this model can be found in references [52][53][54][55][56] Despite its extreme simplicity, the solvent primitive model has at least two virtues. It may be as reasonable a model of a solvent as can be treated conveniently by density functional theory and it does recognize that solvent molecules occupy space.
In this work we restrict our attention to the case of monovalent ions, i.e., Z C 1 = |Z A 1 | = 1 and |Z P j | = 1 for j = 1, . . . , M. The solvent S molecules, however, are uncharged hard spheres of the same diameter, σ. Thus The interactions between all the spherical species (i.e., between the molecules S, A and C and the consecutive segments of the chains P) are given by where η, α = S, C, A, P.
The confined system is in equilibrium with a three component bulk mixture of the components C, A and S. The bulk system is at the same temperature and at the same chemical potentials µ η , η = A, C, S as the confined system. The bulk densities of particular species are ρ η b .
In order to proceed, let us introduce the notation, ρ (P ) (R) and ρ (η) (z), η = S, C, A for the density distribution of chains and of fluid species, respectively. The theory is constructed in terms of the density of particular segments of chains, ρ (P ) s j (z), and the total segment density of chains, ρ (P ) s (z). These densities are introduced via commonly used relations [36][37][38][39]41] ρ (P ) In the system under study, all the local densities are only functions of the distance z from the lower surface that is set at z = 0.
The system is studied in a grand canonical ensemble with the constraint of constancy of the number of tethered chain molecules, i.e., where ρ c is the total number of tethered chain molecules per area of the surface. Since the confined system comprises two surfaces, the surface density of the chains is R C = ρ c /2.
The equilibrium density profiles are obtained by minimizing the thermodynamic potential In the above, F is free energy functional and q(z) is charge density The electrostatic Ψ(z) satisfies the Poisson equation A solution of the differential equation (12) for a fluid confined between two walls is given in reference [51]. It requires denomination of the value of electrostatic potential at the wall, Ψ 0 = Ψ(0) = Ψ(H ).
From the electro-neutrality condition of the system, it follows that The principal task in applying the density functional theory is to derive an expression for the Helmholtz energy as a functional of the density profiles. Following previous works [6][7][8][37][38][39][40], the Helmholtz energy is divided into an ideal term that depends on the bond potentials and the architecture of polymers, and the excesses arising from interactions of various forms. The latter terms are responsible for the thermodynamic nonideality. Within the framework of our model, the Helmholtz energy functional, F , is decomposed into the following contributions, where F id is the ideal contribution, The volume exclusion (the hard-sphere) term, F hs , is calculated according to the Fundamental Measure Theory, cf. references [6][7][8][36][37][38][39], where n i , i = 0, 1, 2, 3 and n V j , j = 1, 2 are, respectively, scalar and vector total weighted densities. The total weighted densities are sums of the weighted densities of individual species. For example, and n (P ) i j . Since the relevant equations defining the weighted densities have been already presented in numerous works [6][7][8][36][37][38][39][40], we have omitted them here.
The contribution F el , arising from the coupling between electrostatic and hard-sphere interactions is written down employing the approach described in detail in references [39,57] We have In the above Γ = 1 + 2κσ − 1 /2σ, T * is the "electrostatic" reduced temperature, T * = kT ǫσ/e 2 , The quantitiesρ s j andρ (η) are the "reference electrostatic system averaged densities", which are calculated according to references [39,57,58].

33801-4 Solvation forces between polyelectrolyte layers
Finally, the excess free energy functional due to the intra-chain correlation is given by where and y j is the contact value of the cavity correlation function.
We are aware of the fact that in our treatment, the calculation of the values of Γ is performed without discriminating the ions belonging to chains and free ions in the confined solution. In the light of the works [1,3,14,59,60], such ions should be distinguished. Therefore, the above expressions (19)- (22), and equation (22) in particular, are approximations. Note that in the case of bulk fluids, the sums in equation (21) are identical.
At equilibrium the density profiles minimize the thermodynamic potential Y , i.e., This condition leads to the equations and where and where the constant C is calculated from the normalization condition (9). The multidimensional density profile equation (24) can be then reduced to the equations for local densities of consecutive segments, ρ (P ) s j (z j ) using the method described in [39].
The solvation force per unit area is calculated from where A is the surface area and p is the pressure of the bulk deference fluid involving the components S, A and C.

Results and discussion
The systems in question are characterized by numerous parameters. In order to reduce their number, all the calculations have been performed assuming constant total bulk density of the fluid, ρ . This value is close to the density of water solutions at standard temperature and pressure (STP) conditions. The composition of the bulk fluid, i.e., the bulk mole fraction of the Also, the value of the reduced electrostatic temperature was kept constant and equal to T * = 0.2. This temperature is lower than the temperature usually used in 33801-5 computer simulations of nonuniform electrolytes, but our calculations have been intentionally carried out at a lower temperature in order to increase the role of electrostatic interactions. The parameters that characterize the brush are the number of segments and the surface brush density, R * C = R C σ 3 . Numerous calculations were carried out assuming that the total charge of the brush is zero, i.e., that the number of positively and negatively charged segments is the same, M + = M − = M/2. However, the distributions of "+" and "−" charges along the chain were different. We have also studied the cases of all the segments of the brush bearing the same charges. Moreover, in some cases, we carried out calculations for uncharged chains and for the chains with non-zero resultant charge, lower than Me. In all cases, the electrostatic potential at the wall, Ψ 0 , was set to be zero. We introduce the following codes to distinguish the systems under study. The symbol m 1 + n 1 − m 2 + n 2 −, . . . abbreviates a chain whose first m 1 segments are positively charged, the next n 1 -negatively charged, etc. When the sequence of the charges along the chain is repeated, we use parentheses to group the repeating units, for example the symbol 5(1 + 1−) means that the chain is built of 10 segments alternately charged with + and −, and whose first (pinned) segment is positively charged. In figure 1 (a), we show how the solvation force depends on the surface density of the grafted chains 10+, whose all segments bear the unit positive charges. For a comparison, we also display here the result for uncharged chains built of M = 10 segments. For a low surface density of the chains, R * C = 0.0025, the course of the solvation force vs. H * reminds us the course for the pore with non-modified walls [6,8]. The solvation force diverges at low wall-to-wall separations H * = H /σ and exhibits oscillations corresponding to attractive and repulsive forces that are well-pronounced at small values of H * and decay to zero at larger plate-to-plate separations. The maxima of the solvation force correspond to the development of consecutive layers of the solvent density profiles, see figure 1 (b). The role of chains and of electrostatic interactions in particular is small and the solvation force is almost entirely determined by the solvent packing effects.
When the surface density of chains increases, the dependence of the solvation force on H * quantitatively changes. For charged chains (10+), the solvation force becomes repulsive for all values of H * . Two factors play a significant role here. The first one is the volume exclusion effect: with an increasing R C , the segments of the chains occupy more and more space, and further compression becomes difficult for small H * . The second factor is electrostatic repulsion between the segments of chains pinned to the opposite walls. For R * C = 0.1, the second effect becomes very important. Indeed, the solvation force for charged and uncharged brush of the surface density 0.1 is very different. For uncharged chains [dotted line in figure 1 (a)], the solvation force becomes even slightly attractive (e.g., for H * ≈ 4). For H * 3 the solvation force becomes repulsive and grows rapidly. Electrostatic forces between the segments cause the appearance of strong repulsion between two walls at larger distances H * . The above results are in quantitative agreement with experiments [42,45]. Figure 1 (b) shows examples of a structure within some selected pores. Upper panels are for R * C = 0.0025 and for the chains 10+. Left upper panel displays the solvent density profiles (solid lines) and the negative total charge profiles (symbols) for two pores with H * = 3 and H * = 4. For H * = 4, three welldeveloped layers of the solvent appear within the pore. This structure is resistant to compression and a local maximum appears on the plot of the solvation force. When H * decreases down to 3, the inner layer is "squeezed out" the pore, and the next structure develops which is resistant to further compression. The total segment density profiles [the right upper panel of figure 1 (b)] behave similarly to the profiles of the solvent.
The average density of a fluid within the pore is much lower than in the bulk reference system (we recall that the latter equals 0.7). A part of the space inside the pore is occupied by the segments of chains [see the upper right panel of figure 1 (b)], but electrostatic interactions cause a lower average density in the confined system that is in chemical equilibrium with the bulk system. The upper left panel shows the negative charge density profile of ionic species of the fluid, −q ′ * (z) = −q ′ (z)σ 3 /e, where q ′ (z) = α=A,C ρ (α) (z).
Despite the fact that the electrostatic potential at the pore walls is zero, Ψ 0 = 0, an effective negative charge on the pore walls is generated. This is the result of the charges on the segments. Unlike the pores with non-modified walls (R * C = 0), the zero surface charge does not correspond to the zero of the electrostatic potential at the wall. It would be of interest to determine the dependence of the so-called "potential of zero charge", PZC, on the pinned chain characteristics, but this problem is out of the scope of the current research.
Lower panels in figure 1 (b) show the structure inside the pore of H * = 10 wide. The grafting density if high, R * C = 0.1. Due to the symmetry, only one half of the profiles is shown. We compare here the situation for uncharged chains built of M = 10 segments and for the chains 10+. Uncharged tethered

33801-7
chains are much more coiled. Indeed, their density is low at the pore center. The stretching of charged chains is caused by electrostatic repulsion between the segments. Positively charged segments attract negative ions, which, in turn, cause the accumulation of positive ions. However, the concentration of solvent molecules inside a pore is very low, because there is simply "no room" for them. In the case of uncharged chains, the concentration of all components of the fluid is very low inside the pore. This is because there is no attraction between the segments and ions, and the volume exclusion effects are responsible for a very low confinement of all the fluid species.
Here we consider the cases of the total charge of a single chain equal to 0. We studied the following distributions of the charges along the chain: (a) 5(1 + 1−), (b) 4(2 + 2−)1 + 1− and (c) 2(5 + 5−). The brush 10+ is treated as a reference. Figure 2 (a) shows the solvation forces for four systems defined above. The surface grafting density was R * C = 0.01. Surprisingly, despite quite big differences between the considered chains (chains with the total net charge, but with different distribution of charges vs. the chain with the total charge equal to 10e), the differences between the solvation force for particular systems are rather small. The solvation force for the chain 10+ is more repulsive. The maxima of f s are higher and the minima are shallower than for all the remaining cases, as expected from the consideration of the role of electrostatic forces. Quite unexpected, however, is a very small effect of different distributions of the charges along the chain on the solvation force. Indeed, the results for the systems (a)-(c) are hardly distinguishable on the figure scale. At a higher grafting density, R * C = 0.1, the solvation forces for the systems (a)-(c) are still very similar. However, the solvation force for the chains 10+ is very different from all the remaining curves and shows that an electrostatic repulsion occurs at larger wall-to-wall separations, cf. figure 2 (b). The pinned chains (a)-(c) can assume coiled configurations at both pore walls, whereas for the chain 10+, electrostatic repulsion between segments inhibits its coiling. Although the latter chains are more stretched than all the remaining chains, the total segment density profile for the chains 10+ is rather low at the pore center. Due to a low surface density, there is enough space for the segments to assume a slightly tilted or parallel to the wall configuration. Right panel in figure 3 (a) shows the profiles of ions, α = A, C. A significant difference between the profiles of anions and cations is observed for the chains 10+ due to electrostatic segmentanion attraction. For electro-neutral chains, the differences between the profiles of anions and cations are small. However, the profiles of ionic species are not identical. This indicates that there should be a small extra charge on the pore walls. In other words, the PZC for electro-neutral chains is not at Ψ 0 = 0.
An interesting question is how PZC depends on the architecture of the chains, but this problem is out of scope of the current research. Figure 3 (b) shows examples of the solvent density profiles (right hand panel) and the charge density profiles (left hand panel). Distribution of the solvent particles inside the pore depends very little on the charges of the segments. The solvent particles "do not feel" the charges, and the volume exclusion effects almost completely determine the structure of the confined solvent. The total charge distribution is sensitive to the architecture of the grafted chains.
Similar investigations of the structure have also been carried out for a higher grafting density, R * C = 0.1. The solvation forces for those systems are shown in figure 3 (b). Figure 4 (a) compares the total segment density profiles for the systems 10+ and for the systems (a)-(c) (left hand panel), as well as the density profiles of fluid components for the systems 10+ and 5(1 + 1−) (right hand panel). The pore is narrow, H * = 5, so the segments of the chains are compressed and form layered structures. The differences between the total segment densities for different systems are not big. Because of a small pore width, there is not much room for different rearrangements of the segments and the observed layered structure permits to effectively pack the segments inside the pore. . Lines without symbols are for anions α = A, with symbolsfor cations, α = C. The latter profile was multiplied by 100 to make it visible. We also shown here the solvent profile for the system 5(1 + 1−) (dotted line, α = S). Part (b). The profiles ρ (P ) s5 (z) (left panel) and the total charge profiles (right panel) for the systems from figure 2 (b). The meaning of the lines is the same as in figure 2 (b). The pore width is H * = 5 and the bulk mole fraction of the electrolyte is x = 0.1.
We have also inspected the density profiles of individual segments, but for the sake of brevity only the profiles of the middle segments, i = 5, are shown in the left hand panel of figure 4 (b). The differences between these profiles are significant. The performed analysis of the shape of the functions ρ (P ) si (z), i = 1, 2, . . . , 10 leads to the conclusion that in the case of 10+, the layers of chains pinned at the opposite pore walls penetrate into one another. However, in the remaining (a)-(c) cases, the chains "turn back" at the middle of the pore and the inter-penetration of the layers attached to the opposite walls is much weaker.
In the case of the pore H * = 5 modified with 10+ chains, the confined fluid contains almost solely anions ( figure 4 (a), left hand panel). The concentration of cations is extremely low; in order to make the cations profile visible on the figure scale, we multiplied it by 100. Similarly, the density of the solvent inside the pore is also extremely low. In the case of a pore modified with 5(1 + 1−) brush, the anions 33801-9 accumulate at the pore walls, while in the pore interior cations prevail. Moreover, for 5(1 + 1−) chains, we also observe a concentration of solvent molecules inside the pore.
The shape of total charge distributions inside the pore exhibits quite large local changes. In all the cases studied, there appears a very high peak at z * = 0.5 on the curves q(z) (we have cut the height of this peak on left hand panel of figure 4 (b) in order to make the figure readable). The height of this peak almost exactly corresponds to the amount of pinned charged segments (in all the cases, the pinned segments are positively charged). For a system 10+, the function q(z) is positive within the region of 0.5 < z < 1.5, while for the same region it is negative for all the remaining systems in question. Similarly to the previously considered cases, the integrals of q(z) over the entire pore are different from zero, i.e., there should be extra charges at the pore walls in order to ensure the total electro-neutrality of the entire system. In other words, the PZC does not correspond to the zero value of electrostatic potential at the wall.  cases, the effect of the bulk concentration of ions on the solvation force is small. It is slightly more pronounced for 8+ chains, especially for t 1.5 < H * < 2. In other words, the most important factors effecting the solvation force are the characteristics of the chain, i.e., its length, grafting density and architecture.
The results presented above have shown that one of the important factors effecting the behavior of the solvation force is the resultant charge of a brush. Therefore, in figure 7, we present the results for the chains of 8 segments with the resultant charge equal to 2e, 4e, 6e and 8e, respectively. Now, the differences between the solvation forces for a particular system are big. For a larger resultant charge a 33801-10  strong repulsion occurs at larger wall-to-wall separations. This is obviously due to electrostatic repulsions between the chains tethered at the opposite walls. When the resultant charge of a chain is 2e, we still observe effective attractions between two pore walls at H * ≈ 2.5 and for H * ≈ 3.5. For the resultant charge of 4e, there appears a marginally attractive force between two walls at z ≈ 2.5. However, for larger resultant charges, the solvation force is repulsive for all H * , although it exhibits several local minima and maxima. When the surface grafting density increases, the local extrema become weaker, and the curve describing the dependence of the solvation force on H * becomes smoother [cf. also figure 2 (a)].

Summary
In this work we applied the density functional theory to the study of the solvation forces between the walls covered with tethered layers of charged chains. We considered the chains possessing different net (resultant) charges and investigated the cases when all the segments had the same charges, as well as polyampholytes with zero and non-zero resultant charge per chain. The results of calculations showed that one of the factors that effects the solvation force very much is the resultant charge of the chain. However, the distribution of the charges along the chain effects the solvation force rather little. Another factor that quantitatively changes the solvation force is the grafting density. On the other hand, the concentration of the solution plays rather a minor role. Of course, the above results are valid for a specific model studied in this work. Moreover, we considered rather short chains and thus cannot exclude that our conclusions regarding the tethered layers built of chains involving hundreds of segments would be different.