text:
# Counterion Penetration and Effective Electrostatic Interactions in Solutions of Polyelectrolyte Stars and Microgels

## I Introduction

Polyelectrolytes (PEs) are ionizable polymers that dissolve in a polar solvent, such as water, by dissociating into polyvalent macroions and small oppositely-charged counterions PE1 . Electrostatic interactions between macroions, mediated by surrounding microions (counterions and salt ions), contribute to the unique macroscopic properties of PE solutions, which are the basis of many industrial applications involving polymer-water systems PE2 . Common synthetic examples of PEs are polyacrylic acid, used in gels and rheology modifiers, and polystyrene sulfonate, a component of reverse osmosis membranes. Naturally occurring examples are biopolymers, such as DNA, proteins, and starches. Colloidal in size, PEs are also routinely added as flocculants and stabilizers to colloidal suspensions, such as foods and water-based paints Evans ; Hunter . Depending on PE concentration, adsorption or grafting of PE chains onto surfaces of colloidal particles can either induce flocculation, by bridging particles, or impart electrosteric stabilization.

Conformations of PE macroions and electrostatic interactions between macroions are strongly influenced by the distribution of microions. If dispersed in solution, microions act to screen the bare Coulomb interactions between ionized monomers. If condensed on the macroion chains, microions may reduce the macroion charge Schiessel . Linear PE chains whose monomers are sufficiently weakly interacting – either because weakly charged or because of strong microion screening or condensation – may form random-walk coils with roughly spherical conformation. With increasing charge and screening length, linear chains stretch into non-spherical conformations because of electrostatic repulsion between ionized monomers deGennes ; Khokhlov . The extent of elongation depends on the chain charge density, salt concentration, and solvent quality. Highly charged chains in good solvents (e.g., DNA in water) often form stiff rod-like macroions, whose effective interactions and complex phase behavior (such as bundling) have been widely studied rods1 ; rods2 . In poor solvents, sufficiently highly charged chains may form necklaces of compact globules joined by narrow threads, as predicted by theory Rubinstein  and confirmed by simulation Kremer .

Although many common PEs are linear, other topologies can be readily synthesized. Examples are stars, microgels, micelles, and brushes. Star polymers Likos1  consist of chains chemically grafted or adsorbed to a common microscopic core. Microgels are mesoscopic polymer networks, synthesized by polymerization in microemulsion microgels . Micelles are formed by association of charged diblock (amphiphilic) copolymers micelles . Brushes are formed by grafting PE chains onto a mesoscopic solid core brushes . Solutions of spherical stars, microgels, micelles, and brushes can be regarded as colloidal suspensions of soft macroions that are permeable to microions.

Electrostatic interactions in charged colloids have received much attention in recent years reviews , motivated largely by anomalous phase behavior that is unexplained by the classic Derjaguin-Landau-Verwey-Overbeek (DLVO) theory DLVO . Most studies have been restricted, however, to hard, impermeable macroions. The objectives of this paper are first, to explore implications of microion penetration for screening of effective electrostatic interactions between spherical macroions, and second, to lay a foundation for future studies of thermodynamic phase behavior of PE solutions. Our approach is based on a recently-proposed theory of effective interactions in charged colloids Denton1 , which we adapt here from hard to penetrable macroions and apply to spherical star-branched and microgel macroions.

The remainder of the paper is organized as follows. Section II describes the assumed model of PE solutions. Section III reviews the theoretical approach, based on second-order perturbation (linear response) theory. Sections IV and V present analytical and numerical results for counterion profiles and effective interactions in bulk solutions of star and microgel macroions. Finally, Sec. VI closes with a summary and conclusions.

## II Model

Adapting the primitive model of ionic liquids HM , the model system comprises $`N_m`$ spherical macroions, of radius $`a`$ (diameter $`\sigma =2a`$) and charge $`Ze`$, and $`N_c`$ point counterions of charge $`ze`$ dispersed in an electrolyte solvent in a volume $`V`$ at temperature $`T`$. Assuming, for simplicity, a symmetric electrolyte and equal salt and counterion valences, the electrolyte contains $`N_s`$ point salt ions of charge $`ze`$ and $`N_s`$ of charge $`ze`$. The microions thus number $`N_+=N_c+N_s`$ positive and $`N_{}=N_s`$ negative, for a total of $`N_\mu =N_c+2N_s`$. Global charge neutrality in a bulk solution constrains average macroion and counterion number densities, $`n_m=N_m/V`$ and $`n_c=N_c/V`$, via $`Zn_m=zn_c`$. The polar solvent is treated as a continuum, characterized by dielectric constant $`ϵ`$ that acts only to reduce Coulomb interactions between ions.

The local number density profiles of counterions, $`\rho _c(r)`$, and of macroion monomers, $`\rho _{\mathrm{mon}}(r)`$, are modelled as spherically symmetric, continuous distributions. Spherical symmetry is a reasonable approximation, considering that equilibrium averaging over macroion orientations tends to smear out any anisotropy. Furthermore, discreteness of the charge distributions can be ignored if we restrict consideration to length scales exceeding the scale of discreteness.

In general, the counterions are distributed over three regions: (1) the immediate vicinity of the PE chains making up the macroions, (2) the region inside of the macroions but away from the chains, and (3) the region outside of the macroions. Counterions in the first two regions are trapped by the macroions, while those in the third region are free. Within the first region, the counterions may be either condensed on a chain or free to move along a tube surrounding a chain. These chain-localized counterions, whether condensed or mobile, tend to distribute uniformly along the chains to favor local charge neutrality. In our model, counterions in region (1) simply renormalize the effective macroion valence $`Z`$.

The detailed form of the monomer density profile depends on the macroion conformation. For star-branched macroions, Coulomb repulsion between charged monomers tends to stiffen and radially stretch the chains into a porcupine conformation Pincus . We assume the ideal case of fully stretched chains and model the monomer density profile by $`\rho _{\mathrm{mon}}(ra)=Z/(4\pi ar^2)`$, where $`r`$ is the radial distance from the star’s center. For microgel macroions, the dense network of chains is well approximated by a uniform monomer distribution, and is modelled here by $`\rho _{\mathrm{mon}}(ra)=3Z/(4\pi a^3)`$. This distribution may also approximate a weakly-charged linear PE chain with a spherical random-coil conformation, although a Gaussian distribution may then be more accurate. For both the star and microgel models, the monomer density profile is cut off sharply at the macroion surface: $`\rho _{\mathrm{mon}}(r>a)=0`$.

## III Theory

For the model PE solutions described above, the theoretical challenge is to predict the distributions of microions inside and outside of the macroions and the effective interactions between macroions. Following the same general strategy as applied previously to charged colloids Denton1 ; Silbert , we reduce the multi-component mixture to an equivalent one-component system governed by effective interactions, which are approximated via perturbation theory. For clarity of presentation, we initially ignore salt ions. The Hamiltonian then decomposes into three terms:

$$H=H_m(\{𝐑\})+H_c(\{𝐫\})+H_{mc}(\{𝐑\},\{𝐫\}),$$
(1)
where $`\{𝐑\}`$ and $`\{𝐫\}`$ denote collective coordinates of macroion centers and counterions, respectively. The first term,

$$H_m=K_m+\frac{1}{2}\underset{ij=1}{\overset{N_m}{}}v_{mm}(|𝐑_i𝐑_j|),$$
(2)
is the bare Hamiltonian for macroions with kinetic energy $`K_m`$ that interact via the bare pair potential $`v_{mm}(r)`$ at center-center separation $`r`$. The form of $`v_{mm}(r)`$ depends on the macroion conformation and is specified in the Appendix. The second term in Eq. (1),

$$H_c=K_c+\frac{1}{2}\underset{ij=1}{\overset{N_c}{}}v_{cc}(|𝐫_i𝐫_j|),$$
(3)
is the Hamiltonian for counterions with kinetic energy $`K_c`$ that interact via the Coulomb pair potential $`v_{cc}(r)=z^2e^2/ϵr`$. The third term in Eq. (1),

$$H_{mc}=\underset{i=1}{\overset{N_m}{}}\underset{j=1}{\overset{N_c}{}}v_{mc}(|𝐑_i𝐫_j|),$$
(4)
is the macroion-counterion interaction. For spherical macroions,

$$v_{mc}(r)=\{\begin{array}{cc}\frac{Zze^2}{ϵr},\hfill & r>a\hfill \\ v_<(r),\hfill & ra,\hfill \end{array}$$
(5)
where the interaction inside a macroion, $`v_<(r)`$, depends on the macroion conformation and is specified in Sec. IV. For later reference, we note that Eq. (4) also may be expressed in the form

$$H_{mc}=d𝐑\rho _m(𝐑)d𝐫\rho _c(𝐫)v_{mc}(|𝐑𝐫|),$$
(6)
where $`\rho _m(𝐑)=_{j=1}^{N_m}\delta (𝐑𝐑_j)`$ and $`\rho _c(𝐫)=_{j=1}^{N_c}\delta (𝐫𝐫_j)`$ are the macroion and counterion number density operators, respectively.

The mixture of macroions and counterions is formally reduced to an equivalent one-component system by tracing over counterion coordinates. Denoting counterion and macroion (classical) traces by $`_c`$ and $`_m`$, respectively, the canonical partition function can be expressed as

$$𝒵=\mathrm{exp}(\beta H)_c_m=\mathrm{exp}(\beta H_{\mathrm{eff}})_m,$$
(7)
where $`H_{\mathrm{eff}}=H_m+F_c`$ is the effective one-component Hamiltonian, $`\beta =1/k_BT`$, and

$$F_c=k_BT\mathrm{ln}\mathrm{exp}[\beta (H_c+H_{mc})]_c$$
(8)
is the free energy of a nonuniform gas of counterions in the presence of the macroions.

At this stage, approximations are necessary for the counterion free energy. It is first convenient to convert the counterion Hamiltonian to the Hamiltonian of a classical one-component plasma (OCP) of counterions by adding to $`H_c`$, and subtracting from $`H_{mc}`$, the energy of a uniform compensating negative background note1 , $`E_b=N_cn_c\widehat{v}_{cc}(0)/2`$, where $`\widehat{v}_{cc}(0)`$ is the $`k0`$ limit of the Fourier transform of $`v_{cc}(r)`$. Now regarding the macroions as an “external” potential for the OCP, we invoke perturbation theory HM ; Denton1 ; Silbert  and write

$$F_c=F_{\mathrm{OCP}}+_0^1d\lambda H_{mc}^{}_\lambda ,$$
(9)
where $`F_{\mathrm{OCP}}=k_BT\mathrm{ln}\mathrm{exp}[\beta (H_c+E_b)]_c`$ is the OCP free energy, the $`\lambda `$-integral charges the macroions, $`H_{mc}^{}=H_{mc}E_b`$ represents the perturbing potential of the macroions acting on the counterions, and $`H_{mc}^{}_\lambda `$ is the mean value of this potential in a solution of macroions charged to a fraction $`\lambda `$ of their full charge. Further progress is facilitated by expressing $`H_{mc}`$ \[Eq. (6)\] in terms of Fourier components:

$`H_{mc}_\lambda ={\displaystyle \frac{1}{V}}{\displaystyle \underset{𝐤0}{}}\widehat{v}_{mc}(k)\widehat{\rho }_m(𝐤)\widehat{\rho }_c(𝐤)_\lambda `$
$`+`$ $`{\displaystyle \frac{1}{V}}\underset{k0}{lim}\left[\widehat{v}_{mc}(k)\widehat{\rho }_m(𝐤)\widehat{\rho }_c(𝐤)_\lambda \right],`$ (10)
where $`\widehat{v}_{mc}(k)`$ is the Fourier transform of Eq. (5) and where $`\widehat{\rho }_m(𝐤)=_{j=1}^{N_m}\mathrm{exp}(i𝐤𝐑_j)`$ and $`\widehat{\rho }_c(𝐤)=_{j=1}^{N_c}\mathrm{exp}(i𝐤𝐫_j)`$ are Fourier components of the macroion and counterion densities.

In first-order perturbation theory, the response of the counterion plasma to the macroions is ignored. Here we apply second-order perturbation (linear response) theory, in which the counterions are assumed to respond linearly to the macroion external potential:

$$\widehat{\rho }_c(𝐤)=\chi (k)\widehat{v}_{mc}(k)\widehat{\rho }_m(𝐤),k0,$$
(11)
where $`\chi (k)`$ is the linear response function of the OCP. Note that the $`k0`$ limit here, and in Eq. (10), must be treated separately, since the average counterion density, $`n_c=\widehat{\rho }_c(0)`$, does not respond to the macroion charge, but rather is fixed by the constraint of global charge neutrality.

Upon combining Eqs. (9)-(11), the effective Hamiltonian can be recast in the form of the Hamiltonian of a pairwise-interacting system:

$$H_{\mathrm{eff}}=K_m+K_c+\frac{1}{2}\underset{ij=1}{\overset{N_m}{}}v_{\mathrm{eff}}(|𝐑_i𝐑_j|)+E_0,$$
(12)
where $`v_{\mathrm{eff}}(r)=v_{mm}(r)+v_{\mathrm{ind}}(r)`$ is an effective macroion pair interaction that combines the bare macroion interaction with a microion-induced interaction

$$\widehat{v}_{\mathrm{ind}}(k)=\chi (k)\left[\widehat{v}_{mc}(k)\right]^2.$$
(13)
The final term in Eq. (12) is the volume energy, formally given by

$`E_0=F_{\mathrm{OCP}}+{\displaystyle \frac{N_m}{2}}\underset{r0}{lim}v_{\mathrm{ind}}(r)`$
$`+`$ $`N_m\underset{k0}{lim}\left[{\displaystyle \frac{n_m}{2}}\widehat{v}_{\mathrm{ind}}(k)+n_c\widehat{v}_{mc}(k)+{\displaystyle \frac{Zn_c}{2z}}\widehat{v}_{cc}(k)\right],`$
which is a natural by-product of the one-component reduction. Although independent of the macroion coordinates, the volume energy depends on the average macroion density and thus can influence thermodynamics.

The linear response function is proportional to the corresponding static structure factor, $`S(k)`$, which may be obtained from liquid-state theory HM . In practice, the OCP is weakly correlated, with coupling parameter $`\mathrm{\Gamma }=\lambda _B/a_c1`$, where $`\lambda _B=\beta e^2/ϵ`$ is the Bjerrum length and $`a_c=(3/4\pi n_c)^{1/3}`$ is the counterion sphere radius. For example, for macroions of diameter $`\sigma =100`$ nm, valence $`Z=100`$, and volume fraction $`\eta =(\pi /6)n_m\sigma ^3=0.01`$, in water at room temperature ($`\lambda _B=0.714`$ nm), we find $`\mathrm{\Gamma }0.014`$. As for charged colloids Denton1 ; Silbert , we adopt the random phase approximation (RPA), which is accurate for weakly-coupled plasmas. The RPA equates the two-particle direct correlation function of the OCP to its exact asymptotic limit: $`c^{(2)}(r)=\beta v_{cc}(r)`$. Using the Ornstein-Zernike relation, $`S(k)=1/[1n_c\widehat{c}^{(2)}(k)]`$, the linear response function then takes the analytical form

$$\chi (k)=\beta n_cS(k)=\frac{\beta n_c}{(1+\kappa ^2/k^2)},$$
(15)
where $`\kappa =\sqrt{4\pi n_cz^2\lambda _B}`$ is the inverse Debye screening length. Note that since permeable macroions do not exclude counterions from their interiors, the excluded-volume corrections required for hard colloidal macroions Denton1  are not relevant here. With $`\chi (k)`$ specified, the counterion density can be explicitly determined from Eqs. (5) and (11) for a given macroion distribution (see Sec. IV). Finally, salt is easily introduced via additional microion response functions. In the process, the pair interaction and volume energy are unchanged, except for a redefinition of the screening constant as $`\kappa =\sqrt{4\pi (n_c+2n_s)z^2\lambda _B}`$, where $`n_s`$ is the average number density of salt ion pairs.

It is worth noting the formal equivalence of the present theory to linearized Poisson-Boltzmann (DLVO) theory. Both are mean-field theories in the sense that they ignore fluctuations in microion distributions. An advantage of linear response theory, however, is that it encompasses the volume energy, which can be important for describing phase behavior Denton1 ; Silbert ; vRH ; Graf ; vRDH ; Warren . Moreover, response theory can be straightforwardly generalized to incorporate nonlinear response, which entails both many-body effective interactions and corrections to the pair potential and volume energy Denton2 . In contrast, nonlinear Poisson-Boltzmann theory is practical only for the simple boundary conditions afforded by cell models. For simplicity, higher-order nonlinear effects are here ignored.

Equations (11)-(LABEL:E0) constitute the main formal expressions from linear response theory. Explicit calculations require specifying the counterion-macroion interaction $`v_<(r)`$ in Eq. (5) for specific macroion models. Below, we apply the theory to obtain analytical and numerical results for counterion profiles and effective interactions in bulk solutions of spherical star-branched and microgel macroions.

## IV Analytical Results

### IV.1 Star Macroions

For our idealized model of a star-branched macroion with $`1/r^2`$ monomer density profile, Gauss’s law gives the electric field as

$$E(r)=\{\begin{array}{cc}\frac{Ze}{ϵr^2},\hfill & r>a\hfill \\ \frac{Ze}{ϵar},\hfill & ra.\hfill \end{array}$$
(16)
Integration over $`r`$ yields the electrostatic potential energy between a star and a counterion:

$$v_{mc}(r)=\{\begin{array}{cc}\frac{Zze^2}{ϵr},\hfill & r>a\hfill \\ \frac{Zze^2}{ϵa}\left[1\mathrm{ln}\left(\frac{r}{a}\right)\right],\hfill & ra,\hfill \end{array}$$
(17)
whose Fourier transform is

$$\widehat{v}_{mc}(k)=\frac{4\pi Zze^2}{ϵk^3a}\mathrm{sinc}(ka),$$
(18)
with $`\mathrm{sinc}(x)_0^xdu\mathrm{sin}(u)/u`$. We can now calculate the counterion number density around a single macroion in the dilute limit, where $`\widehat{\rho }_m(𝐤)=1`$. From Eqs. (11), (15), and (18), the Fourier component of the counterion density profile is given by

$$\widehat{\rho }_c(k)=\frac{Z}{z}\frac{\kappa ^2}{ka(k^2+\kappa ^2)}\mathrm{sinc}(ka),$$
(19)
whose real-space form is

$$\rho _c(r)=\frac{Z}{z}\frac{\kappa }{4\pi ar}\mathrm{sinhc}(\kappa a)e^{\kappa r},r>a$$
(20)
$`\rho _c(r)`$ $`=`$ $`{\displaystyle \frac{Z}{z}}{\displaystyle \frac{\kappa }{8\pi ar}}\left[\mathrm{Ec}(\kappa a,\kappa r)+2\mathrm{sinhc}(\kappa a)\right]e^{\kappa r}`$
$``$ $`\mathrm{Ec}(\kappa a,\kappa r)e^{\kappa r},ra,`$ (21)
where

$`\mathrm{sinhc}(x)`$ $``$ $`{\displaystyle _0^x}du{\displaystyle \frac{\mathrm{sinh}(u)}{u}}`$
$`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{x^{2n+1}}{(2n+1)(2n+1)!}}`$ (22)
and

$$\mathrm{Ec}(x_1,x_2)_{x_1}^{x_2}du\frac{e^u}{u}=\mathrm{ln}\left(\frac{x_2}{x_1}\right)+\underset{n=1}{\overset{\mathrm{}}{}}\frac{x_2^nx_1^n}{nn!},$$
(23)
which can be efficiently computed from the first few terms of the rapidly converging series expansions. Approaching the macroion center, the counterion density profile varies more gradually than the $`1/r^2`$ macroion monomer density profile, diverging logarithmically, according to

$$\underset{r0}{lim}\rho _c(r)=\frac{Z\kappa ^2}{4\pi a}\left[1\mathrm{ln}\left(\frac{r}{a}\right)\right].$$
(24)
Integrating Eq. (21) over the spherical volume of the macroion yields the fraction of counterions inside a star:

$`f_{\mathrm{in}}`$ $`=`$ $`{\displaystyle \frac{z}{Z}}4\pi {\displaystyle _0^a}drr^2\rho _c(r)`$
$`=`$ $`1\left(1+{\displaystyle \frac{1}{\kappa a}}\right)e^{\kappa a}\mathrm{sinhc}(\kappa a).`$ (25)
Note the clear predictions that (1) the counterion distribution is determined entirely by $`\kappa a`$, or the dimensionless ratio of the macroion radius and the Debye screening length, and (2) the fraction of counterions inside increases monotonically with $`\kappa a`$. Thus, for fixed macroion radius, $`f_{\mathrm{in}}`$ increases with increasing macroion valence and concentration. This result is physically sensible: the shorter the screening length, the shorter the range of the counterion response, and thus the tighter the localization of counterions around the macroion centers.

From Eqs. (13) and (18), the induced electrostatic pair interaction is given by

$$\widehat{v}_{\mathrm{ind}}(k)=\frac{4\pi Z^2e^2}{ϵ}\frac{\kappa ^2}{k^4a^2(k^2+\kappa ^2)}\mathrm{sinc}^2(ka).$$
(26)
Fourier transforming, we obtain

$`v_{\mathrm{ind}}(r)`$ $`=`$ $`{\displaystyle \frac{16\pi ^2Z^2e^2\kappa ^2a^2}{ϵr}}{\displaystyle _0^{\mathrm{}}}dx{\displaystyle \frac{\mathrm{sin}(xr/a)}{x^3(x^2+\kappa ^2a^2)}}`$
$`\times `$ $`\mathrm{sinc}^2x.`$ (27)
For nonoverlapping stars, Eq. (27) can be reduced to the analytical form

$$v_{\mathrm{ind}}(r>2a)=\frac{Z^2e^2}{ϵr}+\frac{Z^2e^2}{ϵ}\left[\frac{\mathrm{sinhc}(\kappa a)}{\kappa a}\right]^2\frac{e^{\kappa r}}{r}.$$
(28)
Since nonoverlapping macroions interact via a bare Coulomb potential, $`v_{mm}(r)=Z^2e^2/ϵr`$, the effective pair interaction for this case is

$$v_{\mathrm{eff}}(r>2a)=\frac{Z^2e^2}{ϵ}\left[\frac{\mathrm{sinhc}(\kappa a)}{\kappa a}\right]^2\frac{e^{\kappa r}}{r}.$$
(29)
Thus, at the level of linear response, nonoverlapping star macroions interact via an effective Yukawa (screened-Coulomb) pair potential. The screening constant, $`\kappa `$, in the potential depends on the total density of microions – inside and outside of the macroions – since all microions respond to the macroion charge. Note that the potential has the same $`r`$-dependence as the DLVO potential for hard colloidal macroions Hunter ; DLVO ,

$$v_{\mathrm{DLVO}}(r)=\frac{Z^2e^2}{ϵ}\left[\frac{\mathrm{exp}(\kappa a)}{1+\kappa a}\right]^2\frac{e^{\kappa r}}{r},r>2a,$$
(30)
differing only in the macroion-size-dependent amplitude. For overlapping stars, the bare macroion interaction is somewhat more complex and is relegated to the Appendix.

Finally, from Eqs. (LABEL:E0), (18), (26), and (27), the volume energy is obtained as

$$E_0=F_{\mathrm{OCP}}N_m\frac{8\pi ^2Z^2e^2\kappa ^2a}{ϵ}_0^{\mathrm{}}dx\frac{\mathrm{sinc}^2x}{x^2(x^2+\kappa ^2a^2)}.$$
(31)
For weakly-coupled microion plasmas, the OCP free energy may be approximated by its ideal-gas limit:

$$F_{\mathrm{OCP}}=N_+[\mathrm{ln}(n_+\mathrm{\Lambda }^3)1]+N_{}[\mathrm{ln}(n_{}\mathrm{\Lambda }^3)1],$$
(32)
where $`\mathrm{\Lambda }`$ is the thermal wavelength. Note that if $`Z`$ is allowed to vary (with counterion condensation), then the volume energy per macroion must be augmented by the self (Hartree) energy of a macroion: $`U_H=Z^2e^2/ϵa`$. The first term in Eq. (31), long recognized as important for phase behavior PE1 , represents the entropy of free counterions; the second term accounts for the cohesive electrostatic energy of microion-macroion interactions. The volume energy, analogous to its counterpart for charged colloids Denton1 ; vRDH ; Graf ; vRDH , depends on the average macroion concentration and thus has the potential to influence phase behavior and other thermodynamic properties. Equations (20), (21), (25), (29), and (31) are the main analytical results for star macroions.

It is important to emphasize that the present approach, while including the entropy of the counterions, neglects the configurational entropy of the macroions by assuming rigid (fully stretched) PE chains. Recently, Jusufi et alJusufi  modelled pair interactions between PE stars by both molecular dynamics simulation and a variational free energy that incorporates chain flexibility. An important conclusion of their study is that pair interactions are dominated by counterion entropy. Our approach is complementary: while the macroion model neglects chain flexibility, which is reasonable at least for nonoverlapping stars, the linear response theory refines somewhat the modelling of the counterion distribution.

### IV.2 Microgel Macroions

For our model of microgel macroions, we apply exactly the same procedure as in Sec. IV.1. The electric field of a uniformly-charged sphere is

$$E(r)=\{\begin{array}{cc}\frac{Ze}{ϵr^2},\hfill & r>a\hfill \\ \frac{Zer}{ϵa^3},\hfill & ra,\hfill \end{array}$$
(33)
which integrates to give the macroion-counterion interaction,

$$v_{mc}(r)=\{\begin{array}{cc}\frac{Zze^2}{ϵr},\hfill & r>a\hfill \\ \frac{Zze^2}{2ϵa}\left(3\frac{r^2}{a^2}\right),\hfill & ra.\hfill \end{array}$$
(34)
Equation (34) Fourier transforms to

$$\widehat{v}_{mc}(k)=\frac{12\pi Zze^2}{ϵk^4a^2}\left[\mathrm{cos}(ka)\frac{\mathrm{sin}(ka)}{ka}\right],$$
(35)
which, when substituted into Eq. (11), yields the Fourier transform of the counterion density profile around a single macroion,

$$\widehat{\rho }_c(k)=\frac{Z}{z}\frac{3\kappa ^2}{k^2a^2(k^2+\kappa ^2)}\left[\mathrm{cos}(ka)\frac{\mathrm{sin}(ka)}{ka}\right].$$
(36)
The Fourier transform of Eq. (36) gives the real-space counterion density profile,

$`\rho _c(r)`$ $`=`$ $`{\displaystyle \frac{Z}{z}}{\displaystyle \frac{3}{4\pi a^2r}}\left[\mathrm{cosh}(\kappa a){\displaystyle \frac{\mathrm{sinh}(\kappa a)}{\kappa }}a\right]e^{\kappa r},`$ (37)
$`r>a`$
$`\rho _c(r)`$ $`=`$ $`{\displaystyle \frac{Z}{z}}{\displaystyle \frac{3}{4\pi a^2r}}\left[{\displaystyle \frac{r}{a}}\left(1+{\displaystyle \frac{1}{\kappa a}}\right)e^{\kappa a}\mathrm{sinh}(\kappa r)\right],`$ (38)
$`ra,`$
which approaches a constant as $`r0`$:

$$\rho _c(r=0)=\frac{Z}{z}\frac{3}{4\pi a^3}\left[1(1+\kappa a)e^{\kappa a}\right].$$
(39)
Integration of Eq. (38) yields an analytical result for the internal fraction of counterions:

$$f_{\mathrm{in}}=1\frac{3}{\kappa a}\left(1+\frac{1}{\kappa a}\right)e^{\kappa a}\left[\mathrm{cosh}(\kappa a)\frac{\mathrm{sinh}(\kappa a)}{\kappa a}\right].$$
(40)
Again the theory predicts a counterion distribution depending only on the ratio of macroion radius to screening length and an internal counterion fraction that increases monotonically with this ratio. This prediction may be compared with that of Oosawa’s “two-phase” approximation PE1 , which assumes uniform (but differing) counterion concentrations inside and outside of the macroions. According to the latter approach, for spherical macroions with volume fraction $`\eta `$, the condition for equilibrium between free and bound counterions, in the absence of salt ions, is

$$\mathrm{ln}\left(\frac{f_{\mathrm{in}}}{1f_{\mathrm{in}}}\right)=\mathrm{ln}\left(\frac{\eta }{1\eta }\right)+Z\frac{\lambda _B}{a}(1f_{\mathrm{in}})(1\eta ^{1/3}),$$
(41)
which must be solved numerically for $`f_{\mathrm{in}}`$. The predictions of Eqs. (40) and (41) are compared below in Sec. V.

Next, substituting Eq. (35) into Eq. (13), the induced pair interaction is

$`\widehat{v}_{\mathrm{ind}}(k)`$ $`=`$ $`{\displaystyle \frac{36\pi Z^2e^2}{ϵ}}{\displaystyle \frac{\kappa ^2}{k^6a^4(k^2+\kappa ^2)}}`$
$`\times `$ $`\left[\mathrm{cos}(ka){\displaystyle \frac{\mathrm{sin}(ka)}{ka}}\right]^2.`$ (42)
For nonoverlapping macroions, the Fourier transform of Eq. (42) is straightforward to evaluate and yields an effective pair interaction

$`v_{\mathrm{eff}}(r)`$ $`=`$ $`{\displaystyle \frac{Z^2e^2}{ϵ}}{\displaystyle \frac{9}{\kappa ^4a^4}}\left[\mathrm{cosh}(\kappa a){\displaystyle \frac{\mathrm{sinh}(\kappa a)}{\kappa a}}\right]^2`$
$`\times `$ $`{\displaystyle \frac{e^{\kappa r}}{r}},r>2a.`$ (43)
As for star macroions, a Yukawa form is predicted, but with a different amplitude. The case of overlapping macroions is left to the Appendix.

Finally, from Eqs. (LABEL:E0), (35), and (42), the volume energy is obtained as

$`E_0`$ $`=`$ $`F_{\mathrm{OCP}}N_m{\displaystyle \frac{3Z^2e^2}{ϵa}}\{{\displaystyle \frac{1}{5}}{\displaystyle \frac{1}{2\kappa ^2a^2}}`$
$`+`$ $`{\displaystyle \frac{3}{4\kappa ^3a^3}}[1{\displaystyle \frac{1}{\kappa ^2a^2}}+(1+{\displaystyle \frac{2}{\kappa a}}+{\displaystyle \frac{1}{\kappa ^2a^2}})e^{2\kappa a}]\}`$
$``$ $`(N_+N_{}){\displaystyle \frac{k_BT\kappa ^2a^2}{2}}.`$ (44)
Equations (37), (38), (40), (43), and (44) are the main analytical results for microgel macroions.

## V Numerical Results

The theory developed above can be applied to solutions of arbitrary ionic strength, under the assumption that the macroion PE chains remain stretched. In order to highlight the role of the counterions, we present numerical results for salt-free solutions. Within the model considered, the effect of salt is merely to increase the Debye screening constant. Furthermore, we consider the case of monovalent counterions ($`z=1`$) in aqueous solutions at room temperature ($`\lambda _B=0.714`$ nm). Figure 1 illustrates the form of the counterion number density profiles \[Eqs. (20), (20), and (37), (38)\] inside and outside of a macroion. Inside a star macroion the counterion density diverges logarithmically towards the center, while inside a microgel macroion $`\rho _c(r)`$ remains finite. Evidently, counterions penetrate stars more easily than they do microgels. This property is also reflected in the internal counterion fractions \[Eqs. (25) and (40)\], functions of $`\kappa a`$ only, which are shown in Fig. 2a. In Fig. 2b, we compare predictions of linear response theory \[Eq. (40)\] with those of Oosawa’s two-phase approximation \[Eq. (41)\] for uniformly-charged spherical (microgel) macroions PE1 . Both approaches qualitatively predict an increase in the fraction of bound counterions with increasing macroion concentration. However, linear response theory predicts a considerably more gradual accumulation of bound counterions than does the two-phase approximation.

Counterion penetration strongly influences screening of bare macroion interactions. The effective pair potentials, $`v_{\mathrm{eff}}(r)`$, and corresponding forces, $`F(r)=\mathrm{d}v_{\mathrm{eff}}(r)/\mathrm{d}r`$, are shown in Figs. 3 and 4, respectively. Beyond overlap the effective interaction has Yukawa form, with the amplitude depending on the type of macroion. Figure 5 compares the variation of the macroion-size-dependent amplitude of $`v(r>2a)`$ with the Debye screening constant for the two permeable macroions and for hard macroions. Evidently, the greater the permeability of the macroions to counterions, the weaker the amplitude of long-range repulsion. For overlapping macroions, the bare charge distribution combined with counterion penetration leads to softly repulsive interactions. Note that the interactions are bounded: they do not diverge as the macroions approach complete overlap. It must be emphasized that the effective interactions presented in Figs. 3 and 4 arise physically from electrostatic repulsion and counterion screening, but do not include steric interactions due to compression of overlapping chains Jusufi .

## VI Conclusions

To summarize, we have applied second-order perturbation (linear response) theory to model solutions of spherical polyelectrolyte star-branched and microgel macroions. The theory predicts the counterion density profiles inside and outside of the macroions, effective interactions between pairs of macroions, and a one-body, density-dependent, volume energy that contributes to the total free energy of the system. The main conclusions are: (1) Counterions penetrate stars more easily than they do microgels. (2) Inside a star macroion, the density profile of mobile counterions varies more gradually than the macroion monomer density profile, diverging logarithmically toward the center. (3) The fraction of counterions trapped inside a macroion depends only on the ratio of macroion radius to Debye screening length and increases monotonically with this ratio. (4) Counterion screening significantly weakens the bare electrostatic pair interactions, which remain bounded up to complete overlap of macroions. (5) The effective pair interactions are softly repulsive for overlapping macroions and Yukawa in form for separated macroions, with amplitudes depending on the type of macroion.

It is important to point out some limitations of the theory. First, the linear response approximation limits applicability of the theory to dilute solutions of weakly charged macroions. The quantitative range of validity depends on the relative magnitudes of nonlinear corrections, including three-body and high-order interactions, in the perturbation expansion. The same techniques that have been used to analyze nonlinear response in charged colloids Denton2  can be applied to polyelectrolytes. Second, the mean-field approach taken here ignores fluctuations in the counterion distribution, which may be especially relevant for short-range interactions and multivalent counterions. Third, the neglect of chain flexibility restricts the theory to nonoverlapping macroions. This restriction may be reasonable for dilute solutions of sparsely separated macroions. However, for a sufficient concentration of macroions in a good solvent, chain elasticity and entropy must play a role. A unification of linear response theory and the variational theory of ref. Jusufi  may then prove fruitful.

In principle, the predicted counterion profiles could be probed experimentally, e.g., by neutron scattering, using isotopic labelling to contrast the PE chains and counterions. The macroion-macroion interactions may be less accessible to experiment. Conceivably, the solvent quality might be tuned to minimize the second virial coefficient between neutral monomers of overlapping macroions, effectively highlighting electrostatic interactions by masking any steric interactions. Comparisons of predicted and observed macroscopic properties will provide the most practical, if indirect, tests of the theory. Future applications will examine thermodynamic phase behavior, especially possible implications of the volume energy for the stability and structure of deionized solutions microgels .

###### Acknowledgements.

It is a pleasure to thank C. N. Likos and H. Löwen for helpful discussions and hospitality during a visit to the University of Düsseldorf, where parts of this work were completed. This work was supported by the National Science Foundation under Grant Nos. DMR-0204020 and EPS-0132289.

*

## Appendix A Interactions between Overlapping Macroions

The bare Coulomb interaction between a pair of macroions at center-center separation $`r`$ is given in general by

$$v_{mm}(r)=\frac{e^2}{ϵ}d𝐫^{}d𝐫^{\prime \prime }\frac{\rho _{\mathrm{mon}}(𝐫^{})\rho _{\mathrm{mon}}(𝐫^{\prime \prime })}{|𝐫^{}𝐫^{\prime \prime }𝐫|}.$$
(45)
For nonoverlapping macroions, spherical symmetry reduces the interaction to

$$v_{mm}(r)=\frac{Z^2e^2}{ϵr},r>2a.$$
(46)
For overlapping macroions, the six-dimensional integral in Eq. (45) may be reduced, by exploiting cylindrical symmetry and Gauss’s law, to two-dimensional integrals, which in turn may be evaluated analytically. For star macroions, the result may be expressed piece-wise as follows:

$`v_{mm}(r)={\displaystyle \frac{Z^2e^2}{2ϵa}}\{{\displaystyle \frac{9}{2}}{\displaystyle \frac{7}{4}}{\displaystyle \frac{r}{a}}{\displaystyle \frac{1}{2}}[(3{\displaystyle \frac{a}{r}})(1{\displaystyle \frac{r}{a}})`$ (47)
$`+`$ $`{\displaystyle \frac{r}{a}}\mathrm{ln}\left({\displaystyle \frac{r}{a}}\right)\left]\mathrm{ln}\right({\displaystyle \frac{ar}{a}}){\displaystyle \frac{r}{2a}}{\displaystyle }_1^{a/r1}\mathrm{d}x{\displaystyle \frac{\mathrm{ln}(1+x)}{x}}`$
$`+`$ $`{\displaystyle \frac{r}{2a}}{\displaystyle _1^{a/r}}\mathrm{d}x{\displaystyle \frac{\mathrm{ln}(x1)}{x}}\},0<ra`$
$`v_{mm}(r)={\displaystyle \frac{Z^2e^2}{2ϵa}}\{{\displaystyle \frac{9}{2}}{\displaystyle \frac{7}{4}}{\displaystyle \frac{r}{a}}{\displaystyle \frac{1}{2}}[(3{\displaystyle \frac{a}{r}})(1{\displaystyle \frac{r}{a}})`$ (48)
$`+`$ $`{\displaystyle \frac{r}{a}}\mathrm{ln}\left({\displaystyle \frac{r}{a}}\right)]\mathrm{ln}\left({\displaystyle \frac{ra}{a}}\right)+{\displaystyle \frac{r}{2a}}{\displaystyle _{1a/r}^{a/r}}\mathrm{d}x{\displaystyle \frac{\mathrm{ln}x}{1x}}\},`$
$`a<r2a.`$
For computational purposes, the remaining integrals may be expressed as convergent series:

$`v_{mm}(r)={\displaystyle \frac{Z^2e^2}{2ϵa}}\{{\displaystyle \frac{9}{2}}{\displaystyle \frac{7}{4}}{\displaystyle \frac{r}{a}}{\displaystyle \frac{1}{2}}[(3{\displaystyle \frac{a}{r}})(1{\displaystyle \frac{r}{a}})`$
$`+{\displaystyle \frac{r}{a}}\mathrm{ln}\left({\displaystyle \frac{r}{a}}\right)\left]\mathrm{ln}\right({\displaystyle \frac{ar}{a}})+{\displaystyle \frac{r}{2a}}({\displaystyle \frac{1}{2}}\left[\mathrm{ln}\left({\displaystyle \frac{r}{a}}\right)\right]^2`$
$`{\displaystyle \frac{1}{2}}\left[\mathrm{ln}\left({\displaystyle \frac{r}{ar}}\right)\right]^2+{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(r/a)^n(r/(ar))^n}{n^2}}`$
$`+2{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(r/(ar))^{2n1}2}{(2n1)^2}})\},ra/2`$ (49)
$`v_{mm}(r)={\displaystyle \frac{Z^2e^2}{2ϵa}}\{{\displaystyle \frac{9}{2}}{\displaystyle \frac{7}{4}}{\displaystyle \frac{r}{a}}{\displaystyle \frac{1}{2}}[(3{\displaystyle \frac{a}{r}})(1{\displaystyle \frac{r}{a}})`$
$`+{\displaystyle \frac{r}{a}}\mathrm{ln}\left({\displaystyle \frac{r}{a}}\right)\left]\mathrm{ln}\right({\displaystyle \frac{ar}{a}})+{\displaystyle \frac{r}{2a}}({\displaystyle \frac{1}{2}}\left[\mathrm{ln}\left({\displaystyle \frac{r}{a}}\right)\right]^2`$
$`+{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(r/a)^n+(1a/r)^n2}{n^2}})\},`$
$`a/2<ra`$ (50)
$`v_{mm}(r)`$ $`=`$ $`{\displaystyle \frac{Z^2e^2}{2ϵa}}\{{\displaystyle \frac{9}{2}}{\displaystyle \frac{7}{4}}{\displaystyle \frac{r}{a}}{\displaystyle \frac{1}{2}}[(3{\displaystyle \frac{a}{r}})(1{\displaystyle \frac{r}{a}})`$ (51)
$`+`$ $`{\displaystyle \frac{r}{a}}\mathrm{ln}\left({\displaystyle \frac{r}{a}}\right)\left]\mathrm{ln}\right({\displaystyle \frac{ra}{a}})`$
$``$ $`{\displaystyle \frac{r}{2a}}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(a/r)^n(1a/r)^n}{n^2}}\},`$
$`a<r2a`$
For microgel macroions, the bare interaction may be expressed more compactly as

$`v_{mm}(r)`$ $`=`$ $`{\displaystyle \frac{Z^2e^2}{ϵa}}[{\displaystyle \frac{6}{5}}{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{r}{a}}\right)^2+{\displaystyle \frac{3}{16}}\left({\displaystyle \frac{r}{a}}\right)^3`$ (52)
$``$ $`{\displaystyle \frac{1}{160}}\left({\displaystyle \frac{r}{a}}\right)^5],r2a.`$
Finally, the effective pair interaction between microgels, $`v_{\mathrm{eff}}(r)`$, is the sum of Eq. (52) and the induced interaction, obtained by Fourier transforming Eq. (42):

$`v_{\mathrm{ind}}(r)`$ $`=`$ $`{\displaystyle \frac{9Z^2e^2}{2ϵ\kappa ^4a^4r}}\{(1e^{\kappa r}+{\displaystyle \frac{1}{2}}\kappa ^2r^2+{\displaystyle \frac{1}{24}}\kappa ^4r^4)(1{\displaystyle \frac{1}{\kappa ^2a^2}})+{\displaystyle \frac{2}{\kappa a}}e^{2\kappa a}\mathrm{sinh}(\kappa r)`$
$`+`$ $`\left[e^{2\kappa a}\mathrm{sinh}(\kappa r)+2\kappa ^2ar+{\displaystyle \frac{1}{3}}\kappa ^4(4a^3r+ar^3)\right]\left(1+{\displaystyle \frac{1}{\kappa ^2a^2}}\right){\displaystyle \frac{2r}{a}}\left(1+2\kappa ^2a^2+{\displaystyle \frac{8}{15}}\kappa ^4a^4\right)`$
$``$ $`{\displaystyle \frac{r^3}{3a^3}}(\kappa ^2a^2+{\displaystyle \frac{4}{3}}\kappa ^4a^4){\displaystyle \frac{1}{720}}{\displaystyle \frac{\kappa ^4}{a^2}}r^6\},r2a.`$ (53)


format: text

text:
# Glassy Behavior of Electrons as a Precursor to the Localization Transition

## 1 Introduction

Disorder-driven metal-insulator transitions (MITs) in strongly correlated electronic systems continue to demonstrate fascinating behaviors that remain to be elucidated. The competition of electron-electron interactions and disorder is most pronounced in the insulator, where one generally expects the existence of a large number of low energy configurations of the electronic density, leading to glassy behavior, metastability, and slow relaxation. Far from the transition, the kinetic energy of the electrons is negligible and one can utilize a classical picture, as pioneered in early work by Efros and Shklovskii . The precise role of glassy dynamics is more difficult to assess closer to the MIT, where quantum fluctuations due to electronic mobility become important. In recent work , we have developed a formalism capable of describing such glassy behavior of electrons both in the classical limit and in presence of quantum fluctuations. Here we extend this approach to incorporate Anderson localization effects , which prove to have a dramatic effect on the stability of the glass phase.

The following physical picture emerges from our study: Anderson localization leads to the formation of bound electronic states, thus suppressing quantum fluctuations introduced by mobile electrons. However, we find that glassy behavior is not eliminated as soon as the electrons become delocalized. Very close to the MIT the electronic mobility remains very small, thus limiting the size of quantum fluctuations needed to overcome barriers between metastable states. An intermediate metallic glass phase is predicted, separating the insulator and the normal metal.

## 2 Model

We consider the simplest microscopic model capable of describing the interplay of glassy behavior and Anderson localization effects. It is given by a one-band tight-binding model of spinless electrons at half filling, with nearest-neighbor hopping $`t_{ij}`$, interacting by infinite range inter-site interactions $`V_{ij}`$, in presence of random site energies $`\epsilon _i`$. For technical simplicity, these interactions are chosen also to be random , with zero mean and variance $`<V_{ij}^2>=`$ $`V^2/N`$, where $`N`$ is the number of sites in the systems. In the following, we briefly outline the formalism that we use; the details will be presented elsewhere . We first formally average over random interactions using the replica method , which generates a four-density interaction term. This term is then decoupled by introducing Hubbard-Stratonovich fields, and the resulting action takes the form $`S_{eff}=S_{el}+S_{int}`$, where

$$S_{el}=_o^\beta 𝑑\tau \underset{a}{}\underset{<ij>}{}c_i^a(\tau )[(_\tau +\epsilon _i)\delta _{ij}t_{ij}]c_i^a(\tau ),$$
(1)
$`S_{int}=V{\displaystyle \underset{i}{}}{\displaystyle _o^\beta }𝑑\tau {\displaystyle _o^\beta }𝑑\tau ^{}`$ $`[{\displaystyle \underset{a}{}}\stackrel{~}{Q}_i^{aa}(\tau ,\tau ^{})\delta n_i^a(\tau )\delta n_i^a(\tau ^{})`$ (2)
$`+{\displaystyle \underset{a<b}{}}Q_i^{ab}(\tau ,\tau ^{})\delta n_i^a(\tau )\delta n_i^b(\tau ^{})].`$
We have used standard functional integration over replicated Grassmann fields, where $`a=1,\mathrm{},n(n0)`$ are the replica indices . Here, the operators $`\delta n_i^a(\tau )=(c_i^a(\tau )c_i^a(\tau )1/2)`$ represent the density fluctuations from half filling.

In the classical limit ($`t_{ij}=0`$), these equations can be solved exactly due to the infinite range of the interactions $`V_{ij}`$, but once finite-range hopping is introduced, the resulting electronic degrees of freedom generate additional effective nonlocal interactions between density modes, and we have to resort to approximations. A mean-field treatment of the interactions is then provided by evaluating the $`Q`$-fields in the saddle-point approximation, which also makes it possible to identify the replica symmetry breaking instability used to establish the onset of glassy ordering .

The resulting instability criterion takes the form

$$1=V^2<\frac{1}{N}\underset{ij}{}(\chi _{ij})^2>_{dis}.$$
(3)
Here, the non-local compressibilities $`\chi _{ij}`$ are defined by

$$\chi _{ij}=\frac{1}{\beta }_o^\beta 𝑑\tau _o^\beta 𝑑\tau ^{}\left[<\delta n_i^a(\tau )\delta n_j^a(\tau ^{})>_o<\delta n_i^a(\tau )>_o<\delta n_j^a(\tau ^{})>_o\right],$$
(4)
where the quantum averages $`<\mathrm{}>_o`$ are taken with respect to the full replica-symmetric action $`S_{eff}`$ evaluated at the saddle-point, and $`<\mathrm{}>_{dis}`$ indicates the average with respect to disorder.

Further simplification is obtained  by concentrating on the limit of large disorder, where to leading order one can evaluate the compressibilities by setting $`V=0`$ in $`S_{eff}`$. The calculation then reduces to calculating the compressibilities with respect to a system of noninteracting electrons in presence of disorder. Physically, this calculation amounts to focusing on electrons in the vicinity of the Anderson transition, and examining the leading perturbation introduced by turning on inter-site interactions $`V_{ij}`$. The crucial insight is obtained by observing that in contrast to the average compressibilities $`<\chi _{ij}>_{dis}`$, the quantities $`<(\chi _{ij})^2>_{dis}`$ actually diverge  at $`T=0`$ at the Anderson transition. We immediately conclude that Anderson localization introduces a singular perturbation to the stability of the glass phase, as anticipated in Ref. . In addition, we note that the sum in Eq. (3) runs over positive definite terms, so that a lower bound to the size of the glass phase is obtained by retaining only the local ($`i=j`$) term. We have examined this quantity using several different approaches to Anderson localization, but the simplest closed form expression can be obtained by using the recently developed typical medium theory .

We have used a simple semi-circular model density of states  in presence of a uniform distribution of random site energies of width $`W`$, and have done calculations as a function of the disorder strength $`W`$ and temperature $`T`$. We find that $`<(\chi _{ii})^2>_{dis}(W_cW)^1`$ indeed blows up at $`T=0`$ precisely at the Anderson transition, which for this model occurs at $`W_c1.36B`$, where $`B`$ is the electronic bandwidth. As a result, the instability criterion for the emergence of the glass phase is satisfied before the Anderson transition is reached, for any finite value of the interaction $`V`$. To illustrate our findings, we present the resulting phase diagram for $`W/V=2`$, as a function of $`W`$ and $`T`$ in Fig. 1. Note that the “kink” in the glass transition boundary at $`W=W_c`$ (full line) reflects the crudeness of our treatment of localization, which ignores the fact that inelastic scattering effects will tend to delocalize the electrons at any finite temperature. These inelastic effects can be accounted for by keeping the interaction terms in $`S_{eff}`$, but these corrections will be investigated in more detail elsewhere. To illustrate the qualitative effects of such inelastic scattering, we phenomenologically introduce a temperature-dependent scattering rate of the form $`\eta =AT`$. The resulting modification of the phase boundary for the choice $`A=0.1`$ is presented by a dotted line, leading to the rounding of the kink, as expected.

## 3 Discussion

Our results strongly suggest that the onset of glassy dynamics should be considered as a precursor to disorder-driven metal-insulator transitions. Physically, as the system approaches localization the electrons are barely mobile, allowing for even moderate interactions to induce meta-stability and glassy freezing. We do not expect this glassy freezing to have a direct effect on the average conductivity, but it should be more readily identified by examining the slowing down of the electron dynamics. In particular, as the glassy phase is entered, we expect large enhancements for the low frequency components of the conductivity noise spectra, reflecting slow transitions between emerging metastable states. An interesting question relates to the the size of the proposed metallic glass phase as a function of disorder. This intermediate phase is expected to be extremely small for very strong disorder, since in this limit the effect of the interactions leading to glassy behavior should be very modest. But what should be expected for very weak disorder but strong interactions? In the limit where $`r_s=V/E_F>>1`$, insulating behavior can emerge due to Wigner crystallization, even in absence of disorder. As the Wigner crystal is approached, we expect the local compressibility to be strongly reduced, reflecting strong short range order induced by electronic correlations. As a result, the possibility for many electronic configurations is again impaired, suppressing glassy ordering. We thus expect that the best chance to observe the metallic glass phase is for moderately disordered systems with $`W/V1`$. Remarkably, very recent experiment on a two-dimensional electron gas in silicon seems to support this picture. In this material, clear evidence for an intermediate metallic glass phase has been discovered for low mobility samples , but not in high-mobility devices. In latter systems, the emergence of glassy dynamics seems to coincide with the MIT , as one would expect for systems closer to a (disordered) Wigner crystal.

This work was supported by the NSF grant DMR-9974311 and the National High Magnetic Field Laboratory. We thank S. Bogdanovich, J. Jaroszyński, and D. Popović for useful and stimulating discussions.

## 4 Figure captions

Figure 1

Phase diagram for spinless interacting electrons in presence of disorder, as a function of temperature $`T`$ and the Fermi energy $`E_F`$, expressed in units of the disorder strength $`W`$. Results are presented for moderate interaction strength $`W/V=2`$. The glass transition boundary is shown, in absence of inelastic scattering (full line) and in presence of an inelastic scattering rate $`\eta =0.1T`$ (dotted line). Note that the glassy phase emerges before the metal-insulator transition (MIT) is reached (heavy dot). The (glassy) Anderson insulator is shown with a heavy full line (at $`T=0`$).


text:
# Nonparametric estimation of distribution and density functions in presence of missing data: an IFS approach

## 1 Introduction

Let $`X_1,X_2,\mathrm{},X_n`$ be an i.i.d. sample drawn from a random variable $`X`$ with unknown distribution function $`F`$ with compact support $`[\alpha ,\beta ]`$. The empirical distribution function (e.d.f.)

$$\widehat{F}_n(x)=\frac{1}{n}\underset{i=1}{\overset{n}{}}\chi (X_ix)$$

is one commonly used estimator of the unknown distribution function $`F`$ (here $`\chi `$ is the indicator function). The e.d.f. has an impressive set of good statistical properties such as it is first order efficient in the minimax sense (see , , , , and ). More or less recently, other second order efficient estimators have been proposed in the literature for special classes of distribution functions $`F`$. Golubev and Levit (1996a, b) and  are two of such examples. It is rather curious that a step-wise function can be such a good estimator and, in fact,  shows that, for the class of analytic functions, for small sample sizes, the e.d.f. is not the best estimator. In this paper we study the properties of a new class of distribution function estimators based on iterated function systems (IFSs) introduced by the authors in a previous work . IFSs have been introduced in  and . The main idea on which this method is based consists of thinking the estimation of $`F`$ as the fixed point of a contraction $`T`$ on a complete metric space. The operator $`T`$ is defined in terms of a family of affine maps $`𝒲`$ and a vector of parameters $`p`$. For a given family $`𝒲`$, $`T`$ depends only on the choice $`p`$. The idea, known as inverse approach (see Section 2) is to determine $`p`$ by solving a constrained quadratic optimization problem built in terms of sample moments. In this paper this optimization problem is solved by a penalization method. The nature of affine maps allow to derive easily the Fourier transform of $`F`$ and, when available, an explicit formula for the density of $`F`$ via anti Fourier transform. In this way, given $`𝒲`$ and $`p`$ we have at the same time estimators for the distribution, characteristic and density functions.

The paper is organized as follows. In Section 2 the inverse approach is presented and a penalization method is proposed in order to solve a quadratic optimization problem. We also discuss the choice of the family of maps $`𝒲`$. In Section 3 numerical results and comparisons with classical estimators are shown for small samples via Monte Carlo Analysis.

Finally we show an application of these estimators when the empirical distribution function (or the kernel density estimator for the density) cannot be applied. We will consider situations of missing data when, for example, the data can only be observed on some windows of the support of $`F`$. This can be the case of directional data analysis when, for some reason, instruments are not able for technical or physical reason to collect data in same range of angles say $`A`$ and $`B`$, $`A,B[0,2\pi ]`$. For $`x`$ in $`A`$ or $`B`$ the e.d.f. will be constant and, at the same time, the kernel density estimator will estimate a plurimodal distribution for these data. In this case we will show examples in which the IFS estimator does it job incredibly well.

Tables and figures can be found at the end of the paper after the references.

## 2 An IFS estimator

The theory of distribution function approximation via IFSs we will use to derive estimators is due to . Results from this section, apart from were explicitly mentioned, are from the cited authors. Let $`(X)`$ be the set of probability measures on $`(X)`$, the $`\sigma `$-algebra of Borel subsets of $`X`$ where $`(X,d)`$ is a compact metric space (in our case will be $`X=[\alpha ,\beta ]`$ and $`d`$ the Euclidean metric.)

In the IFSs literature the following Hutchinson metric plays a crucial role

$$d_H(\mu ,\nu )=\underset{f\mathrm{Lip}(X)}{sup}\left\{_Xfd\mu _Xfd\nu \right\},\mu ,\nu (X)$$

where

$$\mathrm{Lip}(X)=\{f:X,|f(x)f(y)|d(x,y),x,yX\}$$

thus $`((X),d_H)`$ is a complete metric space \[\[, see\]\]hutch.

We denote by $`(𝐰,𝐩)`$ an $`N`$-maps contractive IFS on $`X`$ with probabilities or simply an $`N`$-maps IFS, that is, a set of $`N`$ affine contraction maps, $`𝐰=(w_1,w_2,\mathrm{},w_N)`$,

$$w_i=a_i+b_ix,\mathrm{with}|b_i|<1,b_i,a_i,i=1,2,\mathrm{},N,$$

with associated probabilities $`𝐩=(p_1,p_2,\mathrm{},p_N)`$, $`p_i0`$, and $`_{i=1}^Np_i=1`$. The IFS has a contractivity factor defined as

$$c=\underset{1iN}{\mathrm{max}}|b_i|<1$$

Consider the following (usually called Markov) operator $`M:(X)(X)`$ defined as

$$M\mu =\underset{i=1}{\overset{N}{}}p_i\mu w_i^1,\mu (X),$$
(1)
where $`w_i^1`$ is the inverse function of $`w_i`$ and $``$ stands for the composition. In Hutchinson (1981) it was shown that $`M`$ is a contraction mapping on $`((X),d_H)`$ i.e. for all $`\mu ,\nu (X)`$, $`d_H(M\mu ,M\nu )cd_H(\mu ,\nu )`$. Thus, there exists a unique measure $`\overline{\mu }(X)`$, the invariant measure of the IFS, such that $`M\overline{\mu }=\overline{\mu }`$ by Banach theorem. Associated to each measure $`\mu (X)`$, there exists a distribution function $`F`$. In terms of it the previous operator $`M`$ can be rewritten as

$$TF(x)=\{\begin{array}{cc}0\hfill & \text{if }x\alpha \hfill \\ & \\ \underset{i=1}{\overset{N}{}}p_iF(w_i^1(x))\hfill & \text{if }\alpha <x<\beta \hfill \\ & \\ 1\hfill & \text{if }x\beta \hfill \end{array}$$

### 2.1 Minimization approach

For affine IFSs there exists a simple and useful relation between the moments of probability measures on $`(X)`$. Given a $`N`$-maps IFS$`(𝐰,𝐩)`$ with associated Markov operator $`M`$, and given a measure $`\mu (X)`$ then, for any continuous function $`f:X`$,

$$_Xf(x)d\nu (x)=_Xf(x)\mathrm{d}(M\mu )(x)=\underset{i=1}{\overset{N}{}}p_i_X(fw_i)(x)d\mu (x),$$
(2)
where $`\nu =M\mu `$. In our case $`X=[\alpha ,\beta ]`$ so we readly have a relation involving the moments of $`\mu `$ and $`\nu `$. Let

$$g_k=_Xx^kd\mu ,h_k=_Xx^kd\nu ,k=0,1,2,\mathrm{},$$
(3)
be the moments of the two measures, with $`g_0=h_0=1`$. Then, by (2), with $`f(x)=x^k`$, we have

$$h_k=\underset{j=0}{\overset{k}{}}\left(\genfrac{}{}{0pt}{}{k}{j}\right)\left\{\underset{i=1}{\overset{N}{}}p_ib_i^ja_i^{kj}\right\}g_j,k=1,2,\mathrm{},.$$

Set $`X=[\alpha ,\beta ]`$ and let $`\mu `$ and $`\mu ^{(j)}(X)`$, $`j=1,2,\mathrm{}`$ with associated moments of any order $`g_k`$ and

$$g_k^{(j)}=_Xx^kd\mu ^{(j)}.$$

Then, the following statements are equivalent (as $`j\mathrm{}`$ and $`k0`$):

1. $`g_k^{(j)}g_k`$,
2. $`f𝐂(X)`$, $`_Xfd\mu ^{(j)}_Xfd\mu ,`$ (weak\* convergence),
3. $`d_H(\mu ^{(j)},\mu )0`$.

(here $`𝐂(X)`$ is the space of continuous functions on $`X`$). This result gives a way to find and appropriate set of maps and probabilities by solving the so called problem of moment matching. With the solution in hands, given the convergence of the moments, we also have the convergence of the measures and then the stationary measure of $`M`$ approximates with given precision (in a sense specified by the collage theorem below) the target measure $`\mu `$ \[\[, see\]\]bd.

Next result, called the collage theorem is a standard product of the IFS theory and is a consequence of Banach theorem.

(Collage theorem) :  Let $`(Y,d_Y)`$ be a complete metric space. Given an $`yY`$, suppose that there exists a contractive map $`f`$ on $`Y`$ with contractivity factor $`0c<1`$ such that $`d_Y(y,f(y))<\epsilon `$. If $`\overline{y}`$ is the fixed point of $`f`$, i.e. $`f(\overline{y})=\overline{y}`$, then $`d_Y(\overline{y},y)<\frac{\epsilon }{1c}`$.

So if one wishes to approximate a function $`y`$ with the fixed point $`\overline{y}`$ of an unknown contractive map $`f`$, it is only needed to solve the inverse problem of finding $`f`$ which minimizes the collage distance $`d_Y(y,f(y))`$.

The main result in Forte and Vrscay that we will use to build one of the IFS estimators is that the inverse problem can be reduced to minimize a suitable quadratic form in terms of the $`p_i`$ given a set of affine maps $`w_i`$ and the sequence of moments $`g_k`$ of the target measure. Let

$$\mathrm{\Pi }^N=\{𝐩=(p_1,p_2,\mathrm{},p_N):p_i0,\underset{i=1}{\overset{N}{}}p_i=1\}$$

be the simplex of probabilities. Let $`𝐰=(w_1,w_2,\mathrm{},w_N)`$, $`N=1,2,\mathrm{}`$ be subsets of $`𝒲=\{w_1,w_2,\mathrm{}\}`$ the infinite set of affine contractive maps on $`X=[\alpha ,\beta ]`$ and let $`𝐠`$ the set of the moments of any order of $`\mu (X)`$. Denote by $`M`$ the Markov operator of the $`N`$-maps IFS $`(𝐰,𝐩)`$ and by $`\nu _N=M\mu `$, with associated moment vector of any order $`𝐡_N`$. The collage distance between the moment vector of $`\mu `$ and $`\nu _N`$

$$\mathrm{\Delta }(𝐩)=𝐠𝐡_N_{\overline{l}^2}:\mathrm{\Pi }^N$$

is a continuous function and attains an absolute minimum value $`\mathrm{\Delta }_{\mathrm{min}}`$ on $`\mathrm{\Pi }^N`$ where

$$𝐱_{\overline{l}^2}=x_0^2+\underset{k=1}{\overset{\mathrm{}}{}}\frac{x_k^2}{k^2}.$$

Moreover, $`\mathrm{\Delta }_{\mathrm{min}}^N0`$ as $`N\mathrm{}`$. Thus, the collage distance can be made arbitrarily small by choosing a suitable number of maps and probabilities.

The above inverse problem can be posed as a quadratic programming one in the following notation

$$S(𝐩)=(\mathrm{\Delta }(𝐩))^2=\underset{k=1}{\overset{\mathrm{}}{}}\frac{(h_kg_k)^2}{k^2}$$

$$D(X)=\{𝐠=(g_0,g_1,\mathrm{}):g_k=_Xx^kd\mu ,k=0,1,\mathrm{},\mu (X)\}$$

Then by (2) there exists a linear operator $`A:D(X)D(X)`$ associated to $`M`$ such that $`𝐡_N=A𝐠`$. In particular

$$h_k=\underset{i=1}{\overset{N}{}}A_{ki}p_i,k=1,2,\mathrm{}\text{where}A_{ki}=\underset{j=0}{\overset{\mathrm{}}{}}\left(\genfrac{}{}{0pt}{}{k}{j}\right)b_i^ja_i^{kj}g_j$$
(4)
Thus

$$S(𝐩)=𝐩^tQ𝐩+𝐁^t𝐩+C,$$
$`(𝐐)`$
$$\text{where}Q=[q_{ij}],q_{ij}=\underset{k=1}{\overset{\mathrm{}}{}}\frac{A_{ki}A_{kj}}{k^2},i,j=1,2,\mathrm{},N,$$

$$B_i=2\underset{k=1}{\overset{\mathrm{}}{}}\frac{g_k}{k^2}A_{ki},i=1,2,\mathrm{},N\text{and}C=\underset{k=1}{\overset{\mathrm{}}{}}\frac{g_k^2}{k^2}.$$
(5)
The series above are convergent as $`0A_{ni}1`$ and the minimum can be found by minimizing the quadratic form on the simplex $`\mathrm{\Pi }^N`$.

The estimator will then be built by substituting the moments of the target measure with the empirical moments and by truncation of the above series to a finite sum.

### 2.2 Numerical solutions

When practical cases are considered, in particular concerning estimation, the previous series have to be truncated and this implies that the matrix $`Q`$ is assured to be definite positive. Standard numerical procedures for the minimization of constrained quadratic optimization problems involving positive definite quadratic forms cannot be used in this context. To solve this problem an approach is to build the following penalized function $`L_\lambda `$

$$L_\lambda (𝐩)=𝐩^tQ𝐩+𝐁^t𝐩+C+\lambda \left(1\underset{i=1}{\overset{N}{}}p_i\right)^2$$

and then to study the following problem

$$\mathrm{min}L_\lambda (𝐩),0p_i1$$
$`(LOP)`$
It is trivial that an optimizer $`𝐩^{}`$ of (LOP) such that $`_{i=1}^Np_i^{}=1`$ is also an optimizer for the problem

$$\mathrm{min}S(𝐩),𝐩\mathrm{\Pi }^N$$
$`(OP)`$
For solving (LOP) numerically, we have used the method L-BFGS-B due to  which allows to minimize a nonlinear function with box constraints, i.e. when each variable can be given a lower and/or upper bound. The initial value of this procedure must satisfy the constraints. This uses a limited-memory modification of the BFGS quasi-Newton method. The method ‘”BFGS”’ is a quasi-Newton method (also known as a variable metric algorithm).

### 2.3 The choice of affine maps

As we are mostly concerned with estimation, we briefly discuss the problem of choosing the maps. In  the following two sets of wavelet-type maps are proposed. Fixed and index $`i^{}`$, define

$$\gamma _{ij}=\frac{x\alpha +(j1)(\beta \alpha )}{2^i}+\alpha ,i=1,2,\mathrm{},i^{}j=1,2,\mathrm{},2^i$$

and

$$\eta _{ij}=\frac{x\alpha +(j1)(\beta \alpha )}{i},i=2,\mathrm{},i^{}j=2,\mathrm{},i.$$

Then set $`N=\underset{i=1}{\overset{i^{}}{}}2^i`$ or $`N=i^{}(i^{}1)/2`$ respectively. To choose the maps, consider the natural ordering of the maps $`\omega _{ij}`$ and operate as follows

$$𝒲_1=\{w_1=\gamma _{11},w_2=\gamma _{12},w_3=\gamma _{21},\mathrm{},w_6=\gamma _{24},w_7=\gamma _{31},\mathrm{},w_N=\gamma _{i^{}2^i^{}}\}$$

and

$$𝒲_2=\{w_1=\eta _{22},w_2=\eta _{32},w_3=\eta _{33},w_4=\eta _{42},\mathrm{},w_6=\eta _{44},\mathrm{},w_N=\eta _{i^{}i^{}}\}$$

respectively. In  we proposed the following quantile based maps

$$𝒬_1=\{w_i(x)=(q_{i+1}q_i)x+q_i,i=1,2,\mathrm{},N\}$$

where $`q_i=F^1(u_i)`$, and $`0=u_1<u_2<\mathrm{}<u_N<u_{N+1}=1`$ are $`N+1`$ equally spaced points on $`[0,1]`$. With these maps, it has been shown that, there is no need to use a moment matching approach. In particular, given $`p_i=1/N`$, the IFSs turns out to be a smoother of the e.d.f. and so it has nice small sample and asymptotic statistical properties (see cited reference) even for non compact support distribution functions $`F`$. Here we will also mix the quantile information with the moment matching idea. To distinguish the two cases (fixed $`p_i=1/N`$ or $`p`$ solution of $`(\mathrm{𝐐𝐏})`$) we will use the notation $`𝒬_1`$ and $`𝒬_2`$ later on.

### 2.4 Fourier analysis results

We recall, from  results that are rather straight forward to prove but also essential to us since we will use them in density estimation and in particular in presence of missing data. Simplicity is due to affinity of the maps. We assume that the support of the measures is $`X=[0,1]`$ without loss of generality.

Given a measure $`\mu (X)`$, the Fourier transform (FT) $`\varphi :`$, where $``$ is the complex space, is defined by the relation

$$\varphi (t)=_Xe^{itx}d\mu (x),t,$$

with the well known properties $`\varphi (0)=1`$ and $`|\varphi (t)|1`$, $`t`$. It can be shown that the space of characteristic functions $`T(X)`$ can be made metric and complete with an opportune metric. Thus, given a $`N`$-maps affine IFS$`(𝐰,𝐩)`$ it is possibile to define a new linear operator $`B:T(X)T(X)`$ whose unique fixed point reads as

$$\overline{\varphi }(t)=\underset{k=1}{\overset{N}{}}p_ke^{ita_k}\overline{\varphi }(b_kt),t.$$

This $`\overline{\varphi }(t)`$ is the FT of the fixed point of the $`N`$-maps affine IFS$`(𝐰,𝐩)`$. Now \[\[, see e.g.\]\]tarter, suppose that the target distribution $`F`$ admits a density $`f`$. It is possible to write the density $`f`$ via Fourier expansion. In fact,

$$\varphi (t)=_0^1f(x)e^{itx}dx=_0^1e^{itx}dF(x)$$

thus

$$f(x)=\frac{1}{2\pi }\underset{k=\mathrm{}}{\overset{+\mathrm{}}{}}B_ke^{ikx}\text{where}B_k=\varphi (k).$$

## 3 Relative efficiency and estimation in presence of missing data

Suppose to have an i.i.d. sample on $`n`$ observations with common unknown distribution function $`F`$ with compact support on $`[\alpha ,\beta ]`$ which has all the moments up to order $`M`$. An IFS estimator of $`F`$ is the fixed point of the functional $`TF`$ where the $`N`$ maps are choosen in advance and the $`p_i`$ are the solution of the (QP) quadratic programming problem where in the expression on $`A_{ik}`$, $`B_i`$ and $`C`$ we replace, in equations (5) and (4), the true moments $`g_k`$ with the sample moments $`m_k`$, $`k=0,1,\mathrm{},M`$ for a fixed $`M`$ and we consider the first $`M`$ terms of the series involved.

Given the solution of (QP), we have an estimator for $`F`$ and an estimator for the characteristic function of $`F`$, say $`\widehat{\varphi }`$. Suppose that $`F`$ posseses a density $`f`$ then we have further a (Fourier) density estimator for $`f`$

$`\widehat{f}(x)`$ $`={\displaystyle \frac{1}{2\pi }}{\displaystyle \underset{k=m}{\overset{+m}{}}}\widehat{B}_ke^{ikx}`$
$`={\displaystyle \frac{1}{2\pi }}+{\displaystyle \frac{1}{\pi }}{\displaystyle \underset{i=1}{\overset{m}{}}}\left\{\mathrm{Re}(\widehat{B}_k)\mathrm{cos}(kx)\mathrm{Im}(\widehat{B}_k)\mathrm{sin}(kx)\right\}`$
where $`\widehat{B}_k=\widehat{\varphi }(k)`$ and $`m`$, the number of Fourier terms, is choosen in the usual way, i.e.

$$\text{if}\left|\widehat{B}_{m+1}\right|^2\text{and}\left|\widehat{B}_{m+2}\right|^2<\frac{2}{n+1}\text{then use the first }m\text{ coefficents}$$

\[\[, see again\]\]tarter. Tables 1 and 2 show camparisons between the empirical cumulative distribution function $`\widehat{F}_n`$ and the IFS estimator, say $`\widehat{T}_N`$, for some target distributions $`F`$, in terms of average mean square error (AMSE) and sup-norm (SUP) distance. These tables contain Monte Carlo analysis where 100 simulations have been done for each target distribution. Tables report the average ratio of the sup norm (and AMSE) of the IFSs over the corresponding sup norm (respectively AMSE) of the empirical distribution function.

It is possible to notice that the IFS estimator based on maps $`𝒲_1`$ has good properties for symmetric bell-shaped distributions and distributions with not so heavy tails (see also Figure 2). It is also evident the asymptotic equivalence of the IFSs to the e.d.f. when quantile maps are used. Remark that, for $`𝒲_1`$ we have decided to use 62 maps, for $`𝒲_2`$ 28 maps and $`n/2`$ quantiles for the quantiles maps $`𝒬_1`$ and $`𝒬_2`$. So it is evident that for wavelet-type maps an adjustment can be done by choosing a suitable number of maps in terms of the sample size $`n`$.

### 3.1 What if data are missing?

Suppose now that the for some reason, the $`n`$ sample observations from $`F`$ are in fact a subset of a biggest sample, of unknown size. In practice we do not observe the data on the whole support of $`F`$ $`[\alpha ,\beta ]`$ but only on some windows. This sample reduction has happened due to some sort of censoring. So we are in presence of missing data when we do not know how many data are missing and where exactly they were missed, i.e. we are not in a classical censoring setup. A motivation for this scheme of (non)-observation is the following: suppose one wants to estimate the distribution of the angle of the wind registered by some instruments in degrees (0,360). For some reason, data from angles (15,37) and (62,79) are missing for technical failures or physical obstacles. In this case the empirical distribution function will be flat on these windows and a kernel density estimator will probably show a bimodal behaviour.

Heuristically, this is due to the fact that quantile estimation is inappropriate in this context. At the same time, moments estimation tend to be more robust, in particular if the distribution is symmetric. We only report a graphical example of what can happen. Figure 1 is about a sample from a Beta(2,2) distribution when only the observation in $`(.1,.15)(.37,.43)(.7,.8)`$ are available to the observer all the other being truncated by the instrument (we have choose this interval by hazard). The IFS estimator with $`𝒲_1`$ maps seems to be able to reconstruct the underlying distribution and density function, whistle, for obvious reasons both the e.d.f. and the kernel estimators fail. In this example the relative efficiency (IFS/EDF) is 7% for the AMSE and 23% for the SUP-norm which is dramatically better than expected!

### 3.2 Algorithm flow for estimation

1. calculate sample moments
2. choose the family of maps $`𝒲`$
3. build the quadratic form and solve it for $`p`$
4. if you want to estimate $`F`$ at point $`x`$: take any distribution function, for example the uniform over $`[\alpha ,\beta ]`$ and start to iterate $`T`$
5. stop after few iteration (normally 5 is enough)
6. the “fixed point” of $`T`$ evaluated in $`x`$ is the estimate of $`F(x)`$

In case the support of $`F`$ is not known one case use the range of the sample but the resulting IFS estimator will then try to approximate a distribution function which has exactly that support. If any hints on the shape of the distribution $`F`$ is available, use it to choose the maps.

All the examples, tables and graphics have been done by some software developed by the authors. In particular, a package called ifs is freely available for the R environment system  in the CRAN (Comprehensive R Archive Network) http://cran.R-project.org under the contributed section.

## Conclusions

It seems that this kind of approach can be used to make nonparametric inference when data are missing or sample size are small. Remark that with this method it is only possible to work with distributions with compact support. Moreover, a knowledge on the support itself it is needed. Neverthless, it seams a promising approach and the use of different sets of maps merits further investigation.