format: text

text:
# Adsorption of polyelectrolytes from semi-dilute solutions on an oppositely charged surface

## I Introduction

When a charged solid surface (charge density $`\sigma _0`$) is exposed to a solution of polyelectrolytes carrying opposite ions ($`Z`$ charges per chain), a counter-ion exchange occurs such that the chains adsorb and the counter-ions initially retained captive by electrostatic attraction are released in the solution. This is because the system gains the translational entropy of the micro-ions, which is $`Zk_BT`$ (where $`k_BT`$ is the thermal energy), and only looses $`k_BT`$ for each adsorbed chain. Then at sufficiently high charge density, adsorbed polyelectrolytes form a dense layer of coils strongly bounded to the solid. In some favorable cases, because the macro-ions have very long tails dandling in the solution, the adsorbed layer alone may carry a surface charge density larger than $`\sigma _0`$. When the solution has been rinsed away, the decorated surface (solid surface plus attached chains) then behaves as a new charged system surrounded by counter-ions. This feature is the basic principle of a new device to build controlled charged multi-layers of alternated sign Decher .

On the theoretical side, the structure of the adsorbed layer is known in the dilute regime, where the chains in the bulk are isolated coils, and the solution does not affect much the layer Joanny ; DobryninDeshkovski ; Dobrynin ; Dobrynin2 ; varoqui ; stoll1 ; stoll2 . Comparatively, much less is known about adsorption from semi-dilute solutions. Because the dilute regime for polyelectrolytes is found at vanishing monomer concentrations ($`\varphi _b^{}N^2`$ where $`N`$ is the index of polymerization), the semi-dilute regime is certainly relevant from an experimental point of view. Note that the Decher process takes place in the semi-dilute regime in Ref. Decher2 . Borukhov et al. have numerically solved the non-linear mean-field system of equations in the semi-dilute regime for repulsive surfaces Borukhov . They also used scaling arguments to determine the relevant characteristic length scale in terms of the total electrostatic potential drop, but no detailed analytical description of the layer is done. Besides, the case of attractive surfaces does not simply follow from the repulsive case explicitly treated by Borukhov et al. since the concentration profiles are significantly affected by the boundary conditions. Châtellier and Joanny have solved the linearized mean-field system of equations suitable to describe how the charged surface affects the semi-dilute solution Chatellier . However, this approach is restricted to small perturbations. By solving numerically the self-consistent field equations using the ground state dominance approximation, Wang explored the influence of various parameters (surface charge, salt concentration, bulk concentration) on the occurrence of charge inversion wang . The latter was found to be strong at high salt concentration.

In any case, these works do not provide the extension, $`H`$, of the adsorbed layer (which does not identify with the characteristic relaxation length of the concentration), and the amount of material, $`\mathrm{\Gamma }`$, attached to the solid surface (which does not identifies with the excess material driven to the interface $`(\varphi \varphi _b)𝑑z`$, where $`\varphi `$ is the monomer concentration, as defined in Refs. Borukhov ; Andelman ; wang ). There is a technical reason for that : the mean-field approach for adsorbed polyelectrolytes in semi-dilute conditions has to be supplemented by another approach to provide the full picture. This is because the variations of the polymeric concentration and the electrostatic potential does not carry enough information to completely describe the layer. This matters in practical situation since $`\mathrm{\Gamma }`$ is directly related to the charge density carried by the decorated surface after removal from the solution. Accordingly, this quantity deserves a special attention. Moreover, it has been seen experimentally that the extension and the amount of material of the adsorbed layer vary significantly with $`N`$ Chibowski . This effect is not explained using the mean-field theory in the ground-state dominance approximation PGGbook  alone, which is used in Refs. Borukhov ; Chatellier ; Andelman ; wang . Indeed within this approach, the semi-dilute solution is described with $`\varphi `$ and the adsorbed chains cannot be distinguished from free chains in the semi-dilute layer.

In this article, we propose a complete description of the adsorbed layer in equilibrium with the semi-dilute solution. In particular, we estimate $`\mathrm{\Gamma }`$ and $`H`$ as a function of $`N`$, $`\sigma _0`$ and $`\varphi _b`$. We use a mean-field approach supplemented by a scaling type of description. The mean-field theory that we consider is the celebrated Poisson-Boltzmann-Edwards’ (PBE) description of the polyelectrolyte solution (see the reviews Andelman ; Barrat ). The scaling approach that we use is the Scaling Functional Theory (SFT) Manomacromol  that we adapt to treat semi-dilute solutions of polyelectrolytes.

We shall consider a rather standard situation: long ($`N1`$), linear, fully flexible (the persistence length is identified with the monomer size $`a`$), polymer chains, small fraction of charges ($`f1`$), in $`\mathrm{\Theta }`$ solvent conditions ($`v=0`$, and the monomer/monomer interaction is dominated by three body excluded volume interactions), bare surface charge $`\sigma _0`$. We assume no added salt and no specific interaction with the surface, as a first hint into the full problem. For definiteness, we suppose (without loss of generality) that the surface is negatively charged with a bare surface charge density $`\sigma _0`$ (with $`\sigma _0>0`$), and the chain carries $`Z=fN`$ quenched elementary positive charges $`+e`$ where $`f`$ is the fraction of charged monomers.

In view of the number of parameters, this issue is a formidable task. To proceed further, our analysis is based on two major assumptions :
1. That we may find long polyelectrolyte tails belonging to adsorbed chains far from the surface (protruding in the solution), will be our first assumption (see Fig. 1). This hypothesis is inspired by the results found with long neutral chains adsorbed from semi-dilute solution, where we know that the concentration field decreases over a distance comparable to the bulk blob size $`\xi _b\varphi _b^{3/4}`$ (the characteristic length scale of the solution), whereas the layer of adsorbed chains extend over much larger distances $`HR_GN^{1/2}\varphi _b^{7/8}`$ Manomacromol ; Marques ; Daoud . In this example, the essential idea is that the collective response of the chains screens the perturbation induced by the surface at a characteristic length scale which does not depend on the chain size PGG81 . A similar feature is found with polyelectrolytes: in the semi-dilute regime, the Debye-Hückel length of the solution, $`\kappa ^1`$, is independent of the chain length. One may possibly argue that polyelectrolytes lie quite flat on the surface, as with dilute solutions. This is however very unlikely since for such picture extrapolated to very high concentrations of polyelectrolyte, it would mean that the adsorbed chains adopt an almost two-dimensional configuration. Clearly, the system of adsorbed chains would gain large amounts of conformational entropy by allowing more chains to adsorb and recovering an almost Gaussian structure (thus leaving some tails dandling in the solution), without perturbing much the distribution of charges in the vicinity of the solid. In other words, since $`\kappa ^1R_GN^{1/2}`$ in the limit $`N1`$, we expect that the influence of the surface will be small for the chain structure in equilibrium.
2. Looking at the variations of the electrostatic potential, $`\mathrm{\Psi }(z)`$ with the distance to the surface, $`z`$, in the dilute regime, we anticipate a very simple form for $`\mathrm{\Psi }(z)`$ in the semi-dilute regime (see Fig. 2). From the mean-field description of the layer in the dilute regime Joanny ; DobryninDeshkovski , we learn that the electrostatic potential first increases from negative values, then saturates, and eventually relaxes to zero. That there is a regime of parameters in the semi-dilute regime where such profile for $`\mathrm{\Psi }`$ will persist to some extent be our second assumption. This assumption will be justified by the calculation but one can find a simple argument for it. The monomer concentration profile in the dilute regime tells us that the osmotic pressure decreases from the solid surface. Suppose that we increase the concentration of the solution from the dilute regime into the semi-dilute regime. Because the osmotic pressure close to the solid surface is much higher than what is found in the solution, we expect that the presence of the solution will affect the structure of the outer fringe of the layer, not the structure of the inner regions. Accordingly, we do not expect that the form of the electrostatic potential will be profoundly affected, at least in the semi-dilute regime close to the dilute regime.

Following our second hypothesis, we define the plane of vanishing electrostatic field (such that $`d\mathrm{\Psi }/dz=0`$) as the neutral plane. Similarly, we set the plane where the electrostatic potential vanishes (such that $`\mathrm{\Psi }0`$) as the relaxation plane (we implicit assume that we haven chosen $`\mathrm{\Psi }(\mathrm{})=0`$). Then the layer might be divided into three main regions (labeled (1), (2) and (3), see Fig. 2) such that:

1. The compensation region stands between the solid surface and the neutral plane. By definition, the charges provided by the polyelectrolyte material situated in this region exactly compensate the surface charge.
2. The intermediate region stands between the neutral plane and the relaxation plane.
3. The outer region stands between the relaxation plane and the outer border of the layer, as defined by the longest polymeric coil in direct contact with the surface.

It is important to realize that this distinction, although a priori, is not arbitrary. Each of these region is well defined from a physical point of view.

In what follows, we first introduce the Scaling Functional Theory suitable to describe polyelectrolyte semi-dilute layers (Section 2). Then we propose a detailed picture for the adsorbed layer in Section 3. Section 4 presents our concluding remarks focusing on the issue of charge inversion. Scaling laws in polyelectrolyte solutions and the mean-field approach are briefly introduced in Appendices 1 and 2 respectively.

## II Scaling functional theory

This theory was first introduced to describe neutral polydisperse polymeric layers by generalizing the Alexander-de Gennes model for brushes AGR . The basic principle is to describe the layer of adsorbed chains as a thermodynamical ensemble of tails and loops where the entropy associated with the polydispersity in size compete with the elastic energy and the monomer excluded volume. From a formal point of view, the SFT approach is a functional theory where the Hamiltonian is simplified according as follows ManoPRE : a) loops are cut into pseudo-tails, b) the tails and pseudo-tails are described by a single function $`z`$, where $`z(n)`$ is the position of the nth monomer, $`n`$, of each tail. This idea proved successful in describing various situation ranging from polymer brushes, reversible adsorption, irreversible adsorption, whatever the solvent conditions.

Our aim in this article is not to generalize the SFT to describe polyelectrolyte chains. Rather, we take advantage of the fact that in the semi-dilute regime, the chain is both Gaussian and neutral at large scales, and the system of polymers is thus analogous to a neutral solution provided that the blob is renormalized according to the physics at small scales. Essentially, our approach consists in integrating the degrees of freedom of the various species (including electrostatics interactions) into a proper description of the blob.

Within the SFT, the system of loops and tails attached to the solid surface is described by three functions: $`S`$, the local grafting density of tails and pseudo-tails (a function of the curvilinear distance $`n`$), $`z`$, the path, and $`\varphi `$, the volume fraction of monomers. These functions are linked through the conservation of monomers:

$$\varphi (z)=\frac{S(n(z))}{\dot{z}(n(z))}$$
(2.1)
For polyelectrolyte semi-dilute layers, the effective free-energy (per $`cm^2`$) of the system write ($`\beta ^1=k_BT`$):

$$\beta F_0^N𝑑n\left\{\frac{\left(a^2S(n)\right)^{3/2}}{\xi _e^{3/2}\dot{z}^{1/2}(n)}+a^1\xi _e^{1/2}S^{3/2}(n)\dot{z}^{3/2}(n)S^{}(n)\mathrm{ln}\left(\frac{S^{}(n)}{S_0}\right)\right\}$$
(2.2)
where $`S_0`$ is the “grafting density” at the surface and $`\xi _ea(f^2\mathrm{}_B)^{1/3}`$ is the electrostatic blob size (in $`\mathrm{\Theta }`$ solvent conditions, $`\mathrm{}_B`$ is the Bjerrum length in units of $`a`$, cf. Appendix 1). The effective free-energy eq. (2.2) is the sum of three contributions which account respectively for :
a) monomer excluded volume interactions. This contribution is $`\mathrm{\Pi }𝑑z`$, where the polymeric contribution to osmotic pressure, $`\mathrm{\Pi }`$, scales as $`k_BT/\xi ^3`$:

$$\mathrm{\Pi }k_BT\frac{\varphi ^{3/2}}{a^{3/2}\xi _e^{3/2}}$$
(2.3)
With eq. (2.1), it is simple to show that

$$_0^H\mathrm{\Pi }𝑑zk_BT_0^N\frac{(a^2S(n))^{3/2}}{\xi _e^{3/2}\dot{z}^{1/2}(n)}𝑑n$$
(2.4)
where $`H=z(N)`$ is the position of the last monomer and is identified as the layer thickness.
b) the elasticity of the tail at a scale larger than the blob size (after integration by parts)

$$\beta F_{el}=_0^N𝑑nS^{}(n)_0^n\left(\frac{dz/\xi }{dm/g}\right)^2\frac{dm}{g}a^1_0^N\xi _e^{1/2}S^{3/2}(n)\dot{z}^{3/2}(n)𝑑n$$
(2.5)
c) the entropy associated to the polydispersity in size of loops and tails which is similar for charged and neutral monomers as soon as polymers are flexible AGR .

Minimizing the free-energy eq. (2.2) with respect to $`S`$ and $`\dot{z}`$ yields two equations

$`{\displaystyle \frac{3}{2}}a^3\xi _e^{3/2}\left({\displaystyle \frac{S}{\dot{z}}}\right)^{1/2}+{\displaystyle \frac{3}{2a}}\xi _e^{1/2}S^{1/2}\dot{z}^{3/2}+{\displaystyle \frac{S^{\prime \prime }}{S^{}}}`$ $`=`$ $`\mu _b`$ (2.6)
$`{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{a}{\xi _e}}\right)^{3/2}\left({\displaystyle \frac{S}{\dot{z}}}\right)^{3/2}{\displaystyle \frac{3}{2a}}\xi _e^{1/2}S^{3/2}\dot{z}^{1/2}`$ $`=`$ $`\mathrm{\Pi }_b`$ (2.7)
These equations are formally similar to what we find in the neutral case (and their interpretation in physical terms is thus identical), except that the small scale structure introduces different scaling relationship. Equation (2.7) describes the local balance of forces (per unit surface parallel to the solid): elasticity (first term of lhs) competes with osmotic pressure (second term of lhs) and the bulk pressure (rhs). Note that the osmotic pressure of counter-ions is the same in the outer region and in the bulk so that it does not appear in eq. (2.7). Equation (2.6) describes the conservation of the generalized chemical potential. In equation (2.6), $`\mu _b=\mathrm{\Pi }_b/\varphi _b`$ is the chemical potential in the bulk and scales as $`g_b^1`$:

$$\mu _b\varphi _b^{1/2}\left(\frac{a}{\xi _e}\right)^{3/2}$$
(2.8)
The bulk osmotic pressure, $`\mathrm{\Pi }_b1/\xi _b^3`$, in equation (2.7) accounts for the pressure induced by the solution.

We emphasis that the SFT system of equations (2.6)–(2.7) assumes that the properties of the tails and loops are that of neutral strings of blobs, and may only apply to situations where the electrostatics interactions are screened at distances larger than the mesh size $`\xi `$ and especially very close to the substrate.

## III Structure of the layer

### III.1 Compensation region (1).

This region, close to the solid surface, is such that the charges provided by the polyelectrolyte material exactly compensate the surface charge. Hence

$$\mathrm{\Gamma }_1\sigma _0/f$$
(3.9)
and $`\frac{d\mathrm{\Psi }}{dz}0`$ in this region (at the outer boundary of the compensation region: $`\frac{d\mathrm{\Psi }}{dz}=0`$). For this reason, we expect that the counter-ions from polyelectrolytes are not present, $`\rho ^{}0`$. Moreover, at the low salt limit, positive counter-ions (from the charged surface) are negligible since their concentration is proportional to surface/volume and vanishes in the thermodynamic limit.

Each electrostatic blob carries $`fg_e`$ charges and occupies a surface $`\xi _e^2`$. These relation defines a critical value for $`\sigma _0`$ such that a single “carpet” of close packing electrostatic blob is able to fully compensate the surface charge:

$$\sigma _cfa^2$$
(3.10)
The monomer volume fraction inside this region is thus $`\varphi _e`$, given by Eq. (0.41). Moreover, the electrostatic energy of a blob interacting with the surface is of order $`k_BT`$ for the critical value, $`\sigma _0\sigma _c`$, i.e. of the same order of the electrostatic repulsion between monomers. Accordingly, we shall distinguish two cases: $`\sigma _0<\sigma _c`$ and $`\sigma _0>\sigma _c`$ (Fig. 3).

III.A.1 Case $`\sigma _0<\sigma _c`$.

In this regime, we rely on the results found by Dobrynin et al. DobryninDeshkovski  for the case of dilute solutions. The electrostatic attraction of charged monomers to the charged surface is weaker than the electrostatic repulsion between charged monomers. Therefore, the chain statistics is rod-like at scales larger than the blob $`\xi `$ (in the plane parallel to the surface) and Gaussian at smaller scales ($`\xi ag^{1/2}`$). However, the blob size is found by assuming that the electrostatic attraction to the surface is of order $`k_BT`$: $`\frac{fg\sigma _0\xi ^2e^2}{ϵ\xi }k_BT`$ which leads to $`\xi \xi _e\left(\frac{\sigma _c}{\sigma _0}\right)^{1/3}>\xi _e`$. The chain segment is thus confined by the electrostatic attraction to the charged surface. However, contrary to the case of dilute solutions (single chain adsorption), polyelectrolytes do not lie flat on the surface since, as mentioned in the Introduction, the system gains conformational entropy by allowing more chains to adsorb with large loops protruding in the solution. The concentration of blobs in this region is found by assuming that this compensation layer actually compensates the surface charges and that the thickness is given by the confinement blob size $`\xi `$:

$$\mathrm{\Lambda }_1\xi _e\left(\frac{\sigma _c}{\sigma _0}\right)^{1/3}$$
(3.11)
This scaling result has been shown numerically for semi-dilute solutions by Borukhov et al. Borukhov  and Shafir et al. Shafir  using the mean-field equations. The electric field vanishes at $`z=\mathrm{\Lambda }_1`$ in this regime. Of course, we do not have a closed packing of these blobs, $`\xi `$, since the number of blobs per unit surface is $`\sigma _0/fg(\sigma _0/\sigma _c)^{5/3}\xi _e^2<\xi _e^2`$.

This picture is valid only if the average monomer volume fraction close to the surface is higher than the bulk value, i.e. $`\varphi _e(\sigma _0/\sigma _c)^{4/3}>\varphi _b`$ or

$$\sigma _c\left(\frac{\varphi _b}{\varphi _e}\right)^{3/4}\sigma _0\sigma _c$$
(3.12)
III.A.2 Case $`\sigma _0>\sigma _c`$.

In the limit $`\sigma _0\sigma _c`$, the monomer volume fraction close to the surface, $`\varphi _0`$, increases being larger than $`\varphi _e`$, the electrostatic attraction becomes larger than the electrostatic intra-repulsion, and the osmotic pressure tends to swell the compensation region. Assuming that excluded volume interactions dominate over elasticity, as in Ref. DobryninDeshkovski , solving eq. (0.47) of Appendix 2 yields (at the leading order in $`\varphi _e/\varphi _01`$ and assuming $`\mu g_e^1`$ in this compensation region) a parabolic profile: $`\varphi (z)\varphi _e+(\varphi _0\varphi _e)(1z/\mathrm{\Lambda }_1)^2`$, where $`\varphi _0=(3\pi \mathrm{}_B\sigma _0^2)^{1/3}`$ and

$$\mathrm{\Lambda }_1\xi _e\left(\frac{\sigma _0}{\sigma _c}\right)^{1/3}$$
(3.13)
The structure of the compensation region is thus self-similar (see Fig. 3) as already described for the dilute case in Ref. DobryninDeshkovski . For the general case, the polymer elasticity is no more negligible and the full equation (0.47) must be solved without any approximation.

It is important to note that our result Eq. (3.13) is in contradiction with Andelman et al. Borukhov ; Andelman ; Shafir . They predict that the characteristic length scale of the adsorbed layer decreases with $`\sigma _0`$ (Andelman , eq. (71)) following Eq. (3.11). However, we argue that, for strong adsorption ($`\sigma _0\sigma _c`$), the structure of the compensation region is found by balancing the electrostatic adsorption with the osmotic pressure (excluded volume interactions) and that chains are not confined. Instead of blobs of confinement, the structure is described in terms of osmotic blobs, in a similar way of describing neutral polymer layers.

Anyway, as we argue below (Section 4), in the semi-dilute regime and for long polyelectrolytes, we do not expect that this region has a dominant influence on observables such as the thickness of the layer, or the amount of material attached to the surface, and it is enough for our purpose to treat the restricted case where $`\sigma _0\sigma _c`$, and

$$\mathrm{\Lambda }_1\xi _e$$
(3.14)
Anticipating the use of the SFT, it is useful to estimate the amount of monomers (per tail), $`n_1`$, involved in the compensation region. In a simple view where each tail contributes to one electrostatic blob, we find

$$n_1g_e$$
(3.15)
and the local grafting density at the outer limit of the compensation region scales as

$$S(n_1)\frac{1}{\xi _e^2}$$
(3.16)
### III.2 Intermediate region (2).

In this region, the direct electrostatic influence of the surface vanishes, and from this point of view, the structure should be similar to a layer of polyelectrolytes chains “adsorbed” onto a fictive neutral surface positioned at $`z=\mathrm{\Lambda }_1`$. Here, the concentration of monomers decreases from $`\varphi _e`$ to the bulk value, $`\varphi _b`$. Because monomers and counter-ions are ruled by very different equations, the electrostatic potential $`\mathrm{\Psi }(z)`$ does not vanish. As the system is globally neutral and $`\varphi (\mathrm{\Lambda }_1)=\varphi _e>\varphi _b`$ in the semi-dilute regime, charge conservation in the system implies that the volume fraction of counter-ions, $`\rho ^{}(z)`$, will decrease from $`\rho ^{}(\mathrm{\Lambda }_1)>f\varphi _b`$ to the bulk value $`f\varphi _b`$ at the external border of region (2). Knowing that $`\rho ^{}`$ is related to the adimensional electrostatic potential $`\mathrm{\Psi }=\beta eV`$ by $`\rho ^{}=f\varphi _b\mathrm{exp}\mathrm{\Psi }`$ (see Appendix 2), we deduce that $`\mathrm{\Psi }(\mathrm{\Lambda }_1)=\mathrm{\Psi }_1`$ is positive (with the boundary condition $`\mathrm{\Psi }(\mathrm{})=0`$). Thus, we expect that $`\mathrm{\Psi }(z)`$ decreases monotonically from a finite positive value, $`\mathrm{\Psi }_1`$ (to be estimated below), to 0.

When $`\varphi _b`$ is not too small in comparison to $`\varphi _e`$, say $`\varphi _e\varphi _b\varphi _b`$, the linearized PBE system of equations described in Appendix 2 is justified to find the relaxation of the concentration profile. Note that for the case $`\varphi _b\varphi _e`$ the full system should be solved numerically Borukhov ; Shafir . With the appropriate boundary conditions \[$`\varphi (\mathrm{\Lambda }_1)=\varphi _e`$, $`\varphi (\mathrm{})=\varphi _b`$, $`d\mathrm{\Psi }/dz(\mathrm{\Lambda }_1)=0`$ and $`\mathrm{\Psi }(\mathrm{})=0`$\], we find:

$`\phi (z)`$ $`=`$ $`\sqrt{\varphi _b}+\left(\sqrt{\varphi _e}\sqrt{\varphi _b}\right)\left[(1R)e^{(z\mathrm{\Lambda }_1)/\mathrm{}_{}}+Re^{(z\mathrm{\Lambda }_1)/\mathrm{}_+}\right]`$ (3.17)
$`\mathrm{\Psi }(z)`$ $`=`$ $`8\pi \mathrm{}_Bf\sqrt{\varphi _b}(\sqrt{\varphi _e}\sqrt{\varphi _b})[{\displaystyle \frac{1R}{4\pi \mathrm{}_Bf\varphi _ba^2\lambda _{}}}`$ (3.18)
$`\times `$ $`\mathrm{e}^{(z\mathrm{\Lambda }_1)/\mathrm{}_{}}+{\displaystyle \frac{R}{4\pi \mathrm{}_Bf\varphi _ba^2\lambda _+}}\mathrm{e}^{(z\mathrm{\Lambda }_1)/\mathrm{}_+}]`$
$`\rho (z)`$ $`=`$ $`f\varphi _b[\mathrm{\Psi }(z)+1]`$ (3.19)
where

$`\mathrm{}_\pm ^2`$ $`=`$ $`\lambda _\pm ={\displaystyle \frac{2\varphi _b}{a^2}}\left(\pi \mathrm{}_Bf+6\varphi _b\right)[1\pm \sqrt{\mathrm{\Delta }}]`$ (3.20)
$`\mathrm{\Delta }`$ $`=`$ $`1{\displaystyle \frac{12\pi \mathrm{}_Bf(2\varphi _b^2+f)}{\varphi _b(\pi \mathrm{}_Bf+6\varphi _b)^2}}`$ (3.21)
$`R^1`$ $`=`$ $`1\sqrt{{\displaystyle \frac{\lambda _+}{\lambda _{}}}}{\displaystyle \frac{(4\pi \mathrm{}_Bf\varphi _ba^2\lambda _{})}{(4\pi \mathrm{}_Bf\varphi _ba^2\lambda _+)}}`$ (3.22)
In Fig. 4 are plotted $`\delta \phi (z)=\phi (z)\sqrt{\varphi }_b`$ and $`\mathrm{\Psi }(z)`$ for $`\mathrm{}_B=1`$, $`f=0.1`$, and $`\varphi _b=0.26`$ (concentrated semi-dilute solution). Note that the characteristic relaxation lengths $`\mathrm{}_\pm `$ were obtained by Châtellier and Joanny Chatellier . In this reference, the possibility of damped oscillations in the volume fraction profile ($`\mathrm{}_\pm `$ imaginary) was carefully examined for lower values of $`\varphi _b`$. Here, we will only consider the situation where $`\mathrm{}_\pm `$ are real numbers, for simplicity (the same discussion can be carried on with the real parts of $`\mathrm{}_\pm `$ in more general cases). Hence the surface induce an external perturbation to $`\phi `$ of magnitude $`\sqrt{\varphi _e}\sqrt{\varphi _b}`$ which relaxes, at the linear response level, to the bulk value, $`\sqrt{\varphi _b}`$ on the characteristic wave length of the fluctuations in the bulk.

Of course, we identify $`\mathrm{\Lambda }_2`$ with the longest relaxation length, $`\mathrm{}_{}`$:

$$\mathrm{\Lambda }_2\mathrm{}_{}$$
(3.23)
It is instructive to compare the characteristic relaxation length $`\mathrm{}_{}`$ with $`\xi _b`$, which requires that we simplify the expression given by eq. (3.20). Assuming that the intermediate region is dominated by the counter-ions pressure with an effective excluded volume $`v_{\mathrm{eff}}f/\varphi _b`$ (and that contribution arising from the monomers is negligible), the relaxation length $`\mathrm{}_{}`$ writes (with $`\mathrm{\Delta }1`$):

$$\mathrm{}_{}\frac{a}{\sqrt{12f}}$$
(3.24)
In this limit, we have $`\mathrm{}_{}\xi _b`$ for $`\varphi _b(f/\mathrm{}_B)^{1/3}`$, a reasonable case when $`\varphi _b<\varphi _e`$. The amount of monomers is this region is given by

$$\mathrm{\Gamma }_2=a^3_{\mathrm{\Lambda }_1}^{\mathrm{\Lambda }_1+\mathrm{\Lambda }_2}\varphi (z)𝑑z(0.17\varphi _b+0.43\varphi _e+0.40\sqrt{\varphi _b\varphi _e})\frac{\mathrm{}_{}}{a^3}$$

Equation (3.18) also yields an estimate for the extremum value of the electrostatic potential $`\mathrm{\Psi }_1=\mathrm{\Psi }(\mathrm{\Lambda }_1)`$:

$$\mathrm{\Psi }_18\pi \mathrm{}_Bf\sqrt{\varphi _b}\left(\sqrt{\varphi _e}\sqrt{\varphi _b}\right)\left[\frac{1R}{4\pi \mathrm{}_Bf\varphi _ba^2\lambda _{}}+\frac{R}{(4\pi \mathrm{}_Bf\varphi _ba^2\lambda _+)}\right]$$
(3.25)
Equation (3.25) in turn allows to estimate $`\mathrm{\Psi }_0=\mathrm{\Psi }(z=0)`$ through the integrated Poisson equation:

$$\mathrm{\Psi }_0=\mathrm{\Psi }_14\pi \mathrm{}_B\sigma _0\mathrm{\Lambda }_1+4\pi \mathrm{}_Bf_0^{\mathrm{\Lambda }_1}_0^u\varphi (z)𝑑z𝑑u$$
(3.26)
For the moderate adsorption case ($`\sigma _0\sigma _c`$) that we consider in this article, we find:

$$\mathrm{\Psi }_0\mathrm{\Psi }_12\pi \mathrm{}_B\frac{fg_e}{\xi _e}$$
(3.27)
For the case of strong adsorption ($`\sigma _0\sigma _c`$), we obtain:

$$\mathrm{\Psi }_0\mathrm{\Psi }_12\pi \mathrm{}_B\frac{fg_e}{\xi _e}\left(\frac{\sigma _0}{\sigma _c}\right)^{4/3}$$
(3.28)
The SFT approach described in Section 2 is unlikely to describe correctly the intermediate region. At a simple scaling level, we will note that the local grafting density should drop from $`S(n_1)\frac{1}{\xi _e^2}`$ to $`S(n_2)\frac{1}{\xi _b^2}`$, where $`n_2`$ is the curvilinear index at $`z=\mathrm{\Lambda }_1+\mathrm{\Lambda }_2`$. A crude estimate for $`n_2`$ is found with a renormalized linear string approximation for the tails in the intermediate regime:

$$n_2n_1+g_e\frac{\mathrm{\Lambda }_2}{\xi _e}g_b$$
(3.29)
which is equivalent to write

$$\mathrm{\Lambda }_2\xi _b\xi _e\xi _b$$
(3.30)
where the last approximation holds in the limit $`\varphi _e\varphi _b`$. Hence, the thickness of region (2) is then given by equation (3.23) at a mean-field level and for $`\varphi _e\varphi _b\varphi _b`$, and by equation (3.30) at a scaling level.

### III.3 Outer region (3).

Here, the polymer concentration (resp. the counterion concentration, the adimensional electrostatic potential) is the bulk concentration, $`\varphi _b`$ (resp. $`f\varphi _b`$, 0 at scales larger than the bulk blob size), but the tails and loops in this region belong to chains in direct contact with the solid surface. At a scaling level of description, the layer is a close packing of semi-dilute blobs (size $`\xi `$). By analogy with the Rubinstein description of the bulk solution Rubinstein , we expect that electrostatic interactions are screened at distances larger than the mesh size $`\xi `$, and the system is amenable to a description in terms of neutral chains. Hence the SFT approach put forward in Section 2 is now appropriate.

Obviously, free and adsorbed chains interpenetrate, and we expect a slowly vanishing concentration profile for the monomers which belong to chains in direct contact with the surface, as depicted in Fig. 5a. In what follows, however, we assume that the free chains do not penetrate in the layer of tails and loops in direct contact with the surface. Accordingly, we consider a sharp, Heaviside type of profile at the outer fringe (see Fig. 5b). A similar type of simplification is found in various theories of polymeric layers. In the context of polymer brushes in good solvent, e. g., Alexander and de Gennes assume an Heaviside step for the concentration profile in order to find the extension of the layer Alexander ; PGGbrush . Their scaling result was confirmed by more refined theories where this constraint is relaxed (instead a parabolic profile is found) MWC . Similarly, Guiselin assumes no interpenetration in the context of molten solutions of neutral polymers adsorbed at an interface Guiselin . Here again, the scaling result was experimentally confirmed by measuring the amount of material irreversibly adsorbed Liliane . Because the semi-dilute solution is amenable to a description in terms of neutral blobs of strings of electrostatic coils PGGPincus ; Rubinstein , it is likely that the situation is analogous. In any case, we are confident that the scaling results that we find will not be dramatically affected by our simplification.

The SFT equations suitable to describe this region are (2.6)–(2.7). In equilibrium, we expect that the bulk osmotic pressure exerted at the extremity of the layer ($`z=H`$) is balanced by the inner osmotic pressure, and the elasticity can be neglected in eq. (2.7). This yields

$$\frac{S}{\dot{z}}\varphi _b$$
(3.31)
and thus \[eq. (2.1)\], $`\varphi (z)\varphi _b`$. When this relation is introduced into eq. (2.6), this equation simplifies into

$$\frac{\xi _e^{1/2}}{\varphi _b^{3/2}}S^2+\frac{S^{\prime \prime }}{S^{}}0$$
(3.32)
which admits a power law type of solution in the limit $`N\mathrm{}`$,

$$S(n)\frac{\varphi _b^{3/4}}{a^2(\xi _e/a)^{1/4}n^{1/2}}$$
(3.33)
Note that we find $`S(n_2)\xi _b^2`$. Now solving Eqs. (3.31)–(3.33) (with the boundary condition $`z(n_2)\mathrm{\Lambda }_1+\mathrm{\Lambda }_2`$) for $`z`$, gives a Gaussian structure at large scale:

$$z(n)\mathrm{\Lambda }_1+\mathrm{\Lambda }_2+a(\xi _e/a)^{1/4}\varphi _b^{1/4}n^{1/2}$$
(3.34)
which readily yields \[$`H=z(N)`$\]:

$$\mathrm{\Lambda }_3a\frac{N^{1/2}}{(\xi _e/a)^{1/4}\varphi _b^{1/4}}$$
(3.35)
and

$$\mathrm{\Gamma }_3=_{n_2}^NS(n)𝑑n\frac{\varphi _b^{3/4}N^{1/2}}{a^2(\xi _e/a)^{1/4}}$$
(3.36)
Note that the result eqs. (3.35)–(3.36) are precisely what we would find with scaling arguments assuming that the chains are Gaussian at large scale, and the SFT is only a formal way to recover them in this context. Figure 6 shows a scaling representation of the whole adsorbed layer in the case $`\sigma _0\sigma _c`$. On the other hand, $`S`$ is directly related to the loop size distribution $`P`$ through ManoPRE

$$S(n)=_n^NP(u)𝑑u$$
(3.37)
and eq. (3.33) may be experimentally checked with a surface force apparatus senden .

## IV Discussion

The structure that we propose for the adsorbed layer clearly separates the influence of the charged surface which occurs only in the compensation region, from that of the bulk solution which act on the outer region, with an intermediate regime in between. We briefly describe the influence of $`\varphi _b`$ and $`\sigma _0`$ on the layer structure and then discuss the occurrence of charge inversion for layers of very long charged chains.

We expect that as we increase $`\varphi _b`$ from $`\varphi _b^{}`$, we progressively build the outer region at the expense of the intermediate region ($`\mathrm{\Lambda }_2`$ decreases), but without much consequences for the compensation region. This is because the osmotic pressure of the counter-ions (the dominant contribution to the physics in the intermediate region) is much higher near $`z=\mathrm{\Lambda }_1`$ than at $`z=\mathrm{\Lambda }_2`$. Arguably, any increase of the bulk concentration of counter-ions will predominantly affect the low osmotic domain of the intermediate region, and hardly the structure near $`z=\mathrm{\Lambda }_1`$, where the pressure is high. Of course, when $`\varphi _b`$ reaches $`\varphi _e`$, the intermediate region disappears in our simple view, and the bulk solution should then affect directly the compensation region. In this regime where fluctuations are negligible, numerical mean-field approaches are well adapted Borukhov ; wang . On the other hand, if we increase $`\sigma _0`$ from $`\sigma _c`$, at constant bulk concentration, the surface potential $`|\mathrm{\Psi }_0|`$ increases and more and more monomers are involved in compensating the surface charge ($`n_1`$ and $`\mathrm{\Lambda }_1`$ increase). Accordingly, the concentration profile in the compensating region is affected. But for very long chains ($`Nn_2`$), we do not expect that this should influence much the structure of the outer region.

Note that our description of the compensation and intermediate region resembles, but is not quite similar to the description of the adsorbed layer from dilute solution proposed in Ref. DobryninDeshkovski  : the so-called self-similar adsorbed layer regime. Essentially, we keep the parabolic decrease of the monomer volume fraction close to the surface, until we reach $`\varphi =\varphi _e`$ (where neglecting the chain elasticity is no more justified), and then we use the mean-field equations (which account for the elastic contribution and the osmotic pressure of counter-ions) to find the relaxation of $`\varphi `$ below $`\varphi _e`$. In Ref. DobryninDeshkovski , the region where $`\varphi \varphi _e`$ is described in terms of electrostatic blobs. However such a description does not take into account the presence of the counter-ions in this region which is imposed by electroneutrality (the electric potential is positive). In technical terms, we make advantage of the fact that $`\varphi `$ relaxes to a finite value, $`\varphi _b`$, to determine the monomeric profile below $`\varphi _e`$ by linearizing the mean-field equations. Because $`\varphi `$ goes to $`0`$ in dilute solutions, such calculation is not possible and would necessitate a numerical resolution. The difference between the two descriptions is not just a matter of refinement since the relaxation of $`\varphi `$ yields the variations of the electrostatic potential. With our approach, we provide both a coherent and complete picture for the variations of this quantity \[eq. (3.18)\].

A crucial aspect of adsorbed polyelectrolyte layers is their ability to capture a very large number of counter-ions per unit surface, much larger indeed than $`\mathrm{}_B^2`$. In particular, the driving force for the adsorption of oppositely charged chains depends precisely on the amount of potentially free ions. With this in mind, it is interesting to discuss the following situation.

Suppose that we build an adsorbed layer from a semi-dilute solution, as described before. Then we remove the polymeric solution after the layer has equilibrated, but we suppose that the adsorbed chains remain attached to the charged surface, in a fashion similar to the scenario discussed by Guiselin for irreversible adsorption of molten liquids of neutral chains Guiselin . Finally, we put the decorated layer into a solution of pure solvent. The apparent charge (also called ”net charge” in Ref. DobryninDeshkovski ) of the layer in the dilute solution will be zero, because of the counter-ions present within (or in the proximity) of the layer. On the other hand, the decorated layer carries a huge amount per unit surface of potentially free counter-ions. Provided that we disregard the fluffy structure of the layer, the decorated surface effectively behaves as a new charged surface of ”nominal” charge density

$`\mathrm{\Delta }\sigma \sigma \sigma _0`$ $`=`$ $`f(\mathrm{\Gamma }_1+\mathrm{\Gamma }_2+\mathrm{\Gamma }_3)\sigma `$ (4.38)
$`=`$ $`f\mathrm{\Gamma }_2+\mathrm{}_B^{1/12}f^{7/6}\varphi _b^{3/4}N^{1/2}`$
Hence, for long chains, we expect that $`\mathrm{\Delta }\sigma `$ may be much larger than $`\sigma _0`$.

Let us compare: a) a layer of polyelectrolytes reversibly adsorbed from dilute solution in equilibrium, and b) a layer of polyelectrolytes irreversibly adsorbed from semi-dilute solution put into a solution of pure solvent (the solid surfaces have the same $`\sigma _0`$ in both cases). To this aim, we need a precise vocabulary. We define three regimes, depending on the comparison of $`\sigma `$ with $`\sigma _0`$ and $`2\sigma _0`$:

* $`\sigma <\sigma _0`$: undercharging ($`\sigma =\sigma _0`$ corresponding to charge compensation);
* $`\sigma _0<\sigma <2\sigma _0`$: overcharging ($`\sigma =2\sigma _0`$ corresponding to charge inversion);
* $`\sigma >2\sigma _0`$: supercharging.

For case a) the mean-field theory predicts exact charge compensation in the absence of added salt Joanny . With added salt, Joanny finds overcharging, with possibly charge inversion in the limiting case of very high ionic force. In any case, there is no supercharging. It also appears that the surface controls the amount of adsorbed material in case a) (together with the ionic strength of the solution in case of added salt), not the chain length, nor the polymer concentration. A more recent analysis is provided by the work of Dobrynin et al. DobryninDeshkovski  who carefully combine the mean-field theory with scaling arguments. But this study does not change the conclusions concerning the regime of compensation.

In contrast, for case b) our analysis strongly suggests that we have supercharging for sufficiently long chains, $`N>(\sigma _0/f\mathrm{\Gamma }_2)^2/[(\mathrm{}_Bf^2)^{1/6}\varphi _B^{3/2}]`$, even in the absence of salt.

## V Concluding remarks

In this paper, we describe completely the structure of the adsorbed layer from semi-dilute solution by identifying the physical mechanisms in each sublayers : electrostatic attraction by the surface charge vs. three body interactions in the compensation region (1), importance of the polymer elasticity and the counter-ions osmotic pressure in the intermediate region (2), and presence of long tails and loops of adsorbed chains in the distal region (3) at the bulk concentration. We give explicit expressions for the coverage including all monomers belonging to adsorbed polymers, $`\mathrm{\Gamma }`$ and the layer thickness $`H`$ as a function of parameters, bare surface charge, bulk concentration which could be checked experimentally.

Clearly, our description of the intermediate region suffers from several weaknesses. First, we consider a mean-field approach extrapolated outside its domain of validity. Secondly, we use a linearized theory, which is only justified for bulk volume fractions close to $`\varphi _e`$. However, it is very likely that our essential scaling type of conclusions will be robust to any improved description of the intermediate region.

Similarly, the influence of the solvent quality (unless it is bad solvent) is rather limited. This is because the electrostatic interactions dominates over the excluded volume interactions at large scales, and in actual facts, the solvent only matters in the compensation region and in the electrostatic blob. For example, it is simple to show that $`\mathrm{\Gamma }_3a^2\varphi _b^{3/4}(\mathrm{}_Bf^2)^{1/14}N^{1/2}`$ and $`\mathrm{\Lambda }_3a\varphi _b^{1/4}(\mathrm{}_Bf^2)^{1/14}N^{1/2}`$ in good solvent conditions.

This study is limited to the very low salt limit. When the ionic force is very large, electrostatic interactions are screened at length scale larger than the Debye length and lead to an effective excluded volume $`v_{el}=f^2/(2n_s)`$ where $`n_s`$ is the salt concentration Joanny ; Andelman . The structure of the compensation is then modified but far from the surface, the theory of neutral polymer layer Manomacromol  can be applied by introducing the total excluded volume $`v+v_{el}`$ (where $`v`$ is the bare one). One find for instance $`\mathrm{\Gamma }_3a^2\varphi _b^{7/8}[v+f^2/(2n_s)]^{1/8}N^{1/2}`$.

As far as we know, very few experimental measurements of the amount of material attached to the surface as a function of $`N`$ have been done under the conditions we are interested (oppositely charged surface, weakly charged chains, no salt added, semi-dilute solutions). The experiments of Chibowski et al. on adsorption of polyacrylic acid (PAA) and polyacrilamide (PAM) to a $`\mathrm{Fe}_2\mathrm{O}_3`$ solid surface, show a significant enhancement of both the layer thickness and the amount of material when the chain molecular weight is increased Chibowski . For instance, the thickness increases from 4.53 nm at molecular weight $`M_w=`$170 kg/mol to 6.69 nm at $`M_w=`$240 kg/mol for PAA at pH 3 and for a polymer concentration of 500 ppm (semi-dilute regime). This increase is in a qualitative agreement with Eq. (3.35) and is related to the formation of larger loops and tails. Of course, the authors mention the important role of the pH, since dissociation equilibriums are monitored by pH and thus both the number of charges of the chains and of the substrate can vary. Such effect has to be included in the theory to do a quantitative comparison.

Concerning the first step of the layer-by-layer process Decher , it should be pointed out that the adsorbed layer is globally neutral ($`d\mathrm{\Psi }/dz(H)=0`$) because counter-ions which penetrate into the layer compensate the decorated surface charge. Hence it is not the apparent charge which is important but the amount of free counter-ions which can be released during the adsorption of oppositely charged polyelectrolytes. This is a key point since, by assuming that adsorbed chains stay upon rinsing, the driving force of adsorption is the exchange of counter-ions related to the amount of potentially free counter-ions. The good observable is thus the amount of adsorbed polymers $`\mathrm{\Gamma }`$, and not the polymer surface excess relative to the bulk concentration. Hence the structure of the adsorbed layer, described in terms of loops and tails, is necessary to compute the coverage. We thus go beyond standard mean-field approaches Joanny ; varoqui ; stoll1 ; stoll2 ; Borukhov ; wang  which cannot distinguish between adsorbed chains and free chains in the distal region. As a result, we find that the issue of overcharging–supercharging relies more on the chain length $`N`$ and bulk volume fraction $`\varphi _b`$ than salt concentration : charge inversion occurs even at low salt concentration which underlines the role played by the chain length as a tunable experimental parameter. Presumably, this feature opens the way to a both simple and efficient way to tune the size of the successive layers in the Decher process. Obviously, a closer look at the case where two layers of oppositely charged polymers interpenetrate and rearrange ladam  is required before we can drive any conclusion on the Decher process castel .
We are grateful to L. Bocquet, E. Trizac, G. Decher and R.R. Netz for interesting discussions.

## Appendix 1 : scaling laws

In this Appendix, we briefly recall the basic scaling concepts for polyelectrolytes. We also provide the essential equations needed to understand the scaling arguments presented in the article. Our present understanding of weakly charged polyelectrolyte solutions without added salt is based on the concept of electrostatic blob PGGPincus ; Barrat : because the different charges of the backbone repel, the electrostatic interactions will deform the otherwise isotropic coil into an elongated string. At a scaling level, the chain is pictured as a linear string of electrostatic blobs (size $`\xi _e`$ , number of monomers $`g_e`$) such that $`\beta \frac{\left(fg_ee\right)^2}{ϵ\xi _e}1`$. Let us assume $`\mathrm{\Theta }`$ solvent conditions for simplicity (the two-body excluded volume parameter $`v=0`$). Then the Flory scaling law for random walks ($`\xi _eag_e^{1/2}`$) yields

$$\xi _e\frac{a}{\left(f^2\mathrm{}_B\right)^{1/3}}\mathrm{and}g_e\frac{1}{\left(f^2\mathrm{}_B\right)^{2/3}}$$
(0.39)
where $`\mathrm{}_B=\frac{\beta e^2}{4\pi ϵa}`$ is the Bjerrum length, in units of $`a`$. The extension of the polymer is

$$L\frac{N}{g_e}\xi _eaN\left(f^2\mathrm{}_B\right)^{1/3}$$
(0.40)
and the local volume fraction (inside the electrostatic blob),

$$\varphi _e\left(f^2\mathrm{}_B\right)^{1/3}$$
(0.41)
When the bulk concentration of monomers is larger than the overlapping concentration, $`\varphi _b^{}NL^3N^2`$, the widely accepted picture is that electrostatic interactions are screened at distances larger than the mesh size of the solution Rubinstein . Accordingly, the solution is a close packing of semi-dilute blobs (size $`\xi _b`$, number of monomers $`g_b`$, $`\varphi _bg_b/\xi _b^3`$), such that the chain is a Gaussian walk of blobs at a scale larger than $`\xi _b`$, and a linear string of electrostatic blobs at a scale smaller than $`\xi _b`$. We have PGGPincus

$$\xi _ba\frac{\varphi _b^{1/2}}{\left(f^2\mathrm{}_B\right)^{1/6}}\mathrm{and}g_b\frac{\varphi _b^{1/2}}{\left(f^2\mathrm{}_B\right)^{1/2}}$$
(0.42)
Eventually, when the bulk volume fraction exceeds $`\varphi _e`$, the electrostatic blob vanishes, and the chain is Gaussian at all scales.

## Appendix 2 : mean-field theory

In this Appendix, we recall the equations of the mean-field theory used to describe polyelectrolytes (see, e.g. Chatellier ; Andelman ). The basic idea is to describe the monomers as microscopic ions carrying a charge $`fe`$, with the following specific features: a) the translational entropy is negligible (this contribution scales as $`(\varphi /N)\mathrm{log}(\varphi /N)`$ and vanishes for long polymers ($`N1`$)), b) an excluded volume interaction estimated in mean-field, c) an elastic contribution to account for the chain deformation. Then the energy (per cm<sup>2</sup>) writes

$`\beta F`$ $`=`$ $`{\displaystyle \frac{1}{2a^2}}\sigma _0\mathrm{\Psi }(0)+{\displaystyle _0^{\mathrm{}}}\left[\mathrm{\Psi }\rho ^+\mathrm{\Psi }\rho ^{}{\displaystyle \frac{a^2}{8\pi \mathrm{}_B}}(\mathrm{\Psi })^2\right]{\displaystyle \frac{dz}{a^3}}`$
$`+`$ $`{\displaystyle _0^{\mathrm{}}}\left[\rho ^+(\mathrm{ln}\rho ^+1)+\rho ^{}(\mathrm{ln}\rho ^{}1)\right]{\displaystyle \frac{dz}{a^3}}+{\displaystyle _0^{\mathrm{}}}\left[{\displaystyle \frac{1}{3}}\varphi ^3+{\displaystyle \frac{a^2}{12}}(\sqrt{\varphi })^2+f\varphi \mathrm{\Psi }\right]{\displaystyle \frac{dz}{a^3}}`$
where the first term is the self-energy of the surface, the second term is the electrostatic energy of resp. the co-ions and the counter-ions and the self-energy of the electric field, the third term is the translational entropy of micro-ions (in volume fraction $`\rho ^\pm `$), and the fourth term accounts for the energy of the monomers (resp.: three body-interaction with amplitude $`w=a^6`$, Edwards elastic contribution, electrostatic energy of the charged monomers). Minimizing the energy yields the Poisson-Boltzmann-Edwards (PBE) system of equations:

$`{\displaystyle \frac{d^2\mathrm{\Psi }}{dz^2}}(z)`$ $`=`$ $`{\displaystyle \frac{4\pi \mathrm{}_B}{a^2}}\left[f\varphi (z)+n_s\mathrm{e}^{\mathrm{\Psi }(z)}(n_s+f\varphi _b)\mathrm{e}^{\mathrm{\Psi }(z)}\right]`$ (0.44)
$`\rho ^+(z)`$ $`=`$ $`n_s\mathrm{e}^{\mathrm{\Psi }(z)}`$ (0.45)
$`\rho ^{}(z)`$ $`=`$ $`(n_s+f\varphi _b)\mathrm{e}^{\mathrm{\Psi }(z)}`$ (0.46)
$`{\displaystyle \frac{a^2}{6}}{\displaystyle \frac{1}{\sqrt{\varphi (z)}}}{\displaystyle \frac{d^2\sqrt{\varphi (z)}}{dz^2}}`$ $`=`$ $`\varphi (z)^2+f\mathrm{\Psi }(z)\mu `$ (0.47)
where $`n_s`$ is the volume fraction of salt in the bulk. Note that the mean-field theory is strictly valid for concentrations above $`\varphi _e`$ for at least one reason: the average volume fraction, $`\varphi `$, is not the local volume fraction experienced by the monomers, $`\varphi _e`$, in the semi-dilute regime, and the three-body interactions (the term $`\varphi ^3`$ in eq. (Appendix 2 : mean-field theory)) is therefore incorrect. However, in the absence of any other theory to describe these systems, it is quite common to extrapolate the mean-field approach outside its strict domain of validity, and try to solve the PBE system of equation in the semi-dilute regime. Because $`\mu `$ is the chemical potential of the monomers in the solution, and this quantity is sensitive to the local concentration, it is tempting to correct for $`\mu 1/g_e`$ ($`\varphi _e^2`$ in theta solvent), while keeping the mean-field estimate of the chemical potential ($`\varphi (z)^2`$ in our context) in the rhs of eq. (0.47). However, this approximation has a technical drawback since $`\mathrm{\Psi }(z)=0`$, $`\varphi (z)=\varphi _b`$ is then not solution of the Edwards equation, as we expect far away from the surface. Here, we will keep the PBE system of equation with $`\mu \varphi _b^2`$. In Section 3.2, we use the Châtellier and Joanny approximation Chatellier  where both the Boltzmann eq. (0.450.46) and Edwards eq. (0.47) are linearized (resp. around $`\mathrm{\Psi }=0`$, and $`\varphi =\varphi _b`$). In terms of the polymer order parameter $`\phi =\sqrt{\varphi }`$, the PBE system of equation simplifies into

$`{\displaystyle \frac{d^2\mathrm{\Psi }}{dz^2}}(z)`$ $`=`$ $`{\displaystyle \frac{8\pi \mathrm{}_Bf}{a^2}}\sqrt{\varphi _b}\delta \phi (z)+\left(\kappa ^2+{\displaystyle \frac{4\pi \mathrm{}_B}{a^2}}\right)\mathrm{\Psi }(z)`$ (0.48)
$`\rho ^+(z)`$ $`=`$ $`n_s[1\mathrm{\Psi }(z)]`$ (0.49)
$`\rho ^{}(z)`$ $`=`$ $`\left(n_s+f\varphi _b\right)[1+\mathrm{\Psi }(z)]`$ (0.50)
$`{\displaystyle \frac{a^2}{6\sqrt{\varphi _b}}}{\displaystyle \frac{d^2\delta \phi }{dz^2}}(z)`$ $`=`$ $`4\varphi _b^{3/2}\delta \phi (z)+f\mathrm{\Psi }(z)`$ (0.51)
where $`\mathrm{\Psi }1`$, $`\phi =\sqrt{\varphi _b}+\delta \phi `$ (with $`\delta \varphi \sqrt{\varphi _b}`$) and $`\kappa ^2=8\pi \mathrm{}_Bn_s/a^2`$.


text:
# Low and high frequency variability as a function of spectral properties in the bright X–ray binary GX 5–1

## 1 Introduction

Low–mass X–ray binaries (LMXBs) are systems where the compact object, either a neutron star or a black hole, accretes matter from a companion with a mass of less than 1 $`M_{}`$. The neutron–star LMXB systems are subdivided on the basis of the pattern they trace out in an X–ray colour–colour or hardness–intensity diagram (CD or HID, respectively) in atoll and Z sources (?). The X–ray flux measurements combined with the knowledge of the distance have shown that Z sources have a high luminosity ($``$LEdd) and accrete at a rate close to the Eddington accretion rate, whereas atoll sources have typical luminosities and inferred accretion rates 5–100 times lower (e.g. see the compilation of source luminosities by ?). The three branches of the Z traced out in a CD or HID by Z sources are called (from top to bottom); Horizontal Branch, Normal Branch, and Flaring Branch.
Studies of the X–ray variability of Z sources revealed two types of quasi–periodic oscillations (QPOs) with frequencies less than 100 Hz (Horizontal Branch oscillations; HBO and Normal Branch oscillations; NBO), twin kHz QPO peaks, and three types of rapid flickering (“noise”), the very low–frequency noise (VLFN), the low–frequency noise (LFN), and the high–frequency noise (HFN) (see ?; ? for reviews). The properties of noise features and the HBOs, such as the central frequency and the fractional rms amplitude, are strongly correlated with the position of a source along the Horizontal Branch. These correlations, together with the observed increase in the ultra–violet flux in the Z source Sco X–1 (?) when the source moves from the Horizontal Branch via the Normal Branch to the Flaring Branch lend support to the idea that the mass accretion rate increases from the Horizontal Branch via the Normal Branch to the Flaring Branch (see ?). However, recent observations (?; ?; ?), together with problems noted before (most notably secular motion of the Z track; ?) show that the situation may be more complex (e.g. van der Klis 2001).
KHz QPOs have now been seen in just over twenty LMXBs, including all Z sources (?; ?). Potentially they can provide a key to measure the basic properties of neutron stars (spin rates and perhaps magnetic field strengths, radii and masses) and thereby constrain the equation of state of ultra–dense matter, and to verify untested general relativistic effects by tracing spacetime just above the neutron star surface (e.g. ?; ?; ?; ?; ?; ?; ?; ?). The discovery of kHz QPOs in GX 5–1 was reported by ?. The HBOs were discovered by ?, and the NBOs in GX 5–1 by ?. The source was also detected at radio (?) and infrared wavelengths (?). In this paper we present an analysis of all RXTE observations of GX 5–1 obtained before 2001. We show for the first time the kHz QPO separation frequency is not constant in GX 5–1. We show that besides the two harmonically related low–frequency QPOs found previously on the Horizontal Branch, two additional harmonically related broad Lorentzian components are needed to obtain a good fit. We discuss the power spectra and colour–colour diagrams of the black hole candidate XTE J1550–564, the Z source GX 5–1, and the atoll source 4U 1608–52, and conclude that there are some remarkable similarities and differences which have not been appreciated in the past due to differences in analysis conventions between in particular neutron stars and black hole candidates.

## 2 Observations and analysis

GX 5–1 was observed 76 times in the period spanning July 27 1996 to March 3 2000 with the proportional counter array (PCA: ?) on board the Rossi X–ray Timing Explorer (RXTE) satellite (?). A log of the observations is presented in Table 1. The total observing time was $``$564 ksec. During $``$40% of the time only a subset of the 5 detectors was active. Two short observations (2 and 4) were omitted from our analysis due to data overflows. Data were always obtained in a mode providing 16 s time–resolution and a high spectral resolution (129 channels covering the effective 2–60 keV range; the Standard 2 mode). In addition, a variety of high time–resolution modes was used in the various observing campaigns; for all observations data with a time resolution of at least $`2^{11}`$ s were obtained in the energy band spanning 2–60 keV.

From the Standard 2 data we computed HIDs. We only used data obtained with the two detectors (proportional counter units 0 and 2) that were always operational since each detector has a slightly different energy response. The background was subtracted using the standard PCA background subtraction model version 2.1e, where we used the version of the skyvle model appropriate for data from each epoch. The hardness, or hard colour, and intensity are defined as the logarithm of the 10.2–17.5 / 6.5–10.2 keV count rate ratio and as the count rate in the 2.5–17.5 keV band of a 16 s average, respectively (RXTE energy channels 25–44 / 14–24 and 3–44, respectively). The high voltage setting of the PCA detectors was changed on March 22 1999 (i.e., the gain changed). Combined with the finite energy resolution this forced us to use somewhat different energy boundaries in our computation of the HID for data obtained after March 22 1999 (i.e., 10.1–17.6 keV / 6.6–10.1 keV for the hard colour and 2.6–17.6 keV for the intensity for observations 42–76. We used RXTE energy channels 20–36 / 11–19 and 2–36, respectively). Besides the gain change the energy response of the detectors changed slowly with time. These various changes have been partially corrected for using the Crab pulsar as a reference, under the assumption that the spectrum of the Crab pulsar does not vary. This correction is only perfect if the spectral shape of GX 5–1 is equal to that of the Crab (the steps involved in this correction were outlined in ? and ?). So, a small systematic error is introduced both by the correction for the drift in detector response and the slightly different energy boundaries which we used as a result of the gain change. Together with the secular variation of the source (?) this led us to compute 5 separate HIDs combining observations 1–7, 8–18, 19–26, 27–41, 42–76 (see Fig. 1 left panel; Table 1).

In Fig. 1 (right panel) the 5 HIDs are overplotted. The source was found on the Horizontal Branch, Normal Branch, and Flaring Branch in four of the five HIDs; for the HID of observations 42–76 the source was not found on the Horizontal Branch. The position of the source along the Z track is characterized by a parameter called Sz representing curve length along the track (?; ?). The Z tracks were parametrized by fitting a spline through manually selected points along the track in each of the 5 HIDs separately. The Sz values are determined using the logarithm of the colours and count rate (?) and normalized by assigning the hard vertex (the Horizontal Branch–Normal Branch junction) the value ‘1’ and setting the soft vertex (the Normal Branch–Flaring Branch junction) to ‘2’. This defines the measure of curve–length, and therefore, S$`{}_{z}{}^{}=0`$ and S$`{}_{z}{}^{}=3`$ are not special points or vertices. We applied this parameterization to the HID obtained from each of the subsets. The track traced out by the source during observations 42–76 only covers part of the Z. In order to calculate the $`S_z`$ values along this track we assumed that the hard vertex in that data set was at the same position relative to the soft vertex as in the HID of observations 27–41.

Using the high time–resolution data we calculated power–density spectra (2–60 keV) of data stretches of 16 s, up to a Nyquist frequency of 2048 Hz, using a Fast Fourier Transformation algorithm. The power spectra were normalized to fractional rms amplitude squared per Hz, added, and averaged according to their position along the Z track (the selection method will be described in detail below). Due to the relatively low amplitude of the variability in the high frequency part of the power spectrum, simultaneous fits of the entire frequency range would not constrain the properties of the high frequency part. Therefore, we fitted the low (1/16–128 Hz) and high ($``$128–2048 Hz) frequency part of the power spectra separately. The low–frequency part of the average power spectra was fit with a function consisting of an exponentially cut–off power law to describe the noise at low frequencies (LFN) and at most five Lorentzians to fit the QPOs. The component arising in the power spectrum due to the Poisson counting noise was subtracted. For values of Sz $`>`$1.0 a power law was added to the fit function to fit the very low–frequency noise (VLFN). The function used to describe the high–frequency part of the average power spectrum was built up out of two Lorentzians describing the kHz QPOs. Contrary to when we fit the low–frequency part, the continuum due to Poisson noise was not subtracted prior to fitting, therefore we included a constant to account for this Poisson noise. Sometimes also a power law component was added to account for the slope of the underlying continuum, caused by the high frequency tail of one of the low–frequency components. The high frequency tail did not systematically influence our fit–parameters.

From previous studies of GX 5–1 and other Z sources it is known that the frequency of the QPOs on the Horizontal Branch (the HBOs and the kHz QPOs) increases with Sz (?; ?; ?; ?). For all but the HID of observations 42–76 (when the source was not found on the Horizontal Branch or the upper part of the Normal Branch and therefore no HBOs were present in the power spectra; see Fig. 1 left panel) we studied the relation between Sz and HBO frequency. These relations, plotted in Fig. 2, have two notable properties. First, the relation found for observations 27–41 is clearly offset from the other relations. This is an effect of secular motion of the source during the gap in between the observations 27–35 and 36–41. During the first 9 observations (27–35) the source was found on the Horizontal Branch, whilst, after a 9 day gap, during observations 36–41 the source traced a part of the Horizontal Branch close to the Normal Branch, the Normal Branch, and the Flaring Branch. We checked the vertices of the 5 HIDs to search for additional evidence of secular motion but in view of the uncertainties involved in their manual selection, we conclude that the vertices of the 5 HIDs were not significantly different. The second notable property involves the jump in HBO frequency near S$`{}_{z}{}^{}1`$ (see also plots in ?). To investigate the nature of the jump further we plotted the hard colour vs. the frequency of the HBO. A discontinuity in the rate of change in frequency as a function of hard colour was also observed in the hard colour vs. frequency plot at a frequency of $``$50 Hz. Therefore, the jump in frequency when plotted vs. Sz is not an artefact of the Sz parameterization. We investigated whether a similar jump was present in the frequencies of the kHz QPOs, but we lacked the signal–to–noise to conclude on this.

In order to combine all power spectra with similar HBO frequencies, we shifted all Sz–HBO frequency relations to one “parent” relation. To this end, we fitted a polynomial to the relation between Sz and HBO frequency for each of the HIDs separately. We included only measurements for S$`{}_{z}{}^{}<1.0`$ in the fit. The order of the polynomial was determined such that the reduced $`\chi ^2`$ of the fit was $``$1. The parent relation for S$`{}_{z}{}^{}<1.0`$ is the best linear fit to the Sz–HBO frequency relation for all the HIDs combined ($`\nu _{HBO}=9.9+36.8\times S_{z,parent}`$; the drawn line in Fig. 2). The shifting procedure works as follows: given the Sz value (Sz,initial) obtained from the HID corresponding to the observation the power spectrum is calculated from, the HBO will have a certain predicted frequency which follows from the Sz–HBO frequency relation found for that HID. That predicted HBO frequency corresponds, given the Sz–HBO frequency parent relation, to a Sz,parent. This Sz,parent is then assigned to the power spectrum instead of Sz,initial (the change of Sz corresponds to a horizontal shift in Fig. 2). This shifting procedure was done for each 16 s power spectrum separately. Finally, all power spectra were selected according to their Sz,parent values for S$`{}_{z}{}^{}<1.0`$; the selection bin width was 0.1. For S$`{}_{z}{}^{}>1.0`$ the unshifted values obtained from each separate HID were used in the selection, since both the secular motion and the effects of the changes in the response move the Z track in a direction nearly perpendicular to the Normal Branch. Here the selected bin width was 0.05. The selected power spectra were averaged and fitted with the function described above.

The errors on the fit parameters were determined using $`\mathrm{\Delta }\chi ^2=1.0`$ (1$`\sigma `$, single parameter). The error on Sz is the standard deviation of the distribution of Sz values in one selection bin. In cases where components were not significantly detected, 95% confidence upper limits were determined using $`\mathrm{\Delta }\chi ^2=2.71`$. The full–width at half maximum (FWHM) of the Lorentzians was fixed at 10 Hz in case an upper limit was determined to the sub–HBO, its third harmonic or the second harmonic of the HBO, and at 30 Hz in case of the first harmonic (which is the fundamental) of the HBO. The frequency of the component for which an upper limit was determined was not fixed but restricted to a range of values around that expected on the basis of the observed trends. Upper limits to the presence of kHz QPOs were determined using a FWHM of 75 Hz.

## 3 Results

### 3.1 Spectral states

The source was found on the Horizontal Branch, Normal Branch and Flaring Branch. Note the clear “Dipping Flaring Branch” trailing the Flaring Branch (see Fig. 1, see also ?). GX 5–1 was found to reside most often on the Normal Branch during the time of our observations (see Fig. 3, upper panel). The average velocity along the Z track (defined as the average of $`\mathrm{V}_\mathrm{z}(\mathrm{i})=\frac{\mathrm{S}_\mathrm{z}(\mathrm{i}+1)\mathrm{S}_\mathrm{z}(\mathrm{i}1)}{\mathrm{T}(\mathrm{i}+1)\mathrm{T}(\mathrm{i}1)}`$; ?; Fig. 3, lower panel) was approximately constant for $`0.1<`$S$`{}_{z}{}^{}<1`$ but gradually increased for S$`{}_{z}{}^{}>1`$, i.e., in anti–correlation with the time the source spent in each part of the Z diagram, although the product of the two is not exactly constant. Note that the average velocity along the Z track also increased towards the lowest Sz values.

### 3.2 Low–frequency power spectra

We found two new components in the averaged power spectrum of GX 5–1. One of these components was located at frequencies consistent with half the frequency of the HBO, the so–called sub–HBO. This component was reported before in the Z sources Sco X–1, GX 340+0, and GX 17+2 with RXTE (?; ?; ?). The second new component has a frequency which is consistent with three times the frequency of the sub–HBO (or 1.5 times the frequency of the HBO; see Fig. 4). Fitting the average low–frequency (1/16–256 Hz) power spectrum of S$`{}_{z}{}^{}=0.60\pm 0.03`$ with a function consisting of a cut–off power law (LFN), and three Lorentzians (sub–HBO, HBO, and 2nd harmonic of the HBO) a $`\chi ^2`$ of 383 for 209 degrees of freedom was obtained. Adding a fourth Lorentzian component at frequencies $``$3 times the frequency of the sub–HBO gave a $`\chi ^2`$ of 314 for 206 degrees of freedom. An F–test (?) to the $`\chi ^2`$ of the fits with and without the fourth Lorentzian revealed that the probability that the reduction in $`\chi ^2`$ could be achieved by a random process is 6.5$`\times 10^8`$, i.e., the significance of the addition of the fourth Lorentzian component is $`6\sigma `$. For the other power spectra where this component was found a similar value was obtained. Although this component is too broad to qualify formally as a QPO (since the Q–value is less than 2; see Table 2 and Fig. 4, right panel), for now we refer to this component as the third harmonic to the sub–HBO component (see Fig. 4). This classification is rendered strong support by the fact that the Q–values of both the sub–HBO and its third harmonic are consistent with being the same (Fig. 4; right panel).
Similar to our work on the Z–source GX 340+0 (Jonker et al. 2000b) we also experimented with several other fit functions used by various authors to describe the power spectra of other LMXBs. Both using a fit function built up out of a broken power law, to fit the LFN component, and several Lorentzians to fit the QPOs, and a fit function built up out of just Lorentzians resulted in significantly higher $`\chi _{red}^2`$ values than when the fit function described in above was used ($`\chi _{red}^2`$= 2.67 for 207 degrees of freedom (d.o.f.), and $`\chi _{red}^2`$= 1.93 for 205 d.o.f., respectively), while adding the fourth Lorentzian component to these fit functions still improved the $`\chi _{red}^2`$ significantly (to $`\chi _{red}^2`$= 2.48 for 204 d.o.f., and $`\chi _{red}^2`$= 1.73 for 202 d.o.f., respectively).
Besides the two new components, the HBO, its second harmonic and the LFN were also detected (see Table 2, 3). A typical fit (S$`{}_{z}{}^{}=0.60\pm 0.03`$) showing the contributions of the individual components is presented in Fig. 5. In four of the Sz selections the HBO was fit with two Lorentzian peaks in order to obtain a good fit (see Table 2). This is either related to the fact that the HBO moves in frequency within the selection, or the HBO profile itself is asymmetric. Whenever two Lorentzian peaks were used to describe the HBO we obtained the rms and FWHM weighted mean of the two peaks; we used the weighted mean parameters in Fig. 4 and in Table 2. We weighted the frequencies according to one over the square of the FWHM and proportionally to the square of the amplitude. The FWHMs were weighted proportionally to the square of the amplitude. The frequency offsets between the two peaks were ignored in the weighing since one of the two Lorentzian peaks contained several times more power than the other. The powers of the two Lorentzians were added. We would like to remark that the jump in HBO frequency vs. Sz for each HID (Fig. 2) has been reduced due to the shifting to Sz,parent.

The fractional rms amplitude (integrated between frequencies ranging from 0–$`\mathrm{}`$) of all the low–frequency Lorentzian components (excluding the NBO) decreased as a function of Sz (see Fig. 6 upper panels; Table 2). For S$`{}_{z}{}^{}=1.01.5`$ the frequency and the fractional rms amplitude of the HBO was consistent with being constant at $``$50 Hz and $``$2%, respectively. For S$`{}_{z}{}^{}>1.5`$ the HBO was not detected with an upper limit of 1.2%. The FWHM of the Lorentzian components increased as a function of Sz (see Fig. 6 lower panels; Table 2).

The fractional rms amplitude of the LFN (0.1–100 Hz) decreased from 7.4%$`\pm `$0.3% to 3.38%$`\pm `$0.02% as Sz increased from 0.12$`\pm `$0.03 to 1.01$`\pm `$0.02 (Table 3). The LFN component could not be measured for S$`{}_{z}{}^{}>1`$. The LFN power law index increased from 0.07$`\pm `$0.05 to 0.41$`\pm `$0.01 while the cut–off frequency increased from 2.9$`\pm `$0.5 Hz to $``$60 Hz. For S$`{}_{z}{}^{}>1`$ we added a power law to the fit function to represent the power at frequencies below 1 Hz; the very low–frequency noise (VLFN). The VLFN fractional rms amplitude (0.001–1 Hz) increased from less than 1% at the hard vertex to more than 20% at S$`{}_{z}{}^{}>3.0`$. This increase was gradual for S$`{}_{z}{}^{}<2.0`$ but steep for $`S_z>2.0`$ (see Fig. 7, Table 3). The VLFN power law index gradually increased to 1.5 and remained constant for S$`{}_{z}{}^{}<2.0`$, for S$`{}_{z}{}^{}>2.0`$ the power law index was consistent with 2.

The frequency of the NBO decreased from 6.1$`\pm `$0.3 Hz at S$`{}_{z}{}^{}=1.07\pm 0.02`$ to 5.25$`\pm `$0.05 Hz at S$`{}_{z}{}^{}=1.33\pm 0.02`$ before it increased again to 6.1$`\pm `$0.1 Hz at S$`{}_{z}{}^{}=1.82\pm 0.03`$. For S$`{}_{z}{}^{}>1.9`$ its frequency could not be determined with high accuracy since the FWHM increased to 8 Hz. The fractional rms amplitude of the NBO varied between 1.69%$`{}_{0.09}{}^{}{}_{}{}^{+0.05}`$% and 2.38%$`\pm `$0.03% for S$`{}_{z}{}^{}=1.07\pm 0.02`$ and S$`{}_{z}{}^{}=1.52\pm 0.02`$, respectively but it decreased to $``$1.14%$`\pm `$0.05% at S$`{}_{z}{}^{}=1.95\pm 0.05`$. The somewhat erratic behaviour in the fractional rms amplitude of the NBO between S$`{}_{z}{}^{}=1.31.5`$ can be explained by the fact that the HBO was very broad and had power in the same frequency range as the NBO. The fractional rms amplitude is anti–correlated with the frequency and FWHM of the NBO (see Fig. 8). At S$`{}_{z}{}^{}<1.0`$ stringent upper limits on the presence of the NBO were derived.

At S$`{}_{z}{}^{}>2.0`$ the NBO was not detected but instead a cut–off power law was fitted to represent the power at frequencies comparable to the NBO frequencies. This cut–off power law was peaked (power law index of $`0.9\pm 0.2`$) near the soft vertex, with power law index $`0.9\pm 0.2`$ and evolved into a flat–topped noise component at S$`{}_{z}{}^{}=2.15\pm 0.05`$ with a typical power law index of 0.3. The cut–off frequency increased from 4.3$`\pm `$0.6 Hz to 25$`{}_{5}{}^{}{}_{}{}^{+8}`$ Hz. The fractional rms amplitude (0.001-100 Hz) increased from 1.68%$`\pm `$0.07% at S$`{}_{z}{}^{}=2.05\pm 0.05`$ to 2.07%$`\pm `$0.05% at S$`{}_{z}{}^{}=2.15\pm 0.05`$. For S$`{}_{z}{}^{}>2.15`$ the fractional rms amplitude was consistent with 1.6%. The noise component became undetectable at an Sz value of 2.64$`\pm `$0.06 with an upper limit of 1.2%, where the vertex between the Flaring Branch and the “Dipping Flaring Branch” is at S$`{}_{z}{}^{}2.5`$. The fit parameters of the cut–off power law are marked with a b in Table 3.

### 3.3 KHz QPOs

The lower and upper kHz QPO were detected at frequencies ranging from 156$`\pm `$23 Hz to 627$`{}_{13}{}^{}{}_{}{}^{+25}`$ Hz and from 478$`\pm `$15 Hz to 866$`\pm `$23 Hz, respectively (dots and crosses, respectively top panel, Fig. 9 A). The peak separation ($`\mathrm{\Delta }\nu `$) was not constant; the fit of a constant to the peak separation (325$`\pm `$10 Hz) vs. Sz resulted in a $`\chi _{red}^2`$ of 7 for 10 degrees of freedom; an unacceptable fit. We tried various other functions (e.g. see Fig. 9 bottom panel and Table 5). From an F–test it was clear that the use of a broken function or a parabola significantly reduced the $`\chi ^2`$ ($`3.5\sigma `$, $`3\sigma `$, respectively) with respect to a constant, indicating that $`\mathrm{\Delta }\nu `$ is decreasing towards higher Sz at a rate that is not constant but increases. The data indicate a decrease in FWHM of the lower kHz QPO towards larger Sz, but the scatter is large. The FWHM of the upper kHz QPO varies between 100 Hz and 247 Hz (see crosses Fig. 9 B top panel). The fractional rms amplitude of the lower kHz QPO decreased gradually from 3.6%$`\pm `$0.4% to 0.5%$`\pm `$0.1%, that of the upper kHz QPO first increased slightly from 2.0%$`\pm `$0.4% to 2.6%$`\pm `$0.1% to decrease gradually to 0.9%$`\pm `$ 0.1% (see dots and crosses respectively, Fig. 9 C top panel). The Sz changed from $`0.11.1`$ over the plotted interval. All fit results are given in Table 4.
Motivated by the results on GX 17+2 (?) we plotted the frequency of the kHz QPOs and their separation frequency as a function of $`\nu _{HBO}`$ (Fig. 10). Near $`\nu _{HBO}45`$ Hz the peak separation starts to decrease most likely due to an increase in the rate at which the frequency of the lower kHz QPO increases with $`\nu _{HBO}`$. This $`\nu _{HBO}`$ is consistent with the jump at $`\nu _{HBO}44`$$`47`$ Hz in a $`\nu _{HBO}`$ vs. Sz plot (see Section 2).

## 4 Discussion

A detailed analysis of all RXTE observations obtained to date of GX 5–1 was presented. The main result of our high–frequency variability study is that the kHz QPO peak separation is not constant in GX 5–1. The low–frequency power spectra for S$`{}_{z}{}^{}<1.0`$ are complex and can be best described by 4 harmonically related Lorentzians, in addition to the low–frequency noise (LFN). Two of these Lorentzians were not known previously in GX 5–1. At S$`{}_{z}{}^{}>1.9`$ a new component was found in the power spectra. Below we discuss these findings and compare the timing properties and colour–colour diagrams of GX 5–1, the black hole candidate XTE J1550–564, and the atoll source 4U 1608–52.

### 4.1 KHz QPOs

We showed that the kHz QPO peak separation is not constant. The data is consistent with a constant peak separation followed by a steep decrease for S$`{}_{z}{}^{}1.0`$, or a parabolic relation between the peak separation and Sz. A decrease in peak separation was previously found in Sco X–1 (?), 4U 1608–52 (?), 4U 1728–34 (?), 4U 1735–44 (?), 4U 1702–43 (?), and GX 17+2 (?). We found that in GX 5–1 the rate of increase in frequency of the lower kHz QPO changes near $`\nu _{HBO}45`$Hz, causing the peak separation to decrease.
The frequencies of the kHz QPOs of GX 5–1 are low in comparison with the other kHz QPO Z sources (the lowest frequencies of the lower and upper kHz QPO we found are $`156\pm 23`$ Hz, and $`478\pm 15`$ Hz, respectively). The lowest upper kHz QPO frequency found so far in an atoll source (449$`\pm `$20 Hz for 4U 0614+09; ?) is at the low–frequency limit ($``$500 Hz for the upper kHz QPO; ?) imposed by the sonic point beat frequency model. The cause of the lower bound on the Keplerian frequency in the sonic point beat frequency model is that the radiation drag of the luminosity produced near the surface of the neutron star can only remove the required amount of the specific angular momentum of the gas to make it fall to the neutron star when the gas is close to the neutron star (?). Unlike 4U 0614+09, in the case of GX 5–1 QPOs with still lower frequency did not occur because the source was not found at lower Sz values but not because the QPO disappeared. Therefore, observations of kHz QPOs in GX 5–1 at even lower Sz values will provide a strong test of the sonic point beat frequency model.

In the so–called two–oscillator model (?; ?) the angle $`\delta `$ is defined ($`\delta =arcsin[(\nu _{upper}^2\nu _{lower}^2)^{0.5}(\nu _{HBO}\nu _{upper}/\nu _{lower})]`$, where $`\nu _{lower}`$ is the frequency of the lower kHz QPO, $`\nu _{upper}`$ is the frequency of the upper kHz QPO, and $`\nu _{HBO}`$ is the frequency of the HBO; $`\delta `$ is the angle between the magnetosphere equator and the disc plane in the model). The current version of the theory predicts $`\delta `$ to be constant. However, the measurements of ? for GX 17+2 and the our measurements for GX 5–1 show that $`\delta `$ vs. the frequency of the lower kHz QPO changes is not constant (see Fig. 11). A fit of a constant gives $`\delta =6.1\pm 0.2`$ and a $`\chi ^2`$ of 377 for 9 d.o.f. for GX 17+2 and $`\delta =6.3\pm 0.1`$ and a $`\chi ^2`$ of 95.2 for 9 d.o.f. for GX 5–1.

### 4.2 Comparison with other Z sources

We have detected a sub–HBO component similar to the component in the Z–sources Sco X–1 (?), GX 340+0 (?), and GX 17+2 (?). Furthermore, we found a fourth Lorentzian component whose frequency was consistent with 1.5 times the frequency of the HBO. If we interpret the sub–HBO as the fundamental frequency, then the HBO would be the second harmonic, the newly found peak at 1.5 times the HBO frequency the third harmonic, and the second harmonic of the HBO the fourth harmonic. In such a scheme the odd harmonics should either be formed less coherently than the even harmonics, or an additional broadening mechanism has to be invoked which broadens the odd harmonics more than the even harmonics (see Fig. 4; right panel). This harmonic structure could be explained in a scenario where the HBO is caused by a warped accretion disc with a two–fold symmetry. The warp itself is stable to deviations from its two–fold symmetry, i.e., variations in the warp take place on a timescale long compared to the timescale of variations in the difference between the two sides of the warp, leading to broader odd harmonics. This connects to the models proposed by ? and ?. In the magnetic beat–frequency model for the HBO (?; ?) the sub–HBO should reflect differences between the two magnetic poles. Differences in the magnetic field configuration are not likely to vary in time rapidly, so the extra broadening of the sub–HBO with respect to the HBO is unexplained. This necessitates the introduction of an additional broadening mechanism working on the odd harmonics only after their formation.
Prior to the detection of the sub–HBO and its third harmonic, the frequency of the second harmonic of the HBO was not measured to be exactly twice the HBO frequency (?); with the introduction of the third harmonic of the sub–HBO in the fit function this discrepancy has disappeared (see Fig. 4; left panel). Power density spectra calculated from data obtained with the EXOSAT and Ginga satellites (?; ?) were fit using only the even harmonics (the HBO and its 2nd harmonic) and an additional high frequency noise (HFN) component. Now, using the larger collecting area of RXTE we find that the HFN in GX 5–1 is better described by two Lorentzian components which have centroid frequencies of 0.5 and 1.5 times that of the HBO.
The same function we used for our analysis of GX 5–1 was used in the analysis of the similar Z source GX 340+0 (?). There no third harmonic of the sub–HBO was found. Instead we found excess power (a shoulder) close to the HBO peak, which when fit with a Lorentzian was at a frequency entirely inconsistent with 1.5 times that of the HBO. Reanalysis of those data with the insights gained by our work on GX 5–1 did not alter this conclusion. In GX 340+0 both the HBO and the NBO peak became asymmetric when they were strongest. Similar behaviour was found in GX 5–1, where for S$`{}_{z}{}^{}<0.5`$ four of the five HBO measurements were done using two Lorentzians with centroid frequencies up to 2 Hz apart with similar FWHMs but different strengths, indicating that the HBO was asymmetric. However, unlike in our analysis of GX 340+0 where the data were selected according to HBO frequency, in the case of GX 5–1 this asymmetry could be an artifact of the selection method as explained in Section 3.2. An attempt to fit the average power spectra of GX 17+2 (?) with the same fit function as that of GX 5–1 did not lead to the detection of a third harmonic to the sub–HBO. The asymmetry of the HBO found in GX 340+0 and the third harmonic of the sub–HBO found in GX 5–1 could be two different components, one of which is detected in each source because of the special circumstances affecting detections in each case: the larger amplitude of the HBO in GX 340+0 than in GX 5–1 allows to detect a strong asymmetric (shoulder) component, whereas the higher signal–to–noise in GX 5–1 than in GX 340+0 (due to the fact that GX 5–1 is several times brighter than GX 340+0) leads to the detection of the broad third harmonic. New observations with higher signal–to–noise ratios can reveal whether both components are present in both sources or whether they represent power provided by one single component.
The low–frequency power spectra of the three Z sources discussed above (GX 17+2, GX 340+0, GX 5–1) are quantitatively different. In Sco X–1 and GX 17+2 a QPO is found when the source is on the Flaring Branch (the FBO; ?; ?). Such a QPO is absent in GX 340+0, GX 5–1, GX 349+2, and Cyg X–2. The new component we found in the power spectrum of GX 5–1 with frequencies similar to those of the FBO when the source is on the Flaring Branch is much broader than the FBO in Sco X–1 and GX 17+2. The second harmonic of the HBO is relatively strong in GX 17+2 while for GX 5–1 and GX 340+0 only weak harmonics have been detected (see ?; ?). Furthermore, the fractional rms amplitude of the HBO is highest in GX 340+0, and lowest for GX 17+2. The HBO frequency in GX 17+2 was found to decrease on the Normal Branch (?; ?), whereas in GX 340+0 and GX 5–1 the frequency is consistent with being constant (?). Besides the Lorentzians also the LFN properties are different; it is always peaked in GX 17+2 whereas this is never the case in GX 5–1 and GX 340+0 (see ?; ?; ?; ?). The dependence of the HBO, lower and upper kHz QPO Q–values on Sz is shown in Fig. 12. The HBO Q–value decreases gradually for S$`{}_{z}{}^{}>0.5`$, while that of the kHz QPOs increases. This is opposite to what was found for GX 17+2 (?), where the Q–value of the HBO and that of the kHz QPOs increased as a function of Sz.

### 4.3 Comparing GX 5–1 with XTE J1550–564 and 4U 1608–52

In the power spectra of black hole candidates low–frequency QPOs with harmonics are occasionally found in the low, intermediate and very high state (e.g. GS 1124–68, ?; GX 339–4, ?; XTE J1550–564, ?; ?). Occasionally, these black hole candidate QPOs have asymmetric profiles. We found a similar harmonic structure of QPOs in the Z source GX 5–1. The Q–values of the harmonics for XTE J1550–564 are similar to those found for GX 5–1: low Q–values for the odd harmonics and high Q–values for the even harmonics (?). However, we checked the energy dependence of the QPOs in GX 5–1 and there are no indications of changes in the Q–values with energy as were found in XTE J1550–564 (?). In general the energy dependence and time lag behaviour of the QPOs in GX 5–1 (?) is different from what was found in XTE J1550–564 (?). On the Flaring Branch, similar to the High State in black hole candidates, the fractional rms amplitude of the variability is low, and QPOs are weak or absent (in GX 17+2 and Sco X–1, just as in XTE J1550–564 a QPO with a centroid frequency of $``$18 Hz is found, but whether they are related is unclear). The similarities between the Horizontal Branch variability and that of the Intermediate/Low State have been pointed out before (?; ?; ?). Furthermore, the trend of a decrease in variability timescale as the source approaches the Flaring Branch/High State is the same in black hole candidates as in GX 5–1 and 4U 1608–52 (cf. ?; ?).
To further investigate the similarities between a black hole candidate (XTE J1550–564), a Z source (GX 5–1), and an atoll source (4U 1608–52) we created colour–colour diagrams (CDs) for the three sources in two ways. First, we used the colour definitions commonly used for atoll and Z sources (soft and hard colour were defined as the logarithm of the 3.6–6.2 / 2.5–3.6 keV and 9.8–16.0 / 6.2–9.8 keV count rate ratio, respectively). These CDs are shown in the left panels of Fig. 13, 14, and 15. Second, we created CDs in a similar way as ?, an approach often used for black hole candidates (?; ?; although see also ?). We defined the hard and soft colour as the 9.7–16.0/ 2–6.4 keV and 6.4–9.7/ 2–6.4 keV count rate ratio for observations 1–41, respectively and that of the 9.7–15.8 / 2–6.6 keV and 6.6–9.7 / 2–6.6 keV count rate ratio for observations 42–76, respectively (Fig. 13, right panel). In case of the right panels for Fig. 14 and 15 the hard and soft colour were defined as the 16.0–19.4 / 2.2–6.2 keV, and 6.5–15.7 / 2.2–6.2 keV count rate ratio, respectively. Note that in the right panels the soft colour is plotted vs. the hard colour, as the convention, introduced by (?) to plot hard colour vs. soft colour was abandoned in some of the recent black hole work. In the right panels the points have been connected by a line. Lines bridging different parts of the diagrams do not reflect sudden jumps or changes in source state, but are due to gaps introduced by observational windowing.
A similar structure as for XTE J1550–564 was found in the CD of GX 5–1 (Figure 13, 14 right panel). The Flaring Branch is now vertical and located at the left side of the figure, similar to the High State in XTE J1550–564. The Normal Branch connects to the Flaring Branch and points towards the upper right corner of the diagram (both the hard colour and the soft colour increase as the source moves up the Normal Branch). The Horizontal Branch is a continuation of the Normal Branch, with only a slight bend at the Normal Branch–Horizontal Branch junction. (It is interesting to note that there is also no vertex present in GX 17+2 at the Horizontal– Normal Branch junction at energies above 14.8 keV, Figure 3C of ?). The ’Dipping’ Flaring Branch is parallel to the Normal Branch but directed towards the lower left corner. The ’Dipping’ Flaring Branch has been found in the Z sources Cyg X–2, GX 5–1 and GX 340+0 (?; ?; ?). The CD of the atoll source 4U 1608–52 again shows a similar structure as found for XTE J1550–564 (Fig. 14, 15 right panels). The island state is found in the top right corner of the right panel of Fig. 15, whereas the lower banana branch connects to the upper banana branch in the lower left corner. The upper banana branch seems similar to the black hole candidate High State and the Z source Flaring Branch.
The CD of the black hole candidate XTE J1550–564, plotted using colours typically applied for neutron star sources is reminiscent of that of an atoll source, although more structure is observed (Fig. 14 left panel). The source was first found in the upper right part of the CD as the outburst started. The path traced by the source as the outburst progressed is indicated with arrows. For a complete description of the outburst of XTE J1550–564 we refer to ?, ?, and ?.
The similarities between the colour–colour diagrams of the Z source GX 5–1, the atoll source 4U 1608–52, and the black hole candidate XTE J1550–564 strengthen the idea that some of the spectral properties of LMXBs originate in the accretion disc and do not depend on the presence of a solid surface. However, our analysis also indicates that differences between spectral properties of the neutron star Z and atoll sources and the black hole candidates do exist (e.g., the left panel of Fig. 13, 14, and 15) but may not clearly show up in representations of the CD commonly used for black–hole candidates. The explanation for this could be that in the “black–hole colour scheme” only three independent colours are defined whereas in the “neutron star colour scheme” four colours are usually defined.

So, although the CDs of a Z or atoll source and that of a black hole candidate appear similar in shape when plotted using colours typical for black hole studies, their appearance is both dissimilar and more complicated when plotted using colours typical for neutron star studies. Whether the complications can be ascribed to the unique nature of the sources considered can only be checked by comparing more sources but the conclusion that dissimilar CDs can appear more similar than they really are in the “black hole” representation can already be drawn. Conclusions regarding the spectral state or properties of the source based on the colours commonly used in black hole candidates should be regarded with caution since an other choice of colours may lead to a qualitatively different CD.

## Acknowledgments

This work was supported in part by the Netherlands Organization for Scientific Research (NWO). This research has made use of data obtained through the High Energy Astrophysics Science Archive Research Center Online Service, provided by the NASA/Goddard Space Flight Center. This work was supported by NWO Spinoza grant 08-0 to E.P.J.van den Heuvel. P.G.J. would like to thank Rob Fender, Kieran O’Brien, and Marc Klein Wolt for various discussions. RW was supported by NASA through Chandra Postdoctoral Fellowship grant under PF9–10010 awarded by CXC, which is operated by SAO for NASA under contract NAS8–39073.