format: text

text:
# Effects of Quasi-Static Aberrations in Faint Companion Searches

## 1 Introduction

The next generation high order adaptive optics (NGAO) systems that are currently in development will produce high Strehl ratio images in the infrared on 10 m class telescopes. If those NGAO systems are to be used to search for very faint companions (brown dwarfs and exoplanets), it is of great importance to understand the limits for high contrast imaging on current adaptive optics (AO) systems and see how they would affect NGAO systems performances. Static aberrations have long been known to be the limiting factor for space-based telescopes when attempting high contrast imaging for faint companion detection (\[Brown & Burrows 1990\]). If the telescope physical conditions change in time (like point spread function (PSF) breathing for HST), the structure becomes quasi-static and produces PSF structures that evolve slowly in time, making difficult a good PSF calibration when using reference stars. Long AO corrected exposures have shown that ground-based imaging is not free from this problem. When using an AO system that delivers partially diffracted PSFs, the final wavefront is an interference between the AO wavefront residuals and the optical path wavefront aberrations to the detector. If enough random aberrations are removed by the AO system, a coherent PSF starts to appear. This is also true for the quasi-static portion of the wavefront. As more turbulence is removed, the quasi-static wavefront errors not corrected by the AO system will introduce coherent structures that will not average out in time and will mask faint companions. A partially-corrected long exposure AO observation thus suffers from the same problem that limits space telescopes. The current limit for detecting faint companions with existing AO systems is not the atmospheric turbulence correction efficiency but rather a quasi-static PSF calibration problem.

In the past few years, our group has developed a specialized infrared camera, TRIDENT (\[Marois et al. 2000a, Marois et al. 2002\]), to overcome the PSF calibration problem. The main feature of TRIDENT is to acquire three simultaneous images in three distinct narrow spectral bands. It is possible to enhance the star/companion contrast after image combinations by selecting special spectral features that are typical of the companion and not of the star (\[Smith 1987, Rosenthal et al. 1996, Racine et al. 1999, Marois et al. 2000b\]). In TRIDENT, the three wavelengths (1.567 $`\mu `$m, 1.625 $`\mu `$m and 1.680 $`\mu `$m, 1% bandwidth) have been selected across the 1.6 $`\mu `$m methane absorption bandhead that is only present in the spectrum of cold ($`T_{\mathrm{eff}}`$ $`<`$ 1470 K, \[Fegley and Lodders 1996\]) substellar objects. The 1.567 $`\mu `$m image ($`I_{\lambda _1}`$) shows the star and companion while the 1.625 and 1.680 $`\mu `$m images ($`I_{\lambda _2}`$ and $`I_{\lambda _3}`$) show the star with a fainter companion due to the methane absorption bandhead. The 1.625 and 1.680 $`\mu `$m images can thus be used as reference PSFs to correct the first and second order PSF structure evolution with wavelength. Since the three wavelengths are simultaneous, each observed wavelength sees the same PSF atmospheric distortion, a good PSF correlation between the three wavelengths can thus be achieved. This concept would remove atmospheric speckles and problems associated with the PSF evolution with telescope pointing, thermal changes and atmospheric r0 variations. PSF subtraction is achieved by calculating simple and double image differences (SD and DD) obtained with the following algorithms (\[Marois et al. 2000b\]):

$`sd_{12}=I_{\lambda _1}I_{\lambda _2}`$ and $`dd=sd_{12}`$k$`sd_{13}`$

The DD can bring the atmospheric speckle noise below the level of the photon noise and remove the PSF quasi-static structure to its second order evolution with wavelength.

## 2 Data acquisition and reduction

Data were obtained on 2001 July 8-12 and on 2001 November 21-24, at the f/20 focus of the 3.6 m CFHT AO bonnette PUEO (\[Rigaut et al. 1998\]) with TRIDENT. Flat fields and darks were obtained during the day. In July, data well inside the linear regime of the detector were acquired to test the PSF stability. Total integration time was typically 1 h per target. Seeing conditions for this run were good (Strehl of 0.5 in the $`H`$ band). In November, reference stars were acquired and instrument rotations ($`+/`$ 90 degrees in steps of 2 to 5 degrees) were done to calibrate or smooth common and non-common path aberrations discovered to be the limiting factor in the July mission data. Saturated and non-saturated images were acquired to minimize readout noise. Total integration time was 1h30 with medium seeing conditions (Strehl from 0.15 to 0.4 in the $`H`$ band). In total, 35 stars have been observed during these two missions with spectral type ranging from B to M. Some stars were observed on different nights to increase their signal to noise ratios and study PSF stability. More details on data reduction and analysis are given in Marois et al. (2003a) and Marois et al. (2003b).

## 3 Simulations

A long AO-corrected exposure is a superposition of a coherent PSF with a residual atmospheric halo. The coherent PSF does not look like a perfect diffraction pattern, but is distorted with some structures. Figure 1 shows the observed PUEO-TRIDENT PSF compared to a simulated PSF of 100 independent atmospheric speckle patterns (approximately 1.5 s exposure).

The observed Ups And AO corrected PSF shows a series of broken rings with a diffuse atmospheric halo while the simulated PSF shows a series of perfect diffraction rings with the same diffuse halo. Structures seen in long exposure AO corrected PSFs are not due to turbulence but to a wide range of low to high frequency telescope and instrument aberrations. Increasing the integration time beyond a few seconds increases the signal to noise ratio of the coherent PSF structure. The total estimated quasi-static aberrations residual is 110 nm RMS or a $`H`$ band Strehl ratio of 0.82. This result is in good agreement with the measured value of the PUEO AO system (approximately 0.84 in the $`H`$ band, \[Rigaut et al. 1998\]), comfirming that most of the aberrations are coming from the PUEO AO system. The estimated power spectrum has the following shape: a decreasing spectrum for the first 12 Zernike terms that correspond to 65% of the total aberrations, a constant power spectrum for 150 Zernike terms that have 20% of the total aberrations, a mid-frequency component that have 5% of the total aberrations and a high frequency component with 10% of the total aberrations. Figure 2 shows the same simulation as figure 1 (100 coadditions of independent atmospheric speckle patterns) but including, for each pattern, a static phase error that interferes with the atmospheric turbulence. The simulated PSF now shows some broken diffracted rings as observed in the Ups And PSF.

The simulations were done to determine the SD and DD subtraction performances. Images were magnified to the same number of pixels per Full Width Half Max to correct for the PSF chromatic dependence. The $`3\sigma `$ detection limits are presented in figure 4 and correspond to the theoretical limits of the TRIDENT differential simultaneous imaging technique with the current estimated static aberrations power spectrum.

Simulated aberrations in figures 2 and 3 reproduce the power spectrum of the observed data. The asymetries that can be seen in figure 3 (three upper frames) come from the aberrations randomly generated in the simulation. A companion 10 magnitudes fainter than the primary at $`0.5^{\prime \prime }`$ separation would be detected with the SD alone, while a companion 14 magnitudes fainter than the primary would be detected at the same separation with the DD.

## 4 Observations

Figure 5 show the observed PSF attenuation SD and DD for Ups And. Corresponding $`3\sigma `$ detection thresholds are presented in figure 6.

As seen in figure 6, a companion 7 magnitudes fainter than the primary would be detected at $`0.5^{\prime \prime }`$ separation. The Ups And SD and DD do not improve the detection limit contrary to simulated predictions (figure 4). Non-common path aberrations between the three optical paths (not included in the simulation presented in figure 4) are believed to interfere with common path aberrations (prior to TRIDENT beam splitter) and decorrelate the three PSFs. This shows the importance of minimizing non-common path aberrations when designing multiple optical path instruments.

## 5 PSF calibration techniques

PSF calibration techniques are needed to further reduce the PSF residual. A usual calibration technique consists of using other observations made during the same night but with different telescope orientations. This approach was found to be ineffective since the PSF structure varied slowly in time with telescope pointing and evolving atmospheric conditions. Figure 7 shows a set of July 2001 simple differences for a range of telescope pointings.

The PSF is stable for a given telescope pointing (Gl526 data set from July 8 and 9 2001, see two upper frame in figure 7) but it is evolving when changing the telescope orientation (figure 7, three lower frames). This PSF evolution is probably generated by optical flexions with changing telescope orientation and the differential atmospheric refraction between the optical wavefront sensor and the infrared detector, making the three optical paths change with time. However, a much more accurate PSF calibration is obtained if a reference star is acquired at approximately the same telescope pointing orientation within a few minutes of the target star observations.

### 5.1 Reference star subtraction technique

To minimize the effect of PSF evolution with telescope pointing, the reference stars were selected at approximately the same celestial declination and 15 minutes East or West of the target star.

A typical gain of 3 magnitudes in detection sensitivity is obtained with the combination of reference star subtraction and simultaneous imaging (figure 8) as compared with simultaneous imaging alone (figure 6), thus confirming the quasi-static nature of the PSF structure since the reference and object were acquired within a few minutes. It can be seen that the PSF subtraction residuals are still 3 magnitudes away from being photon noise limited at 0.5′′. Difficulty to calibrate optical aberrations due to PSF evolution with telescope pointing and atmospheric variations is the main limitation. A companion 10 magnitudes fainter than the primary would be detected at 0.5′′ separation.

### 5.2 Aberration smoothing by instrument rotation

Another way to minimize PSF structure effects generated after the instrument rotator is to smooth them by rotating the instrument (PUEO AO system and TRIDENT camera). The PSF attenuation obtained by PSF subtraction after smoothing, with a Strehl ratio of 0.15, is similar to that obtained with reference star PSF subtraction for data having a Strehl ratio of 0.5, showing that aberrations were effectively smoothed by instrument rotation and thus confirming that they are coming mostly from PUEO and TRIDENT (Figure 9).

The aberration smoothing technique has increased the correlation between the three optical paths by averaging non-common path aberration after the TRIDENT beam-splitter. The SD and DD with the PSF rotation technique (figure 9) are thus more efficient to subtract the PSF than differential imaging alone (figure 6). A companion 8.5 magnitudes fainter than the primary would be detected at 0.5′′ separation.

Figure 10 shows the $`\frac{3\sigma (I)}{I}`$ PSF suppression efficiency (also known as the figure of merit, or Q factor) for the differential simultaneous imaging technique with reference star subtraction or aberration smoothing by instrument rotation.

Both the reference subtraction and aberration smoothing techniques deliver similar PSF subtraction limits. Reference subtraction is slightly better for separation smaller than 0.6′′ while the aberration smoothing technique is slightly better for separations greater than 0.6′′. The next step would be to combine both calibration techniques.

## 6 TRIDENT new developments

Photon noise limited subtraction requires the minimization of optical aberrations. To increase PSF correlation a new optical design is presently under evaluation that mainly consists of a better beam splitter and filters placed near the detector to minimize wavefront degradation and non-common path aberrations. The goal is photon noise limited performance for a few hours integrations.

## 7 Conclusion - Quasi-static aberrations, an important factor

Quasi-static aberrations in AO corrected images are the main problem when attempting high contrast imaging for faint companion searches. If quasi-static aberrations cannot be minimized or removed in NGAO systems, an optimized camera will be essential to calibrate and subtract them. If a multiple optical path camera is considered, non-common path aberrations need to be minimized to avoid PSF decorrelation. Reference star subtraction and aberration smoothing by instrument rotation improve by one order of magnitude the faint companion detection threshold. Typical attenuation in good seeing conditions are 10 magnitudes in $`H`$ band at $`0.5^{\prime \prime }`$ on a 3.6 m telescope.

This work is supported in part through grants from NSERC, Canada and from Fonds FQRNT, Québec.


text:
# Variational charge renormalization in charged systems

## 1 Introduction

The behavior of charged systems has attracted renewed interest in the last few yearsHolm . One of the topics of interest is centered around the effects of multi-valent counterions, in which case correlations between ions become important and mean-field type theories break down. In such situations it is known experimentally and theoretically that equally charged bodies can attract each other. But even for monovalent ions and moderately to highly charged surfaces many open questions remain: In such situations the coupling parameter (which measures to which extent corrections to the mean-field or saddle-point solution are important) is not very high, but the surface potential can exceed thermal energy by far such that non-linear effects become important. The complication in this case is that the mean-field or Poisson-Boltzmann equation is a non-linear differential equation which in the presence of salt ions can be solved in closed form only in the planar geometrygouy . In the cylindrical and spherical geometry a number of different approximations have been proposed, which all more or less agree on the fact that counterions condense in the vicinity of the charged surface such as to reduce its effective surface chargeManning2 ; Pincus . The original Manning condensation argument states that the effective charge of a charged cylinder is maintained at a constant level of one charge per Bjerrum lengthManning . The prescriptions for calculating the effective charge are numerous and so are the predictions for the precise value of the effective charge as a function of external parameters such as salt, temperature, pH and so onBelloni ; Levin ; Deserno .

In this paper we show how to approach the problem using a field-theoretic formulation of a variational theory, within which the effective surface charge is treated as a variational parameter. There are in principle two different approaches, one is based on the standard density functional theory, within which the counterion distribution is variationally determined. This approach is clumsy for more complicated geometries and in principle is difficult to generalize such as to include strong-coupling (i.e. deviations from mean-field) effects. We therefore developed an alternative approach based on an exact field-theoretic formulation using the fluctuating imaginary electric potential. The surface charge is used as a variational parameter and the variational field-theoretic action is Gaussian. This is probably the most unambiguous definition of the renormalized or effective surface charge, since no additional assumptions (other than that the variational Hamiltonian is Gaussian) is used. As a trivial byproduct, the effective interaction between charges is at large distance of the Debye-Hückel type with renormalized effective charges. Our scheme in principle allows to go beyond mean-field, but we do not pursue such studies in this paper. As we show in the appendix, the standard Gibbs variational principle can be reformulated as a perturbative variational scheme which can be systematically improved by a perturbation series (which we also leave as a problem for the future). We test our approach for the planar charged wall, for which we essentially exactly reproduce the Poisson-Boltzmann result. We then tackle the charged cylinder. The effective cylinder charge interpolates continuously between the unrenormalized limit (at high salt concentrations) and the Manning charge-renormalized limit (at low salt concentrations), in agreement with previous resultsLevin . However, this crossover occurs for fully charged polyelectrolytes over 4 orders of magnitude of the salt concentration, and therefore it is important to take this slow crossover into account. For the controversial debate on the electrostatic persistence length it means that the standard Odijk-Skolnick-Fixman resultOdijk ; Skolnick ; Joanny , obtained on the linear (Debye-Hückel) level, is modified by the salt dependent effective polymer charge, which changes the salt dependency from $`\mathrm{}_{\mathrm{OSF}}\kappa ^2`$ to $`\mathrm{}_{\mathrm{OSF}}\kappa ^{1.4}`$, closer to experimental valuesForster . For the charged sphere the resultant behavior is different: Here charge renormalization takes place only at intermediate salt concentrations and both for high and low salt concentrations the standard Debye-Hückel potential with unrenormalized charges is validBelloni . Our method can be easily generalized to more complicated geometries and situations.

## 2 The model

The partition function of a system consisting of $`N_+`$ positively charged and $`N_{}`$ negatively charged ions of the same valency $`q`$, interacting with an arbitrary charge distribution $`\sigma (𝐫)`$, can be written as

$`𝒵={\displaystyle \frac{1}{N_+!}}{\displaystyle \underset{j=1}{\overset{N_+}{}}}\left[{\displaystyle \frac{\mathrm{d}𝐫_j^+}{\lambda _t^3}\mathrm{\Omega }(𝐫_j^+)}\right]{\displaystyle \frac{1}{N_{}!}}{\displaystyle \underset{j=1}{\overset{N_{}}{}}}\left[{\displaystyle \frac{\mathrm{d}𝐫_j^{}}{\lambda _t^3}\mathrm{\Omega }(𝐫_j^{})}\right]`$
$`\mathrm{exp}\{{\displaystyle \frac{\mathrm{}_B}{2}}{\displaystyle }\mathrm{d}𝐫\mathrm{d}𝐫^{}[q\widehat{\rho }_+(𝐫)q\widehat{\rho }_{}(𝐫)+\sigma (𝐫)]v(𝐫𝐫^{})`$
$`[q\widehat{\rho }_+(𝐫^{})q\widehat{\rho }_{}(𝐫^{})+\sigma (𝐫^{})]+{\displaystyle \frac{q^2\mathrm{}_B(N_++N_{})}{2}}v(0)\}.`$ (1)
The pair potential $`v(𝐫)`$ denotes the Coulomb interaction, $`v(𝐫)=1/r`$, but the first part of our discussion is rather general and independent of the specific form of $`v(𝐫)`$. The self-energy term in Eq.(2), proportional to $`q^2\mathrm{}_Bv(0)/2`$, simply subtracts the diagonal term in the double integral. It will be later shown to renormalize the chemical potential (or fugacity) in a rather trivial way. The Bjerrum length $`\mathrm{}_B=e^2/4\pi \epsilon k_BT`$ is the distance at which two unit charges interact with thermal energy $`k_BT`$. The length $`\lambda _t`$ is the thermal wavelength, which results from integrating out kinetic degrees of freedom, the precise value of which is unimportant. The function $`\mathrm{\Omega }(𝐫)`$ takes into account the presence of hard walls and restricts the positions of ions to an appropriate region in space, as will become more clear when we consider explicit examples in the following sections. The counterion density operators are

$$\widehat{\rho }_+(𝐫)=\underset{j=1}{\overset{N_+}{}}\delta (𝐫𝐫_j^+)$$

for positively charged ions with a similar expression for the negative ions. At this point, there is not much one can do with the partition function, mainly because the density operators $`\widehat{\rho }_+`$ and $`\widehat{\rho }_{}`$ enter in a quadratic fashion. The main step in deriving a field theory consists in introducing a unit operator which couples to the density operator $`\widehat{\rho }_+`$,

$`1`$ $`=`$ $`{\displaystyle 𝒟\rho _+\delta (\rho _+\widehat{\rho }_+)}`$ (2)
$`=`$ $`{\displaystyle 𝒟\rho _+𝒟\psi _+\mathrm{exp}\left\{ıd𝐫\psi _+(𝐫)[\rho _+(𝐫)\widehat{\rho }_+(𝐫)]\right\}},`$
(and a similar expression for the density of negative ions), where we used the integral representation of the delta function. This unit operator can be introduced into the partition function. Once introduced, it allows us to replace density operators $`\widehat{\rho }_+(𝐫)`$ and $`\widehat{\rho }_{}(𝐫)`$ by the corresponding fluctuating density fields $`\rho _+(𝐫)`$ and $`\rho _{}(𝐫)`$ and the configurational integrals over the ion positions can be performed. As a result, the partition function now reads

$`Z={\displaystyle 𝒟\rho _+𝒟\rho _{}𝒟\psi _+𝒟\psi _{}}`$
$`\mathrm{exp}\{{\displaystyle \frac{\mathrm{}_B}{2}}{\displaystyle }\mathrm{d}𝐫\mathrm{d}𝐫^{}[q\rho _+(𝐫)q\rho _{}(𝐫)+\sigma (𝐫)]v(𝐫𝐫^{})`$
$`[q\rho _+(𝐫^{})q\rho _{}(𝐫^{})+\sigma (𝐫^{})]`$
$`+ı{\displaystyle }\mathrm{d}𝐫[\psi _+(𝐫)\rho _+(𝐫)+\psi _{}(𝐫)\rho _{}(𝐫)]\}`$
$`{\displaystyle \frac{1}{N_+!}}\left[\mathrm{e}^{q^2\mathrm{}_Bv(0)/2}{\displaystyle \frac{\mathrm{d}𝐫}{\lambda _t^3}\mathrm{\Omega }(𝐫)\mathrm{e}^{ı\psi _+(𝐫)}}\right]^{N_+}`$
$`{\displaystyle \frac{1}{N_{}!}}\left[\mathrm{e}^{q^2\mathrm{}_Bv(0)/2}{\displaystyle \frac{\mathrm{d}𝐫}{\lambda _t^3}\mathrm{\Omega }(𝐫)\mathrm{e}^{ı\psi _{}(𝐫)}}\right]^N_{}.`$ (3)
The partition function becomes simpler upon transformation to the grand-canonical ensemble according to

$$Z_\lambda =\underset{N_+,N_{}=0}{\overset{\mathrm{}}{}}\lambda _0^{N_++N_{}}Z$$
(4)
where $`\lambda _0`$ is the bare fugacity which is related to the particle chemical potential $`\mu `$ by $`\lambda _0=\mathrm{e}^\mu `$. We assume here that the chemical potential of plus and minus ions is the same. Using the definition of the exponential function, $`\mathrm{e}^x=_{N=0}^{\mathrm{}}x^N/N!`$, the grand-canonical partition function can be written as

$$Z_\lambda =𝒟\rho _+𝒟\rho _{}𝒟\psi _+𝒟\psi _{}\mathrm{exp}\left\{H[\rho _+,\rho _{},\psi _+,\psi _{}]\right\}$$
(5)
with the field-theoretic action given by

$`H[\rho _+,\rho _{},\psi _+,\psi _{}]=`$
$`{\displaystyle \frac{\mathrm{}_B}{2}}{\displaystyle d𝐫d𝐫^{}[q\rho _+(𝐫)q\rho _{}(𝐫)+\sigma (𝐫)]v(𝐫𝐫^{})}`$
$`[q\rho _+(𝐫^{})q\rho _{}(𝐫^{})+\sigma (𝐫^{})]`$
$`ı{\displaystyle d𝐫\left[\psi _+(𝐫)\rho _+(𝐫)+\psi _{}(𝐫)\rho _{}(𝐫)\right]}`$
$`\lambda {\displaystyle d𝐫\mathrm{\Omega }(𝐫)\left[\mathrm{e}^{ı\psi _+(𝐫)}+\mathrm{e}^{ı\psi _{}(𝐫)}\right]}`$ (6)
where we defined a rescaled fugacity as $`\lambda =\lambda _0\mathrm{e}^{q^2\mathrm{}_Bv(0)/2}/\lambda _t^3`$. As follows from the definition of the grand-canonical partition function, Eq.(4), the expectation value of the total particle number is given by

$`N_++N_{}=\lambda _0{\displaystyle \frac{\mathrm{ln}Z_\lambda }{\lambda _0}}=\lambda {\displaystyle \frac{\mathrm{ln}Z_\lambda }{\lambda }}=`$
$`\lambda {\displaystyle d𝐫\mathrm{\Omega }(𝐫)\mathrm{e}^{ı\psi _+(𝐫)}+\mathrm{e}^{ı\psi _{}(𝐫)}}.`$ (7)
As can be be seen, the expectation value of the particle number, $`N`$, is independent of any multiplicative factor of the fugacity. The standard way of saying this is that the chemical potential is only defined up to a constant.

### 2.1 Density description

One type of approximations starts by eliminating the fluctuating fields $`\psi _+`$ and $`\psi _{}`$ from the problem. The simplest way of doing this is to replace the functional integral over these fields by the value of the integrand at the extremum or optimal value. The saddle point equations which determine the optimal values of the fluctuating fields, $`\psi _+^{SP}`$ and $`\psi _{}^{SP}`$, are given by

$$\frac{\delta H[\rho _+,\rho _{},\psi _+,\psi _{}]}{\delta \psi _+(𝐫)}=0$$
(8)
and with a similar equation for $`\psi _{}`$. Using $`H`$ as defined in Eq.(2), the solutions are

$$\psi _+^{SP}(𝐫)=ı\mathrm{ln}(\rho _+(𝐫)/\lambda ))$$
(9)
(and a similar solution for $`\psi _{}`$) which, upon insertion into the Hamiltonian Eq.(2), leads to the saddle-point Hamiltonian

$`H_\lambda [\rho _+,\rho _{},\psi _+^{SP},\psi _{}^{SP}]=`$
$`{\displaystyle \frac{\mathrm{}_B}{2}}{\displaystyle d𝐫d𝐫^{}[q\rho _+(𝐫)q\rho _{}(𝐫)+\sigma (𝐫)]v(𝐫𝐫^{})}`$
$`[q\rho _+(𝐫^{})q\rho _{}(𝐫^{})+\sigma (𝐫^{})]`$
$`+{\displaystyle d𝐫\rho _+(𝐫)\left(\mathrm{ln}[\rho _+(𝐫)/\lambda ]1\right)}.`$
$`+{\displaystyle d𝐫\rho _{}(𝐫)\left(\mathrm{ln}[\rho _{}(𝐫)/\lambda ]1\right)}.`$ (10)
This expression is the starting point for local density functional theoryStevens . The extremum of Hamiltonian Eq.(2.1), determined by the equation

$$\frac{\delta H[\rho _+,\rho _{},\psi _+^{SP},\psi _{}^{SP}]}{\delta \rho _+(𝐫)}=0$$
(11)
(and with a similar equation for $`\rho _{}`$), is the mean-field or Poisson-Boltzmann density distribution. As is well-known, the mean-field equation (11) can be solved in closed form only for a few cases, and for more complicated situations one has to rely on additional approximations. One popular way of doing this is to minimize the Hamiltonian Eq.(2.1) for a restricted functional form of the ion distribution, for example a box distribution. In the appendix A1 we show how this can be done for a planar charged wall (where comparison with the closed-form mean-field solution is possible). The draw back of box models is that the ion distribution is discontinuous and thus rather artificial and for more complicated charge distributions the box geometry becomes complicated. Another, more basic drawback of the density functional Eq.(2.1) is that correlation effects are not included on the local level. Non-local terms can in principle be taken into account by expanding around the saddle-point in Eq.(9), but the resulting functional is complicated. In the next section we show that a variational approach is more useful when applied in the potential description.

### 2.2 Potential description

A more powerful variational formulation is possible by choosing a description in terms of the fluctuating potentials $`\psi _+`$ and $`\psi _{}`$. The starting point is the observation that the integral over the density fields $`\rho _+`$ and $`\rho _{}`$ in Eq. (5) can be performed exactly. Introducing the new field $`\varphi (𝐫)=[\psi _+(𝐫)\psi _{}(𝐫)]/2`$, the partition function can be written exactly as

$$Z_\lambda =\frac{𝒟\varphi }{Z_v}\mathrm{exp}\left\{H[\varphi ]\right\}$$
(12)
with

$$H[\varphi ]=d𝐫\left\{\frac{[\varphi (𝐫)]^2}{8\pi \mathrm{}_Bq^2}+\frac{ı\sigma (𝐫)\varphi (𝐫)}{q}2\lambda \mathrm{\Omega }(𝐫)\mathrm{cos}\varphi (𝐫)\right\}$$
(13)
where we used the inverse of the Coulomb potential $`v^1(𝐫)=^2\delta (𝐫)/4\pi `$. The symbol $`𝒵_v`$ denotes the determinant of the Gaussian integral

$$Z_v=𝒟\varphi \mathrm{exp}\left\{\frac{1}{2q^2\mathrm{}_B}d𝐫d𝐫^{}\varphi (𝐫)v^1(𝐫𝐫^{})\varphi (𝐫^{})\right\}$$
(14)
and is a measure of the free energy of vacuum fluctuations.

Next, we rescale all lengths by the Gouy-Chapman length $`\mu =1/(2\pi q\mathrm{}_B\sigma _s)`$ according to $`𝐫=\mu \stackrel{~}{𝐫}`$, where $`\sigma _s`$ denotes the surface charge density of the charged object under consideration. We obtain the modified partition function

$$Z_\lambda =\frac{𝒟\varphi }{Z_v}\mathrm{exp}\left\{\stackrel{~}{H}[\varphi ]/\mathrm{\Xi }\right\}$$
(15)
with the rescaled action

$$\stackrel{~}{H}[\varphi ]=\frac{\mathrm{d}\stackrel{~}{𝐫}}{2\pi }\left[\frac{1}{4}\left(\varphi (\stackrel{~}{𝐫})\right)^2+ı\varphi (\stackrel{~}{𝐫})\stackrel{~}{\sigma }(\stackrel{~}{𝐫})\frac{\mathrm{\Lambda }\stackrel{~}{\mathrm{\Omega }}(\stackrel{~}{𝐫})\mathrm{cos}\varphi (\stackrel{~}{𝐫})}{2}\right]$$
(16)
where we used the rescaled charge distribution $`\stackrel{~}{\sigma }(𝐫/\mu )=\mu \sigma (𝐫)/\sigma _s`$ and ion exclusion function $`\stackrel{~}{\mathrm{\Omega }}(𝐫/\mu )=\mathrm{\Omega }(𝐫)`$. The rescaled fugacity $`\mathrm{\Lambda }`$ is defined by

$$\mathrm{\Lambda }=8\pi \lambda \mu ^3\mathrm{\Xi }=\frac{2\lambda }{\pi \sigma _s^2\mathrm{}_B}$$
(17)
and is determined by the normalization condition (in the case of fixed particle number)

$`N_++N_{}=\mathrm{\Lambda }{\displaystyle \frac{\mathrm{ln}Z_\mathrm{\Lambda }}{\mathrm{\Lambda }}}=`$
$`{\displaystyle \frac{\mathrm{\Lambda }}{4\pi \mathrm{\Xi }}}{\displaystyle d\stackrel{~}{𝐫}\stackrel{~}{\mathrm{\Omega }}(\stackrel{~}{𝐫})\mathrm{cos}\varphi (\stackrel{~}{𝐫})}.`$ (18)
As one can see, there are two parameters left in the field-theoretic formulation, namely the coupling parameter

$$\mathrm{\Xi }=2\pi \mathrm{}_B^2\sigma _sq^3,$$
(19)
which measures to which extent deviations from the saddle-point are important, and the ion fugacity $`\mathrm{\Lambda }`$. Infinitely far from any charged object the potential is constant and the concentration of minus and plus ions is equal and given by

$$c_s=\frac{N_+}{V}=\frac{N_{}}{V}=\frac{\mathrm{\Lambda }}{8\pi \mathrm{\Xi }\mu ^3}\mathrm{cos}\varphi $$
(20)
which, using the definition of the screening length, $`\kappa ^1=(8\pi q^2\mathrm{}_Bc_s)^{1/2}`$ and the Gouy Chapman length can be rewritten as

$$\kappa ^2\mu ^2=\stackrel{~}{\kappa }^2=\mathrm{\Lambda }\mathrm{cos}\varphi .$$
(21)
The rescaled fugacity $`\mathrm{\Lambda }`$ thus measures the ratio between the screening length and the Gouy-Chapman length: For $`\stackrel{~}{\kappa }1`$ the screening by salt is dominant, non-linear effects are unimportant and Debye-Hückel theory is valid, for $`\stackrel{~}{\kappa }1`$ on the other hand non-linear effects are important and Debye-Hückel theory fails. It is the latter case where our variational theory becomes important.

## 3 Variational Theory

The action Eq.(16) is too complicated to be solved exactly. Previously, perturbative theories have been proposed for the weak-coupling (small $`\mathrm{\Xi }`$) regime in terms of a loop-expansion around the saddle-point solution (which corresponds to the solution of the non-linear Poisson-Boltzmann equation)Netz1  and for the strong-coupling (large $`\mathrm{\Xi }`$) regime in terms of a virial expansionNetz2 ; Moreira1 ; Moreira2 . Both expansions are of limited usefulness because of their bad convergence properties as one moves into the intermediate regime between the two asymptotic limitsMoreira3 . It is this point where variational theories come in as a robust means to obtain reliable, useful results.

The standard Gibbs variational procedure consists in minimizing the following variational free energy

$$_{Gibbs}=_0\stackrel{~}{H}_0\stackrel{~}{H}_0/\mathrm{\Xi }$$
(22)
with respect to the variational parameters included in the variational Hamiltonian $`\stackrel{~}{H}_0`$. As discussed in detail in Appendix B the Gibbs variational principle is equivalent to the first-order perturbational variation on the observable which is conjugate to the variational parameter. The perturbational variation scheme can be systematically improved by going to higher order, the Gibbs variation cannot. Since we stay at the first-order level throughout this paper, we will use the Gibbs variational notation, though one should keep in mind that this can be interpreted as the first order result of a perturbational variational calculation.

### 3.1 Most general Gaussian variation

The most general, yet exactly tractable variational Hamiltonian, is a Gaussian one

$$\stackrel{~}{H}_0[\varphi ]=\frac{1}{2}_{𝐫,𝐫^{}}[\varphi (𝐫)ı\varphi _0(𝐫)]v_0^1(𝐫,𝐫^{})[\varphi (𝐫^{})ı\varphi _0(𝐫^{})]$$
(23)
where the Gaussian kernel $`v_0^1`$ and the mean potential $`\varphi _0`$ are variational parameters. The variational free energy according to Eq. (22) reads

$`_{Gibbs}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{Tr}\mathrm{log}v_0{\displaystyle _𝐫}{\displaystyle \frac{(\varphi _0)^2}{8\pi \mathrm{\Xi }}}`$ (24)
$`+{\displaystyle _{𝐫,𝐫^{}}}\delta (𝐫𝐫^{}){\displaystyle \frac{_r_r^{}v_0(𝐫,𝐫^{})}{8\pi \mathrm{\Xi }}}`$
$`{\displaystyle _𝐫}{\displaystyle \frac{\stackrel{~}{\sigma }(𝐫)\varphi _0(𝐫)}{2\pi \mathrm{\Xi }}}`$
$`{\displaystyle \frac{\mathrm{\Lambda }}{4\pi \mathrm{\Xi }}}{\displaystyle _𝐫}\stackrel{~}{\mathrm{\Omega }}\mathrm{e}^{\mathrm{\Xi }v_0(𝐫,𝐫)/2}\mathrm{cosh}(\varphi _0).`$
The variational equations are given by $`\delta _{Gibbs}/\delta \varphi _0(𝐫)=0`$ and $`\delta _{Gibbs}/\delta v_0^1(𝐫,𝐫^{})=0`$ and read

$`{\displaystyle \frac{^2\varphi _0}{4\pi }}={\displaystyle \frac{\stackrel{~}{\sigma }}{2\pi }}+{\displaystyle \frac{\mathrm{\Lambda }}{4\pi }}\mathrm{\Omega }e^{\mathrm{\Xi }v_0(𝐫,𝐫)/2}\mathrm{sinh}\varphi _0`$ (25)
$`{\displaystyle \frac{^2v_0(𝐫,𝐫^{})}{4\pi \mathrm{\Xi }}}+{\displaystyle \frac{\mathrm{\Lambda }}{4\pi }}e^{\mathrm{\Xi }v_0(𝐫,𝐫)/2}\mathrm{cosh}\varphi _0v_0(𝐫,𝐫^{})`$
$`=\delta (𝐫𝐫^{})`$ (26)
These self-consistent equations are quite complicated even for the simplest charge distribution. In the following we will devise a simpler, one-parameter variational scheme which will lead -as a byproduct- to a very simple and clear definition of the renormalized or effective charge of an object.

### 3.2 Variational Charge Renormalization

As we will show in the following, non-linear effects are captured in an essentially exact way already by a much simpler variational form, which as the only variational parameter contains a charge-renormalization factor. The variational Hamiltonian reads

$$\stackrel{~}{H}_0[\varphi ]=\frac{\mathrm{d}\stackrel{~}{𝐫}}{2\pi }\left[\frac{1}{4}\left(\varphi (\stackrel{~}{𝐫})\right)^2+ı\eta (\stackrel{~}{𝐫})\varphi (\stackrel{~}{𝐫})\stackrel{~}{\sigma }(\stackrel{~}{𝐫})+\frac{\stackrel{~}{\kappa }^2\varphi ^2(\stackrel{~}{𝐫})}{4}\right]$$
(27)
where the variational function $`\eta (𝐫)`$ multiplies the charge distribution and therefore corresponds to a spatially dependent charge renormalization factor. The Hamiltonian is quadratic in the fluctuating field and therefore all expectation values can be calculated by using Green’s function formalism. Also, the geometry-dependent factor $`\mathrm{\Omega }(𝐫)`$ is for the present formulation omitted and the salt concentration dependent term proportional to $`\stackrel{~}{\kappa }^2`$ has the same strength throughout the space. This is certainly a simplification, which however allows to use a simple translationally invariant Green’s function (future elaborations should probe variational forms with a spatially dependent mass term). Using the Debye-Hückel potential

$$v_{\mathrm{DH}}(𝐫)=\frac{\mathrm{e}^{\stackrel{~}{\kappa }r}}{r}$$
(28)
and the potential distribution

$$\varphi _0(𝐫)=\frac{1}{2\pi }_𝐫^{}\eta (𝐫^{})\stackrel{~}{\sigma }(𝐫^{})v_{\mathrm{DH}}(𝐫𝐫^{})$$
(29)
the variational Hamiltonian Eq.(27) can be rewritten as

$`\stackrel{~}{H}_0[\varphi ]`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle _{𝐫,𝐫^{}}}[\varphi (𝐫)+ı\varphi _0(𝐫)]v_{\mathrm{DH}}^1(𝐫,𝐫^{})[\varphi (𝐫^{})+ı\varphi _0(𝐫^{})]`$ (30)
$`+{\displaystyle \frac{1}{4\pi }}{\displaystyle _𝐫}\eta (𝐫)\stackrel{~}{\sigma }(𝐫)\varphi _0(𝐫).`$
The variational free energy according to Eq.(22) is easily calculated and reads

$`_{Gibbs}`$ $`=`$ $`_{\mathrm{DH}}+{\displaystyle \frac{1}{4\pi \mathrm{\Xi }}}{\displaystyle _𝐫}(2\eta (𝐫))\stackrel{~}{\sigma }(𝐫)\varphi _0(𝐫)`$
$`+{\displaystyle \frac{\stackrel{~}{\kappa }^2}{8\pi \mathrm{\Xi }}}{\displaystyle _𝐫}[\varphi _0^2(𝐫)\mathrm{\Xi }v_{\mathrm{DH}}(0)]`$
$`{\displaystyle \frac{\stackrel{~}{\kappa }^2}{4\pi \mathrm{\Xi }}}{\displaystyle _𝐫}\stackrel{~}{\mathrm{\Omega }}(𝐫)\mathrm{cosh}[\varphi _0(𝐫)]`$
where we used the short-hand notation

$`_{\mathrm{DH}}`$ $`=`$ $`\mathrm{ln}{\displaystyle 𝒟\varphi \mathrm{e}^{_{𝐫,𝐫^{}}\varphi (𝐫)v_{\mathrm{DH}}^1(𝐫𝐫^{})\varphi (𝐫^{})/2\mathrm{\Xi }}}`$ (32)
$`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{d^3q}{(2\pi )^3}\mathrm{log}(𝐪^2+\stackrel{~}{\kappa }^2)}`$
and fixed the fugacity $`\mathrm{\Lambda }`$ according to Eq.(21). The variational equation follows by taking a functional derivative of the free energy with respect to the charge renormalization function, $`\delta _{Gibbs}/\delta \eta (𝐫)=0`$, and reads after a few algebraic manipulations

$`{\displaystyle _𝐫}(1\eta (𝐫))\stackrel{~}{\sigma }(𝐫)\varphi _0(𝐫)=`$
$`{\displaystyle \frac{\stackrel{~}{\kappa }^2}{2}}{\displaystyle _𝐫}\varphi _0(𝐫)\left(\stackrel{~}{\mathrm{\Omega }}(𝐫)\mathrm{sinh}[\varphi _0(𝐫)]\varphi _0(𝐫)\right).`$ (33)
This integral equation is to be solved by a suitably chosen function $`\eta (𝐫)`$. Since in this paper we will be only concerned with situation where translational invariance holds, we will in the following section specialize the variational approach to the homogeneous case.

### 3.3 Homogeneous variational approach

The homogeneous variational approach is obtained by choosing a charge renormalization which is uniform in space, i.e. $`\eta (𝐫)=\eta _0`$. We redefine the bare potential as

$$\varphi _0(𝐫)=\frac{1}{2\pi }_𝐫^{}\stackrel{~}{\sigma }(𝐫^{})v_{\mathrm{DH}}(𝐫𝐫^{})$$
(34)
and the Hamiltonian can be rewritten as

$`\stackrel{~}{H}_0[\varphi ]`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle _{𝐫,𝐫^{}}}[\varphi (𝐫)+ı\eta _0\varphi _0(𝐫)]v_{\mathrm{DH}}^1(𝐫,𝐫^{})[\varphi (𝐫^{})+ı\eta _0\varphi _0(𝐫^{})]`$ (35)
$`+{\displaystyle \frac{\eta _0^2}{4\pi }}{\displaystyle _𝐫}\stackrel{~}{\sigma }(𝐫)\varphi _0(𝐫).`$
The variational free energy according to Eq.(22) reads

$`_{Gibbs}`$ $`=`$ $`_{\mathrm{DH}}+{\displaystyle \frac{1}{4\pi \mathrm{\Xi }}}{\displaystyle _𝐫}(2\eta _0)\eta _0\stackrel{~}{\sigma }(𝐫)\varphi _0(𝐫)`$
$`+{\displaystyle \frac{\stackrel{~}{\kappa }^2}{8\pi \mathrm{\Xi }}}{\displaystyle _𝐫}[\eta _0^2\varphi _0^2(𝐫)\mathrm{\Xi }v_{\mathrm{DH}}(0)]`$
$`{\displaystyle \frac{\stackrel{~}{\kappa }^2}{4\pi \mathrm{\Xi }}}{\displaystyle _𝐫}\stackrel{~}{\mathrm{\Omega }}(𝐫)\mathrm{cosh}[\eta _0\varphi _0(𝐫)]`$
The variational equation follows by taking a derivative with respect to the charge renormalization factor, $`_{Gibbs}/\eta _0=0`$, and is given by

$`2(1\eta _0){\displaystyle _𝐫}\stackrel{~}{\sigma }(𝐫)\varphi _0(𝐫)`$ (37)
$`=\stackrel{~}{\kappa }^2{\displaystyle _𝐫}\varphi _0(𝐫)\left(\stackrel{~}{\mathrm{\Omega }}(𝐫)\mathrm{sinh}[\eta _0\varphi _0(𝐫)]\eta _0\varphi _0(𝐫)\right).`$
In the following we will solve this self-consistent equation for the charge-renormalization factor $`\eta _0`$ in the planar, cylindrical, and spherical geometries.

#### 3.3.1 Planar case

For the calculation of the uniformly charged planar wall we assume that salt is present on both sides of the wall, the geometry factor therefore equals unity $`\stackrel{~}{\mathrm{\Omega }}(𝐫)=1`$ (it is an easy exercise to generalize the calculation to the case where salt is present only on one side of the charged wall). With the bare charge distribution $`\sigma (𝐫)=\sigma _s\delta (z)`$ the rescaled charge distribution becomes $`\stackrel{~}{\sigma }(\stackrel{~}{𝐫})=\delta (\stackrel{~}{z})`$. The bare potential is according to Eq.(34) given by $`\varphi _0(\stackrel{~}{z})=\mathrm{e}^{\stackrel{~}{\kappa }\stackrel{~}{z}}/\stackrel{~}{\kappa }`$. The integrals in Eq.(37) can all be performed exactly, and the resulting equation becomes

$$1\frac{\eta _0}{2}\frac{\stackrel{~}{\kappa }^2}{\eta _0}\left(\mathrm{cosh}[\eta _0/\stackrel{~}{\kappa }]1\right)=0.$$
(38)
In the limit $`\stackrel{~}{\kappa }1`$ the solution is

$$\eta _01\frac{1}{24\stackrel{~}{\kappa }^2}$$
(39)
and therefore only weak charge-regulation takes place, as expected for a weakly charged surface or high salt concentrations. In the opposite limit of a highly charged surface or at low salt concentrations, for $`\stackrel{~}{\kappa }1`$, the solution is

$$\eta _0\stackrel{~}{\kappa }\mathrm{ln}(\stackrel{~}{\kappa }).$$
(40)
These results can be directly compared with the Poisson-Boltzmann results derived explicitly in Appendix A2. The quantity one wants to compare is the surface potential, since this is conjugate to the surface charge density, i.e. the variational parameter, and one can reinterpret the Gibbs variational principle as a first-order perturbational variation on the expectation value of the surface potential (as is explained in more detail in Appendix B). Within our scheme, the surface potential is defined by

$$\psi _0=\eta _0\varphi _0(\stackrel{~}{z}=0)$$

and is given in units of $`k_BT`$ and valency $`q`$. For the planar case, we obtain $`\psi _0=\eta _0/\stackrel{~}{\kappa }`$, and therefore find for the surface potential $`\eta _0/\stackrel{~}{\kappa }=1/\stackrel{~}{\kappa }`$ at high salt ($`\stackrel{~}{\kappa }1`$) and $`\eta _0/\stackrel{~}{\kappa }=\mathrm{ln}(\stackrel{~}{\kappa })`$ at low salt ($`\stackrel{~}{\kappa }1`$). These are the surface potentials one obtains from the full non-linear Poisson-Boltzmann theory (except a factor of two for the low-salt case). This shows that our variational procedure captures non-linear effects accurately and justifies to consider more complicated geometries, where the Poisson-Boltzmann equation can only be solved numerically.

#### 3.3.2 Cylindrical case

We assume a charged cylinder with radius $`R`$ and linear charge density $`\tau `$. The surface charge density follows as $`\sigma _s=\tau /(2\pi R)`$. Defining the Gouy-Chapman length (like in the planar case) as $`\mu =1/(2\pi q\sigma _s)`$, we obtain $`\mu =R/(q\tau \mathrm{}_B)`$. The rescaled cylinder radius, $`\stackrel{~}{R}=R/\mu =q\tau \mathrm{}_B`$ is thus equal to the Manning parameterManning : For $`\stackrel{~}{R}>1`$ one expects counterions to be bound to the cylinder even in the limit of vanishing salt concentration, for $`\stackrel{~}{R}<1`$ no counterions stay bound in that limit. In the following, we assume the cylinder to be penetrable for ions, the geometry function is $`\stackrel{~}{\mathrm{\Omega }}(\stackrel{~}{r})=1`$ for all values of $`\stackrel{~}{r}`$. Using cylindrical coordinates, the potential $`\varphi _0`$ as defined in Eq.(34) satisfies the differential equation

$$\frac{1}{\stackrel{~}{r}}\frac{}{\stackrel{~}{r}}\stackrel{~}{r}\frac{}{\stackrel{~}{r}}\varphi _0(\stackrel{~}{r})+2\delta (\stackrel{~}{r}\stackrel{~}{R})+\stackrel{~}{\kappa }^2\varphi _0(\stackrel{~}{r})=0$$

with the boundary conditions $`\varphi _0(\stackrel{~}{R}0)=\varphi _0(\stackrel{~}{R}+0)`$ and $`\varphi _0^{}(\stackrel{~}{R}+0)\varphi _0^{}(\stackrel{~}{R}0)=2`$. The complete solution is constructed from the two solutions of the homogeneous differential equation, $`\varphi _0(\stackrel{~}{r})=K_0(\stackrel{~}{\kappa }\stackrel{~}{r})`$ and $`\varphi _0(\stackrel{~}{r})=I_0(\stackrel{~}{\kappa }\stackrel{~}{r})`$, and reads

$$\varphi _0(\stackrel{~}{r})=2\stackrel{~}{R}I_0(\stackrel{~}{\kappa }\stackrel{~}{R})K_0(\stackrel{~}{\kappa }\stackrel{~}{r})$$
(41)
for $`\stackrel{~}{r}>\stackrel{~}{R}`$ and

$$\varphi _0(\stackrel{~}{r})=2\stackrel{~}{R}K_0(\stackrel{~}{\kappa }\stackrel{~}{R})I_0(\stackrel{~}{\kappa }\stackrel{~}{r})$$
(42)
for $`\stackrel{~}{r}<\stackrel{~}{R}`$. The solution of the integral equation (37) can be easily obtained graphically by calculating $`\eta _0`$ as a function of $`\stackrel{~}{\kappa }\stackrel{~}{R}`$ and $`\eta _0/\stackrel{~}{\kappa }`$ via numerical integration followed by a parametric plot. In Fig. 1a we show the renormalization factor $`\eta _0`$ as a function of the Manning parameter $`\stackrel{~}{R}`$ for fixed values of the rescaled inverse screening length $`\kappa R=5,1,0.2,0.01`$ (thin solid lines from top to bottom). The thick solid line denotes the Manning limit

$$\eta _0=1/\stackrel{~}{R}=1/(\tau \mathrm{}_Bq)$$
(43)
which is very slowly approached as the salt concentration decreasesLevin . Indeed, an asymptotic analysis of the self-consistent equation (37) with the cylinder potential Eqs. (41,42) gives for $`\kappa R1`$

$$\eta _0\frac{1}{\stackrel{~}{R}}\left(1+\frac{\mathrm{ln}\stackrel{~}{R}}{2\mathrm{ln}(1/\kappa R)}\right)$$
(44)
which shows that the approach to the Manning limit is logarithmically slow. In the unrescaled units, the effective cylinder charge density is

$$\tau _{\mathrm{eff}}=\tau \eta _0\frac{1}{\mathrm{}_Bq}\left(1+\frac{\mathrm{ln}\tau \mathrm{}_Bq}{2\mathrm{ln}(1/\kappa R)}\right)$$
(45)
In Fig.1b we show the data for $`\eta _0`$ for fixed Manning parameter $`\stackrel{~}{R}=3`$ (which roughly corresponds to a fully charged vinyl-based polyelectrolyte with a distance of $`0.254`$ $`nm`$ between charges along the backbone) as a function of $`\kappa R`$ together with the asymptotic prediction Eq.(44) as a broken line in a semi-logarithmic plot and in Fig.1c as a function of $`1/\mathrm{ln}(\kappa R)`$. It is seen that the asymptotic law is valid for low salt concentrations satisfying $`\kappa R<0.1`$. For large salt concentrations, asymptotic analysis predicts that the charge-renormalization vanishes according to

$$\eta _01\frac{1}{24\stackrel{~}{\kappa }^2}$$
(46)
(identical to the result for the planar case) which is shown in Fig.1b as a dotted line.

#### 3.3.3 Nonlinear electrostatic persistence length

The electrostatic contribution to the bending rigidity has been calculated on the linear levelOdijk ; Skolnick  but also on the non-linear level using the full Poisson-Boltzmann equation for a bent cylinderLeBret ; Fixman . Here we re-address the relevance of non-linear effects for the persistence length of highly charged polyelectrolytes. In Fig.1d we plot the charge renormalization factor on double logarithmic scales. Over roughly a decade in $`\kappa R`$, which corresponds to two decades in the salt concentration, the charge renormalization factor can be approximately described by a power law $`\eta _0\kappa ^\alpha `$ with $`\alpha 0.3`$ (shown as a solid line). This range of salt concentrations is quite relevant for experimentsForster , since it has been reported that the experimentally measured salt-concentration dependence of the persistence length deviates from the Odijk-Skolnick-Fixman prediction, according to which

$$\mathrm{}_{\mathrm{OSF}}=\frac{\mathrm{}_B\tau ^2}{4\kappa ^\beta }$$
(47)
with $`\beta =2`$Odijk ; Skolnick . Replacing the bare line charge density $`\tau `$ in Eq.(47) by the effective charge density

$$\tau _{\mathrm{eff}}=\eta _0\tau \kappa ^\alpha \tau $$
(48)
(which is permissible since to leading order in an expansion for large bending radii modifications of the charge renormalization due to polymer bending are irrelevant) one obtains a modified exponent $`\beta =22\alpha 1.4`$ which is closer to the experimentally measured exponentsForster . Other effects that have been shown to reduce the exponent $`\beta `$ are relevant for weakly charge chains which are crumpled at small length scalesJoanny ; Sim2 ; Micka ; Thirumalai  although recent simulations using linear Debye-Hückel potentialsRalf ; Boris  as well as general Gaussian variational calculationsNetz4  tend to confirm the original Odijk-Skolnick-Fixman prediction with $`\beta =2`$. It might therefore turn out that the true explanation for the discrepancy between experimental and theoretical predictions for the electrostatic persistence length is due to non-linear effects which had been neglected in most theoretical treatments but are of course omnipresent in experiments. Needless to say, we expect modification of a whole number of other scaling relations for polyelectrolytes which depend on the line charge density $`\tau `$.

#### 3.3.4 Spherical Case

We now treat a sphere of charge $`Z`$ and radius $`R`$, leading to a surface charge density $`\sigma _s=Z/(4\pi R^2)`$ and Gouy-Chapman length $`\mu =1/(2\pi q\mathrm{}_B\sigma _s)`$. The rescaled radius is given by $`\stackrel{~}{R}=R/\mu =q\mathrm{}_BZ/(2R)`$.

Using spherical coordinates, the potential $`\varphi _0`$ as defined in Eq.(34) satisfies the differential equation

$$\frac{1}{\stackrel{~}{r}^2}\frac{}{\stackrel{~}{r}}\stackrel{~}{r}^2\frac{}{\stackrel{~}{r}}\varphi _0(\stackrel{~}{r})+2\delta (\stackrel{~}{r}\stackrel{~}{R})+\stackrel{~}{\kappa }^2\varphi _0(\stackrel{~}{r})=0$$

with the boundary conditions $`\varphi _0(\stackrel{~}{R}0)=\varphi _0(\stackrel{~}{R}+0)`$ and $`\varphi _0^{}(\stackrel{~}{R}+0)\varphi _0^{}(\stackrel{~}{R}0)=2`$. The two solutions of the homogeneous differential equation are $`\varphi _0(\stackrel{~}{r})=\mathrm{e}^{\stackrel{~}{\kappa }\stackrel{~}{r}}/\stackrel{~}{r}`$ and $`\varphi _0(\stackrel{~}{r})=\mathrm{e}^{\stackrel{~}{\kappa }\stackrel{~}{r}}/\stackrel{~}{r}`$, from which we obtain the solutions

$$\varphi _0(\stackrel{~}{r})=2\stackrel{~}{R}\mathrm{sinh}(\stackrel{~}{\kappa }\stackrel{~}{R})\mathrm{e}^{\stackrel{~}{\kappa }\stackrel{~}{r}}/(\stackrel{~}{\kappa }\stackrel{~}{r})$$
(49)
for $`\stackrel{~}{r}>\stackrel{~}{R}`$ and

$$\varphi _0(\stackrel{~}{r})=2\stackrel{~}{R}\mathrm{e}^{\stackrel{~}{\kappa }\stackrel{~}{R}}\mathrm{sinh}(\stackrel{~}{\kappa }\stackrel{~}{r})/(\stackrel{~}{\kappa }\stackrel{~}{r})$$
(50)
for $`\stackrel{~}{r}<\stackrel{~}{R}`$ where we have again assumed that salt ions are present everywhere, i.e. $`\stackrel{~}{\mathrm{\Omega }}(\stackrel{~}{r})=1`$ throughout space.

In Fig.2a we show parametric plots of the charge renormalization factor for $`\kappa R=100,10,1,0.4`$ (from top to bottom) and in Fig.2b for $`\kappa R=0.4,,0.1,0.01,0.001`$ (from bottom to top). In contrast to the cylindrical case, charge renormalization is most pronounced at intermediate salt concentrations corresponding to a screening length which roughly equals the sphere radiusBelloni . For large salt concentrations, asymptotic analysis predicts that the charge renormalization vanishes according to

$$\eta _01\frac{1}{24\stackrel{~}{\kappa }^2}$$
(51)
in agreement with the planar and cylindrical geometries. For low salt concentrations, an asymptotic analysis of the self-consistent equation (37) with the spherical potential Eqs. (49,50) gives for $`\kappa R1`$

$$\eta _0\frac{1}{\stackrel{~}{2R}}\mathrm{ln}\left(\frac{1}{\stackrel{~}{R}\stackrel{~}{\kappa }^2}\right).$$
(52)
In unrescaled units this leads to an effective sphere charge

$$Z_{\mathrm{eff}}=\eta _0Z\frac{R}{q\mathrm{}_B}\mathrm{ln}\left(\frac{q\mathrm{}_BZ}{R^3\kappa ^2}\right)$$
(53)
which only depends very weakly on the bare sphere charge (in agreement with previous qualitative argumentsPincus ). In Fig.2c we compare data for a sphere with rescaled radius $`\stackrel{~}{R}=R/\mu =10`$ with the asymptotic low-salt prediction Eq.(52) (broken line) and the high-salt prediction Eq.(51) (dotted line). Note that the asymptotic low-salt prediction is only valid in a finite salt-concentration limit: for very low salt concentration, the counterions are entropically driven away from the sphere and the bare charge is recovered. This might be detectable in light-scattering experiments on de-salted charged-colloid suspensions.

## 4 Conclusions

We have introduced the framework for a general variational technique which allows to capture non-linear effects of counterion distributions around highly charged objects in an essentially exact way. The scheme builds on the field-theoretic formulation of fluctuating counterions and uses the effective surface charge as a variational parameter. Since we use a Gaussian variational Hamiltonian, the resulting interactions between charged particles are of the Debye-Hückel type, but with the bare charges replaced by effective ones. This scheme is thus very close to the treatment by Kjellander et al.Kjellander , where a similar procedure was used within an integral-equation scheme. We show how in principle to systematically improve the leading term results by using perturbational variational methods. We tested our approach for the planar geometry, where comparison with the analytical solution of the non-linear Poisson-Boltzmann equation is possible, and then presented a detailed analysis of the cylindrical and spherical case. In the future, we will apply these methods to the interaction between charged bodies (including the experimentally relevant case of a charged cylinder adsorbing on a charge planefleck ), soft (fluctuating) charged polymers and membranes. We also plan to investigate higher-order perturbation variational terms which will allow to include effects at high coupling parameter $`\mathrm{\Xi }`$ where deviations from the mean-field like behavior are important.

## Appendix A Appendix: Charged planar wall

In this section we will treat the single charged wall immersed in a salt solution using a box model and using the Poisson Boltzmann equation.

### A.1 Box model

The mean-field Hamiltonian in terms of the ionic densities, Eq.(2.1), can be used as a starting point for a simplified variational treatment using a box model, where the true ionic density distribution is replaced by a step profile. The fugacity can be fixed by the ionic density in the bulk where we assume the electrostatic interactions to vanish. Minimization Eq.(2.1) with respect to the density of plus ions yields

$$\frac{\delta H}{\delta \rho _+}=\mathrm{ln}(\rho _+/\lambda )=0$$
(54)
and thus $`\rho _+=\lambda `$. We conclude that the fugacity is given by the ionic density in the bulk reservoir, i.e. $`\lambda =\rho _+^b`$.

As a box model we assume a negatively charged planar wall of surface charge density $`\sigma _s`$ which is neutralized by a uniformly charged box of width $`d`$. We assume that the wall charge is neutralized by an excess of positive ions, i.e. the density of negative ions stays constant. The electric field distribution thus is

$$\frac{eE}{k_BT}=4\pi \mathrm{}_B\sigma _s\frac{dz}{d}$$
(55)
and the electrostatic energy (per unit area and per $`k_BT`$) follows as

$`U`$ $`=`$ $`{\displaystyle \frac{\epsilon }{2k_BT}}{\displaystyle dzE^2}={\displaystyle \frac{1}{8\pi \mathrm{}_B}}{\displaystyle dz\left(\frac{eE}{k_BT}\right)^2}`$ (56)
$`=2\pi \mathrm{}_B\sigma _s^2d/3.`$
The excess entropy of the layer of condensed counterions is

$$S=d\left\{\left(\rho _+^b+\frac{\sigma }{qd}\right)\left(\mathrm{log}\left[1+\frac{\sigma _s}{qd\rho _+^b}\right]1\right)+\rho _+^b\right\}.$$
(57)
Minimization of the free energy $`F=US`$ leads for high salt concentrations (for $`\kappa \mu =\stackrel{~}{\kappa }1`$) to the optimal layer thickness

$$d^2\frac{3}{4\pi \mathrm{}_Bq^2\rho _+^b}=6\kappa ^2$$
(58)
which is proportional to the screening length $`\kappa ^1`$ (as expected) and for low salt concentrations (for $`\kappa \mu 1`$) to

$$d\frac{3}{2\pi \mathrm{}_Bq\sigma _s}=3\mu $$
(59)
and thus to a layer thickness which differs from the Gouy Chapman length $`\mu `$ only by a numerical factor. The surface free energies are given by

$$F\frac{2\sqrt{6}\sigma _s}{3q\kappa \mu }$$
(60)
for high salt and

$$F\frac{2\sigma _s}{q}\mathrm{log}\left(\frac{1}{\kappa \mu }\right).$$
(61)
for low salt.

### A.2 Poisson-Boltzmann equation

The Poisson-Boltzmann equation follows from minimizing the Hamiltonian Eq.(16) with respect to the potential distribution. Defining the real electrostatic potential $`\psi =ı\varphi `$ the saddle point equation can be written as

$$\psi ^{\prime \prime }(z)+2\delta (z)\stackrel{~}{\kappa }^2\mathrm{sinh}\psi =0.$$
(62)
The solution is

$$\psi (z)=2\mathrm{log}\left(\frac{1\gamma \mathrm{e}^{\kappa z}}{1+\gamma \mathrm{e}^{\kappa z}}\right)$$
(63)
where the constant $`\gamma `$ is determined by the boundary conditions. In the case when the charged wall is impenetrable to ions, the boundary condition is $`\psi ^{}(0)=2`$, which leads to the equation $`\gamma ^21=2\gamma \stackrel{~}{\kappa }`$. For low salt concentrations the constant follows as $`\gamma 1+\kappa `$ and the surface potential is given by

$$\psi _02\mathrm{ln}(1/\stackrel{~}{\kappa }).$$
(64)
The surface free energy in this limit is $`F=\psi _0\sigma _s/(2q)=(\sigma _s/q)\mathrm{ln}(1/\stackrel{~}{\kappa })`$ and thus differs from the box-model calculation only by a prefactor. For high salt concentrations the constant is $`\gamma =1/(2\stackrel{~}{\kappa })`$ and the surface potential is

$$\psi _02/\kappa .$$
(65)
The surface free energy is $`F=\sigma _s/(\stackrel{~}{\kappa }q)`$ and is also reproduced by the box-model calculation. In the case when the surface is penetrable to ions, as is used in our variational calculations, the surface potential in the low salt limit is not modified to leading order,

$$\psi _02\mathrm{ln}(1/\stackrel{~}{\kappa }),$$
(66)
while in the high salt limit one obtains

$$\psi _01/\kappa .$$
(67)
## Appendix B Appendix: Connection between the Gibbs Variational Principle and Perturbative Variation

Let us start by considering the free energy

$$=\mathrm{ln}𝒟\varphi \mathrm{e}^{}$$
(68)
where $``$ depends on the fluctuating field $`\varphi `$. We introduce a variational Hamiltonian $`_0(𝐱)`$ which depends on the fluctuating field $`\varphi `$ and on the vector $`𝐱`$, the components of which may act as variational parameters and/or as generating fields to calculate observables. Adding and subtracting the variational Hamiltonian to the action, the free energy can be written as

$$=\mathrm{ln}𝒟\varphi \mathrm{e}^{_0(𝐱)}\mathrm{ln}\left[\frac{𝒟\varphi \mathrm{e}^{_0(𝐱)+\zeta Y}}{𝒟\varphi \mathrm{e}^{_0(𝐱)}}\right]$$
(69)
where we defined

$$Y=_0(𝐱)$$
(70)
and introduced a dummy parameter $`\zeta `$ which is used to count orders in a perturbation series (and will eventually be set to unity). The free energy can now be expanded in cumulants as

$$_n=_0(𝐱)\underset{m=1}{\overset{n}{}}\frac{\zeta ^m}{m!}Y^m_0^C$$
(71)
where $`_0(𝐱)`$ denotes the free energy of the variational Hamiltonian

$$_0(𝐱)=\mathrm{ln}𝒟\varphi \mathrm{e}^{_0(𝐱)}$$
(72)
and $`\mathrm{}_0^C`$ denotes the cumulant average with respect to the variational Hamiltonian. Clearly, the exact free energy is approached in the infinite limit of the series, i.e. $`=_n\mathrm{}`$. The standard perturbation technique is often slowly converging for problems of interest. Variational techniques typically perform much better. The standard Gibbs-variational technique exploits the fact that the first order perturbation result $`_1`$ is always larger than the true free energy, $`_1`$. The approximate free energy is therefore given by minimizing $`_1`$ with respect to the set of variational parameters,

$$_{\mathrm{Gibbs}}=\underset{𝐱}{\mathrm{min}}_1(𝐱).$$
(73)
However, it is not easy to systematically improve this result (since the property of the perturbative free energy to be bound by the true free energy is lost at higher order), other than by choosing a more refined variational Hamiltonian, which is often not possible. The perturbational variation, which has been introduced by Edwards for the calculation of the scaling behavior of self-avoiding polymersEdwards , consists in requiring the perturbative free energy of any order to be equal to the zeroth order free energy,

$$_0(𝐱)=_n(𝐱)$$
(74)
which is equivalent to the condition

$$\underset{m=1}{\overset{n}{}}\frac{\zeta ^m}{m!}Y^m_0^C=0.$$
(75)
This equation has to be solved by fixing the set of variational parameters $`𝐱`$. To do this, it is convenient to expand the parameters in a power series $`𝐱=_{i=1}^n\zeta ^{i1}𝐱^{(i)}`$ and to solve Eq.(75) term by term. (If the variational parameter is not expanded in a series but Eq.(75) is solved directly, one recovers the naive perturbation result. The hope is therefore that the expansion of the variational parameters improves the convergence behavior of the series for the free energy.) One notes that the prediction for the free energy from the Gibbs variational approach, Eq.(73), and the perturbational variation, Eq.(74), are not equal, not even at the leading order. It is a priori not clear which approximation closer to the exact free energy.

Often one is not interested in the free energy itself but in a certain observable. We define an arbitrary observable $`X_j`$ as being conjugate to the field $`x_j`$ according to

$$X_j=\frac{_0(𝐱)}{x_j}.$$
(76)
The first few terms of the expansion of the expectation value of $`X_j`$ in powers of $`\zeta `$ are

$`X_j=X_j_0+\zeta \left[X_jY_0X_j_0Y_0\right]`$
$`+`$ $`{\displaystyle \frac{\zeta ^2}{2}}\left[X_jY^2_02X_jY_0Y_0X_j_0Y^2_0+2X_j_0Y_0^2\right].`$
As before, the perturbative variation consists of demanding that the corrections to the leading contribution of the expectation value vanish, i.e.

$$X_j=X_j_0,$$
(78)
which can be solved by adjusting one or more of the variational parameters $`𝐱`$.

A quite instructive connection exists between the Gibbs variational method and the perturbative variation: In the case when the Gibbs variational parameter is the field $`x_j`$ conjugate to the observable $`X_j`$, the Gibbs variational equation and the first order perturbative variation Eq.(78) are identical. This means that the higher order perturbative variational method on the expectation value can be viewed as a generalization of the Gibbs variational technique. One might hope to approximate the true free energy better, though the true free energy ceases to be a strict lower bound for the perturbational free energy expression.

Some comments are in order:

i) In principle, the calculated observable in equation (B) does not need to be conjugate to the variational parameter. In any case, it is clear that in principle the observable can be calculated to any wanted precision by taking enough terms in the perturbation-variational treatment.

ii) Several expectation values can be calculated, for which an equal number of variational parameters are needed to solve the uniform expansion condition.

v) It is not clear that equation (78) always has a solution. In that case one could try to find the best estimate to the solution.