text:
# 1 Introduction

## 1 Introduction

Polyelectrolytes represent a broad and interesting class of materials  that enjoy an increasing attention in the scientific community. For example, in technical applications polyelectrolytes are wildly used as viscosity modifiers, precipitation agents, superabsorbers, or leak protectors. In biochemistry and molecular biology they are of interest because virtually all proteins, as well as DNA, are polyelectrolytes.

In contrast to the theory of neutral polymer systems, which is well developed, the theory of polyelectrolytes faces several difficulties. Simple scaling theories, which have been proven so successfully in neutral polymer theory, have to deal with additional length scales set by the long range Coulomb interaction. Furthermore there is a delicate interplay between the electrostatic interaction of the distribution of the counterions and the conformational degrees of freedom, which in turn are governed by a host of short range interactions, which renders the problem difficult. There are only two limiting cases that are easy to solve. These are the case of high salt excess, effectively screening out the electrostatic interaction (which in turn allows one to treat it as a perturbation), or the case of an overwhelming dominance of the Coulomb force, which results in a strongly elongated chain. Unfortunately it is often just the intermediate case, which proves to be the most interesting regime in terms of application, experiment and theory.

Computational simulations provide some unique ways to elucidate the properties of charged systems. We first give a more general introduction to the relevant simulation methods, and focus then on some recently obtained results.

## 2 Simulations techniques

### 2.1 Some Basics on Simulations

There exist a number of nice reviews and books  which deal extensively with various aspects of computer simulations of complex, not necessarily charged, systems. Thus, the present introduction can also be viewed as a guide to the literature.

There are two basic concepts, which are used in computer simulations of complex systems. The conceptionally most direct approach is the molecular dynamics (MD) method. One numerically solves Newton’s equation of motion for a collection of particles, which interact via a suitable interaction potential $`U(\stackrel{}{r}_i)`$, where $`\stackrel{}{r}_i`$ are the positions of the particles. Through the equation of motion a natural time scale is built in, though, this might not be the physically realistic time scale (e. g. if the solvent is replaced by a dielectric background). Running such a simulation samples phase space for the considered system deterministically. Though this sounds very simple, there are many technical and conceptual complications, which one encounters on the way. The second approach, the Monte Carlo (MC) method, samples phase space stochastically. Monte Carlo is intrinsically stable but has no natural time scale built in. This can be reinterpreted, however, by an adjustment of suitable ”time amplitudes”. The MD and MC approaches are the basic simulation methods for exploring the statistical properties of complex fluids. At present, many applications employ variants thereof, or even hybrid methods, where combinations of both are used. Before going into detail we ask when is either kind of model appropriate?

At first sight it is tempting to perform a computer simulation of a polyelectrolyte solution where all details of the chemical structure of the monomers are included. For instance, the chain diffusion constant D could be measured by monitoring the mean square displacement of the monomers of the chains. This, however, is tempting only at the very first glance. Even for the fastest computers one would need an exceeding amount of computer time. As for all disordered, complex, macromolecular materials, (charged) polymers are characterized by a hierarchy of different length and time scales, and these length and especially the time scales span an extremely wide range . On the atomistic level the properties are dominated by the local oscillations of bond angles and lengths.111For reactions or to study exited states, the electronic structure is treated explicitely. Such methods (Carr-Parrinello simulations, quantum chemistry etc.) are beyond the scope of the present paper. The typical time constant of about $`10^{13}`$ sec results in a simulation time step of $`10^{15}`$ sec. On the semi-macroscopic level the behavior is dominated by the overall relaxations of conformation of the objects or even larger units (domains etc). The times, depending on chain length, temperature and density, can easily reach seconds. To cover that many decades in time within a conventional computer simulation is certainly impossible at present. On the other hand, it is important to relate the chemical structure of a system to its properties.

### 2.2 Molecular Dynamics

MD simulations date back to the early fifties. For a rather complete overview about simulations in condensed matter we refer to . Consider a cubic box of Volume $`V=L^3`$, containing $`N`$ identical particles. In order to avoid surface effects and (as much as possible) finite size effects, one typically uses periodic boundary conditions. The particle number density is given by $`\rho =N/L^3`$. The first simulations employed hard spheres of radius $`R_o`$, leading to a volume fraction $`\rho _v=4/3\pi R_o^3\rho `$. Though still used extensively for some studies on the glass transition of colloidal systems we focus here on soft potentials.

A key thermodynamic quantity, the temperature, is imposed via the equipartition theorem $`\frac{m}{2}\stackrel{}{\dot{r}}_i^2=\frac{3}{2}kT`$, $`m`$ being the particle masses. Note that in hard sphere systems temperature only defines a time scale, but is otherwise irrelevant. One can find many different soft potentials in the literature. However, most widely used is the Lennard-Jones potential $`U^{LJ}(r_{ij})`$, derived originally for interactions of noble gases (Ar, Kr …), $`r_{ij}`$ being the distance between particle $`i`$ and $`j`$. In its simplest form for two identical particles it reads

$`U^{LJ}(r_{ij})=4ϵ\left[({\displaystyle \frac{\sigma }{r_{ij}}})^{12}({\displaystyle \frac{\sigma }{r_{ij}}})^6\right]`$ (1)
Usually, a cutoff $`r_c`$ is introduced for the range of the interaction. This typically varies between 2.5 $`\sigma `$ (classical LJ interaction with an attractive well of depth $`ϵ`$, used for the poor solvent chains in Sec. 4.2) and 2$`{}_{}{}^{1/6}\sigma `$ (the potential is cut off at the minimum, leading to a pure repulsive interaction, as is typically done for good solvent chains). For chain molecules a bonding interaction for $`r<R_0`$

$`U_{\mathrm{FENE}}(r)={\displaystyle \frac{1}{2}}kR_0^2\mathrm{ln}(1{\displaystyle \frac{r^2}{R_0^2}})`$ (2)
is added, which keeps the bond length below a maximum of $`R_0`$. The spring constant $`k`$ varies between 5 and 30 $`ϵ/\sigma ^2`$. Electrostatics is included via the Coulomb interaction

$`U_c(r_{ij})={\displaystyle \frac{l_Bk_BT}{r_{ij}}}q_iq_j`$ (3)
Here $`l_B:=\frac{e_0^2}{4\pi \epsilon _0\epsilon _rk_BT}`$, where $`e_0`$ is the elementary unit charge, $`k_B`$ is the Boltzmann constant, $`T`$ denotes temperature, $`\epsilon _0`$ and $`\epsilon _r`$ are the vacuum and relative dielectric permeability of the solvent, respectively, and $`q_{i,j}`$ are charges measured in units of $`e_0`$. Two monovalent charges separated by the Bjerrum length $`l_B`$ have an interaction energy equal to $`k_BT`$. The Bjerrum length thus is a measure of the interaction strength. It is equal to 7.14 Å for water at room temperature. The computational aspects of the long range potential will be discussed shortly in Sec.2.4.

The unit of energy is $`ϵ`$, of length $`\sigma `$ and of mass $`m`$. This defines the ”LJ-units” for temperature $`[T]=ϵ/(k_BT)`$, time $`[t]=\sqrt{\sigma ^2m/ϵ}`$ and number density $`[\rho ]=\sigma ^3`$. In most practical programs $`\sigma `$, $`m`$, $`ϵ`$ are used as the basic units and set to one. The straight forward simulation technique is to integrate Newton’s equations of motion for the particles:

$`m_i\ddot{\stackrel{}{r}}_i=\stackrel{}{}{\displaystyle \underset{j,ji}{}}U(r_{ij})`$ (4)
Since energy in such a simulation is conserved we have a microcanonical ensemble. Presently other thermodynamic ensembles are commonly used for practical applications (NPT: isobaric-isothermal, NVT: isothermal (canonic) …). Because we often employ a stochastic MD method, known as the Langevin thermostat, we will briefly describe its main ingredients. Instead of integrating Newton’s equations of motion, one solves a set of Langevin equations

$`m_i\ddot{\stackrel{}{r}}_i=\stackrel{}{}{\displaystyle \underset{j,ji}{}}U(r_{ij})\mathrm{\Gamma }\dot{\stackrel{}{r}}_i+\stackrel{}{\xi }_i(t)`$ (5)
with $`\stackrel{}{\xi }_i(t)`$ being a $`\delta `$-correlated Gaussian noise source with its first and second moments given by

$`\stackrel{}{\xi }_i(t)=0\text{ and }\stackrel{}{\xi }_i(t)\stackrel{}{\xi }_j(t^{})=6k_BT\mathrm{\Gamma }\delta _{ij}\delta (tt^{}).`$ (6)
The friction term $`\mathrm{\Gamma }\dot{\stackrel{}{r}}_i`$ and the noise $`\stackrel{}{\xi }_i(t)`$ is thought of as imitating the presence of a surrounding viscous medium responsible for a drag force and random collisions, respectively. The second moment of $`\stackrel{}{\xi }_i(t)`$ is adjusted via an Einstein relation in order to reach the canonical state in the limit $`t\mathrm{}`$. The dynamics generated by the Langevin equation can alternatively be written as a general Fokker-Planck process. This permits a transparent proof of two important facts: (i) the stationary state of the process is the Boltzmann distribution and (ii) the system will evolve and converge to the Boltzmann distribution . Since for small times the stochastic part is more important than the deterministic one ($`\sqrt{t}t`$ for small $`t`$) it is actually not necessary to use Gaussian random variables in the simulation . It suffices to use equidistributed random variables with first and second moment being identical to the Gaussian deviate.

A simple but very efficient and stable integration scheme is the Verlet algorithm (more complicated methods, which do not have time inversion symmetry, do not, in general, perform significantly better). With a simulation time step $`\delta t`$ , where $`\delta t<<2\pi /\omega _{max}`$ and $`\omega _{max}`$ is the typical highest frequency of the system (for crystals the Einstein frequency), we have in one dimension

$`r_i(t+\delta t)`$ $`=`$ $`r_i(t)+\delta tv_i(t)+{\displaystyle \frac{\delta t^2}{2}}a_i(t)+{\displaystyle \frac{\delta t^3}{6}}\dot{a_i}(t)+𝒪(\delta t^4)`$
$`r_i(t\delta t)`$ $`=`$ $`r_i(t)\delta tv_i(t)+{\displaystyle \frac{\delta t^2}{2}}a_i(t){\displaystyle \frac{\delta t^3}{6}}\dot{a}_i(t)+𝒪(\delta t^4),`$
where $`v_i(t)=\dot{r}_i(t)`$ and $`a_i(t)=\dot{v}_i(t)`$. An addition of the two lines yields

$`r_i(t+\delta t)=2r_i(t)r_i(t\delta t)+\delta t^2a_i(t)+𝒪(\delta t^4)`$ (7)
Therefore the position calculations have an algorithmic error of $`𝒪(\delta t^4)`$. Subtraction of the lines yields

$`v_i(t)={\displaystyle \frac{1}{2\delta t}}[r_i(t+\delta t)r_i(t\delta t)]+𝒪(\delta t^3)`$ (8)
leading to errors of $`𝒪(\delta t^3)`$. There are many variants of this basic integrating scheme used throughout the literature . One can follow the realistic time evolution of a system, as long as the forces/potentials are realistic for the modeled system and as long as classical mechanics is sufficient. In a purely deterministic simulation the accumulation of small errors can cause significant deviations from the real trajectory. If the system is ergodic, which requires mixing of normal modes (recall the well-known Fermi-Pasta-Ulam problem, where one asks how anharmonic a potential has to be in order to equilibrate a one dimensional chain of particles) one can determine ensemble averages from time averages

$`<A>={\displaystyle \frac{1}{M}}{\displaystyle \underset{i=1}{\overset{M}{}}}A(t_i)`$ (9)
of any physical quantity $`A`$ of interest. This describes the most elementary Ansatz for a microcanonical simulation . Here all extensive thermodynamic variables of the system, namely $`N`$, $`V`$, $`E`$ are kept constant. Sometimes this is also called $`NVE`$ ensemble. As mentioned before, most applications employ other ensembles such as the canonic ($`NVT`$), the isobaric-isothermal ($`NPT`$) or even the grand canonical $`(\mu ,P,T)`$ ensemble, $`\mu `$ being the chemical potential. As a general rule, in cases such as two phase coexistence or calculations of transport properties, one should choose an ensemble with many intensive variables kept constant as possible. For charged systems, however, it is rather difficult to perform efficient simulations in the ($`N`$, $`P`$, $`T`$) or ($`\mu `$, $`P`$, $`T`$) ensembles. Therefore the most common ensemble is the NVT since it is easy to use the deterministic equations with additional stochastic terms to constrain the temperature (Eq. 5).

### 2.3 Monte Carlo Method

The classical Monte Carlo approach goes to the other extreme, namely to purely stochastic sampling. Starting from a particular configuration, randomly a particle (or a number of particles) is selected and displaced by a random jump. For hard sphere systems the move is accepted if the new configuration complies with the excluded volume; if not the old configuration is retained. The approach is also called simple sampling. This cycle is repeated over and over. Once every particle on average has a chance to move, one Monte Carlo step is completed. This is the most basic Monte Carlo simulation (see e.g.). Since there is no energy involved, it trivially fulfills detailed balance

$`W(\{x\}\{y\})P_{eq}(\{x\})=W(\{y\}\{x\})P_{eq}(\{y\})`$ (10)
where $`W(\{x\}\{y\})`$ is the probability to jump from state $`\{x\}`$ to state $`\{y\}`$ and $`P_{eq}(\{x\})`$ the equilibrium probability of state $`\{x\}`$. All states have exactly the same probability. Detailed balance is a sufficient condition for a MC simulation to relax into thermal equilibrium, though this may take a very long time. Special cases of algorithms without detailed balance will not be discussed here.

It is useful to compare the basic aspects of MC simulation to the examples discussed above for molecular dynamics simulations. The Hamiltonian depends for simplicity only on the positions of all particles $`\{r_i\}`$ and is denoted by $`H(\{r_i\}`$. The expectation value of any observable A is given by

$`<A>={\displaystyle \underset{\{r_i\}}{}}A(\{r_i\})P_{eq}(\{r_i\})`$ (11)
with

$`P_{eq}(\{r_i\})`$ $`=`$ $`\mathrm{exp}(H(\{r_i\})/k_BT)/Z`$
$`Z`$ $`=`$ $`{\displaystyle \underset{\{r_i\}}{}}\mathrm{exp}(H/k_BT)`$ (12)
An exact way would be to sample all possible states, which in all but the most trivial cases is impossible. Thus we sample phase space stochastically. Taking a particle at random, calculating its energy, one moves it and calculates the new energy. With $`P(\{x\})`$ being the Boltzmann-probability of the original state and $`P(\{y\})`$ being that of the new state, detailed balance is obeyed if

$`{\displaystyle \frac{W(\{x\}\{y\})}{W(\{y\}\{y\})}}=\mathrm{exp}\{(H(\{x\})H(\{y\})/k_BT)\}`$ (13)
Under this condition the algorithm is ergodic and the system relaxes into equilibrium. The Metropolis method is the most frequently used prescription one to accept or reject a move:

$`W(\{x\}\{y\})=\mathrm{\Gamma }\{\begin{array}{ccc}\mathrm{exp}(H(\{x\})H(\{y\})/k_BT)& ,& \mathrm{}H>0\\ 1& ,& \mathrm{}H<0.\end{array}`$ (16)
Since only the ratio of the W’s is relevant, $`\mathrm{\Gamma }`$ is an arbitrary constant between zero and one, usually $`\mathrm{\Gamma }`$ = 1. A random number $`x`$, equally distributed between 0 and 1, is used to decide upon the acceptance of a move. If $`x<W(\{x\}\{y\})`$ the move is accepted, otherwise rejected. (For $`\mathrm{\Gamma }`$ = 1 any move, which lowers the energy is accepted.) This is the basic MC procedure used for sampling phase space in statistical physics.

In many cases, however, one also would like to gain information on the dynamics of a system or even better, of a model system. How can a MC simulation, with no intrinsic time scale, be used to obtain information on the dynamics? In the method described above the system evolves from one state to another by a local move. Through these local stochastic moves the configurations of particles change with ”time”. This is a dynamic MC method based on a Markov process, where subsequent configurations $`\{x\}\{x^i\}\{x^{ii}\}\mathrm{}`$ are generated with a transition probability $`W(\{x^i\}\{x^{ii}\})`$. To a large extent the choice of the move is arbitrary, as long as one can interpret it as a local elementary unit of motion. The prefactor $`\mathrm{\Gamma }`$ can actually be interpreted as an attempt rate $`\mathrm{\Gamma }=\tau _o^1`$ for the moves and introduces a timescale. This ”changes” the purely statistical transition probability into a transition probability per unit time. To compare simulated (overdamped) dynamics with an experiment, it essentially requires determination of $`\tau _o`$ (e. g. diffusion constants). It is obvious, that this simulation does not include any hydrodynamic effects since there is no momentum involved. There are very interesting, more advanced methods like DPD (dissipative particle dynamics) and Lattice-Boltzmann methods, currently under development in order to include this efficiently ). Using the interpretation of a MC step as a time step, ensemble averages can be written as time averages:

$`<A>={\displaystyle \frac{1}{MM_o}}{\displaystyle \underset{i=M_o+1}{\overset{M}{}}}A(\{x^i\}){\displaystyle \frac{1}{tt_o}}{\displaystyle _{t_o}^t}𝑑t^{}A(t^{}).`$ (17)
We view one attempted move per system particle as one time-step.

The first configurations in a simulation are usually not yet equilibrium configurations. One first has to ”relax” the system into equilibrium, meaning the data for the first $`M_o`$ steps are omitted. In this interpretation the dynamic Monte Carlo procedure is nothing but a numerical realization of a Markov process described by a Markovian master equation

$`{\displaystyle \frac{d}{dt}}P(\{x\},t)`$ $`=`$ $`{\displaystyle \underset{\{x^i\}}{}}W(\{x\}\{x^i\})P(\{x\},t)`$ (18)
$`+{\displaystyle \underset{\{x^i\}}{}}W(\{x^i\}\{x\})P(\{x^i\},t)`$
with $`P(\{x\},t)`$ the time dependent probability of state $`\{x\}`$. The condition of detailed balance is sufficient that $`P_{eq}(\{x\})`$ is the steady-state solution of the master equation. If all states are mutually accessible $`P(\{x\},t)`$ must relax towards $`P_{eq}(\{x\})`$ as $`t\mathrm{}`$ irrespective of the starting state. Note however that the choice of a ”good” starting state can save enormous amounts of CPU time.

So far, the two extreme cases for classical, particle based computer simulations were discussed, microcanonical MD and canonical MC. There are many approaches ”in between” which are used depending on the problem under consideration. The techniques range from pure MD, where Newton’s equations of motion are solved ($`\ddot{x}=U)`$, MD coupled to a heat bath and added friction (”Langevin MD”, ”Noisy MD”), ($`\ddot{x}=U\zeta \dot{x}+f(t))`$, $`\zeta `$ friction, f (t) random force), Brownian Dynamics (BD) ($`\dot{x}=U`$ \+ random displacement), force biased MC (attempted moves are selected from the very beginning according to local forces), to plain MC as described above.

For the application to polymers one should keep in mind, that the conformational entropy of the chains add additional complications, which make proper equilibration especially difficult or time consuming. Thus, wherever possible, methods should be used, which are faster than the slow intrinsic dynamics of the chains .

### 2.4 Methods for long range interactions

One of the biggest problems for the simulations of charged systems is the long range nature of the Coulomb interactions. In principle, each charge interacts which all others, leading to a computational effort of $`𝒪(N^2)`$ already within the central simulation box. For many physical investigations one wants to simulate bulk properties and therefore introduces periodic boundary conditions to avoid boundary effects. The standard method to compute the merely conditionally convergent Coulomb sum

$`E={\displaystyle \frac{1}{2}}{\displaystyle \underset{\stackrel{}{n}}{\overset{}{}}}{\displaystyle \underset{ij}{}}{\displaystyle \frac{q_iq_j}{|\stackrel{}{r}_{ij}+\stackrel{}{n}L|}},`$ (19)
where the prime denotes that for $`\stackrel{}{n}=\stackrel{}{0}`$ the term $`i=j`$ has to be omitted, is the traditional Ewald summation . The basic idea is to split the original sum via a simple transformation into two exponentially convergent parts, where the first one, $`\varphi _r`$, is short ranged and evaluated in real space, the other one, $`\varphi _k`$, is long ranged and can be analytically Fourier transformed and evaluated in Fourier space:

$`{\displaystyle \frac{1}{r}}={\displaystyle \frac{1f(r)}{r}}+{\displaystyle \frac{f(r)}{r}}\varphi _r(r_c,\alpha )+\varphi _k(k_c,\alpha )`$ (20)
Traditionally, one uses for $`f`$ the error function $`\mathrm{erf}(\alpha r):=2\pi ^{1/2}_0^{\alpha r}\mathrm{exp}t^2\mathrm{d}t`$, though other choices are possible and sometimes more advantageous. For any choice of the Ewald parameter $`\alpha `$ and no truncation in the sums the formula yields the exact result. In practice one wants to cut off the infinite sum at some finite values $`r_c`$ and $`k_c`$ to obtain $`E`$ to a user controlled accuracy, which is possible by using error estimates . The aformentionend procedure results in the well known Ewald formula for the energy of the box

$`E=E^{(r)}+E^{(k)}+E^{(s)}+E^{(d)},`$ (21)
where the contributions from left to right are the real space, Fourier space, self , and dipole-correction energy terms. These are given by

$`E^{(r)}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{i,j}{}}{\displaystyle \underset{𝐦^3}{\overset{}{}}}q_iq_j{\displaystyle \frac{\text{erfc}(\alpha |𝐫_{ij}+𝐦L|)}{|𝐫_{ij}+𝐦L|}}`$ (22)
$`E^{(k)}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{L^3}}{\displaystyle \underset{𝐤0}{}}{\displaystyle \frac{4\pi }{k^2}}e^{k^2/4\alpha ^2}|\stackrel{~}{\rho }(𝐤)|^2`$ (23)
$`E^{(s)}`$ $`=`$ $`{\displaystyle \frac{\alpha }{\sqrt{\pi }}}{\displaystyle \underset{i}{}}q_i^2`$ (24)
$`E^{(d)}`$ $`=`$ $`{\displaystyle \frac{2\pi }{(1+2ϵ^{})L^3}}\left({\displaystyle \underset{i}{}}q_i𝐫_i\right)^2,`$ (25)
and the Fourier transformed charge density $`\stackrel{~}{\rho }(𝐤)`$ is defined as

$`\stackrel{~}{\rho }(𝐤)={\displaystyle _{V_b}}\text{d}^3r\rho (𝐫)e^{i𝐤𝐫}={\displaystyle \underset{j=1}{\overset{N}{}}}q_je^{i𝐤𝐫_j}.`$ (26)
The dipole term depends not on $`\alpha `$, hence is independent of the splitting function. It reflects the way Eq. (19) is summed up, here in a spherical way towards infinity. It also includes a correction for the dielectric constant $`ϵ^{}`$ outside the summed up sphere volume. For metallic boundary conditions, $`ϵ^{}=\mathrm{}`$, and the dipole term vanishes. In principle the thermodynamic properties should be independent of the choice of boundary conditions. The Ewald sum has complexity $`𝒪(N^{3/2})`$ in its optimal implementation , and therefore is not suitable for the study of large systems ($`N>𝒪(1000)`$). Implementing a fast Fourier transformation (FFT) for the Fourier part results in the so-called particle-mesh-Ewald formulations, which improve the efficiency to $`𝒪(N\mathrm{log}N)`$ , and which can also be efficiently be parallelized. The most versatile and accurate method of all mesh-methods is the oldest P3M algorithm, for which also precise error estimates exist. Another way of computing Eq.(19) is via a convergence factor

$`E=\underset{\beta 0}{lim}{\displaystyle \frac{1}{2}}{\displaystyle \underset{\stackrel{}{n}}{\overset{}{}}}{\displaystyle \underset{ij}{}}{\displaystyle \frac{q_iq_j\mathrm{exp}(\beta |\stackrel{}{r}_{ij}+\stackrel{}{n}L|)}{|\stackrel{}{r}_{ij}+\stackrel{}{n}L|}}.`$ (27)
This approach is used in the Lekner  and Sperb  methods to efficiently sum up the 3D Coulomb sum. Although the method in its original versions has $`𝒪(N^2)`$ complexity, Sperb et al. have developed a factorization approach which yields an $`𝒪(N\mathrm{log}N)`$ algorithm .

Other advanced methods of $`𝒪(N\mathrm{log}N)`$ are tree algorithms, which are the first order approximation of even better, so-called fast multipole methods . These can reach a linear complexity, but at the expense of a heavy computational overhead which makes these methods advantageous only for a very large number of charges ($`N\mathrm{100\hspace{0.17em}000}`$) .

For thin polyelectrolytes films or membrane interactions one is also interested in summations where only 2 dimensions are periodically replicated and the third one is of finite thickness $`h`$ ($`2D+h`$ geometry). For this geometry Ewald based formulas are only slowly convergent, have mostly $`𝒪(N^2)`$ scalings and no “a priori” error estimates exist . Recently Arnold  developed a method which is based on convergent factors, whose errors are well controlled, and which uses a factorization approach resulting in an $`𝒪(N^{5/3})`$ scaling (MMM2D). In two dimensions the convergence factor based methods and the Ewald sum methods yield exactly the same results, there is no dipolar correction term needed. This is in contrast to the 3D methods. However, an even better scaling can be achieved, if one returns to the 3D Ewald formula, and allows for a large empty space between the unwanted replicas in the third dimension . However, so far the method has been only checked on a trial and error basis. We recently improved this situation by computing an analytic error term which accounts for the contributions of the unwanted replicas. By simply subtracting this term from the 3D sum, which is a linear operation in $`N`$, one can in principle come as close as desired to the real $`2D+h`$ sum, by allowing for just some arbitrary small amount of empty space between the layers. Using then again the P3M method we obtain an $`N\mathrm{log}N`$ scaling with well controlled errors also for the $`2D+h`$ geometry, which up to now seems to be the optimal choice.

To simulate the structure of water (or other dipolar solvents), one also needs to treat the dipolar interactions in a similar fashion. Also here the Ewald method is applicable, and error estimates exist.

## 3 Mean-field models: Debye-Hückel chains

While we deal here mainly with systems where the charges are explicitly taken into account, historically and even up to now many studies consider the ions solely in a mean field approximation. In the first step all non bonded charges are considered as a smeared continuous charge density. Such a situation is described by the Poisson-Boltzmann (PB) equation. While this in most cases is not exactly solvable, many studies employ the Debye-Hückel approximation, which is the solution of the linearized PB equation. The resulting potential between charges is the screened Coulomb potential

$`V_{DH}={\displaystyle \frac{l_B}{e_0}}k_BT{\displaystyle \frac{\mathrm{exp}(\kappa r)}{r}},`$ (28)
with$`1/\kappa `$ being the Debye screening length. The polymer can now easily be modeled as a random walk of $`N`$ monomers from which a fraction $`f`$ is monovalently charged. If one in addition introduces the stiffness along the backbone of the chain by a cosine-potential, the total Hamiltonian reads

$`{\displaystyle \frac{H}{k_BT}}=A{\displaystyle \underset{i=1}{\overset{N1}{}}}(\stackrel{}{b}_i\stackrel{}{b}_{i+1})+{\displaystyle \underset{i=2}{\overset{N}{}}}{\displaystyle \underset{j=1}{\overset{i1}{}}}\theta (q_iq_j)l_B{\displaystyle \frac{\mathrm{exp}(\kappa r_{ij})}{r_{ij}k_BT}}.`$ (29)
The positive amplitude $`A`$ defines the strength of the angular potential and the Heavyside $`\theta `$-function is one if the monomers $`i`$ and $`j`$ are charged, and zero otherwise. The symbol $`\stackrel{}{b}_i`$ is the bond vector between monomer number $`i`$ and $`i+1`$. For the present simulations the bond length $`|\stackrel{}{b}|`$ and the Bjerrum length $`l_B`$ are fixed to one $`\sigma `$, and thus set the basic length scale. One bond mimics several neutral monomers as usual for coarse grained simulation models. The longest chains considered contained up to 2049 repeat units with a charge fraction of $`f=\frac{1}{16}`$. Mapping this onto a PS-NaPSS copolymer for a more flexible case and keeping in mind that for these polymers roughly three repeat units fit into one Bjerrum length one arrives at a molecular weight of more than 600 000 g/mol. Thus, for these kind of questions computer simulations are quite capable of covering the experimentally interesting regime. This allows us to systematically vary not only the chain length and the screening length but also the chain bending stiffness. This is only of limited experimental relevance since very long isolated chains in dilute solutions cannot be experimentally analyzed so far. However one of the central questions in the theory of polyelectrolytes is whether the characteristic electrostatic length is a linear or quadratic function of the Debye length (i.e. proportional to the square of charge density or to the charge density itself). Analytic results mainly predict a $`\kappa ^2`$ asymptotic dependency of the persistence length. However, an estimate of the required the chain length to reach the asymptotic regime for flexible weakly charged polyelectrolytes shows that the chains are so long that all experimentally realizable concentrations are in the semi-dilute regime. Thus, this is a typical question of interest which is of no direct importance for experiments. Nevertheless, it can add significantly to our understanding which in turn can also influence the interpretation of experiments in the semi-dilute regime. Therefore, it is worthwhile to undertake some efforts in computer simulation to investigate this question. There are a number of attempts to do this, which so far do not lead to a clear-cut answer. A typical result from such a simulation is given in Fig. 1.

It shows the dependence of the electrostatic persistence length on the screening length. The two dotted lines indicate the variational (Var) and the Odijk-Skolnik-Fixman (OSF) results ( $`L_e\kappa ^1,`$ or $`\kappa ^2`$ respectively). As one can see for the regime covered there seems to be a continuous transition from a very weak dependency on $`\kappa `$ towards the asymptotically expected $`\kappa ^2`$ behavior for the semi-flexible polyelectrolyte. In the limit of very stiff chains the OSF result becomes exact and the data has to agree with that. Even though the simulated chains lie very well within the experimentally relevant regime, there is, at least for typical experimental flexibilities, no clear sign of the crossover into the asymptotic OSF regime. This allows for two important conclusions. First of all, the chains are certainly not long enough to display asymptotic behavior. Secondly, typical experimental chains are also not long enough to display the asymptotic behavior as predicted by mean-field theories. Unlike neutral polymers, polyelectrolytes normally are not in an asymptotic limit where the predicted scaling laws can be cleanly observed. Another consequence is revealed in a closer analysis of the simulation data. For semiflexible polyelectrolyte chains there is no unique persistence length anymore, as all theoretical pictures assume. Over short distances the intrinsic stiffness dominates and gives a clear signal in the bond direction correlation function. Only over large distances along the backbone of the chain does the electrostatic effect show up and introduces a second characteristic length scale into the system. This can lead to the special situation that the electrostatic contribution to the persistence length for stiff chains is smaller than for flexible chains. The reason simply lies in the fact that for the stiff chains the charges along the backbone are already much further apart because of the intrinsic stiffness as opposed to flexible chains. These larger distances then experience a dramatically weaker electrostatic repulsion due to the screening of the electrostatic interaction resulting in a weaker effect on the persistence length. For details we refer to Ref. .

Altogether, the main conclusion from these simulations of simplified mean-field like models is on the one hand that polyelectrolytes have to be extremely long to be in the asymptotic regime, and that one has to be very careful in deriving general statements from either simulations, which typically are not asymptotic just as experiments, and analytic theory, which is usually in an asymptotic limit.

Another conclusion is that the concept of persistence length for polyelectrolytes certainly is not well defined in terms of the properties of the chains. As soon as intrinsic stiffness is included there is no longer a unique length scale that describes the internal structure of the chain. On the other hand this is the basis of all classical models used in analytic theory.

## 4 Solutions of flexible polyelectrolytes

### 4.1 Good solvent chains with explicit counterions

The first investigation of totally flexible many chain polyelectrolyte systems in good solvent with explicit monovalent counterions was performed several years ago. The simulations were carried out mostly with systems of 8 or 16 chains with $`N_m`$ = 16, 32 and 64. Instead of the P3M algorithm a spherical approximation in a truncated octahedral simulation box was used which, for values smaller than $`N_{total}500`$ is faster than the PME method. More details of the whole study can be found in . All beads and counterions interacted with the truncated Lennard-Jones potential plus the full Coulomb interaction.

In this work experimental values of the osmotic pressure and the maximum position in the interchain structure facture were successfully reproduced. One of the important findings was that the chains essentially are never rodlike. Counterion-chain correlations can dramatically shrink the polyelectrolyte chain. The end-to-end distance shortens significantly as the density increases from dilute towards the overlap density. The chain structure is highly anisotropic in the very dilute limit, and the scaling with respect to $`N_m`$ is asymmetric; but as the overlap density is approached, the structural anisotropy dissipates and the scaling becomes approximately symmetric. On long length scales the chain structure continuously changes from very elongated to neutral-like coils. Yet, on short length scales, the chain structure is density independent and elongated more than neutral chains.

It was found that in the dilute limit the scaling for the extension perpendicular to the chain was $`R_{}N^{0.650.70}`$, and for the extension parallel $`R_{}N^{0.901.00}`$. Near the density, where the rodlike chains in disordered solution overlap, $`\rho N^2,R_{}`$ grows at the expense of $`R_{}`$ until at the overlap density $`\rho ^{}`$ the effective exponent is about 0.82. The transition regime ranges from $`\rho N^2`$ to about $`\rho N^{1.4}`$ where the coils start to overlap and one eventually reaches $`\nu =1/2`$ in the semidilute regime. The exponents reported should not necessarily be taken as asymptotic ($`N\mathrm{}`$), however they should be relevant for many experimental systems.

### 4.2 Poor solvent chains with explicit counterions

Many polyelectrolytes possess a carbon based backbone for which water is a poor solvent. Therefore, in aqueous solution, there is a competition between the solvent quality, the Coulombic repulsion, and the entropic degrees of freedom. The conformation in these systems can under certain conditions assume pearl-necklace like structures. These also exist for strongly charged polyelectrolytes at finite densities in the presence of counterions. The simulations in Ref. used 16 chains of length $`N_m`$ = 94, with a charge fraction of $`f`$ = 1/3, and monovalent counterions. The hydrophobic interaction strength was tuned by means of the Lennard-Jones parameter $`ϵ`$. There we showed that the polymer density $`\rho `$ can be used as a very simple parameter to separate different conformation regimes. This can already be seen in the plots of the end-to-end distance $`R_e`$ and $`r=\frac{R_E^2}{R_G^2}`$ versus $`\rho `$ in Fig. 2. At very high densities the electrostatic interaction is highly screened, so that the hydrophobic interaction wins, and the chains collapse to dense globules. If one slightly decreases the density, the chains can even contract further, because there are no more steric hinderences from the other chains or counterions, and the screening is smaller. The collapsed globules, however, have still a net charge, and repel each other, so that this phase resembles a charged stabilized colloid or microgel phase.

With decreasing density the electrostatic interaction will dominate over the hydrophobic one. The chains will tend to elongate, assuming pearl-necklace conformations, like in Fig. 3, as they have been predicted for weakly charged polyelectrolytes in Ref. . The more the chain stretches, the smaller become the locally compact regions. Note that in contrast to the analytical theories, the pearls are stable, even though there are counterions localized near and/or inside the pearls.

Experimentally there are some hints for the existence of pearl-necklace chains. One of the obstacles to observing them in scattering experiments could be related to the strong fluctuations of the pearl number. Even in equilibrium we have found coexistence of several pearl states. In Fig. 4 we see the time evolution of one single chain composed of 382 monomers with a charge fraction $`f=1/3`$, $`l_B=1.5`$, and $`ϵ=1.75`$ in a many chain system at density $`\rho =1.48\times 10^5`$. One observes jumps between a five and four pearl configuration. Also the position of the pearls move quite vividly.

The different length scales appearing in a chain can be analyzed by looking at the spherically averaged form factor $`S_1(q)`$ of the chain.

The maximum seen at $`q=6`$ comes from the monomer extension. In the range $`1<q<2`$ we observe a sharp decrease in $`S_1`$, which comes from the scattering from the pearls, because it shows the typical Porod scattering of $`S_1(q)q^4`$. The kink at $`q1.66`$ appears at the position expected from the pearl size, but is broadly smeared out due to large size fluctuations. The shoulder which can be seen at $`q0.5`$ does not come from the intra-pearl scattering but is due to the scattering of neighboring pearls along the chain (inter-pearl contribution), which have a mean distance of $`r_{PP}=13.3.\sigma `$. It is also smeared out due to large distribution of inter-pearl distances. We conclude that the signatures of the pearl-necklaces are weak already for monodisperse samples. A possible improvement could be achieved for chains of very large molecular weights and only few pearl numbers, which could lead to stable and large signatures. Many more interesting results on poor solvent polyelectrolytes can be found in Ref. and will be published soon.

### 4.3 Counterion distribution around finite polyelectrolytes

We recently completed a study of the spatial distribution of the counterions around strongly charged, flexible polyelectrolytes in good and poor solvent. There we demonstrated that by partially neutralizing the quenched charged distribution on the chain backbone the inhomogeneous distribution of counterions leads to the same qualitative effects that are observed in weakly charged polyelectrolytes with an annealed charge distribution. This is due to the presence of the mobile partially neutralizing counterions, which results in an annealed backbone charge distribution. The common underlying physical mechanism for the end-effect is the differences in the electrostatic field of the chain along its backbone. The strength of the end-effect depends on parameters like chain length, charge fraction and ionic strength, and those dependencies were found in agreement with the scaling predictions. We found a saturation of the end-effect for long chains, when the chain extension, namely $`R_e`$, is at least twice as large as the Debye screening length. A simple Debye-length criterion appeared to be sufficient to explain the penetration depth of the end-effect. However, looking at the amplitude dependency on density and ionic strength of the solution, we found that both parameters, the number of annealing ions and the ionic strength of the solution, influence the end-effect and that the first one dominated. The amplitude of the end-effect was shown to depend strongly on the charge parameter $`\xi :=l_B/b`$, where $`b`$ is the distance of the bare charges on the backbone of the chain. The definition of such an end-effect via close mobile counterions can not be made for an effective charge $`\xi <<1`$, because under dilute conditions there are almost no counterions close to the chain.

Even though the chain conformation is very different in the poor solvent case the end-effect was found to be qualitatively the same, namely the counterions are more likely to be found at the middle of the chain than at the ends. We could also clearly see the necklace structure by looking at the effective charge along the contour length. However, the string length of our simulated pearl-necklaces was too short to show any charge difference between the pearls and the strings, as has been predicted in Ref. .

We also obtained a fairly good agreement of the simulated ion distribution with the PB solution of the cell model of an infinitely extended charged rod. This supports the idea that the description of polyelectrolytes as rodlike objects in mean-field theory is valid in the dilute regime. Further improvements could probably be achieved along the lines of Ref., where a combination of a cylindrical and spherical cell model is used to describe the solution properties of polyelectrolytes.

### 4.4 Rodlike polyelectrolytes

Stiff linear polyelectrolytes can be approximated by charged cylinders. This is a relevant special case, applying to quite a few biologically important polyelectrolytes with a large persistence length, like DNA, actin filaments or microtubules. Within PB theory  and on the level of a cell model the cylindrical geometry can be treated exactly in the salt-free case , providing for instance new insights into the phenomenon of the Manning condensation . For low line charges, the agreement between PB theory and the simulations of the full interaction system is rather nice. However, PB theory fails quantitatively (underestimated condensation) and qualitatively (overcharging, charge oscillations and attractive interactions); see, e.g. Ref..

Recently the osmotic coefficient of a synthetic stiff polyelectrolyte, a poly(para-phenylene), was measured in a salt-free environment. We have compared this data to predictions of PB theory, and a local density functional theory which includes a correlation correction of the basis of a recently proposed Debye-Hückel-Hole-cavity theory (DHHC) , and simulational results within the cell model. We find that correlation effects enhance condensation and lower the osmotic pressure, yet are not fully able to explain the discrepancy with the experimental data. Here the approach of working within the “primitive model” breaks down. In our opinion, specific interactions between the counterions, the macroion, and the solvent particles are needed to explain the discrepancy. Other theoretical approaches beyond the cell model which try to incorporate finite-size effects and interactions of the macroion itself will in general lead to a higher osmotic coefficient which is in contrast to the experimental data.

Attractive interactions have been observed and predicted between like-charged macromolecules. However, there are nice rigorous results which prove that these effects cannot be described by mean-field theories. Especially in the community of biological inspired physics, these interactions are thought to be important for the clarification of the mechanism behind DNA compactification in viral heads, the chromatin structure, and novel methods for gene delivery, to name just the most prominent examples. There are numerous simulations which show similar attractions on a distance of few counterion diameters .

The mechanism which is driving the observed attractions for rod-like systems has been speculated to be correlations between the counterion layers around the macroion. However, until now, no unique theoretical picture has emerged that can clarify the detailed mechanism behind the attractions. There is the low temperature Wigner crystal theory, initiated by Refs. , which postulates an ordered ground state of the counterions. Then there are theories which are based on Van der Waals type correlated fluctuations , that are in principle hight $`T`$ theories. There are also theories which are fluctuation based, but are valid at low $`T`$. Integral equation  theories on various approximation levels have been demonstrating the existence of these attractions for a long time, but from these theories it is difficult to extract the detailed mechanism behind the observed correlations. Here also simulations can be helpful, because they have in principal access to all correlations. More details of our results in rod-like geometries can be found in Refs..

## 5 The energetic path to understand overcharging

There has been a recent interest in the study of systems which are strongly coupled by Coulomb interactions. These systems show a variety of, at first sight, surprising behaviors, which can not be accounted for by the mean-field PB theory. For example, there are attractions between like charged objects and a charge reversal of macroions occurs when viewed from some distance. This means that there are more ions of the opposite charge within a certain radius around the macroion then necessary to charge neutralize it. This overcompensation is called “overcharging”.

In this section we want to demonstrate that there are situations for charged colloidal objects in which one can understand the phenomenon of overcharging by very simple energetic arguments. By overcharging, in general, we mean that the bare charge of the macroion is overcompensated at some distance by oppositely charged “microions”. To achieve this in nature we have to add salt to the system. For the sake of simplicity, however, we will consider non-neutral systems, because they can on a very simple basis explain why colloids prefer to be overcharged.

### 5.1 The Model

Our model is solely based on electrostatic energy considerations, meaning that we only look at the ground state of a system of charges. We consider a colloid of radius $`a`$ with a central charge $`Z`$. In the ground state the counterions of this colloid are located on the surface, because there they are closest to the central charge. On the other hand they want to be in such a configuration that they minimize their mutual repulsion. For two, three, and four counterions these configurations correspond to a line, an equilateral triangle, and a tetrahedron, respectively, regardless of the central charge magnitude. The problem of the minimal energy configuration of electrons disposed on the surface of a sphere dates back to Thomson , and is actually unsolved for large $`N`$. The reason is, that there are many metastable states which differ only minimally in energy, and their number seems to grow exponentially with $`N`$. Also chemists developed the valence-shell electron-pair repulsion (VSEPR) theory  which uses similar arguments to predict the molecular geometry in covalent compounds, also known as the Gillespie rule.

A simple illustration of energetically driven overcharging is depicted in Fig. 6. The central charge is +2e, and the neutral system has two counterions of valence 1. If we add successively more counterions of the same valence, and put them on the surface such that their mutual repulsion is minimized, we can compute the total electrostatic energy according to

$`E(n)=k_BT(l_B/a)\left[nZ_m+f(\theta _i)\right],`$ (30)
where $`f(\theta _i)`$ is the repulsive energy part which is only a function of the ground state configuration. We surprisingly find, that actually the minimal energy is obtained when four counterions are present, hence we overcharged the colloid by two counterions, or by 100 %! That is, the excess counterions gain more energy by assuming a energetically favorable configuration around the macroion than by escaping to infinity, the simple reason behind overcharging.

In our example, the minimum is reached when four counterions are present. The colloid radius and the Bjerrum length enter as prefactors and change only the energy difference between neighboring states.

The spatial correlations of the counterions are fundamental to obtain overcharging. Indeed, if we apply the same procedure and smear $`Z`$ counterions onto the surface of the colloid of radius $`a`$, we obtain for the energy

$`E=l_B\left[{\displaystyle \frac{1}{2}}{\displaystyle \frac{Z^2}{a}}{\displaystyle \frac{Z_mZ}{a}}\right].`$ (31)
The minimum is reached for $`Z=Z_m`$, hence no overcharging can occur.

The important message to be learned is that, from an energetic point of view, a colloid always tends to be overcharged by discrete charges. Other important geometries like infinite rods or infinitely extended plates cannot be treated in such a simple fashion because they are not finite in all directions. One needs therefore enough screening charges in the environment to limit the range of the interactions in the infinite directions, there is a need for a minimal amount of salt present to allow for overcharging, which is not the case for a colloid.

Obviously, for a large number of counterions the direct computation of the electrostatic energy by using the exact equation (30) becomes unfeasible. Therefore we resort to simulations for highly charged spheres.

### 5.2 One colloid

The electrostatic energy as a function of the number of overcharging counterions $`n`$ is displayed in Fig. 7. We note that the maximal (critical) acceptance of $`n`$ (4, 6 and 8) increases with the macroionic charge $`Z_m`$ (50, 90 and 180 respectively). Furthermore for fixed $`n`$, the gain in energy is always increasing with $`Z_m`$. Also, for a given macroionic charge, the gain in energy between two successive overcharged states is decreasing with $`n`$.

In the ground state the counterions are highly ordered. Rouzina and Bloomfield  first stressed the special importance of these crystalline arrays for interactions of multivalent ions with DNA strands, and later Shklovskii  showed that the Wigner crystal (WC) theory can be applied to determine the interactions in strongly correlated systems. In two recent short contributions  we showed that the overcharging curves obtained by simulations of the ground state, like Fig. 7, can be simply explained by assuming that the energy $`\epsilon `$ per counterion on the surface of a macroion scales as $`\sqrt{c}`$, where $`c`$ denotes the counterion concentration $`c=N/A`$, $`N`$ is the total number of counterions on the surface and $`A`$ the total macroion area. This can be justified by a simple argument, where each ion interacts in first approximation only with the oppositely charged background of its Wigner-Seitz (WS) cell, which can be approximated by a disk of radius $`h`$, yielding the same WS cell area.

For fixed macroion area we can write the energy per counterion as

$`\epsilon ^{(h)}(N)={\displaystyle \frac{\overline{\alpha }^{(h)}\mathrm{}}{\sqrt{A}}}\sqrt{N}=\overline{\alpha }^{(h)}\mathrm{}\sqrt{c},`$ (32)
where $`\mathrm{}=l_BZ_c^2`$ and the simple hole theory gives $`\overline{\alpha }^{(h)}=2\sqrt{\pi }3.54`$.

For an infinite plane, where the counterions form an exact triangular lattice, one obtains the same functional form as in Eq. (32), but the prefactor $`\overline{\alpha }^{(h)}`$ gets replaced by the numerical value $`\overline{\alpha }^{WC}=1.96`$.

Not knowing the precise value of $`\overline{\alpha }`$ we can still use the simple scaling behavior with $`c`$ to set up an equation to quantify the energy gain $`\mathrm{\Delta }E_1`$ by adding the first overcharging counterion to the colloid. To keep the OCP neutral we imagine adding a homogeneous surface charge density of opposite charge ($`\frac{Z_ce}{A}`$) to the colloid. This ensures that the background still neutralizes the incoming overcharging counterion and we can apply Eq. (32). To cancel our surface charge addition we add another homogeneous surface charge density of opposite sign $`\frac{Z_ce}{A}`$. This surface charge does not interact with the now neutral OCP, but adds a self-energy term of magnitude $`\frac{1}{2}\frac{\mathrm{}}{a}`$, so that the total energy difference for the first overcharging counterion reads as

$`\mathrm{\Delta }E_1=(N_c+1)\epsilon (N_c+1)N_c\epsilon (N_c)+{\displaystyle \frac{\mathrm{}}{2a}}.`$ (33)
By using Eq. (32) this can be rewritten as

$`\mathrm{\Delta }E_1={\displaystyle \frac{\overline{\alpha }\mathrm{}}{\sqrt{A}}}\left[(N_c+1)^{3/2}N_c^{3/2}\right]+{\displaystyle \frac{\mathrm{}}{2a}}.`$ (34)
Completely analogously one derives for the energy gain $`\mathrm{\Delta }E_n`$ for $`n`$ overcharging counterions

$`\mathrm{\Delta }E_n={\displaystyle \frac{\overline{\alpha }\mathrm{}}{\sqrt{A}}}\left[(N_c+n)^{3/2}N_c^{3/2}\right]+{\displaystyle \frac{\mathrm{}}{a}}{\displaystyle \frac{n^2}{2}}.`$ (35)
Using Eq. (35), where we determined the unknown $`\overline{\alpha }`$ from the simulation data for $`\mathrm{\Delta }E_1`$, we obtain a curve that matches the simulation data almost perfectly (Fig. 7). The second term in Equation (35) also shows why the overcharging curves of Fig. 7 are shaped parabolically upwards for larger values of $`n`$. If one successively removes each of $`n`$ counterions from a neutral colloid, one can derive in a similar fashion the ionization energy cost

$`\mathrm{\Delta }E_n^{ion}={\displaystyle \frac{\overline{\alpha }\mathrm{}}{\sqrt{A}}}\left[(N_cn)^{3/2}N_c^{3/2}\right]+{\displaystyle \frac{\mathrm{}}{a}}{\displaystyle \frac{n^2}{2}}.`$ (36)
Using the measured value of $`\overline{\alpha }`$ we can simply determine the maximally obtainable number $`n_{max}`$ of overcharging counterions by finding the stationary point of Eq. (35) with respect to $`n`$:

$`n_{max}={\displaystyle \frac{9\overline{\alpha }^2}{32\pi }}+{\displaystyle \frac{3\overline{\alpha }}{4\sqrt{\pi }}}\sqrt{N_c}\left[1+{\displaystyle \frac{9\overline{\alpha }^2}{64\pi N_c}}\right]^{1/2}.`$ (37)
The value of $`n_{max}`$ depends only on the number of counterions $`N_c`$ and $`\overline{\alpha }`$. For large $`N_c`$ Eq. (37) reduces to $`n_{max}\frac{3\overline{\alpha }}{4\sqrt{\pi }}\sqrt{N_c}`$ which was derived in Ref.  as the low temperature limit of a a neutral system in the presence of salt. What we have shown is that the overcharging in this limit has a pure electrostatic origin, namely it originates from the energetically favorable arrangement of the ions around a central charge. We also showed in Ref.  that $`\overline{\alpha }`$ reaches the perfect WC value of $`1.96`$ if the colloid radius $`a`$ gets very large at fixed $`c`$, or when $`c`$ becomes large at fixed $`a`$.

If instead of a central charge scheme one uses discrete charge centers distributed randomly over the colloidal surface we find counterion structures which are quite far away from the WC array, especially when the counterions are pinned to their counter charges. This depends on the interaction energy at contact, which depends of course on $`l`$ and distance of closest approach. However, we still find overcharging, although reduced in value, of the form given by Eq. 35

#### Macroion-counterion interaction profile at $`T=0K`$

The interaction profile between a completely neutralized macroion and one excess counterion is obtained by displacing adiabatically the excess counterion from infinity towards the macroion. From far away the counterion sees only a neutral object and has no measurable interaction, whereas upon approach to the macroion the WC hole gets created in the counterion layer, and we observe a distance dependant attraction towards the macroion. We investigated cases of $`Z_m=2\mathrm{}288`$. All curves can be nicely fitted with an exponential fit of the form

$`E_1(r)=\mathrm{\Delta }E_1e^{\tau (ra)},`$ (38)
where $`\mathrm{\Delta }E_1`$ is the measured value for the first overcharging counterion, and $`\tau `$ is the only fit parameter. In all our results for $`\tau `$ versus $`\sqrt{N_c}`$ we observe a linear dependence for a wide range of values for $`N_c`$, $`\tau \sqrt{N_c}`$, which again can be explained by applying the WC hole picture .

### 5.3 Two Colloids

Now we apply what we have learned about a single colloid to two equal-sized, fixed charged spheres of bare charge $`Q_A`$ and $`Q_B`$ separated by a center-center separation $`R`$ and surrounded by their neutralizing counterions, which give concentrations $`c_A`$ and $`c_B`$, respectively.

All these ions making up the system are immersed in a cubic box of length $`L=80\sigma `$, and the two macroions are held fixed and disposed symmetrically along the axis passing through the centers of opposite faces. This leads to a colloid volume fraction $`f_m=2\frac{4}{3}\pi (a/L)^38.4\times 10^3`$. For finite colloidal volume fraction $`f_m`$ and temperature, we know from the study carried out above that in the strong Coulomb coupling regime all counterions are located in a spherical “monolayer” in contact with the macroion. Here, we investigate the mechanism of strong, long range attraction stemming from monopole contributions; that is, one colloid is overcharged and the other one undercharged.

#### Observation of metastable ionized states

For the charge symmetrical situation we have $`c_A=c_B`$. When we brought this system to room temperature $`T_0`$ and generated initially the counterions randomly inside the box we observed in some cases that one of the colloids remained undercharged, and the other one was overcharged, and these configurations turned out to be extremely long lived in the course of our MD simulations( more than $`10^8`$ MD time steps). However it is clear that such a state is “metastable” because by symmetry arguments it cannot be the lowest energy state. The observed barrier is the result of the WC attraction, because close to the macroion surface the energy is reduced. For very distant macroions the barrier height for the first overcharged state has to equal $`\mathrm{\Delta }E_1`$ from Eq. 35. The barrier profile at $`T=0`$ can also be extremely well approximated by an application of Eqs. (35) and (36), plus taking into account the distant dependent monopole contribution. This leads to a barrier height which scales as $`\sqrt{c}`$ for large separations. For smaller separations one has to take into account also the effect of strong mutual polarization of both macroions, which leads effectively to a sharing of their proximal counterion layer into a superlattice. This can be taken into account by a higher effective counterion density close to the surface, leading to an almost linear scaling of the barrier height with $`c`$ .

#### Asymmetrically charged colloids

The most interesting phenomenon, however, appears when the two colloids have different counterion concentrations, here $`c_A>c_B`$, since then stable ionized states can appear. The physical reason is that a counterion can gain more energy by overcharging the colloid with $`c_A`$ then it loses by ionizing colloid $`B`$. A straight forward application of the procedure outlined for the barrier calculation  yields a simple criterion (more specifically a sufficient condition), valid for large macroionic separations, for the charge asymmetry $`\sqrt{N_A}\sqrt{N_B}`$ to produce an ionized ground state of two unlike charged colloids with the same size:

$`\left(\sqrt{N_A}\sqrt{N_B}\right)>{\displaystyle \frac{4\sqrt{\pi }}{3\overline{\alpha }^A}}1.2.`$ (39)
#### Finite temperature analysis

We have also demonstrated that the ground state phenomena survive for finite temperatures, i.e. an ionized state can also exist at room temperature $`T_0`$. The left part of Figure 8 shows the time evolution of the electrostatic energy of a system $`Z_A=180`$ with $`Z_B=30`$, $`R/a=2.4`$ and a colloidal volume fraction of $`710^3`$, where the starting configuration is the neutral state (DI = 0). One clearly observes two jumps in energy, $`\mathrm{\Delta }E_1=19.5`$ and $`\mathrm{\Delta }E_2=17.4`$, which corresponds each to a counterion transfer from colloid B to colloid A. These values are consistent with the ones obtained for the ground state, which are$`20.1`$ and $`16.3`$ respectively. Note that this ionized state (DI = 2) is more stable than the neutral but is expected to be metastable, since it was shown previously that the most stable ground state corresponds to DI = 5. The other stable ionized states for higher DI are not accessible with reasonable computer time because of the high energy barrier made up of the correlational term and the monopole term which increase with DI. In the right part of Fig. 8 we display a typical snapshot of the ionized state (DI = 2) of this system at room temperature.

Obviously, these results are not expected by the DLVO theory even in the asymmetric case (see e. g. ). Previous simulations of asymmetric (charge and size) spherical macroions  were also unable to predict such a phenomenon since the Coulomb coupling was weak (water, monovalent counterions). Note that the appearance of (meta-)stable ionized states can alter the effective interactions between charged colloids in solution. The monopole attraction will lead to attraction between like charged colloids, flocculation, and related phenomena.

At this stage, we would like to stress again, that the appearance of a stable ionized ground state is due merely to correlation. An analogous consideration with smeared out counterion distributions along the lines of Eq. (31) will again always lead to two colloids exactly neutralized by their counterions . Our energetic arguments are quite different from the situation encountered at finite temperatures, because in this case even a PB description would lead to an asymmetric counterion distribution. However, in the latter case this happens due to purely entropic reasons, namely in the limit of high temperatures, the counterions want to be evenly distributed in space, leading to an effective charge asymmetry.

Note also, that there can exist parameter regions, such as high molar electrolytes, where the overcharging of a single macroion is due to mainly entropic effects , whose exact mechanism is currently under investigation.

## 6 Acknowledgments

We gratefully acknowledge collaborations at various stages with A. Arnold, J. DeJoannis, M. Deserno, H.J. Limbach, R. Messina, U. Micka, and Z. Wang.