text:
# Correlated Trapped Bosons and the Many-Body Efimov Effect

## Introduction.

The novel theoretical formulation in sor01  was constructed to describe correlations in boson systems and applied for a realistic short-range repulsive interaction to Bose-Einstein condensates. The method goes beyond the mean-field approximation dal99 ; bay96  and is as well applicable to attractive finite range potentials with very large scattering lengths where the Efimov effect occurs efi70 ; fed93 . Experimental properties are available for various condensates from the first observations of both effectively repulsive and95  and attractive bra95  interatomic interactions to the recent observations cor00 ; rob01  of condensates in a magnetic field used to tune the effective interaction via a Feshbach resonance to almost any scattering length.

The condensate dal99  being an excited state of the full many-body system is clearly unstable. The three-body recombination into bound dimers is a dominating decay channel nie99 ; esr99 ; bra02b . In the condensate this occurs independently for neighbouring pairs and much more frequently in a coherent process best described as a collective or macroscopic collapse pit96 ; adh02 . Two-body correlations therefore must be crucial for this collapse, which becomes more likely and eventually inevitable as the scattering length is increased. Unfortunately a theoretical description is hindered by the difficulties, especially pronounced for large scattering lengths, of finding the decisive correlated structure dal99 ; blu01 ; cow01 .

A promising form of a correlated wave function suggested for nucleons rip84  was recently extended to more general systems bar99a . Another related formulation uses generalized hyperspherical coordinates and an adiabatic expansion with the hyperradius as the adiabatic coordinate boh98 . It was applied for many identical bosons, still only with a zero-range interaction and the lowest (constant and thus non-correlated) hyperspherical angular wave function. These crude approximations are removed in a novel method still using hyperspherical coordinates and adiabatic expansion, but now the wave function is determined from a variationally established equation sor01 . The purpose of this letter is to investigate the structure of boson systems as a function of (large) scattering length for attractive finite range potentials with emphasis on an emerging novel many-body Efimov effect.

## Theory.

A system of $`N`$ identical particles of masses $`m`$, trapped in an external field approximated by a harmonic oscillator potential of angular frequency $`\omega `$, can in the center of mass frame be described by hyperspherical coordinates, i.e. $`3N4`$ hyperangles bar99a ; bar99b  denoted collectively by $`\mathrm{\Omega }`$ and one hyperradius, $`\rho `$, given by sor01

$`\rho ^2`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{i<j=1}{\overset{N}{}}}r_{ij}^2={\displaystyle \underset{i=1}{\overset{N}{}}}r_i^2NR^2,`$ (1)
where $`\stackrel{}{r}_i`$ are single-particle coordinates, $`\stackrel{}{R}`$ is the center of mass coordinate, and $`r_{ij}=|\stackrel{}{r}_j\stackrel{}{r}_i|\sqrt{2}\rho \mathrm{sin}\alpha _{ij}`$ with $`\alpha _{ij}`$ varying between $`0`$ and $`\pi /2`$. The hamiltonian is

$`\widehat{H}`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{N}{}}}\left({\displaystyle \frac{\widehat{p}_i^2}{2m}}+{\displaystyle \frac{1}{2}}m\omega ^2r_i^2\right)+{\displaystyle \underset{i<j=1}{\overset{N}{}}}V(r_{ij}),`$ (2)
where $`V`$ is the two-body interaction. It separates into a center of mass part, a radial part, and an angular part depending respectively on $`\stackrel{}{R}`$, $`\rho `$, and $`\mathrm{\Omega }`$ sor01 :

$`\widehat{H}=\widehat{H}_{\text{cm}}+\widehat{H}_\rho +{\displaystyle \frac{\mathrm{}^2\widehat{h}_\mathrm{\Omega }}{2m\rho ^2}},`$ (3)
with $`\widehat{H}_{\text{cm}}=\widehat{p}_R^2/(2Nm)+Nm\omega ^2R^2/2`$, $`\widehat{H}_\rho =\widehat{T}_\rho +m\omega ^2\rho ^2/2`$ and $`\mathrm{}^2\widehat{h}_\mathrm{\Omega }/(2m\rho ^2)=\widehat{T}_\mathrm{\Omega }+_{i<j}V_{ij}`$, where $`\widehat{T}_\rho `$ and $`\widehat{T}_\mathrm{\Omega }`$ are related kinetic energy operators. A suitable expansion of the wave function is

$`\mathrm{\Psi }`$ $`=`$ $`\rho ^{(3N4)/2}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}f_n(\rho )\mathrm{\Phi }_n(\rho ,\mathrm{\Omega }),`$ (4)
where $`\mathrm{\Phi }_n`$ is an eigenfunction of the angular part of the hamiltonian with an eigenvalue $`\mathrm{}^2\lambda _n(\rho )/(2m\rho ^2)`$:

$`\widehat{h}_\mathrm{\Omega }\mathrm{\Phi }_n(\rho ,\mathrm{\Omega })`$ $`=`$ $`\lambda _n(\rho )\mathrm{\Phi }_n(\rho ,\mathrm{\Omega }).`$ (5)
Neglecting couplings between the different $`n`$-channels yields the radial eigenvalue equation for the energy $`E_n`$:

$`\left({\displaystyle \frac{\mathrm{}^2}{2m}}{\displaystyle \frac{d^2}{d\rho ^2}}+U_n(\rho )E_n\right)f_n(\rho )`$ $`=`$ $`0,`$ (6)
$`U_n(\rho )={\displaystyle \frac{m\omega ^2\rho ^2}{2}}+{\displaystyle \frac{\mathrm{}^2(3N4)(3N6)}{8m\rho ^2}}`$ $`+`$ $`{\displaystyle \frac{\mathrm{}^2\lambda _n}{2m\rho ^2}}.`$ (7)
The second term in the radial potential $`U`$ is a generalized centrifugal barrier. We now decompose the angular wave function $`\mathrm{\Phi }`$ in the symmetric Faddeev components $`\varphi `$

$`\mathrm{\Phi }(\rho ,\mathrm{\Omega })={\displaystyle \underset{i<j=1}{\overset{N}{}}}\varphi _{ij}(\rho ,\mathrm{\Omega }){\displaystyle \underset{i<j=1}{\overset{N}{}}}\varphi (\rho ,r_{ij}),`$ (8)
where the restricting assumption is that the interparticle potentials only act in $`s`$-waves leaving only the dependence on the distance $`r_{ij}=\sqrt{2}\rho \mathrm{sin}\alpha _{ij}`$. The capability of this decomposition for large scattering length has been demonstrated for $`N=3`$ by the proper description of the intricate Efimov effect fed93 ; jen97 .

The eigenvalue equation eq. (5) can by a variational technique be rewritten as a second order integro-differential equation in the variable $`\alpha _{12}`$ sor01 . For atomic condensates the interaction range is very short compared to the spatial extension of the $`N`$-body system. Then this equation simplifies even further to contain at most one-dimensional integrals. The principal interaction dependence is through the parameter $`a_\text{B}m_0^{\mathrm{}}𝑑rr^2V(r)/\mathrm{}^2`$, which is the Born-approximation to the $`s`$-wave scattering length $`a_\text{s}`$. The validity of our approximations only relies on a small *range* of the potential whereas the scattering length can be as large as desired.

## Angular potentials.

The crucial ingredient is the angular eigenvalue $`\lambda `$ obtained from eq. (5). We choose gaussian two-body potentials $`V(r)=V_0\mathrm{exp}(r^2/b^2)`$ of range $`b`$ and strength $`V_0`$. We measure lengths in units of $`b`$ and energies in units of $`\mathrm{}^2/(2mb^2)`$. We solve the equation by choosing discrete sets of mesh points and thus constructing a matrix equation ready for diagonalization sor01 . This computational scheme increases strongly with $`N`$ and we therefore illustrate by relatively small values of $`N`$.

Fig. 1 shows the angular eigenvalues for different strengths expressed in terms of $`a_\text{B}`$ and $`a_\text{s}`$. For short range potentials finite at the origin all $`\lambda `$ approach their values at $`\rho =0`$ as $`\rho ^2`$. This is a model dependent region and therefore not considered in the following. When the potential is too weak to support a two-body bound state, ($`a_\text{B}/b>1.18934`$), the lowest $`\lambda _0`$ is negative and approaches zero as $`a_\text{s}/\rho `$ with increasing $`\rho `$. In fig. 1 we also show the eigenvalue $`\lambda _1`$ corresponding to the first excited angular state. It is positive and approaches the value $`4(4+3N5)`$ as $`a_\text{s}/\rho `$ for large $`\rho `$ jen97 ; boh98 .

Increasing the attraction of the potential a two-body bound state of energy $`B_2`$ appears (positive $`a_\text{s}`$) and the corresponding lowest angular eigenvalue diverges proportional to $`2mB_2\rho ^2/\mathrm{}^2`$ when $`\rho \mathrm{}`$. The many-body system would obviously then prefer to recombine into bound dimer states. The first excited angular eigenvalue, similar to the dotted curve in fig. 1, is then responsible for the unstable state referred to as the condensate.

At the two-body threshold the scattering length is infinite. Then the lowest $`\lambda `$ approaches a negative constant $`\lambda _{\mathrm{}}`$ for large $`\rho `$ ($`\rho \rho _b45\sqrt{N}b`$), see fig. 1. Using $`N=3,10,20`$ we find numerically $`\lambda _{\mathrm{}}(N)24(N2)5N(N1)`$, where the form reflects our expectations of $`N`$-dependence of kinetic and potential energies.

These angular eigenvalues behave qualitatively different from those obtained in boh98 , where the expectation value of a $`\delta `$-interaction of strength proportional to $`a_\text{s}`$ is computed with a constant angular wave function without any correlations. The eigenvalues are then necessarily proportional to $`a_\text{s}/\rho `$ for all $`\rho `$, i.e. always diverging when $`\rho 0`$ and converging to zero for $`\rho \mathrm{}`$. The differences are especially pronounced for large scattering lengths and when a bound two-body state is present. However, even in the qualitatively similar case in fig. 1 of $`a_\text{s}/b=15.7`$ our $`\lambda `$ is lower than that of boh98  by a factor between $`1.5`$ and $`2`$ for $`\rho /b[10^3,10^5]`$.

At large $`\rho `$, $`\mathrm{\Phi }`$ approaches a non-correlated angular wave function. Using a constant as in boh98  the expectation value of our gaussian interaction gives for large $`\rho `$ and $`N`$ that $`\lambda =\lambda _{\text{n-c}}(\rho )kN^{7/2}a_\text{B}/\rho `$ with $`k=(3/2)^{3/2}/\sqrt{2\pi }`$, i.e. we get $`a_\text{B}`$ instead of $`a_\text{s}`$ as assumed in boh98 . For large $`a_\text{s}`$ our computed $`\lambda `$ bends upwards from the plateau of $`\lambda _{\mathrm{}}`$ around a point $`\rho \rho _{\text{th}}`$ determined by the large but finite value of $`a_\text{s}`$. If $`\rho >\rho _{\text{th}}`$ then $`\lambda `$ can be estimated as $`\lambda _{\text{n-c}}(\rho )`$ with $`a_\text{s}`$ substituted for $`a_\text{B}`$. The transition point is then given by $`\lambda _{\mathrm{}}=\lambda _{\text{n-c}}(\rho _{\text{th}})`$ which for large $`N`$ yields

$`\rho _{\text{th}}={\displaystyle \frac{1}{5}}k|a_\text{B}|N^{3/2}{\displaystyle \frac{1}{7}}|a_\text{B}|N^{3/2}{\displaystyle \frac{1}{7}}|a_\text{s}|N^{3/2},`$ (9)
where the more realistic estimate for attractive potentials is indicated by exchanging $`a_\text{B}a_\text{s}`$, see jen97 ; boh98 .

## Radial potentials.

In eq. (7) the external harmonic trap of angular frequency $`\omega `$ corresponds to a trap length $`b_\text{t}\sqrt{\mathrm{}/(m\omega )}`$. We model the experimentally studied systems and95 ; rob01  of 85Rb and 87Rb-atoms with oscillator frequencies $`\nu =\omega /(2\pi )=205`$ Hz and $`200`$ Hz and interaction range $`b=10`$ a.u., thus yielding $`b_\text{t}/b=1442`$.

We show in fig. 2 the radial potentials for the angular eigenvalues corresponding to unbound two-body systems ($`a_\text{s}/b=15.7`$ in fig. 1). The external field in eq. (7) is negligible for small $`\rho `$ and therefore the radial potential is negative when $`\lambda +(3N4)(3N6)/4<0`$. Then genuinely bound states of negative energy are possible in our model without the confinement from the trap.

The radial potential for the lowest angular eigenvalue has a global minimum with $`U_{\text{min}}<0`$ at small $`\rho `$ ($`\rho _b`$) separated by a barrier at intermediate $`\rho `$ from a local minimum at $`\rho \rho _\text{t}\sqrt{3N/2}b_\text{t}`$. At large $`\rho `$ the radial potential diverges due to the harmonic term, at small $`\rho `$ due to the centrifugal term. The radial potential corresponding to the first excited angular state, also in fig. 2, coincides with the lowest radial potential at large $`\rho `$, but does not contain the attraction at small $`\rho `$ (see fig. 1) and diverges therefore to $`+`$ infinity for small $`\rho `$ without going through another minimum.

The radial equation in eq. (6) has only one solution with negative energy $`E`$ and the corresponding wave function, shown in fig. 2, is located in the global minimum. The first of the infinitely many excited states in this potential is located in the local minimum at larger $`\rho `$ created by the competition between the centrifugal barrier and the external harmonic oscillator potential. This excited state is usually referred to as the condensate, but the moderate attraction for the relatively small $`a_\text{s}=15.7b`$ already results in a lower lying many-body state about $`\rho _\text{t}/\rho _b37`$ times smaller than the condensate.

## Large scattering length.

At the threshold for binding of the two-body system the angular eigenvalue outside the interaction range must approach a negative constant $`\lambda \lambda _{\mathrm{}}`$, see fig. 1. The scattering length is then infinitely large and the radial potential has the form nie01

$`U(\rho )`$ $``$ $`{\displaystyle \frac{\mathrm{}^2}{2m}}\left({\displaystyle \frac{\xi ^21/4}{\rho ^2}}+{\displaystyle \frac{\rho ^2}{b_\text{t}^4}}\right),`$ (10)
where $`\xi ^2\lambda _{\mathrm{}}(3N4)(3N6)/41/411N^2/4`$ for large $`N`$, where $`\lambda _{\mathrm{}}5N^2`$. Without the external $`\rho ^2`$ potential in eq. (10) the infinitely many radial solutions to eq. (7) all behave like

$`f_{\mathrm{}}(\rho )=\sqrt{\rho }\mathrm{sin}\left(|\xi |\mathrm{ln}(\rho /\rho _{\text{sc}})\right),`$ (11)
with some hyperradius scale $`\rho _{\text{sc}}`$. The energies and sizes of the eigenstates labeled $`j`$ are related by efi70 ; nie01

$`{\displaystyle \frac{E_j}{E_{j+1}}}`$ $`=`$ $`{\displaystyle \frac{\rho ^2_{j+1}}{\rho ^2_j}}=e^{2\pi /|\xi |}.`$ (12)
With increasing quantum number these states become exponentially larger with exponentially smaller energies approaching zero. They originate from the constant $`\lambda `$ and the generic $`1/\rho ^2`$ potential in eq. (10). They are many-body states, but not the embedded three-body Efimov cluster states which also arise from a $`1/\rho ^2`$ potential. However, both multi- and three-body Efimov states appear for very large two-body scattering lengths.

The external harmonic oscillator limits the possible number of these new states with $`E<0`$. These negative energy states have to be located inside the trap and outside the two-body potential. The number of nodes of $`f_{\mathrm{}}`$ allowed in this region equals the number of bound states $`N_E|\xi |\pi ^1\mathrm{ln}(\rho _{\text{max}}/\rho _{\text{min}})`$, where the first and last zero points then are at the end points of the interval. We obtain $`\rho _{\text{max}}/\rho _{\text{min}}\rho _\text{t}/\rho _b`$ (see also fig. 3)

$`N_E{\displaystyle \frac{|\xi |}{\pi }}\mathrm{ln}\left({\displaystyle \frac{b_\text{t}}{37b}}\right),\text{for }b_\text{t}N|a_\text{s}|,`$ (13)
where the quoted condition of validity is found from $`\rho _\text{t}\rho _{\text{th}}`$, i.e. the plateau value must extend beyond the trap. The number of these available Efimov states therefore scales proportional to $`N`$ ($`\xi N`$) and logarithmically with the ratio of trap length and interaction range.

We show in fig. 3 the radial potential for $`N=20`$ and infinite scattering length corresponding to $`\lambda _{\mathrm{}}1340`$ or $`\xi ^2584`$, see fig. 1. Using the estimate in eq. (13) we get $`N_E28`$ in good agreement with the computed 30. The lower curve in the inset of fig. 3 is according to eq. (12) a straight line of slope $`1`$. The very lowest states deviate, since they “feel” non-constant $`\lambda `$ at small $`\rho `$, and the states close to $`E=0`$ deviate due to the influence of the external potential. The energy spectrum becomes even denser above zero energy. For large positive energies in the upper part of the inset the harmonic potential dominates and a straight line with slope $`+1`$ is obtained. The small positive energies are influenced by both external trap and interaction potentials.

Using eq. (1) we get $`2\rho ^2=(N1)r_{12}^22(N1)r_i^2`$. Even the most bound state, $`\rho ^2^{1/2}136b`$, has then a root mean square distance between two particles, $`r_{12}^2^{1/2}44b`$, much larger than the interaction range. Also the root mean square radius $`r_i^2^{1/2}31b`$ is large.

At the plateau, $`\lambda \lambda _{\mathrm{}}`$, $`\xi ^2>0`$ and the radial potential has no intermediate barrier. The centrifugal barrier can only be at $`\rho >\rho _{\text{th}}`$, i.e. in the region where $`\lambda \lambda _{\text{n-c}}`$, which provides the estimate of the barrier position $`\rho 2kN^{3/2}|a_\text{s}|/3>\rho _{\text{th}}`$, see eq. (9). A criterion of stability is then obtained by using $`\lambda _{\text{n-c}}`$ to estimate when the corresponding radial barrier disappears. This yields in analogy to boh98  that the condensate will be (meta-) stable if $`N|a_\text{s}|/b_\text{t}<0.671`$. This critical value is from the Gross-Pitaevskii equation found to be about $`0.6`$ dal99 ; rup95 , and recently measured to $`0.46`$ rob01  for a 85Rb-condensate.

Stability is determined by decay through a barrier, like in fig. 2, when negative energy states and the condensate state in the second minimum are present. The tunneling process then populates the states in the first minimum, i.e. many-body states of smaller extension and correspondingly larger density than the initial condensate. Competing with this macroscopic (all particles together) tunneling are the two- and three-particle recombination processes, which can occur from the initial condensate, or more likely from the collapsed many-body states at larger density and smaller $`\rho `$. Thus we have established a decay mechanism, i.e. macroscopic tunneling followed by recombination processes from the intermediate model independent states.

An even more dramatic collapse could be initiated experimentally by first creating a condensate corresponding to fig. 2, and then suddenly increasing the scattering length to that of fig. 3 rob01 ; adh02 ; don01 . The multi-particle Efimov states of fig. 3 would then quickly be populated, since the barrier is removed, and these states would subsequently leak out of the trap due to recombination.

## Conclusions.

We use a new formulation designed to investigate two-body correlations in boson systems. The qualitative properties of Bose-Einstein condensates are explained by using the hyperspherical adiabatic expansion. Finite range and attractive potentials introduce new features compared to the mean-field approximation. Large scattering lengths are specifically treated. Many-body Efimov states appear for scattering lengths of the order of the trap length as model-independent states of much smaller extension than the condensate. Recombination through these intermediate states is suggested. This opens the possibility for realistic calculations of few-body recombination processes and for macroscopic collapse processes.


format: text

text:
# Kilo-Second Quasi-Periodic Oscillations in the Cataclysmic Variable DW Cancri

## 1 Introduction

Cataclysmic variables (CVs) are semi-detached binaries which contain a white dwarf and a late-type secondary star (Warner (1995)). The overflowing gas from the Roche lobe-filling secondary star forms an accretion disk, which plays a key role in various activities of CVs. From the photometric and spectroscopic characteristics, the systems are classified with some sub-groups, for example, classical novae, dwarf novae (DNe), nova-like variables (NLs), and magnetic CVs of polars and intermediate polars (IPs).

Quasi-periodic oscillations (QPOs) have been detected in CVs on a variety of time scales. Rapid QPOs with a period of a few seconds have been observed in some magnetic CVs. They are considered to originate from the shock region between a white dwarf and accretion streams (Langer et al. 1981, 1982). During DN outbursts, we know three types of temporary short-term modulations, that is, dwarf nova oscillations (DNOs), QPOs, and super-QPOs. DNOs are highly coherent modulations with typical periods of a few tens of seconds. The period of QPOs during DN outbursts is generally longer than that of DNOs ($`P_{\mathrm{QPO}}=30`$$`1000`$ s). DNOs and QPOs have been reported also in some NLs and classical novae (Patterson (1981)). The origin of DNOs has been proposed to be a white dwarf (Warner, Robinson (1972); Papaloizou, Pringle (1978)), an inner edge of the disk (Bath (1973)) or a boundary layer (Okuda et al. (1992); Collins et al. (2000)); however, it is still an issue. On the other hand, the oscillation of a disk or a boundary layer has been suggested for QPOs (Okuda et al. (1992); Collins et al. (2000)). While the amplitudes are quite small ($`0.05\mathrm{mag}`$) in DNOs and QPOs, super-QPOs have higher amplitudes, reaching $`0.2\mathrm{mag}`$. This modulation, whose period is typically a few minutes, has been detected only in the early and late phase of superoutbursts in some SU UMa-type DNe (Kato et al. (1992); Nogami et al. (1998); Kato (2002)).

Stepanian (1982) discovered new variable stars in the course of their search for galaxies with a low-dispersion objective prism survey. According to Stepanian (1982), the “No. 2” object showed brightness variations of $`V=15`$$`17.5`$. In conjunction with strong Balmer emission lines and a weak He ii 4686 Å emission line, Stepanian (1982) suggested that it is a DN. This object was called FBS 0756+164 by Kopylov et al. (1988) and given a general GCVS name of DW Cnc (Kholopov et al. (1981)). Kopylov et al. (1988) supported this classification based on a new spectrum whose characteristics are analogous to that reported by Stepanian (1982).

Here, we report on our time-series photometric monitoring of DW Cnc, in which we discovered unprecedentedly large-amplitude and long-period QPOs. In the next section, we describe our observation method. Our observation results are presented in section 3; we then propose possible interpretations of them in section 4. We summarize our findings in the final section.

## 2 Observation

We performed CCD photometric observations at Kyoto, Ouda Observatory, Nicholas Copernicus Observatory and Planetarium, and Erftstadt. Our equipment and observation log are listed in tables 1 and 2, respectively. After dark subtraction and flat fielding on our obtained images, the differential magnitudes of the object were calculated with point spread function (PSF) photometry for images taken at Kyoto. Aperture photometry was performed for images taken at the other sites. The magnitude scales of each observatory were adjusted to that of the Kyoto system, in which we used a comparison star, GSC 1363.2124. The constancy of the comparison brightness was checked by GSC 1363.2014. We could obtain magnitudes almost equal to the $`R_\mathrm{c}`$ system with an unfiltered ST-7 camera, since its sensitivity peak is almost equal to that of $`R_\mathrm{c}`$ system and the object has $`BV0`$. Heliocentric corrections to the observed times were applied before a following analysis.

## 3 Result

### 3.1 Long-Term Light Curve and Re-Classification of DW Cnc

Figure 1 shows long-term light curves of DW Cnc. The crosses with error bars represent the average magnitudes with standard errors each night. The light curve in the upper panel includes all visual observations reported to Variable Star Network (VSNET111$``$http://www.kusastro.kyoto-u.ac.jp/vsnet/$``$), which are depicted with the open circles and the triangles (upper-limits). The lower panel is a light curve focused on our CCD observations. As can be seen in the upper panel, the object varied between 14–15 mag most of the time. During our observations, shown in the lower panel, the object showed no significant fading or brightening trend. We calculated the average magnitude to be $`R_\mathrm{c}=14.68`$.

Stepanian (1982) proposed that DW Cnc is a DN with the variation range of $`V=15`$$`17.5`$. The characteristic which we observed, as shown in figure 1, is, however, inconsistent with that of the DN light curves. Z Cam stars, which form a sub-group of DNe, are known to experience a phase called “standstill”, which is a phase with constant brightness interrupting ordinary outburst cycles (Warner (1995)). It is, however, less likely that the period of constant brightness of DW Cnc is a standstill, because the long monitoring of visual observations shows no ordinary outburst phase over six years. It may be possible that the long bright state is indeed a quiescent state of a dwarf nova having a long recurrence time, such as WZ Sge stars. This scenario, however, apparently contradicts the presence of the faint phase, which is indicated by a few observations fainter than 15 mag in figure 1 and the originally proposed amplitude of $`V=15`$$`17.5`$ (Stepanian (1982)).

The long-term variations of DW Cnc are analogous to NLs or magnetic CVs, rather than DNe. Figure 1 indicates that the duration of the possible faint state is probably shorter than that of the high state. Although we could obtain no estimation of the magnitude at the time when the spectroscopy was performed by Stepanian (1982) and Kopylov et al. (1988), it is less likely that both of the spectra were taken during the temporal short-duration faint states. The case of a magnetic CV is thus apparently inconsistent with the lack of a strong He ii emission line (Stepanian (1982); Kopylov et al. (1988)). We hence suggest that DW Cnc is neither a DN nor a magnetic system, but a NL.

We cannot conclude whether DW Cnc has changed its nature as a DN or it was just a mis-classification, since no light curve is available in Stepanian (1982). On the other hand, the hints of possible faint states are reminiscent of the long-term variations of VY Scl-type stars, a sub-group of NLs showing temporary low states (Leach et al. (1999)). The variations between 14–15 mag seen in the visual observations can be explained by short-term variations with relatively large amplitudes, which we describe in the following section.

### 3.2 Discovery of Kilo-Second Quasi-Periodic Oscillations

During our CCD monitoring, we detected oscillations whose typical light curves are shown in figure 2. As shown in this figure, they appear to have no constant period, but variable profiles and a quasi-periodic nature. They have a coherence time of at least a few cycles, as can be seen in figure 2. We cannot find any firm evidence of a long-term coherence of these variations within our limited available data due to the short base-lines of each observation. Through a period analysis, we confirmed that no clear coherent peak exists in the power spectrum between periods of 1–200 min. We show a wide-range power spectrum in figure 4. In this figure, the abscissa and the ordinate denote the frequency in cycle day-1 and the power in arbitrary units, respectively. We can find two prominent peaks with wide widths, which demonstrate the QPO nature of the oscillations. We determined their center periods to be $`37.5\pm 0.1`$ and $`73.4\pm 0.4`$ min by fitting Gaussian functions. The shorter QPO period is probably the first harmonic of the longer one; however, due to the small difference between them, we cannot completely exclude the possibility that they are periodicity independent. We show the average light curve of modulations folded with the 73.4-min period in figure 3. We can see the main peak around the phase $`0.8`$ and the sub peak around the phase $`0.4`$.

The QPO in DW Cnc is unique regarding its long lifetime, large amplitude and long period of kilo-second. QPOs in CVs generally have periods of over one order of magnitude shorter than that in DW Cnc (Patterson (1981)). Although kilo-second QPOs have been reported in some CVs (e.g. Udalski (1988); Kraicheva et al. (1999)), their amplitudes are quite small ($``$ a few percent) compared with those in DW Cnc. The amplitude of QPOs in DW Cnc is typically about 0.3 mag and, furthermore, we detected a large flare-like peak whose amplitude reached to 0.6 mag, as can be seen in the upper-right panel of figure 2. We clearly detected the oscillation, or a part of it, in our all observations with a long time-coverage, as shown in figure 2. Although the light curves in the lower-middle and right panel in figure 2 are more noisy due to the bad seasonal condition, they show a clear presence of oscillations with amplitudes and period similar to those in the upper panels. It is thus probable that the QPO has appeared throughout our observations for 61 days.

Except for the two QPOs, the power spectrum shows other noteworthy characteristics. It indicates the presence of a significant power for frequencies lower than 10 $`\mathrm{cycle}\mathrm{d}^1`$. The power suddenly weakens at frequencies higher than the 37.5-min QPO. These features of the power spectrum are qualitatively analogous to that in black hole binaries (Rutledge et al. (1999)). The poisson noise dominates for frequencies higher than 100 cycle d-1, where the spectrum becomes flat.

## 4 Discussion

In this section, we compare DW Cnc with the other CVs which also show similar oscillations to evaluate whether their proposed models can be applied to the QPO in DW Cnc.

### 4.1 Spin of White Dwarf

Some IPs show periodic variations on the order of kilo-second (e.g. Schrijver et al. (1985); Patterson, Price (1981)), which are interpreted as being the rotational period of a white dwarf. In the case of GK Per, the long-period QPO ($`5000\mathrm{s}`$) reported in Watson et al. (1985) is proposed to be a beat phenomenon between the spin of a white dwarf and blobs orbiting at a disk-overflow impact site (Hellier, Livio (1994)) or between the accreting curtain and ordinary DN QPOs (Morales-Rueda et al. (1999)). The long period seen in DW Cnc may prefer a beat scenario rather than a direct spin variation. The most essential objection against this model is that DW Cnc shows no direct evidence for an IP. The two spectra obtained by Stepanian (1982) and Kopylov et al. (1988) show no strong He ii emission and, moreover, we detected no firm evidence of a short periodicity originating from the white-dwarf spin. It may be possible that the shorter periodicity of 37.5-min is not the first harmonic of the longer periodicity, but is associated with the spin period. In this case, the kilo-second QPO can be interpreted with the beat between the white-dwarf spin and $`25`$-min QPO. Within our available data, the broad profile of the 37.5-min peak in the power spectrum, however, supports that this also has the QPO nature, and probably corresponds to the first harmonic of the 73.4-min QPO.

### 4.2 Disk Oscillation

Thin-disk oscillations have been proposed to explain the characteristics of QPOs both in DNe and X-ray binaries (Kato (2001)). In a standard disk during DN outbursts, the oscillation wave propagating inward is reflected at the radius where the wave frequency is equal to the epicyclic frequency. Since the wave cannot propagate in an outer cool disk, the wave is trapped in a narrow region. The wave is thus excited enough to be detected as QPOs in DNe (Yamasaki et al. (1995); Yamasaki, Kato (1996)). Although DW Cnc is definitely not a DN, we can expect the condition of an accretion disk analogous to that during DN outbursts, since the object is most likely to be a NL, which is generally believed to have a sufficiently high mass-accretion rate to form a fully-ionized disk (Osaki (1996)). The Balmer emission lines, furthermore, indicate the presence of an optically thin plasma around the disk. In the case of DN outbursts, the QPO frequency is theoretically expected to become higher with time due to an inward propagation of the transition layer (Yamasaki et al. (1995)). On the other hand, the model can explain the long lifetime of the QPO in DW Cnc with no apparent period shortening, since we can expect a quasi-steady disk system in the case of NLs. Kraicheva et al. (1999) reported 50-min QPOs in a NL, MV Lyr and also suggest applying of the trapped disk oscillation model.

According to Yamasaki et al. (1995), the trapped oscillation model provides the expected QPO period, which is the inverse of the Keplerian frequency at the transition layer, that is, $`P_{\mathrm{QPO}}600\mathrm{s}\times (r_{\mathrm{front}}/10^{10}\mathrm{cm})^{3/2}(M/M_{})^{1/2}`$. With the $`73\mathrm{min}`$ QPO, we estimate $`r_{\mathrm{front}}=3.8\times 10^{10}\mathrm{cm}`$, assuming a white-dwarf mass of $`1M_{\mathrm{solar}}`$. This $`r_{\mathrm{front}}`$ is reasonable in CVs with an orbital period longer than about three hours. The disk-oscillation scenario can thus be an acceptable model for the QPO in DW Cnc if its orbital period is above the period gap. This scenario, however, was originally proposed to explain the ordinary DN QPO, which has amplitudes typically smaller than $`0.1\mathrm{mag}`$ (e.g. Yamasaki et al. (1995)). While Kraicheva et al. (1999) proposed this model for $`0.2\mathrm{mag}`$ QPOs observed in MV Lyr, the amplitude reaching $`0.3\mathrm{mag}`$ in DW Cnc may be too large for this model to reproduce.

### 4.3 Superhump

The large amplitude and the relatively long period are rather analogous to superhumps which are observed during superoutbursts in SU UMa-type DNe (Warner (1985)). Superhumps are modulations which generally have a period a few percent longer than an orbital period. According to the tidal instability model, this phenomenon is caused by the gradual precession of an eccentric disk which is developed by an enhanced tidal effect at the 3:1 resonance (Osaki (1989)). As expected from this model, superhumps have been reported not only in SU UMa-type DNe, but also in NLs and AM CVn stars (Patterson et al. (1997); Taylor et al. (1998)). If QPOs in DW Cnc are superhumps, its period is probably the longer one, i.e., 73.4 min. This is because the 37.5-min period indicates that the system is an AM CVn object; however, it contradicts the Balmer emissions observed in the optical spectra. According to Patterson (2001), a system with a shorter orbital period generally has a smaller superhump excess, $`\epsilon =|P_{\mathrm{superhump}}P_{\mathrm{orb}}|/P_{\mathrm{orb}}`$. Assuming $`0.02<\epsilon <0.03`$ (Patterson (2001)), we can roughly estimate the orbital period of DW Cnc to be 71.3–72.0 min. This is the shortest orbital period in hydrogen-rich CVs, except for the SU UMa-type DNe, V485 Cen ($`P59\mathrm{min}`$; Augusteijn et al. (1996)) and 1RXS J232953.9+062814 ($`P64\mathrm{min}`$; Uemura et al. (2001)). DW Cnc may thus be a permanent superhumper below the period minimum if its orbital period is 70–80 min.

## 5 Summary

Our photometric monitoring of the cataclysmic variable DW Cnc revealed that the object remained at $`R_\mathrm{c}=14.68\mathrm{mag}`$ for 61 days. This characteristic of the long term light curve indicates that it is not a dwarf nova, which was originally proposed in Stepanian (1982), but a nova-like object. It is less likely that the object includes a white dwarf with a strong magnetic field, since both optical spectra obtained by Stepanian (1982) and Kopylov et al. (1988) show only weak He ii emission. The possible presence of faint states implies the VY Scl-type nature of DW Cnc. We discovered kilo-second quasi-periodic oscillations (QPOs) in DW Cnc, whose center periods were calculated to be $`37.5\pm 0.1`$ and $`73.4\pm 0.4\mathrm{min}`$. Compared with typical dwarf-nova QPOs, the oscillation of DW Cnc is unique regarding the large amplitude, long period, and long lifetime. We propose possible interpretations for DW Cnc QPOs: (i) The trapped disk oscillation scenario, which can be acceptable if the orbital period is above the period gap. (ii) The superhump scenario, which can be acceptable if it is below the period minimum.

We are grateful to many amateur observers for supplying their vital visual CCD estimates via VSNET. Part of this work was supported by a Research Fellowship of the Japan Society for the Promotion of Science for Young Scientists (MU). We would like to thank to Masaryk University for giving an observing time.


text:
# 1 Introduction

## 1 Introduction

“Polyelectrolytes are polymers bearing ionizable groups, which, in polar solvents, can dissociate into charged polymer chains (macroions) and small counterions” . The combination of macromolecular properties and long-range electrostatic interactions produces an impressive variety of phenomena . It makes these systems interesting from a fundamental as well as a technological point of view.

A thorough understanding of polyelectrolytes has become increasingly important in biochemistry and molecular biology. This is due to the fact that virtually all proteins, as well as the DNA, are polyelectrolytes. Their interactions with each other and with the charged cell-membrane are still far from being understood. For instance, a puzzling question is why two equally charged objects should attract each other in the first place .

In this article we focus on stiff linear polyelectrolytes, which we subsequently approximate by charged cylinders. This is a relevant special case, applying to quite a few important polyelectrolytes, like DNA, actin filaments or microtubules. We treat the solvent in dielectric approximation and explicitely describe only the small ions. Within Poisson-Boltzmann (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 .

The paper is organized as follows: First we review some essential features of the PB mean-field solution of the cell model. Then we discuss the applicability of PB theory to the ion distribution functions and show under which circumstances PB theory fails quantitatively (underestimated condensation) and qualitatively (overcharging, charge oscillations and attractive interactions). Next we present measurements of the osmotic pressure for the nematic phase of hexagonally packed polyelectrolytes as a function of density and compare it to PB predictions and the Manning limiting law. In the next section we study the particular case of poly(para-phenylene) by means of PB theory, including a correlation correction of the basis of a recently proposed Debye-Hückel-Hole-cavity theory (DHHC) , and simulational results. The results are compared to recent experimental data on this system . We find that correlation effects enhance condensation and lower the osmotic pressure, yet are not fully able to explain the discrepancy to the experimental data. At the end we try to shed light onto the role of specific ionic correlations. Two-dimensional correlations on the surface of the rods are found to be present, but weakly developed. No hexatic order of the ions is observed.

## 2 Simulation method and model system

### 2.1 The Langevin thermostat

We utilize molecular dynamics (MD) simulations using a Langevin thermostat  to study the equilibrium properties of our model system within the canonical ensemble. Technically this is achieved by integrating the stochastic differential equation

$$m_i\ddot{𝒓}_i=U(\{𝒓_i\})\mathrm{\Gamma }\dot{𝒓}_i+𝝃_i(t)$$
(1)
on the computer, where $`m_i`$ are the masses of particles with coordinates $`𝒓_i`$ subject to a potential energy function $`U`$, a friction $`\mathrm{\Gamma }`$ linear in the velocity and a stochastic white noise $`𝝃_i(t)`$. The latter two can be thought of as imitating the presence of a surrounding viscous medium responsible for a drag force and random collisions, respectively. Since both effects have the same origin, they are related to each other by a simple version of the fluctuation-dissipation-theorem. This can be exploited to choose their strength such as to converge towards the canonical state:

$$𝝃_i(t)=0\text{and}𝝃_i(t)𝝃_j(t^{})=6k_\mathrm{B}T\mathrm{\Gamma }\delta _{ij}\delta (tt^{}),$$
(2)
where $`k_\mathrm{B}T`$ is the thermal energy.

We use a standard Verlet integrator to integrate (1). Time step $`\delta t`$ and friction coefficient $`\mathrm{\Gamma }`$ were set to 0.01 and 1.0 in Lennard-Jones units (see below).

### 2.2 Interaction potentials

We use the purely repulsive Lennard-Jones (LJ) potential to give the counterions an excluded volume:

$$V_{\mathrm{LJ}}(r)=\{\begin{array}{cc}\hfill 4ϵ\left[\left(\frac{\sigma }{r}\right)^{12}\left(\frac{\sigma }{r}\right)^6+\frac{1}{4}\right]:& 0<rr_{\mathrm{cut}}2^{1/6}\sigma \hfill \\ \hfill 0:& r_{\mathrm{cut}}<r.\hfill \end{array}$$
(3)
The advantage of including the $`r^6`$ contribution instead of merely using the purely repulsive $`r^{12}`$ is that Eqn. (3) is exactly zero beyond $`r_{\mathrm{cut}}`$ and merges smoothly to this value at $`r_{\mathrm{cut}}`$, allowing a larger time step.

The Coulomb potential of a charge $`Q`$ is written as

$$\beta e_0V_\mathrm{C}(r)=(Q/e_0)\frac{\mathrm{}_\mathrm{B}}{r},$$
(4)
where $`\mathrm{}_\mathrm{B}=\beta e_0^2/4\pi \epsilon `$ is the Bjerrum length, the distance at which two unit charges have an interaction energy equal to $`1/\beta :=k_\mathrm{B}T`$, $`e_0`$ is the positive elementary charge, and $`\epsilon `$ is the product of the dielectric constants of vacuum, $`\epsilon _0`$, and the medium, $`\epsilon _\mathrm{r}`$, respectively. The Lennard-Jones unit $`\sigma `$ is always used as unit of length and $`ϵ`$ is always set to the thermal energy $`k_\mathrm{B}T`$. Temperature is implemented via the Bjerrum length. Mass is irrelevant and set equal to one – it would only be needed to translate the Lennard-Jones time $`\tau _{\mathrm{LJ}}=\sigma \sqrt{m/ϵ}`$ into “real” time. We intend to model an aqueous environment at room temperature, which implies $`\mathrm{}_\mathrm{B}7.14\text{Å}`$.

Within the periodic boundary conditions employed during the simulations, the presence of such long-range interactions poses both mathematical and technical difficulties. We use the most efficient FFT accelerated Ewald sum, the P3M code, which scales almost linearly with the number of charges .

### 2.3 Generating a cell-geometry

Compared to the spherical cell model, the cylindrical one presents one additional but crucial complication: the charged rod is infinitely long. Several methods have been proposed in the literature to handle this problem . They essentially all use as a unit cell a hexagonal prism with a certain height. This approximates the cylindrical cell by a space-filling object. In the present work we take a cube of side length $`L_\mathrm{b}`$ and place a rod along the main diagonal. Upon periodically replicating this system the diagonal rod becomes infinitely long and an infinite triangular array of such replicated rods emerges. The resulting Wigner-Seitz cell of this lattice is a regular hexagon, which can alternatively be viewed as the unit cell. This is illustrated in Fig. 2. Observe that the symmetry of the replicated system is still cubic. The radius $`R`$ of a circle with the same area $`A`$ is then given by $`R=L_\mathrm{b}/\sqrt{\pi \sqrt{3}}`$. This value is most appropriately used for comparing results between the hexagonal and the cylindrical cell model.

If the line charge density of the rod is $`\lambda `$, electroneutrality requires the number of $`v`$-valent counterions to be $`N=\sqrt{3}L_\mathrm{b}\lambda /ve_0`$. Hence the average counterion density is given by $`n=N/L^3=\sqrt{3}\lambda /ve_0L^2`$ and is thus inversely proportional to $`L^2`$ instead of $`L^3`$. The number of counterions can therefore be written as $`N=(\sqrt{3}\lambda /ve_0)^{3/2}n^{1/2}`$, which implies that a smaller density requires more particles. While this makes the simulation of dilute systems rather expensive, it gives at the same time rather small dense systems. The latter problem can be circumvented by combining blocks of $`2\times 2\times 2`$, $`3\times 3\times 3`$ or even more elementary cubes to a big cube and using the latter as the unit box for the periodic boundary conditions.

The ratio between $`\mathrm{}_\mathrm{D}=(4\pi \mathrm{}_\mathrm{B}v^2n)^{1/2}`$ (the average Debye length of the counterions) and the rod separation $`d_{\mathrm{rod}}`$ can be written as

$$\frac{\mathrm{}_\mathrm{D}}{d_{\mathrm{rod}}}=(8\pi \xi v/\sqrt{3})^{1/2}\mathrm{\hspace{0.33em}0.2625}/\sqrt{\xi v},$$
(5)
where the definition of the Manning parameter $`\xi :=\lambda \mathrm{}_\mathrm{B}/e_0`$ was used, see also Sec. 3. Obviously, for strongly charged rods the Debye length is smaller than the separation of the two rods. Even for $`\xi v=1`$ it is only half as large as the distance between rod axis and Wigner-Seitz boundary, and the neighboring cells effectively decouple. This statement is independent of density.

## 3 Poisson-Boltzmann Essentials

### 3.1 The Poisson-Boltzmann Solution

At this stage we want to briefly recall the necessary knowledge about the PB solution of the cylindrical cell model . Consider an infinitely long cylinder of radius $`r_0`$ and line charge density $`\lambda >0`$, coaxially enclosed in a cylindrical cell of radius $`R`$. Global charge neutrality of the system is ensured by adding an appropriate amount of oppositely charged $`v`$-valent counterions. In the following only the case of no extra salt will be discussed.

Within PB theory the individual counterions are replaced by a cylindrical counterion density $`n(r)`$, where $`r`$ denotes the radial distance from the cylinder axis. Defining the reduced electrostatic potential $`y`$ and and the screening constant $`\kappa >0`$ as

$$y(r)=\beta e_0v\psi (r)\text{and}\kappa ^2=\mathrm{\hspace{0.33em}4}\pi \mathrm{}_\mathrm{B}v^2n(R),$$
(6)
the PB equation can be written as

$$y^{\prime \prime }+\frac{y^{}}{r}=\kappa ^2\mathrm{exp}(y).$$
(7)
It will be useful to define the dimensionless *Manning parameter* $`\xi =\lambda \mathrm{}_\mathrm{B}/e_0`$ which counts the number of unit charges on the rod per Bjerrum length. In the following the main focus will be on the strongly charged case characterized by $`\xi v>1`$. Fixing the two boundary conditions for the electric field,

$$y^{}(r_0)=2\xi v/r_0\text{and}y^{}(R)=0,$$
(8)
the correctly normalized solution to Eqns. (7,8) can be written as

$$y(r)=2\mathrm{ln}\left\{\frac{r}{R}\sqrt{1+\gamma ^2}\mathrm{cos}\left(\gamma \mathrm{ln}\frac{r}{R_\mathrm{M}}\right)\right\}.$$
(9)
Insertion of the general solution from Eqn. (9) into the boundary conditions from Eqn. (8) yields two coupled transcendental equations for the two integration constants $`\gamma `$ and $`R_\mathrm{M}`$:

$$\gamma \mathrm{ln}\frac{r_0}{R_\mathrm{M}}=\mathrm{arctan}\frac{1\xi v}{\gamma }\text{and}\gamma \mathrm{ln}\frac{R}{R_\mathrm{M}}=\mathrm{arctan}\frac{1}{\gamma }.$$
(10)
Subtracting the left part from the right part in Eqn. (10) eliminates $`R_\mathrm{M}`$ and provides a single equation

$$\gamma \mathrm{ln}\frac{R}{r_0}=\mathrm{arctan}\frac{1}{\gamma }+\mathrm{arctan}\frac{\xi v1}{\gamma },$$
(11)
from which $`\gamma `$ can be obtained numerically. The second integration constant $`R_\mathrm{M}`$, which will be referred to as the Manning radius, can be obtained from either one of the Eqns. in (10) as soon as $`\gamma `$ is known. Note finally that $`\kappa `$ and $`\gamma `$ are connected via

$$\kappa ^2R^2=2(1+\gamma ^2).$$
(12)
Since the $`\xi `$ and $`v`$ only enter in the combination $`\xi v`$, changing valence or electrostatic interaction strength is equivalent on Poisson-Boltzmann level. In particular, at given cell geometry $`\{r_0,R\}`$ the Manning radius $`R_\mathrm{M}`$ is a unique function of $`\xi v`$. Eqn. (11) implies the sequence of inequalities $`\mathrm{ln}(R/r_0)\pi /\gamma \mathrm{ln}(R/r_0)+\xi v/(\xi v1)`$. The resulting asymptotic for $`\gamma `$ in the dilute limit $`R\mathrm{}`$ gives rise to what is often called the “Manning limiting laws”.

The integrated probability distribution of finding a mobile ion at distance $`r`$ within a cylinder of radius $`r[r_0;R]`$ can be determined analytically by integrating the charge density:

$$P(r)=\frac{1}{\lambda }_{r_0}^r\text{d}r^{}\mathrm{\hspace{0.17em}2}\pi r^{}e_0vn(r^{})=\left(1\frac{1}{\xi v}\right)+\frac{\gamma }{\xi v}\mathrm{tan}\left(\gamma \mathrm{ln}\frac{r}{R_\mathrm{M}}\right).$$
(13)
This is the fraction of counterions found within a cylinder of radius $`r`$. $`P(r)`$ will be referred to as counterion distribution function. At $`r=R_\mathrm{M}`$ the last term in $`P`$ vanishes, giving the Manning fraction $`f_\mathrm{M}:=11/\xi v`$ of ions within $`R_\mathrm{M}`$. These are the ions which can not be diluted away, if the cylinder radius is sent to infinity. This phenomenon is referred to as Manning condensation .

### 3.2 How to quantify counterion condensation

If the counterion distribution function $`P`$ is known, the condensed counterion fraction can be characterized in the following “geometric” way: Eq. (13) shows that $`P`$ viewed as a function of $`\mathrm{ln}(r)`$ is merely a shifted tangent-function with its center of symmetry at $`(\mathrm{ln}(R_\mathrm{M});f_\mathrm{M})`$. Since, however, $`\mathrm{tan}^{\prime \prime }(0)=0`$, Manning radius and Manning fraction can be found by plotting $`P`$ as a function of $`\mathrm{ln}(r)`$ and localizing the point of inflection.

This property of $`P`$, derived within the framework of PB theory, can in turn be used to define the condensed fraction . It provides a suitable way to quantify counterion condensation beyond the scope of PB theory, and it is exact in the salt-free PB limit. From here on this method will be referred to as the inflection point criterion. It has the advantages of (i) not fixing by definition the amount of condensed counterions ($`f_\mathrm{M}`$ and $`R_\mathrm{M}`$ can be determined independently of each other), (ii) reproducing the salt-free PB limit, and (iii) quantifying the breakdown of the coexistence of condensed and uncondensed counterions in the high salt limit.

## 4 Generic ion distribution functions

On the plain PB level the radius $`r_0`$ of the rod is not a completely independent variable, since it enters only in the combination $`R/r_0`$, see Eqn. (11). The rod in our simulation is built up from a sequence of small Lennard-Jones particles lined up along the main diagonal of the simulation cell. The distance of closest approach of two Lennard-Jones particles can then be identified with the radius $`r_0`$ of the rod and yields an effective rod radius of $`\sigma `$.

Upon leaving the PB level, the ratio between counterion diameter $`\sigma `$ and rod radius $`r_0`$ becomes relevant, which has for instance been investigated in Ref. . Here we do not intend to perform a systematic investigation of such effects, but instead present later in Sec. 6 results for systems mapped to physically relevant parameters, DNA and two kinds of poly(p-phenylenes), in which the rod has a considerably larger diameter than the counterions.

### 4.1 Density dependence within monovalent systems

At fixed rod radius of $`r_0`$ the relevant variable $`R/r_0`$ is changed by varying the cell size $`R`$ and therefore the density. Here we present results for such a density scan for systems with $`\mathrm{}_\mathrm{B}/\sigma \{2,3\}`$. The monovalent and positively charged particles forming the rod were placed along the main diagonal with a center-center distance of $`1.04245\sigma `$, giving a line charge density of $`\lambda =0.959279e_0/\sigma `$. In connection with the two presented values for the Bjerrum length this yields Manning parameters of $`\xi =1.919`$, $`2.878`$ and corresponding Manning fractions of $`f_\mathrm{M}=11/\xi =0.4788`$, $`0.6525`$ in the monovalent case. The cell radius $`R`$ has been varied between $`2.06\sigma `$ and $`124\sigma `$, which for monovalent counterions corresponds to average ion densities of $`7.2\times 10^2\sigma ^3`$ and $`2.0\times 10^5\sigma ^3`$, respectively. The inflection point criterion from Sec. 3.2 has been used to determine the radial extension of the condensed layer, $`R_\mathrm{M}^{}`$, and the fraction of ions within it, $`f_\mathrm{M}^{}=P(R_\mathrm{M}^{})`$, where the star denotes the measured quantities. A graphical illustration of the distribution functions is given in Fig. 3.

The important things to observe are the following. The measured condensed fraction $`f_\mathrm{M}^{}`$ is always larger than the fraction predicted by the PB theory. However, the fraction from the simulation decreases monotonically with decreasing density towards the Manning limit $`f_\mathrm{M}`$. Such a deviation is to be expected, since PB is essentially a low density theory .

In contrast to the clear tendency of the measured condensed fraction to decrease upon dilution, the behavior of the condensation radius $`R_\mathrm{M}^{}`$ appears to be more complicated. There does not seem to exist a simple monotonic convergence of $`R_\mathrm{M}^{}`$ towards $`R_\mathrm{M}`$. Rather, for high densities the measured condensation distance is larger than the Manning radius, while for the investigated low densities it is smaller. Unfortunately, a clear-cut statement is difficult since the localization of the point of inflection in $`P`$ as a function of $`\mathrm{ln}(r)`$ is only possible with an error estimated to be of the order of $`1\%`$.

For Bjerrum length $`\mathrm{}_\mathrm{B}=1\sigma `$, corresponding to $`\xi =0.9593<1`$, the radial distribution functions are compared to the PB predictions in Fig. 4. It can be seen that the measured and PB predicted distributions almost coincide for all cell sizes under investigation. While this is to be expected for low densities, the remarkable agreement at high densities is somewhat surprising. Apparently, PB theory is a fairly good description of weakly charged systems, i.e., ones which are below the Manning threshold.

### 4.2 Fitting to a generalized Poisson-Boltzmann distribution

Since the measured Manning radius $`R_\mathrm{M}^{}`$ can be smaller than the PB prediction, but $`f_\mathrm{M}^{}>f_\mathrm{M}`$ and $`R_\mathrm{M}^{}<R_\mathrm{M}`$ is incompatible with PB theory, it is impossible to describe such a measured $`P(r)`$ with a PB distribution function in which a somewhat enlarged effective Manning parameter $`\xi ^{}`$ is used. Nevertheless, it would be desirable to use at least some of the knowledge about the ion distribution function from PB theory for evaluating relevant observables in a simulation or in a real experiment. To the lowest order, e.g., within a two state model, the increase of counterion condensation can be ascribed to such an effective $`\xi ^{}v=1/(1f_\mathrm{M}^{})`$. To break the monotonic connection between $`f_\mathrm{M}`$ and $`R_\mathrm{M}`$ one assumes that the functional form of $`P`$ is given by the PB form in Eqn. (13), but neglects the relation in Eqn. (11). Using Eqn. (10i), this suggests fitting the measured distribution in a region around the inflection point to the form

$$P_{\mathrm{fit}}(r)=\mathrm{\hspace{0.33em}1}\frac{1}{\xi ^{}v}+\frac{\gamma ^{}}{\xi ^{}v}\mathrm{tan}\left(\gamma ^{}\mathrm{ln}\frac{r}{r_0^{}}+\mathrm{arctan}\frac{1\xi ^{}v}{\gamma ^{}}\right)$$
(14)
with the three fit parameters $`\xi ^{}`$, $`\gamma ^{}`$ and $`r_0^{}`$.333It might seem strange to regard the rod radius as a free parameter, but there are two reasons for this. First, this does not force the fit to coincide with the PB form outside the chosen fitting-region. Second, not even in a real experiment is the radius of a macroion always known in advance. Rather, it is often necessary to obtain this information within the same measurement . The condensed fraction is then given by $`f_\mathrm{M}^{}=11/\xi ^{}v`$ and the condensed radius by $`R_\mathrm{M}^{}=r_0^{}\mathrm{exp}\{\mathrm{arctan}[(\xi ^{}v1)/\gamma ^{}]/\gamma ^{}\}`$. Fig. 5 illustrates this procedure for one system. While this approach might seem to be very powerful, it suffers from the drawback that the actual numbers depend on the chosen fitting region. Nevertheless, it provides an independent way of quantifying condensation and can even be applied to very noisy data, for which the localization of an inflection point in $`P`$ as a function of $`\mathrm{ln}(r)`$ is otherwise virtually impossible.

### 4.3 Multivalent ions

The universal dependence of the ion distribution on the product $`\xi v`$ is an artifact of Poisson Boltzmann theory. Fig. 6 shows examples of systems with different Manning parameters but the same value of $`\xi v`$. Not only is the condensation enhanced as compared to PB theory, but the enhancement is stronger for the case involving multivalent ions. For instance, in the dense trivalent system the condensation is enhanced by 30% with respect to the PB prediction.

The reason is that at given charge density multivalent systems have a lower number density and thus fewer particles. This lowers any kind of excluded volume interactions, and also makes the entropic contribution less pronounced. Moreover, the relevant variable for describing the strength of correlation effects is $`\mathrm{}_\mathrm{B}v^2`$ . Hence, an increase in valence is more important than an increase in Bjerrum length.

Within PB theory the expression for the contact density $`n(r_0)`$ in the limit of infinite dilution is given by

$$n(r_0)\stackrel{R\mathrm{}}{=}\mathrm{\hspace{0.33em}2}\pi \mathrm{}_\mathrm{B}\stackrel{~}{\varsigma }^2\left(1\frac{1}{\xi v}\right)^2,$$
(15)
where $`e_0\stackrel{~}{\varsigma }`$ is the surface charge density of the rod. For high rod charge this expression becomes independent of valence. Hence, replacing monovalent counterions by multivalent ones reduces their total *number* in the cell, but for a highly charged rod not their *density* at the rod surface. This result is a consequence of the contact value theorem and thus holds beyond the Poisson-Boltzmann approximation.

### 4.4 Addition of salt

The influence of salt on the distribution functions has been discussed within the PB theory in detail in Ref. . The general finding was that a low salt content leaves the picture of Manning condensation qualitatively unchanged, while at increasing salt concentration a crossover between Manning condensation and simple salt-screening occurs.

If one starts adding $`N`$ salt molecules to the salt free solution one finds that the inflection point gets shifted to smaller values of $`r`$, hence the layer of condensed counterions contracts. The second finding is that the amount of condensation is only marginally increased. From a certain $`N`$ on two more inflection points appear close to the cell boundary. This happens typically for a corresponding Debye length being of the order of the cell size itself, indicating the appearance of a characteristic, salt induced, change in the convexity of $`P`$ (as a function of $`\mathrm{ln}r`$). Upon a further increase in $`N`$ one of the two new inflection points shifts towards smaller values of $`r`$, finally fusing with the Manning inflection point and “annihilating” it. Roughly speaking, the inflection points vanish if the Debye length characterizing the salt content becomes smaller than the radius of the condensed layer. In this case it is no longer meaningful to distinguish between condensed and uncondensed counterions. Indeed, for a very high salt content, where the Debye length is much smaller than the radius of the rod, the solution of the PB equation would be the one of a charged plane and one may consider all excess counterions being condensed – no matter what the charge density of the rod is.

For three systems with numbers of salt molecules $`N\{0,104,3070\}`$ a simulation has been performed and compared to the PB prediction in Fig. 7. As in the salt-free case the computer simulations show a more pronounced condensation effect towards the rod. Nevertheless, the shape of the distribution functions remains qualitatively the same. Note in particular that the appearance and disappearance of two points of inflection at $`N=104`$ and $`N=3070`$ respectively, which leads to extremely small curvatures in the PB distribution functions, also leads to very straight regions in the measured distribution functions. The crossover from Manning condensation to screening, as described within the PB theory, can be expected to be essentially correct. It should not be overlooked, however, that the addition of salt ions will affect the distribution function much more dramatically if their valence is larger than the valence of the counterions. In this case, the salt counterions will accumulate in the vicinity of the rod at the expense of the lower valent “real” counterions as has previously been demonstrated in a system with a valence mixture of counterions .

The PB approach fails to describe the physical situation if one or more of the following conditions apply: (i) the electrostatic interactions are strong, (ii) the counterions are multivalent or (iii) the density is high. Results of a simulation under such conditions can be seen in Figure 8. Here, a system with box length $`L_\mathrm{b}/\sigma =36.112`$, i.e., 60 monovalent counterions and a cell size of $`R/\sigma =15.481`$, and $`\mathrm{}_\mathrm{B}/\sigma =4.1698`$, i.e., $`\xi =4`$, has been investigated after adding 1000 molecules of a 2:2 salt. Since this gives almost 17 times as much salt as counterions and a salt Debye length of $`\mathrm{}_\mathrm{D}/\sigma =0.33R/\sigma `$, this can essentially be viewed as a charged rod in a bulk 2:2 electrolyte. The most characteristic feature of the charge distribution function is that it overshoots unity, showing a charge reversal of the rod at distances around $`r1.5\sigma `$, while the simple PB prediction is clearly qualitatively off. This phenomenon is usually referred to as overcharging and has been predicted for the rod geometry first from hypernetted chain calculations  and later by a modified PB approach , and has also been observed in planar geometries. Since $`P(R)=1`$ for the reason of global electroneutrality, the overshooting above 1 at small distances implies the existence of a range of $`r`$-values at which the mobile ion system is locally positively charged, i.e., with the same charge as the rod, such that $`P(r)`$ can eventually decay to 1. This is seen in the right frame of Fig. 8, which shows that $`n_{+2}(r)>n_2(r)`$ at $`r2\sigma `$. Since $`P(1.5\sigma )1.45`$, the rod and its innermost layer of condensed ions could be viewed as an effective rod of radius $`1.5\sigma `$ which is negatively charged with Manning parameter $`\xi =1.8`$. Since this value is again larger than 1, it entails ion condensation, but this time of positive ions. In fact, it even leads to a second overcharging, as can clearly be seen in Fig. 8, where $`P(r)`$ – in decaying from 1.45 – overshoots the value of 1 again. Overcharging can thus give rise to layering. In the presented example no less than three layers can clearly be made out. These local charge oscillations also reflect themselves in oscillations of the electrostatic potential, as demonstrated in the inset in the right frame of Fig. 8. Notice that these oscillating potentials will also have pronounced effects on the interaction between such rigid polyelectrolytes. HNC/MSA theory can qualitatively and quantitatively describe such effects to a very good degree, as has been demonstrated theoretically a long time ago . A recent comparison with simulations can be found in Ref. .

## 5 Pressure

### 5.1 Pressure measurements within the generic model

Using the results from the Appendix A we now turn towards pressure measurements for the generic systems described above. Because of the periodic boundary conditions, we are actually looking at the pressure of a nematic, hexagonally packed solution of charged rods. We specifically compute the dimensionless osmotic coefficient $`\widehat{p}`$, which is simply the pressure normalized by its ideal gas contribution

$$\widehat{p}=\frac{p}{p^{\mathrm{ig}}}=\frac{p}{nk_\mathrm{B}T}.$$
(16)
For an isolated cell the pressure is given by the particle density at the outer cell boundary. This is a rigorous statement, true for the spherical, cylindrical and planar cell model . It is a merit of the PB equation that it retains the validity of this exact relation. For the extended density functional theories this no longer holds and additional terms appear .

Within Poisson-Boltzmann theory the osmotic coefficient as determined from the boundary density is

$$\widehat{p}=\frac{1+\gamma ^2}{2\xi v}\stackrel{R\mathrm{}}{=}\frac{1}{2\xi v},$$
(17)
where $`\gamma `$ is the density dependent integration constant from Eqn. (11),

Fig. 9 illustrates these measurements for monovalent counterions and three values of $`\mathrm{}_\mathrm{B}/\sigma =1,2,3`$. Several things may be noted: The osmotic coefficient from the simulations is always smaller than the PB prediction, but for low density both values converge. This also illustrates that the Manning limiting law from the r.h.s of Eqn. (17) becomes asymptotically correct for dilute systems. Upon increasing the density, the osmotic coefficient rises weaker than the PB prediction. This is more pronounced for systems with higher Bjerrum length, and consequently, higher Manning parameter, and is due to enhanced counterion condensation which has been observed in Sec 4. Notice that this has a very remarkable side-effect. Over a considerable range of densities the measured osmotic coefficient is much closer to the limiting law than to the actual PB prediction. This makes the Manning limit look much more accurate than it really ought to be. However, the surprising effect should not be over-interpreted, since the underlying reason is nothing but a fortunate cancellation of two contributions of approximately the same size which are not contained in the limiting laws.

Also for the pressure it is interesting to investigate complementary systems in which the values of Bjerrum length $`\mathrm{}_\mathrm{B}/\sigma `$ and valence $`v`$ have been interchanged to keep the product $`\xi v`$ fixed. However, at given density the cell radius depends on the valence, so $`\widehat{p}(n)`$ does no longer universally depend on this product. Since in the limit of high density $`\widehat{p}`$ is easily seen to increase with density, while at low density it decreases, the functions $`\widehat{p}(n)`$ for different valences intersect at some point.

Fig. 10 summarizes the results of measurements on the multivalent systems with $`v=1,2,3`$, which yield the same values of $`\xi v`$ as the monovalent ones investigated before. The most striking feature is the appearance of a negative osmotic coefficient in a certain density region of the trivalent case. If the constraint of a fixed rod-separation were to be replaced, the system would phase separate, hence, attractive interactions must be present between the rods. Similar observations have been reported in Refs. . In order to trace back the origin of those attractions, Fig. 11 displays the osmotic coefficient in the divalent and trivalent cases, split into two contributions: (i) the non-electrostatic part comprising ideal gas contribution and the short-range virial and (ii) the electrostatic part, which for visual convenience is plotted with a reversed sign. Hence, the difference between those two curves gives the curves in Fig. 10. We observe that the electrostatic part leads to a monotonically increasing attraction with increasing density. It is almost twice as strong in the trivalent case. The very strong increase of the osmotic coefficient at large densities is due to the virial, i.e., due to repulsive ion-rod and ion-ion interactions. The measured negative pressure in the trivalent case is the result of a “sudden” drop in the virial contribution. Nothing particular can be observed in the electrostatic part. We suggest the explanation that in an intermediate range of densities the presence of surrounding rods lifts condensed ions from the surface they are pushing on, thereby reducing the short range excluded volume contribution. In too dilute systems this effect is not operative, while in too dense systems the effect of the other rods is actually to push ions more closely to the surface. Hence, there is a density range in which the virial is reduced by the presence of other rods.

Contrary to the simulations, the osmotic coefficient from the PB theory is always positive. This is the consequence of the rigorous proof that such attractive interactions are absent on the PB level . An extension of this statement to ions of finite size and to a wider class of boundary conditions and density functionals can be found in Ref. .

Finally it should be noted that the above measurements can not be used to infer that attractive forces between charged rods require the counterions to be at least trivalent. The reason is twofold: First, at given valence one can vary Bjerrum length and line charge density. Increasing the Manning parameter will lead to negative pressure in the divalent system. Second, keeping all interaction potentials fixed, the radius $`r_0`$ of the charged rod is a relevant observable, as has been demonstrated in Ref. . Hence, a general statement about presence or absence of attractive interactions is difficult, since a five-dimensional parameter space is involved: $`\{\lambda ,\mathrm{}_\mathrm{B},v,r_0,n\}`$.

## 6 Examples: poly(p-phenylene)

This section discusses simulations in which the parameters are explicitly mapped to an experimental system. This affects the values of ion diameter $`\sigma `$, rod radius $`r_0`$, Bjerrum length $`\mathrm{}_\mathrm{B}`$ and line charge density $`\lambda `$. In order to have $`r_0`$ different from $`\sigma `$, we use as the ion-rod-potential a LJ like expression, in which $`r`$ is replaced by $`rr_\mathrm{s}`$ and the shift $`r_\mathrm{s}`$ is related to the desired rod radius by $`r_0=r_\mathrm{s}+\sigma /2`$.

Since these systems have been investigated either in the presence of salt or at rather high Manning parameter, the phenomena happening in the vicinity of the rod are always well decoupled from the cell boundary. Therefore, an even simpler geometry has been chosen: One rod placed parallel to one of the edges of a cubic simulation cell.

### 6.1 Distribution functions and osmotic coefficient for salt free poly(p-phenylene) solutions

The Manning limiting laws are claimed to apply in the limit of low density. However, an experimental verification of this effect using DNA as rod-like polyelectrolytes is difficult, since the double helix starts to unwinds at low ionic strength. It would therefore be desirable to have a stiff model polyelectrolyte which does not suffer from this problem. One such system belonging to the class of poly(para-phenylenes) (ppp) has recently been investigated in Refs. . The constitution formula is given in Fig. 1. The fully aromatic backbone exhibits an excellent chemical stability and has a persistence length of approximately $`20\text{nm}`$. The degree of polymerization used in the abovementioned studies was between 20 and 40, so that the contour length equals approximately one persistence length at most. The dominant contribution to the signal in small angle X-ray scattering experiments stems from the iodine ions, since the excess electron density of the backbone is very low. Table 1 lists four systems which have been simulated based on a mapping to ppp.

In a first step the simulated ion distribution functions shall be compared with the PB prediction as well as with the Debye-Hückel-hole-cavity theory from Ref. , which tries to incorporate correlations missing in the PB approach. Figure 12 shows the corresponding distribution functions for the systems from Table 1. Qualitatively the PB prediction is already a good description of those systems, if one graciously disregards the inevitable problems at the cell boundary. However, the presence of correlations shifts the distribution functions up in all four cases. Since this shift is rather small, i.e., the correlations are weak, the Debye-Hückel-hole-cavity theory is in fact an excellent description of those systems. Observe, for instance, that in the weakly charged dilute case its prediction can no longer be distinguished from the simulated curve on the chosen scale.

With the help of the condensation criterion from Sec. 3.2 the extent of the correlation-enhanced condensation can be further quantified. The Manning fraction from the molecular dynamics simulation increases only by a fairly small amount. It is at most 4% larger than the PB value. This translates to an effective Manning parameter $`\xi _{\mathrm{eff}}=1/(1f_\mathrm{M})`$ being $`5`$$`10`$% larger than the bare one. This increase is very accurately captured by the Debye-Hückel-hole-cavity theory. Its prediction for $`f_\mathrm{M}`$ is at most 1% smaller than the value obtained in the simulation. Interestingly, it is even independent of density. This, however, is not a feature to be generally expected and should therefore not be over-interpreted.

In a second step the osmotic coefficients of the four model systems have been computed. Two distinct approaches have been used for analyzing the simulation results: The first computes the pressure from the stress tensor, as described in detail in Appendix A. The second approach exploits the connection from Eqn. (17) between osmotic coefficient, Manning parameter and the integration constant $`\gamma `$ of the PB equation. The idea is to use the fitting procedure described in Sec. 4.2 and from that obtain values $`\xi ^{}`$ and $`\gamma ^{}`$, leading to the osmotic coefficient $`(1+\gamma _{}^{}{}_{}{}^{2})/2\xi ^{}`$. The most straightforward approach of measuring the boundary density is less recommendable, since in the simulation the boundary has a quadratic and not a circular cross section, i.e., the merely approximate representation of the cell model becomes most visible and questionable here. In contrast to that, the two other approaches “look” at regions of the cell, which are away from the boundary. For comparison, also the predictions from the Manning limiting law and the full PB expression are calculated. For the two dilute systems there also exist measurements in Ref.  of $`\widehat{p}`$ via osmometry, while for the dense systems this was unfortunately impossible due to counter diffusion problems. Table 2 collects the coefficients for comparison.

The first thing which should be noted is that PB theory gives again a surprisingly good description of the measured coefficients, including the experimentally determined one. A closer look reveals that – as expected – it overestimates $`\widehat{p}`$, and the deviations are larger for the dense systems. Still, PB theory is off by only 7% or 18% for the dilute or dense systems, respectively, thereby confirming the observation already made when looking at the distribution functions. A further indication of this point is that the Manning limit $`1/2\xi `$ appears to be a lower boundary, which is also in accord with the pressure measurements on the generic systems with monovalent ions from Sec. 5.1. Since there it has also been found that the Manning limit need no longer act as a boundary in the multivalent case, see Fig. 10, it would be very interesting to perform experiments on those systems with, e.g., divalent ions. Observe that the experimental value of $`\widehat{p}`$ for the highly charged dilute system 2 is already at the Manning limit.

Since the experimentally determined osmotic coefficient appears to be smaller even than the molecular dynamics results, this indicates effects to be relevant which go beyond the model used for simulation. Most obvious candidates for this are the neglect of additional chemical interactions between the ions and the polyelectrolyte as well as solvation effects, i.e., interactions between the ions or the polyelectrolyte with the water molecules from the solution. It is for instance demonstrated in Ref.  that the osmotic coefficient also depends on whether on uses chlorine or iodine counterions. While one could certainly account for the different radii of these ions when computing the distance of closest approach entering the PB equation, the implications of the different hydration energies is much less obvious to incorporate and in principle requires very expensive all-atom simulations. See also the discussion in Ref. , which includes the effects of small excess salt to the osmotic coefficient.

Finally it should be noted that the two different methods of analyzing the molecular dynamics data lead to very similar results. While this fact is not particularly useful in a simulation, its consequences for the analysis of experimental data seem promising. In Ref.  poly(p-phenylene) solutions are investigated by means of small angle X-ray scattering, which due to the rod-like geometry turns out to be sensitive to the radial ion distribution function. Therefore, the measured structure factors can be related to PB distribution functions with effective values of $`\xi `$, $`\gamma `$ etc. Fitting the scattering intensity and thereby determining those values permits in principle a measurement of the osmotic coefficient along the lines already described above. What makes this approach so attractive is that it could ideally complement osmometry. If the density becomes too large, the latter suffers from severe problems with counter diffusion, but it is exactly this high density which meets the requirements of good contrast necessary for scattering.

A similar fitting procedure has been applied to the measured data from small angle X-ray scattering experiments on the systems 3 and 4; see also Ref. . From a fit to the structure factor the radius $`r_0^{}`$ of the rod has been determined. The integration constant $`\gamma ^{}`$ and the Manning parameter $`R_\mathrm{M}^{}`$ the follow from the PB theory. The obtained values are listed in Table 3. The resulting Manning radius is surprisingly close to the one determined by simulation, particularly for the system with $`\xi =3.32`$. The osmotic coefficients constructed from the PB formula are less accurate, 10-20% too large.

## 7 The importance of correlations

We have seen that the nonlinear PB equation suffers from systematic deviations in strongly coupled or dense systems. It underestimates the extent of counterion condensation and at the same time overestimates the osmotic coefficient. The common reason for both problems is the neglect of correlations: More elaborate theories  which try to incorporate correlations reproduce the trends of the simulations, because they favor an increase in density close to the macroion surface, thereby leading to a stronger condensation and a concomitant drop in the boundary density and thus pressure.

On the one hand there is general agreement about the presence of correlations being liable for the failure of PB theory. On the other hand, this insight does not shed any light onto the question, which kind of correlations are important. In the following we investigate a simple pair-correlation function based on a one-rod property.

If the surface charge density of the rod is high, a fairly large number of counterions will stay within a condensed layer of small radial extent. Addition of salt will further increase the ionic density in this layer. It is clear that beyond a certain point (high $`\lambda `$, high $`\mathrm{}_\mathrm{B}`$ or enough salt) the ions will no longer distribute independently of each other but get locally correlated. This effect will now be measured for a DNA-like system, because for those systems correlation effects are assumed to play a prominent role. The rod radius was chosen to be 7.86 Å($`1.85\sigma `$) and the ion diameter is correspondingly 4.25 Å($`1\sigma `$). In addition 0.5 mol/l of a 2:2 salt has been added in excess to the divalent counterions. More information about theses systems can be found in Table 4.

The first issue is to define the innermost layer. For this, a distance from the rod is chosen which contains many counterions but virtually no co-ions. A distance of roughly $`11.5\text{Å}`$ from the rod axis turned out to be suitable. This is about a third ion diameter farther out than the distance of closest approach. To avoid difficulties with remaining co-ions, only the counterions within this distance are taken into account in all what follows. In a second step the coordinates of those ions are radially projected onto the surface of the cylinder of closest approach, and this surface is then rolled out to a flat plane, see Fig. 13 for an illustration of this procedure. Finally, the two-dimensional pair correlation function $`g(r)`$ of these projected points is computed.

One might object that the unrolled flat plane leads one to believe that the pair correlation function does only depend on distance, while in reality such a rotational symmetry cannot be expected, due to the rod curvature. However, an investigation of $`g(𝒓)`$ revealed that the pair correlation function is in fact rotationally symmetric up to essentially a distance which corresponds to half the circumference of the cylinder of closest approach, which is roughly $`30\text{Å}`$. Larger distances can of course only be realized along the rod instead of around it. A possible reason for this surprising symmetry is that at short distances the curvature is not yet perceptible while at large distances the particles are already uncorrelated.

Figure 14 shows the measured $`g(r)`$ for four differently charged systems with Manning parameter 10.5, 8.4, 6.3, and 4.2. The last value corresponds to DNA in aqueous solution. The most important thing to observe is that for the strongly charged systems $`g(r)`$ shows definite signs of correlations. Apart from the trivial correlation hole at small $`r`$ there is a distance $`r_{\mathrm{max}}`$ at which ions are more likely to be found than under the assumption of independent distribution. Notice that the maximum in $`g(r)`$ is not located at $`r\sigma `$ but much farther out. Hence, its existence is not merely an artifact of close packing of repulsive Lennard-Jones spheres. However, if the condensed ions are assumed to form a triangular lattice on the surface in order to maximize their mutual repulsion, the resulting distance $`d_{\mathrm{}}`$ is 25 - 35 % larger than $`r_{\mathrm{max}}`$. Together with the only weakly pronounced oscillations in $`g(r)`$ this proves the correlation induced interactions to be rather short-ranged, yet less local than a pure hard core. In any case, the range is several times larger than the average salt Debye length $`\mathrm{}_\mathrm{D}0.5\sigma `$.

Such two-dimensional pair correlation functions for mobile ions adsorbed on charged planes have recently been investigated theoretically in Ref. . The mean separation $`d_\varsigma =(ve_0/\varsigma )^{1/2}`$ of $`v`$-valent ions on the completely neutralized plane of surface charge density $`\varsigma `$ is identified as the important scaling length. The main predictions for the strongly coupled regime $`\mathrm{\Gamma }_2=\mathrm{}_\mathrm{B}v^2/d_\varsigma >1`$ are: (i), $`g(r)`$ should have a single first peak at a distance about the size $`d_\varsigma `$ of the correlation hole. (ii), the breadth of this peak should decrease with $`\mathrm{\Gamma }_2`$ while its maximum should increase. This trend is indeed seen to be true for the functions in Fig. 14, only $`d_\varsigma `$ is roughly 10–20% larger than the actual peak position. One might want to argue that the presence of much salt entails a further increase of the layer density on top of the usual Gouy-Chapman prediction, but the results obtained in Ref.  are claimed to be independent of the bulk ionic strength, if the latter is smaller than the layer density. This is here the case, even though not always by an order of magnitude. An alternative explanation may be based on the presence of co-ions: Although there reside very few of them within the innermost condensed layer, there will be an appreciable amount beyond and possibly very close to $`r_\mathrm{l}`$. Those ions can act as “bridges”, they attract counterions and thereby reduce the average closest distance between them.

Although $`r_{\mathrm{max}}<d_{\mathrm{}}`$, this does not exclude a local hexatic ordering of the counterions. But since $`g(𝒓)`$ is on average rotationally symmetric, establishing signs for this requires the investigation of a suitable 3-point correlation function. The trick is to break the rotational symmetry, which suggests the following procedure: $`g(r)`$ is proportional to the probability of finding a particle at a distance $`r`$, given that there is also a particle at the origin. Now define $`g_{}(𝒓)`$ to be proportional to the probability of finding a particle at position $`𝒓`$, given that there is also a particle at the origin and given that the particle which is closest to the origin is situated to the right side. This definition is further restricted to be sensitive only for next nearest neighbors and not arbitrary other particles. Thence, it answers the question: “If the nearest neighbor is on the right side, where is the next nearest neighbor?” A scatter plot of this observable is shown in the left part of Fig. 15. For this, the most highly charged DNA-sized system with Manning parameter $`\xi =10.5`$ has been used. Clearly visible are the two correlation holes around the origin and the nearest neighbor as well as a tendency of the next nearest neighbor to be on the side opposing the nearest neighbor. All those effects are trivially explained. However, on top of that no preferential ion accumulations in the “hexatic directions” indicated as dotted lines is perceptible. The right part of Fig. 15 shows the probability distribution of finding the next nearest neighbor at an angle $`\varphi `$ with respect to the line joining nearest neighbor and origin. The correlation hole generated by the nearest neighbor is visible, but no increase in $`p(\varphi )`$ at $`60^{}`$, $`120^{}`$ or $`180^{}`$. Rather, the measured data are compatible with a fairly structureless distribution, as indicated by the solid line.

Since the most strongly charged system does not show hexatic order, it will certainly be absent in the other ones. However, the coupling constant $`\mathrm{\Gamma }_2`$ introduced above has a value only about $`3`$ in the system with Manning parameter $`\xi =10.5`$. For larger values of $`\mathrm{\Gamma }_2`$ one would not only expect much stronger correlations and hexatic ordering but even crystallization of the adsorbed counterions into a two-dimensional Wigner crystal . More extensive investigations to clarify the role of correlations are presently under way

## 8 Summary

Theoretical and numerical studies of stiff linear polyelectrolytes within the framework of a cell model have been presented. We outlined methods to quantify condensation for general measurements of ion distribution functions via an inflection point criterion and by using a generalized PB fit function. We compared our simulational results to predictions from mean-field PB theory, and found excellent agreement for weakly charged systems. The enhanced counterion condensation caused by ionic correlations and its dependence on parameters such as density, Bjerrum length, valence and ionic strength was measured, and effects which qualitatively go beyond mean-field theory, e.g., charge reversal and attractive interactions between like-charged macroions have been found. For systems with a high salt content, integral equation theories have proved to give a good qualitative and quantitative description of such systems

For the osmotic pressure measurements we displayed for a large density range that the measured pressure stays close to the Manning limiting law prediction due to a fortunate error cancellation of neglected ionic correlations and density effects which point in different directions. We showed for the specific example of a ppp-polyelectrolyte what differences to a PB description can be expected for both, the ionic distribution and the osmotic pressure. We finally analyzed some correlations for a DNA-like system, and showed that they show signs of a strongly correlated liquid, as has been advocated recently .

To assess these theoretical ideas, computer simulations will certainly be indispensable, since they are presently the only practicable way of obtaining sufficiently detailed information on ionic distributions and correlation functions. The numerical investigations presented in this work and the correlation analysis of the obtained data are a further step in this direction, however a more detailed analysis is presently under way .

## 9 Acknowledegments

We wish to thank M. Ballauff, M. Barbosa J. Blaul, B. Guilleaume, F. Jimenez-Angeles, K. Kremer, M. Lozada-Cassou and S. May for various contributions to this work. In addition we acknowledge a large computer time grant hkf06 from NIC Jülich and financial support by the German Science foundation.

Appendix:

## Appendix A Defining and computing the pressure

The pressure for the primitive cell model is a nontrivial variable to compute for two reasons: First, the long-range electrostatics has to be properly taken into account, and secondly, the system is inherently anisotropic. Hence, the relevant observable is the stress tensor.

For isotropic systems, the thermodynamic definition of the pressure as the derivative of the free energy with respect to volume leads to the general equation

$$pV=Nk_\mathrm{B}TV\frac{U}{V},$$
(18)
where the angular brackets $`\mathrm{}`$ denote a canonical average. For the case of short-range interactions the contribution from $`U`$ can be further simplified by using

$$\frac{U}{V}=\underset{i}{}\frac{U}{𝒓_i}\frac{𝒓_i}{V}=\underset{i}{}(𝑭_i)\frac{𝒓_i}{3V}.$$
(19)
Substituting this into the pressure equation (18) and using $`𝑭_{ij}=𝑭_{ji}`$ gives

$$pV=Nk_\mathrm{B}T\frac{1}{3}\underset{i<j}{}𝒓_{ij}𝑭_{ij}$$
(20)
The electrostatic contribution to the pressure can be derived  by realizing that the electrostatic energy $`U^\mathrm{C}`$ is a homogeneous function of volume with degree $`1/3`$, and hence by Euler’s theorem $`U^\mathrm{C}/V=U^\mathrm{C}/3V`$. Adding this to Eqn. (20) finally gives the desired pressure equation

$$pV=Nk_\mathrm{B}T\frac{1}{3}\underset{i<j}{}𝒓_{ij}𝑭_{ij}^{\mathrm{sr}}+\frac{1}{3}U^\mathrm{C},$$
(21)
where $`𝑭_{ij}^{\mathrm{sr}}`$ are the short-range pair forces. Note that for a repulsive hard core the virial contribution is positive, while the electrostatic contribution is one third of the energy density, which for neutral systems is usually negative.

For anisotropic systems the stress tensor is the relevant observable to compute. Whereas the ideal gas contribution to the pressure still remains isotropic, the scalar product in the virial is replaced by the tensor product “$``$”. For the case of electrostatic interactions the derivation is more complicated and will not be presented here. The reader is referred to Ref. . The result is that within the framework of Ewald techniques the electrostatic contribution to the stress tensor, $`\text{p}^\mathrm{C}`$, can be decomposed additively into a real space contribution $`\text{p}^{(r)}`$ and a Fourier space contribution $`\text{p}^{(k)}`$. They are given by

$`\text{p}^{(r)}V`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{i,j}{}}q_iq_j{\displaystyle \underset{𝒎^3}{\overset{}{}}}\left[{\displaystyle \frac{2\alpha }{\sqrt{\pi }}}\mathrm{e}^{\alpha ^2𝒓_{ij𝒎}^2}+{\displaystyle \frac{\text{erfc}(\alpha |𝒓_{ij𝒎}|)}{|𝒓_{ij𝒎}|}}\right]{\displaystyle \frac{𝒓_{ij𝒎}𝒓_{ij𝒎}}{|𝒓_{ij𝒎}|^2}}`$ (22)
$`\text{p}^{(k)}V`$ $`=`$ $`{\displaystyle \frac{1}{2V}}{\displaystyle \underset{𝒌\mathrm{𝟎}}{}}{\displaystyle \frac{4\pi }{k^2}}\mathrm{e}^{k^2/4\alpha ^2}\left|\stackrel{~}{\rho }(𝒌)\right|^2\left[𝕀2\left(1+{\displaystyle \frac{k^2}{4\alpha ^2}}\right){\displaystyle \frac{𝒌𝒌}{k^2}}\right].`$ (23)
Here, $`𝒓_{ij𝒎}=𝒓_i𝒓_j+𝒎L_\mathrm{b}`$, and the canonical average has not been denoted explicitly. For details of the notation and the Ewald techniques in general, the reader is referred to the literature . Note that the connection to the isotropic case requires the trace of the pressure tensor to be equal to the mean electrostatic energy density, i.e., $`\mathrm{Tr}\text{p}^\mathrm{C}=U^\mathrm{C}/V`$.

After computing the electrostatic stress tensor, it can be rotated such that the rod points along the $`z`$-axis. The $`xx`$\- and $`yy`$-components of the new tensor then give the pressure perpendicular to the rod, while the $`zz`$-component is the contribution parallel to the rod.

List of Tables:

List of Figures