text:
# Josephson spectroscopy of a dilute Bose-Einstein condensate in a double-well potential

## I Introduction

The behaviour of dilute Bose-Einstein condensates in multi-well potentials is remarkably rich. The tunnelling of a condensate through an optical lattice potential  provides an atomic physics analogue of a Josephson junction array. In principle, the analogue of a single Josephson junction can be realised experimentally by a condensate in a double-well potential , and this system has attracted considerable theoretical interest . However, the experimental observation of the transition between the dc and ac Josephson effects is challenging because the small energy splitting associated with Josephson oscillations means that thermal or quantum fluctuations will tend to destroy the effect even at the lowest achievable temperatures . While the energy splitting can be increased by lowering the barrier height, it then becomes comparable with that of higher-lying states, and the accuracy of a two-mode approximation to describe the system  becomes questionable.

In this paper, we investigate the Josephson dynamics of a dilute Bose-Einstein condensate in a double-well potential, by considering transitions between energy eigenstates. We show that in the regime of most interest for experiments where the energy splittings are large, the influence of higher-lying states cannot be ignored. The paper is organised as follows. First, we calculate the eigenenergies of the double-well potential as a function of the barrier position. We show that in the vicinity of a level crossing the non-linearity leads to triangular structures in the eigenenergy curves. Similar structures have been shown to occur in a non-linear Landau-Zener model , and near the zone boundary in optical lattices . Next we review the two-state model , and show that as long as the influence of higher-lying states is small, the level structure can be reproduced using three parameters determined by numerical solution of the time-independent Gross-Pitaevskii equation. Finally in Section V we consider the time-dependent evolution as the barrier is moved, and analyse the transition between dc and ac tunnelling in terms of transitions between eigenstates.

## II Eigenenergies of the double-well potential

The eigenenergies of the double-well potential, $`V_{x_0}`$, are found by numerical solution of the one-dimensional (1D) Gross-Pitaevskii (GP) equation

$$i_t\psi (t)=\frac{1}{2}^2\psi (t)+V_{x_0}\psi (t)+g|\psi (t)|^2\psi (t).$$
(1)
Here $`g=4\pi na`$ is the self interaction parameter, $`a`$ is the $`s`$-wave scattering length, and $`n`$ is the number of atoms per unit area in the $`yz`$ plane. Distance, time and energy are measured in terms of the standard harmonic oscillator units (h.o.u.), $`(\mathrm{}/m\omega )^{1/2}`$, $`\omega ^1`$ and $`\mathrm{}\omega `$, respectively. In particular, we have concentrated on a system confined by the potential

$`V_{x_0}=\frac{1}{2}x^2+h\mathrm{exp}[(xx_0)^2],`$ (2)
describing a harmonic trap with a Gaussian potential barrier with height $`h`$ and unit width. We have also carried out three-dimensional simulations, and verified that the essential physics is described by a 1D model. This model facilitates an exploration of a wide range of values of the self-interaction, $`g`$, barrier height, $`h`$, and barrier velocity.

Time-independent states of the form, $`\psi (x,t)=\mathrm{e}^{i\mu t}\phi (x)`$, where $`\mu `$ is the chemical potential, are found by numerical solution of the discretized system of equations,

$`f_i`$ $``$ $`\left(\phi _{i+1}+\phi _{i1}2\phi _i\right)/2\mathrm{\Delta }^2+V_i\phi _i+g\phi _i^3\mu \phi _i`$ (3)
$`=`$ $`0,`$ (4)
where $`\phi _i=\phi (x_i)`$ is a real function, and $`\mathrm{\Delta }=1/75`$ is the grid spacing. Note that in general, one must solve for both the real and imaginary parts of $`\phi (x_i)`$ . Since $`\mu `$ is also a variable we have $`n`$ equations and $`n+1`$ unknowns. An additional equation is provided by fixing the normalisation of the wave function,

$$f_{n+1}\underset{i=1}{\overset{n}{}}\phi _i^21=0.$$
(5)
A solution is found iteratively by linearising equations (4) and (5) (with $`\mu \phi _{n+1}`$)

$$f_i(\phi _i^{(p)})+\underset{i=1}{\overset{n+1}{}}\left(\phi _i^{(p+1)}\phi _i^{(p)}\right)\left[\frac{f_i}{\phi _i}\right]^{(p)}0,$$
(6)
where $`\phi ^{(p+1)}`$ is determined from the approximation $`\phi ^{(p)}`$ by solving (6) using the bi-conjugate gradient method . The iterative solution depends on the symmetry of the initial trial wave function $`\phi ^{(0)}`$. For example, to find the lowest energy excited state we begin with a ground state solution at $`x_0=0`$ and change the parity. Higher lying excited states can be found by beginning with an excited state solution of the harmonic trap and moving the barrier in from the side.

In Fig. 1 we show the eigenenergies of the ground level $`g`$ and the first and second excited levels, $`e_1`$ and $`e_2`$, as a function of the barrier position $`x_0`$ with $`g=0.5`$ and $`h=12`$. A triangular structure appears in the upper level at each level crossing. For states $`g`$ and $`e_1`$ this structure is essentially the same as that discussed by Wu and Niu in the context of a non-linear two-state model . The three states in the triangular region correspond to an anti-symmetric state with the population balanced between the two wells, and two higher energy states with most of the population in either the lower or the upper well. The higher energy states are referred to as self-trapping states . Any dissipative mechanism that damps the system towards the ground state will tend to destroy self-trapping, as pointed out in Refs. .

For small $`x_0`$ the energy of the ground state is almost independent of the barrier position. This is because the transfer of atoms through the barrier exactly compensates for the change in potential energy. However, this cannot continue when all the atoms reach one side and the energy becomes strongly position dependent again at the critical displacement, $`x_\mathrm{c}`$ ($`x_\mathrm{c}=0.07`$ for the parameters in Fig. 1).

At $`x_00.5`$ there is a second level crossing where it becomes energetically favourable for the first excited state $`e_1`$ to move from the upper to the lower well. The energy level structure is similar to that at $`x_0=0`$, where it is the ground state population that moves from the upper to the lower well. In both cases, following the eigenenergy curve produces a large change or ‘step’ in the well populations. When considering a moving barrier these steps give rise to a resonance in the tunnelling current. In Section V, we suggest the optimal scheme for observing such tunnelling resonances.

The energy splitting between the ground and lowest excited state is extremely small for $`h=12`$, see Fig. 1(inset). From an experimental viewpoint, lower barriers resulting in larger energy splittings are more interesting. In Fig. 2 we illustrate how the triangular structures evolve as a function of barrier height with $`g=0.5`$. The appearance of the loop structure, which coincides with the threshold for self-trapping, occurs at critical height $`h=2.916`$ (see Section IV). The structure becomes more triangular as $`h`$ increases.

In Fig. 3 we show the eigenenergies for $`g=0.5`$ and $`h=4`$. In this case, the loop structure appears for the first but not the second excited state. The appearance of a loop structure results in the breakdown of adiabatic following of an eigenenergy curve . In Section V, we will see that this breakdown is associated with a discontinuity in the population difference between the two wells.

Finally in Fig. 4, we illustrate the effect of increasing the self-interaction parameter to $`g=5`$. In this case the energy splittings are a significant fraction of the harmonic oscillator energy level spacing and the critical displacement coincides with the position of the crossings between states $`e_1`$ and $`e_2`$. In this regime, the influence of state $`e_2`$ cannot be neglected, and a two-state approximation cannot be used to describe the system.

We stress that by solving the Gross-Pitaevskii equation in three dimensions we have obtained eigenenergy curves qualitatively similar to those shown in Figs. 1 to 4.

## III Non-linear two-state model

In this section we review the two-state model , and apply it to the specific case of a Bose-Einstein condensate in an asymmetric double-well. The wave function, $`\psi (t)`$, is written as

$`\psi (t)`$ $`=`$ $`\varphi _1(t)\psi _1+\varphi _2(t)\psi _2,`$ (7)
where the zeroth-order (real) mode functions $`\psi _i`$, $`i=1,2`$, are eigenfunctions of the time-independent GP equations

$`\mu \psi _i`$ $`=`$ $`\frac{1}{2}^2\psi _i+V_i\psi _i+g|\psi _i|^2\psi _i.`$ (8)
The mode potentials $`V_1`$ and $`V_2`$ are single-well potentials displaced to the left and the right of $`x=0`$, respectively, and satisfy $`V_1(x)=V_2(x)`$ so that $`\psi _1(x)=\psi _2(x)`$. The time-dependent coefficients $`\varphi _i(t)`$ can, in turn, be expressed as

$$\varphi _i(t)=\sqrt{N_i(t)}\mathrm{e}^{i\theta _i(t)},$$
(9)
where $`N_i(t)`$ and $`\theta _i(t)`$ are the population and phase of state $`i`$, with $`N_1(t)+N_2(t)=1`$.

Recalling that the two mode functions satisfy $`\psi _1V_{x_0=0}\psi _1=\psi _2V_{x_0=0}\psi _2`$ and redefining the zero of energy, we obtain the following non-linear equations of motion for $`\varphi (t)=(\varphi _1(t),\varphi _2(t))^\mathrm{T}`$:

$`i\dot{\varphi }(t)`$ $`=`$ $`H(\varphi (t))\varphi (t).`$ (10)
The non-linear Hamiltonian matrix, $`H(\varphi (t))`$, is given by

$$H(\varphi (t))=\frac{1}{2}\left[\begin{array}{cc}\mathrm{\Delta }+NE_\mathrm{C}& E_\mathrm{J}\\ E_\mathrm{J}& \mathrm{\Delta }NE_\mathrm{C}\end{array}\right]$$
(11)
where $`N=N_1N_2`$ is the population difference between the left and right sides of the well, and $`E_\mathrm{J}`$ and $`E_\mathrm{C}`$ are the coupling or Josephson energy and the self-interaction energy, respectively. The self-interaction energy of each mode is

$`E_\mathrm{C}=g\psi _1|\psi _1|^2\psi _1=g\psi _2|\psi _2|^2\psi _2,`$ (12)
and we define the shift in energy of the mode states due to the displacement of the potential barrier to be

$`\frac{1}{2}\mathrm{\Delta }`$ $`=`$ $`\psi _1\left(V_{x_0=0}V_{x_0}\right)\psi _1`$ (13)
$`=`$ $`\psi _2\left(V_{x_0}V_{x_0=0}\right)\psi _2.`$ (14)
If states $`\psi _1`$ and $`\psi _2`$ correspond to the population on the left and right of the barrier, respectively, then linearising yields

$`\frac{1}{2}\mathrm{\Delta }`$ $``$ $`x_0\psi _1\left({\displaystyle \frac{dV_{x_0}}{dx_0}}\right)_{x_0=0}\psi _1=\alpha x_0,`$ (15)
where $`\alpha <0`$ is the rate of change of the energy with displacement. As will be discussed below, this linear approximation works well when the influence of higher-lying states is negligible. Finally, the coupling energy of the two-modes is

$`\frac{1}{2}E_\mathrm{J}`$ $``$ $`\psi _1\left(\frac{1}{2}^2+V_{x_0=0}\right)\psi _2,`$ (16)
where use has been made of the fact that

$`\psi _1\left({\displaystyle \frac{dV_{x_0}}{dx_0}}\right)_{x_0=0}\psi _2`$ $`=`$ $`0.`$ (17)
The Hamiltonian, equation (11), is similar to that considered by Wu and Niu in their study of non-linear Landau-Zener tunnelling . In the next section we discuss the eigenstates of $`H`$ in relation to the numerical solutions discussed in Section II.

## IV Eigenenergies of the two-state model

Within the two-state model, one may find the eigenenergies of the system by substituting $`\theta _1(t)=\theta _2(t)=ϵt`$. The eigenenergies are given by the roots of the quartic equation

$`ϵ^4E_\mathrm{C}ϵ^3`$ $`+`$ $`\frac{1}{4}(E_\mathrm{C}^2E_\mathrm{J}^2\mathrm{\Delta }^2)ϵ^2`$ (18)
$`+`$ $`\frac{1}{4}E_\mathrm{C}E_\mathrm{J}^2ϵ\frac{1}{16}E_\mathrm{C}^2E_\mathrm{J}^2=0.`$ (19)
The population difference between the left and right wells is obtained from the eigenenergies via

$`N`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Delta }}{(E_\mathrm{C}2ϵ)}}.`$ (20)
An important feature of equation (19) is that there are four real roots when the coefficient of the quadratic term $`E_\mathrm{C}^2E_\mathrm{J}^2\mathrm{\Delta }^2`$ is positive and only two when it is negative. For a symmetric double-well ($`\mathrm{\Delta }=0`$), the additional roots appear at the critical point where the self-interaction energy satisfies $`E_\mathrm{C}=E_\mathrm{J}`$. For $`E_\mathrm{C}>E_\mathrm{J}`$, the eigenenergies of the stationary states are $`\pm E_\mathrm{J}/2`$ and $`E_\mathrm{C}/2`$, with the latter being doubly degenerate. The additional roots disappear if the displacement of the barrier is such that $`|\mathrm{\Delta }|=|2\alpha x_0|>\sqrt{E_\mathrm{C}^2E_\mathrm{J}^2}`$.

To test the applicability of the two-state model, we have compared the eigenenergies with those determined from the numerical solution of the GP equation. The energy curves are parametrised by three numbers: the splitting between the two lower levels, which is equal to $`E_\mathrm{J}`$; the energy of the self-trapping states, $`E_\mathrm{C}/2`$; and an energy gradient, $`\alpha `$. In a two-state model, $`\alpha `$ is both the energy gradient of the self-trapping state, and that of the ground state for $`x>x_\mathrm{c}`$. For small $`g`$ and high barriers, such that $`E_\mathrm{C}E_\mathrm{J}`$, $`\alpha E_\mathrm{C}/2x_\mathrm{c}`$.

Fig. 5 shows a comparison between the exact eigenenergies and the model curves for $`h=4`$, $`g=0.5`$. The values of $`E_\mathrm{C}`$, $`E_\mathrm{J}`$, and $`\alpha `$ are taken from the exact solutions. The energy gradient is matched to that of the self-trapping states at $`x_0=0`$ ($`\alpha =3.187`$). This value gives the best agreement when comparing population dynamics (see Section V). However, due to the slight curvature of the eigenenergy curves a smaller value ($`\alpha =2.663`$) gives a better fit to the triangular structure. The agreement between the two-state model and the exact eigenenergies becomes less good as the influence of higher-lying states increases. For example, the model does not predict the almost flat position dependence of the upper levels in $`e_1`$ for $`h=4`$, $`g=5`$, see Fig. 4.

To apply the two-state model one needs to known the value of $`E_\mathrm{J}`$, $`E_\mathrm{C}`$ and $`\alpha `$ for any particular barrier height or non-linearity. In Fig. 6 we show how the energy splittings vary with the non-linearity and the barrier height. The critical point, $`E_\mathrm{C}=E_\mathrm{J}`$, appears as either a critical non-linearity or a critical height depending on which parameter is varied. Note that for large non-linearity, the existence of the second excited state, $`e_2`$, effectively puts an upper limit on the value of $`E_\mathrm{C}`$.

In Fig. 7, we show the left-right population difference of the self-trapping state. A population asymmetry appears at the critical point, $`E_\mathrm{C}=E_\mathrm{J}`$, corresponding to the critical barrier height, $`h=2.916`$ for $`g=0.5`$. The asymmetry increases with increasing barrier height. As a function of the non-linearity, the self-trapping population first increases, then saturates, and finally decreases at large $`g`$ due to the influence of $`e_2`$.

## V Josephson dynamics

We now apply the eigenstate picture to analyse the population dynamics when the barrier is moved uniformly through the condensate at velocity $`v`$. The time-dependent GP equation is integrated using a Crank-Nicholson algorithm and the time evolution analysed in terms of transitions between eigenstates. A similar approach was used previously to study vortex nucleation .

In Fig. 8 we show both the eigenstate population differences and the calculated evolution of the population difference as the barrier is moved from the canter towards the right at a speed $`v=4\times 10^5`$ (bold line). This example is similar to the case studied in . The eigenstate population differences appear as straight diagonal lines. For low velocities, the evolution is adiabatic, and the time-dependent solution follows the ground state population difference curve resulting in a dc tunnelling current. Above a critical velocity, there is a transition to a superposition of ground and excited states such that the population difference remains constant, apart from a decaying oscillatory component (see inset). Subsequently, the excited state component encounters a level crossing with higher-lying states which leads to further ‘steps’ in the population difference or resonances in the tunnelling current. For the example of $`v=4\times 10^5`$ shown in Fig. 8, the speed is just below the critical velocity for the transition from $`e_2`$ to $`e_3`$ leading to a large tunnelling current as the object moves past $`x_00.9`$.

### A Critical velocity

By defining $`\theta (t)=\theta _1(t)\theta _2(t)`$, and rearranging equations (911), one obtains the coupled equations

$$\dot{N}=E_\mathrm{J}\sqrt{1N^2}\mathrm{sin}\theta ,$$
(21)
$$\dot{\theta }=\mathrm{\Delta }NE_\mathrm{C}\frac{E_\mathrm{J}N}{\sqrt{1N^2}}\mathrm{cos}\theta .$$
(22)
Differentiating again and substituting for $`\dot{N}`$ and $`\dot{\theta }`$ gives

$$\ddot{\theta }=\dot{\mathrm{\Delta }}\frac{E_\mathrm{J}(E_\mathrm{C}N\mathrm{\Delta })}{\sqrt{1N^2}}\mathrm{sin}\theta \frac{E_\mathrm{J}^2(1+N^2)}{2(1N^2)}\mathrm{sin}2\theta .$$
(23)
If $`N`$ is roughly constant, this equation may be re-written as $`\ddot{\theta }=_\theta U`$, where $`U(\theta )`$ describes a ‘tilted-washboard’ potential. The critical velocity is determined by setting the minimum gradient of $`U(\theta )`$ equal to zero. For $`N0`$ at $`t=0`$ and using $`\mathrm{\Delta }=2\alpha x_0=2\alpha vt`$, one finds

$$v_\mathrm{c}=\frac{E_\mathrm{J}E_\mathrm{C}}{2\alpha }\frac{E_\mathrm{J}^2}{4\alpha }.$$
(24)
For $`E_\mathrm{C}E_\mathrm{J}`$, one can use the linear approximation, $`\alpha =E_\mathrm{C}/2x_\mathrm{c}`$, where $`x_\mathrm{c}`$ is the critical displacement, which yields $`v_\mathrm{c}E_\mathrm{J}x_\mathrm{c}`$. Taking the values from Fig. 1, $`E_\mathrm{J}=7\times 10^4`$ and $`x_\mathrm{c}=0.07`$, gives $`v_\mathrm{c}5\times 10^5`$. According to equations (21) and (22) the maximum population difference occurs at approximately $`v_\mathrm{c}/\sqrt{2}=3.4\times 10^5`$, which is in excellent agreement with the value determined by numerical integration of the GP equation.

### B Asymmetric initial condition

As stated above, the small critical velocity for high barriers (typically less than a micron per second) makes experimental verification challenging. Although higher critical velocities are obtained for lower barriers, the transition between the dc and ac regime becomes less sharp. This problem can be partially circumvented by noting that one can induce larger population changes by starting with an asymmetric well and moving the barrier back through the origin. This is illustrated in Fig. 9 with a lower barrier height, $`h=4`$. Again for low barrier speeds the population follows the ground state distribution, whereas for faster speeds a transition to an oscillatory current is observed. The critical velocity is a factor of 200 times larger than for $`h=12`$.

For $`x_0<0.6`$ the population difference follows that of the excited state, $`e_1`$, therefore to observe a large population difference, it is important to stop the motion before the barrier reaches $`x_00.4`$. This restricts the evolution to short times or low velocities. In Fig. 10 we show the population difference as a function of the barrier speed for both a symmetric and asymmetric starting condition. Fig. 10(a) shows that the two-state model can no longer predict the correct result if the second level crossing is reached, which for $`t=30`$ occurs when $`v>0.015`$.

Comparison between Fig. 10(a) and (b) illustrates that a larger population difference and a better demarcation of the critical velocity is obtained by moving the barrier from the edge of the condensate inwards (asymmetric initial condition). The numerical results indicate that the critical velocity for this case is similar to that for the symmetric initial state. However, it is not straightforward to predict the critical velocity for the asymmetric initial condition because the relative importance of the coefficients in equation (23) depends on the instantaneous value of $`N(t)`$. For $`h=12`$, $`g=0.5`$, the critical velocity is a factor of three smaller than for the symmetric initial state.

Finally, we consider the interesting case where the influence of higher-lying states severely limits the applicability of the two-state model. In Fig. 11 we show the population change as a barrier with height $`h=4`$ is moved through the centre of a condensate with $`g=5`$. For these parameters the critical velocity is sufficiently high that the level crossing to a higher-lying state is reached before the transition to the ac regime is completed. However, a large population difference for different barrier velocities can still be observed if the motion is halted when the barrier reaches $`x_00.5`$.

The large critical velocity, $`v_\mathrm{c}0.08`$, offers the best potential for the experimental observation of Josephson-like tunnelling resonances. For example, taking a typical sodium condensate with a trap frequency of $`20`$ Hz, the black and grey curves in Fig. 11 correspond to a barrier formed by a blue detuned laser sheet with waist 3 microns, moved by 5 microns in 1.5 s and 1 s, respectively. For these parameters, a model involving higher-lying excited states is essential to give a clear picture of the Josephson population dynamics.

To restate the richness of the double-well system, it is worth noting that as one increases the non-linearity further, the transition from the ground to first excited state corresponds to the formation of a soliton. Consequently by varying only two parameters, the barrier height and the interaction strength, for example using a Feshbach resonance , one can explore the complete parameter space between Josephson tunnelling and soliton formation.

## VI Conclusions

In summary, we have studied the eigenstates of a dilute Bose-Einstein condensates in an asymmetric double-well potential. We have shown that in the regime of interest for experiments, i.e., with a large non-linearity and a low barrier, a two-state description is insufficient to describe the tunnelling dynamics. By determining the influence of higher-lying states we have demonstrated that large population differences can be observed if the barrier motion is halted before the upper tunnelling resonance is reached.

## ACKNOWLEDGMENTS

We thank Andrew MacDonald whose preliminary study stimulated our interest in this problem. We would also like to thank the University of Durham and the Engineering and Physical Sciences Research Council (EPSRC) for financial support.


format: text

text:
# Untitled Document

3rd.April 2002

J. Superconductivity: Incorporating Novel Magnetism, to be published

Interpretation of large room-temperature diamagnetism at low magnetic fields in films of oxidised atactic polypropylene in terms of superconducting current loops

D.M. Eagles1

Abstract

> A simple model is used to analyse published results on large room-temperature diamagnetism for two films of oxidised atactic polypropylene (OAPP) at low magnetic fields. The model involves induced currents expected in circular closed loops of superconductors in fields below the lower critical field $`H_{c1}`$ at which flux penetration would first occur if a metamagnetic transition did not intervene as in OAPP, and the assumption that resistance would be restored at $`H_{c1}`$ (negligible pinning). Fits to the data for the more strongly magnetic sample with the model, allowing two different types of loops with different loop radii $`b_1`$ and $`b_2`$, but with the same cross section $`a`$ of loop material yield $`H_{c1}5260`$ Oe, and fits to the data for the less strongly magnetic sample with two loop sizes and with the same value of $`H_{c1}`$, combined with the knowledge that the minimum number of closed loops of any type is one, requires that the radius $`a`$ of the cross section of the material should be less than about $`0.8\mu `$m, in fair agreement with a maximum radius of 1 $`\mu `$m obtained previously from other data.

Keywords: Room-temperature superconductivity, Oxidised atactic polypropylene, Magnetization, Lower critical field $`H_{c1}`$.

1 19 Holt Road, Harold Hill, Romford, Essex RM3 8PN, England. e-mail: d.eagles@ic.ac.uk

1. INTRODUCTION

There is evidence from at least three different types of experiments that narrow channels through films of oxidised atactic polypropylene (OAPP) are superconducting at room temperature. The three types of evidence are (i) lower limits for conductivity several orders of magnitude greater than that of copper, found by direct  and indirect  methods, (ii) non-thermal destruction of ultra-high conductivity by high pulsed currents, with critical current densities greater than $`10^9`$ A cm-2 , and (iii) negligible electronic contribution to the thermal conductivity . Further support for high-temperature superconductivity in channels of a different polymer is that the thermopower between 87 and 233 K in films of poly(octylmethacrylate) is zero to within estimated errors of measurements .

The magnetic properties of films of OAPP are also unusual. The unusual properties observed in samples in which highly conducting channels occur are (i) a metamagnetic transition in fields of a few kilooersted \[5-7\], (ii) large diamagnetism observed at low fields in about 10% of the samples showing a metamagnetic transition , and (iii) spontaneous forces occurring in some field range tending to push samples to lower magnetic field regions in inhomogeneous magnetic fields .

At least three different models have been suggested for the ultra-high conductivity or the superconductivity \[9-14\], two for the high critical current densities in channels \[10-14\], and two for the unusual magnetic properties . Authors of papers on both types of models for the magnetic properties are agreed that the large diamagnetism in some samples and the occasional spontaneous forces pushing samples out to low-field regions in inhomogeneous fields are associated with superconducting channels which form closed loops, but differ as to what is happening in these closed loops. In  I suggested that large spontaneous currents occur when closed loops form, and gave support for this hypothesis on the basis that some data in a figure in  appear to show that the susceptibility is approximately proportional to (1/field), indicating a constant moment, independent of field, whereas the authors of  appear to think the large diamagnetism (of the order of a percent of a complete Meissner effect at low fields in one sample) is associated with the percentage of the film which is superconducting. We think that this hypothesis is not compatible with the assumption that the diamagnetism is associated with closed loops, since the total fraction of material occupied by conducting channels is typically at most a few percent (channels separated by 7-8 $`\mu `$m estimated in one very thin film, thickness 0.3 $`\mu `$m, showing conducting channels , and channel diameters always less than 2 $`\mu `$m ), the fraction of material in closed loops can be expected to be considerably less than this, and also the Meissner effect may be further reduced if the magnetic-field penetration depth is comparable to channel radii. However, recent experimental results  do not appear to support the spontaneous moment hypothesis well either, as there is at best only a small range of fields for which the diamagnetic moment is approximately constant.

In this paper we explore the hypothesis that the diamagnetic moment is associated with induced supercurrents at fields below which the loop becomes resistant. We make use of a theory for induced currents in circular loops of Type I superconductors below the critical field, discussed in Shoenberg’s book (1952 version) . Looking at the derivation of his results, it appears that they will apply to Type II superconductors below $`H_{c1}`$ if there is negligible flux line pinning, so that resistance is restored as soon as fields reach $`H_{c1}`$. A metamagnetic transition may occur before $`H_{c1}`$ is reached, but we assume that the metamagnetic transition, being a cooperative phenomenon, will not occur until the average field throughout the material of the loop reaches a certain value $`H_M`$, and so supercurrents can continue to flow even if the field near the outer surface of the loop is greater than $`H_M`$, provided that the average field in the material is smaller than this value.

Since only a small fraction of the conducting channels form closed loops, we assume that the observed diamagnetism at low fields is a superposition of a diamagnetism associated with closed loops, and a smaller positive magnetisation as found by a combination of results of observations made on the next day after closed loops have been destroyed by large fields and extrapolation. Thus we fit the larger diamagnetism found after correcting for this effect.

2. FITTING DATA WITH TWO TYPES OF SUPERCONDUCTING LOOPS

Let us introduce the notation $`H_R`$ for the field at which resistance appears in the superconducting channels of which the loops are made, either $`H_c`$ for a Type I superconductor, or $`H_{c1}`$ for a Type II superconductor with no pinning. Let us define a field

$$H_B=0.5H_R.$$
(1)
Then Shoenberg’s type of theory  shows that, for a circular loop of superconductor of radius $`b`$ with a circular cross section of material of radius $`a`$, initially in zero field, the magnetisation commences to vary linearly with applied field $`H_e`$ perpendicular to the plane of the loop to a value

$$m_A=\frac{\pi ab^2H_B}{1+La/\pi b^2}$$
(2)
when the field reaches a value $`H_A`$ given by

$$H_A=\frac{H_B(La/\pi b^2)}{1+La/\pi b^2},$$
(3)
where $`L`$ is the inductance of the loop. The inductance is given by

$$L=4\pi bX,$$
(4)
where, for supercurrents confined to the surface of the material of the loop,

$$X=[\mathrm{ln}(8\mathrm{b}/\mathrm{a})2],$$
(5)
whereas for supercurrents through the bulk of the material of the loop,

$$X=[\mathrm{ln}(8\mathrm{b}/\mathrm{a})7/4].$$
(6)
For, e.g. $`b/a`$ = 50, the difference in the two values of $`X`$ is 6%. For definiteness we shall do our calculations with the first value of $`X`$, which will be correct if the magnetic-field penetration depth is small compared with $`a`$. The initial susceptibility $`\chi `$ associated with the loop is

$$\chi =\pi b^3/4X.$$
(7)
As the applied field is increased above $`H_A`$, the magnitude of the magnetic moment decreases linearly with field until the field $`H_B`$ is reached, at which field the magnetic moment is

$$m_B=\pi ba^2H_B.$$
(8)
At fields above $`H_B`$, ignoring a small discontinuity in moment at $`H_B`$ mentioned by Shoenberg but not discussed in detail, there is a linear decrease of the magnitude of the magnetic moment from the value given by Eq. (7) to zero at a field approximately equal \[to lowest order in $`(a/4b)`$\] to $`2H_B=H_R`$.

In order to fit the observed moments starting from zero applied field and continuing up to the field at which the metamagnetic transition occurs for the sample of  with the largest moment, we use a model with two types of loops with different $`b/a`$ ratios,

$$r_1=b_1/a>b_2/a=r_2,$$
(9)
with $`a`$ the same for both loop types, and suppose that there are $`N_1`$ and $`N_2`$ loops of each type. Such a model is used as a first approximation to a model with a semicontinuous distribution of loop parameters. We also assume that the planes of the loops are parallel to the film surfaces, i.e. perpendicular to the applied field. For the larger loops, geometrical consraints for a film of thickness 10 $`\mu `$m  force any closed loops to have approximately this orientation. If some of the smaller loops have other orientations, then probably the theory will still hold approximately with the reinterpretation of loop areas as the average projections of their areas on to planes parallel to the film surfaces.

The model with two types of loop parameters has four straight-line segments on the moment versus field curve, with discontinuities in slope at fields $`H_{A1}`$, $`H_{A2}`$, and $`H_B`$, where $`H_{A1}`$ and $`H_{A2}`$ are given by Eq. (3) for the two different values of $`b/a`$. The magnetic moments at the points of discontinuity of slope are given by

$$m(H_{A1})=N_1m_{A1}+N_2m_{A2}(H_{A1}/H_{A2}),$$
(10)
$$m(H_{A2})=N_1m_{A1}N_1(m_{A1}m_{B1})(\frac{H_{A2}H_{A1}}{H_BH_{A1}})+N_2m_{A2},$$
(11)
$$m(H_B)=N_1m_{B1}+N_2m_{B2},$$
(12)
with the magnetisation going to zero at $`H_R=2H_B`$. In Eqs. (10) to (12), the second suffices 1 and 2 on $`m_A`$ and $`m_B`$ refer to the two types of loops.

After correcting for a probable effect of superposition of diamagnetism associated with closed loops with a smaller positive moment associated with the majority of conducting channels which do not form closed loops, as discussed in the introduction, we fit the data for the more strongly magnetised sample with five adjustable parameters, viz. $`(b_1/a)`$, $`(b_2/a)`$ (assuming $`b_1>b_2`$), $`N_1m_{A1}`$, $`N_2m_{A2}`$, and $`H_B`$. We find that $`b_1/a=133`$, $`b_2/a=7.9`$, $`H_B=2631`$ Oe, implying $`H_{c1}=5262`$ Oe, and, with use of the fact that the volume of the film $`10^3`$ cm3 , that $`N_1m_{A1}=1.69\times 10^4`$ emu, and $`N_2m_{A2}=1.25\times 10^4`$ emu. Using Eqs. (2), (4) and (5) we deduce that $`N_2/N_1380`$, corresponding to a ratio of area covered by the smaller loops to that covered by the larger loops of about 1.3. We used a program AMOEBA given in a book  to perform the least squares fitting. The fit to the data up to the field of the metamagnetic transition ($`3360`$ Oe for this sample) is shown in Fig. 1. The rms accuracy of the fit is 2.2% of the mean value of the inferred diamagnetism, or 3.3% of the mean value of the net observed magnetisation.

For the second sample we assume that the corrections for a positive moment from the channels which do not form closed loops have values half of those for the first sample, based on the estimate in  that the average electron concentration in the second sample is about half that of the first sample. We use use the same type of model with two different values of $`(b/a)`$, but keep $`H_B`$ as before, and so we have a four-parameter fit. We find $`b_1/a=106`$, $`b_2/a=14.4`$, and, with the volume of the film as $`10^3`$ cm3 , $`N_1m_{A1}=0.35\times 10^4`$ emu, $`N_2m_{A2}=0.25\times 10^4`$ emu. From Eqs. (2), (4) and (5) we deduce that $`N_2/N_158`$, corresponding to a ratio of the area covered by the smaller loops to that covered by the larger loops of about 1.1. The rms accuracy of the fit is 2.4% of the mean value of the inferred diamagnetism, or 6.3% of mean value of the net observed magnetisation. The fit to the data up to the field of about 1390 Oe for the metamagnetic transition for this sample is also shown in Fig. 1. Since $`N_1`$ cannot be less than 1, we deduce from Eqs. (2), (4) and (5) and the value of $`N_1m_{A1}`$ that the radius $`a`$ of the cross sections of the channels forming the loops is less than 0.76 $`\mu `$m. This is in fair agreement with an upper limit of 1 $`\mu `$m estimated from other data . A lower limit for channel radii of 0.1 $`\mu `$m is estimated in , and a stricter upper limit of 0.35 $`\mu `$m is mentioned in , but with only a reference to a future publication (which has not appeared yet as far as I know) for an explanation of how this limit is obtained.

From the parameters obtained we find that the fraction of the first sample occupied by the material of the closed loops is $`7.9\times 10^5`$ for the set of smaller loops and $`0.35\times 10^5`$ for the set of larger loops, independent of the value of $`a`$. The total is sufficiently small compared with the probable total fractional volume of the order of a percent occupied by all conducting channels that we are justified in using approximately the same correction for positive magnetisation as that for all channels after the closed loops are broken, inferred from magnetisation measurements on the day after the original measurements . The contributions to the initial susceptibility from the larger and smaller loops are $`4.9\times 10^4`$ emu cm-3 and $`0.9\times 10^4`$ emu cm-3 respectively, corresponding to a total of 0.7% of that for a complete Meissner effect. The appreciably smaller fraction of superconducting material compared with the fraction of the full Meissner effect for the initial susceptibility arises because the susceptibility due to induced currents for a closed loop is larger by a factor of $`(1/2X)(b/a)^2`$ than that due to a Meissner effect keeping the flux completely out of the material of the loop. Although, for the model used, most of the material of the superconducting loops is associated with the smaller loops, the dominant contribution to the initial susceptibility comes from the larger loops because of the factor $`(b/a)^2`$.

For $`a=0.76\mu `$m, the maximum induced currents in the loops, at the fields $`H_{Ai}`$ (i=1,2), vary between 1.0 A and 1.7 A, depending on the loop size. These currents are considerably smaller than the critical currents of about 60 A through films in pulsed measurements with microprobe contacts of diameter 10 $`\mu `$m on the top surface of the films . Probably there was contact with only one conducting channel in these measurements, but in any case there could not have been many channels involved. We presume that the reason for the larger critical currents in the pulsed measurements is that any magnetic flux associated with the current does not have time to enter the channel, and so critical currents in this case are determined by other factors, and may equal depairing currents.

3. DISCUSSION

The use of two different values of $`(b/a)`$ is a simplification of a model with a quasi continuous distribution of $`(b/a)`$’s. With a quasi continuous distribution, rounding out of corners of the magnetisation curves would occur.

Since a transition temperature greater than room temperature is very high for a superconductor, we expect a small coherence length and Type II superconductivity. Also, because of the high temperature and softness of the material, pinning of flux may be difficult. Thus our model may be appropriate. A value of about 5260 Oe is inferred from the data for the lower critical field $`H_{c1}`$ which would exist for fields perpendicular to the superconducting channels if a metamagnetic transition did not intervene at a lower field . At 4.2 K, high conductivity does not disappear for fields up to 9 tesla  (probably approximately parallel to the channel lengths), and so $`H_{c2}`$ at this temperature for some orientation of the field is greater than 9 T. Since $`T_c`$ has been estimated indirectly to be greater than 700 K in , the low-temperature critical fields may be close to those at room temperature. For the high-temperature oxide superconductor YBa2Cu3O7, $`H_{c1}`$ for fields perpendicular to the film planes is about 700 Oe at low $`T`$ . It would not be surprising to find much higher values of $`H_{c1}`$ in oxidised atactic polypropylene in view of the much higher $`T_c`$.

We have ignored interaction between current loops, and between current loops and the larger amount of positively magnetised material in channels not forming closed loops. Such interactions can on average be taken into account by demagnetisation fields. Since the maximum susceptibility is always below 1% of that corresponding to a full Meissner effect, such corrections can be expected to be represented by changes in the internal field with respect to the applied field by amounts equal to a fraction of a percent of the applied field.

We know of no non-superconducting material which can show diamagnetism as large as a few tenths of a percent of a full Meissner effect. For fields along the $`c`$-axis, graphite has a susceptibility of $`21.1\times 10^6`$ emu g-1 , which corresponds to 0.06 % of a complete Meissner effect, about an order of magnitude smaller than the initial susceptibility for the more magnetic of the two films discussed here.

Although ballistic transport of electrons in arrays of small current loops may give some magnetic properties similar to those of superconductors , the characteristic temperature $`T^{}`$ below which ballistic effects occur for a loop with $`N`$ electrons and radius $`R`$ is

$$k_BT^{}\mathrm{}^2N/2mR^2,$$
(13)
with $`m`$ the electron mass. Taking $`NR/a_0`$, where $`a_0`$ is the period of an assumed periodic system, this reduces to $`k_BT^{}(\mathrm{}^2/2m)(1/a_0R)`$. Assuming e.g $`R=1\mu `$m and $`a_0=0.5`$ nm, we find $`T^{}0.9`$ K. To move $`T^{}`$ up to room temperature, we would thus need to reduce $`R`$ to 3 nm, which would imply a radius of the material of the loop of $``$ 1 nm or less. Although we could obtain an initial diamagnetic susceptibility as high as 0.7% of a Meissner effect by postulating very large numbers of such tiny loops and using Eq.(7), such small sizes would be incompatible with what is known about conducting channels through the films, which have cross sections with radii of the order of 1 $`\mu `$m . Thus it appears that an interpretation of the diamagnetism in terms of ballistic effects in mesoscopic systems at room temperature is very unlikely.

While our model for the diamagnetism associated with superconducting loops has several arbitrary parameters, and may not be the only type of superconductor model which can explain the diamagnetism, we have shown that a plausible model to interpret the diamagnetic properties can be found based on concepts similar to those used to explain the strong evidence for superconductivity in channels from the three types non-magnetic experiments mentioned at the beginning of the paper.

One way to test the model discussed here would be to cycle the magnetic field in a suitable sample, keeping within a field range below that at which the metamagnetic transition occurs, and to compare the predictions for the complete cycle with Shoenberg’s theory . However, since only about 10% of the samples with conducting channels are reported to show large diamagnetism , finding a suitable film for such measurements may not be easy. The only recent reports of work on oxidised atactic polypropylene come from Professor Grigorov and members or ex members of his group, several of whom are now working for commercial companies in the USA, and from Shlimak and Martchenkov , of Bar-Ilan University, who have done more work on polydimethylsiloxane (PDMS), (which also shows some of the unusual properties found in OAPP ), but also some on oxidised atactic polypropylene. Besides having some room-temperature properties similar to those of OAPP, the Josephson effect has been observed for PDMS sandwiched between two superconductors at temperatures below $`T_c`$ of the superconducting contacts .

Another type of polymer which shows narrow channels through films with fairly high conductivity is poly(3,3’-phthalidylene-4,4-biphenylylene) (PPB) (see e.g. ). Although no suggestions that these channels are superconducting at high temperature have been published, resistance too low to be detected has been found through PPB films with Sn contacts at low temperatures . Josephson effects have been reported in polyimide . Further studies of both these materials, especially of magnetic properties, may be worthwhile.

Two other examples of what may be quasi one-dimensional systems with superconductivity at room temperature are (i) carbon deposits containing multiwalled nanotubes \[31-33\], and (ii) powdered mixtures of PbCO3.2PbO $`+`$ Ag2O . The structure of the superconducting components of the system studied in  has been suggested to contain well separated Ag-O chains which are thought to be the main channels for possible superconductivity in this system.

4. CONCLUSIONS

Fair fits to the magnetisation curves showing large diamagnetism at low fields at room temperature in two samples of films of oxidised atactic polypropylene have been obtained using a model involving superconducting current loops in applied magnetic fields below those at which metamagnetic transitions occur, and also below the field $`H_{c1}`$ at which resistance would appear assuming no pinning if a metamagnetic transition did not set in first. Two sizes of loops with the same radius $`a`$ of the cross section of the superconducting channels of the loops are assumed. For the second sample, $`a`$ has been estimated from the fit to be less than about $`0.8\mu `$m, in fair agreement with $`a<1\mu `$m estimated from other data .

ACKNOWLEDGMENT

I wish to thank Professor L.N. Grigorov for correspondence about the results given in .

REFERENCES

1. V.M. Arkhangorodskiĭ, A.N. Ionov, V.M. Tuchkevich, and I.S. Shlimak, Pis’ma Zh. Eksp. Teor. Fiz. 51, 56 (1990) \[JETP Lett. 51, 67 (1990)\].

2. O.V. Demicheva, D.N. Rogachev, S.G. Smirnova, E.I. Shklyarova, M.Yu. Yablokov, V.M. Andreev, L.N. Grigorov, Pis’ma Zh. Eksp. Teor. Fiz. 51, 228 (1990) \[JETP Lett. 51, 258 (1990)\].

3. L.N. Grigorov, O.V. Demicheva, S.G. Smirnova, Sverkhprovodimost’ (KIAE) 4, 399 (1991) \[Superconductivity, Phys. Chem. Tech. 4, 345 (1991)\].

4. A.V. Krayev, T.V. Dorofeeva, E.I. Shklyarova, L.N. Grigorov, 9th. CIMTEC - World Forum on New Materials, Florence, Italy, 14-19 June 1998; Advances in Science and Engineering Technology, Vol. 23: Science and Engineering of HTC Superconductivity, edited by P. Vincennzini (Faenza Techna., Srl., 1999), pp. 459-466.

5. S.G. Smirnova, O.V. Demicheva, L.N. Grigorov, Pis’ma Zh. Eksp. Teor. Fiz. 48, 212 (1988) \[JETP Lett. 48, 231 (1988)\].

6. N.S. Enikolopyan, L.N. Grigorov, S.G. Smirnova, Pis’ma Zh. Eksp. Teor. Fiz. 49, 326 (1989) \[JETP Lett. 49, 371 (1989)\].

7. D.N. Rogachev, L.N. Grigorov, J. Supercond. 13, 947 (2000). The straight dashed lines through the origin shown in Fig. 1 of this paper have smaller slopes than indicated by the labelling for $`\chi `$. The values for $`\chi `$ given are based on an assumed form of extrapolation of the observed moments to lower fields than those to which they were measured (L.N. Grigorov, personal communication).

8. L.N. Grigorov, D.N. Rogachev, A.V. Kraev, Vysokomol. Soedin. B 35, 1921 (1993) \[Polymer Science 35, 1625 (1993)\].

9. A.M. Elyashevich, A.A. Kiselev, A.V. Liapzev, G.P. Miroshnichenko, Phys. Lett. A 156, 111 (1991).

10. D.M. Eagles, Physica C 225, 222 (1994); erratum ibid. 280, 335 (1997).

11. D.M. Eagles, J. Supercond. 11, 189 (1998). On p. 192 of this paper, the magnetic moments of two samples inferred from published data were too large by a factor of $`4\pi `$, besides having some uncertainty because of lack of precise knowledge by me at that time of the volumes of the samples. This implies that the radii of loops required to fit the moments on the assumption of spontaneous currents close to the critical current can be smaller by a factor of the order of $`(1/4\pi )^{\frac{1}{2}}`$ than given in the reference. In fact, judging by the magnetisations given in the recent paper for a sample previously reported on in , the moments estimated in  for that sample are too large by another factor of two, which appears to have arisen because the low-field susceptibility given in  was, I infer from correspondence from Professor Grigorov about reference , based partly on an assumed extrapolation of the measured moment curves to lower fields than those for which measurements were made. The smaller moments required from those estimated in  give rise to the possibility of alternative explanations in terms of smaller induced currents as discussed in the present paper for the results of .

12. D.M. Eagles, Revista Mexicana de Fisica 45, Suplemento 1, 118-121 (1999).

13. L.N. Grigorov, Phil. Mag. B 78, 353 (1998).

14. L.N. Grigorov, 9th. CIMTEC - World Forum on New Materials, Florence, Italy, 14-19 June 1998; Advances in Science and Engineering Technology, Vol. 23: Science and Engineering of HTC Superconductivity, edited by P.Vincennzini (Faenza Techna. Srl., 1999), pp. 675-684.

15. L.N. Grigorov, D.N. Rogachev, Molec. Cryst. Liquid Cryst. 230, 625 (1993).

16. S.G. Smirnova, E.I. Shklyarov, L.N. Grigorov, Vyosokomol. Soedin. B 31, 667 (1989).

17. V.M. Arkhangorodskiĭ, E.G. Guk, A.M. El’yashevich, A.N. Ionov, V.M. Tuchkevich, I.S. Shlimak, Dokl. Akad. Nauk. SSSR 309, 603 (1989) \[Sov. Phys. Doklady 34, 1016 (1989)\].

18. D. Shoenberg, Superconductivity (Cambridge University Press, 1952) Sec.2.6. Note that in Fig. 12 of this reference, the arbitrary units of magnetic moment used appear to be negative, i.e. opposed to the applied field, as expected. This can be seen from Shoenberg’s equations.

19. W.H. Press, B.P. Flannery, S.A. Teukolsky, W.T. Vetterling, Numerical Recipes (Cambridge University Press, 1986).

20. A.V. Kraev, S.G. Smirnova, L.N. Grigorov, Vysokomol. Soedin. A 35, 1308 (1993). \[Polymer Science 35, 1308 (1993)\].

21. I.S. Shlimak, L.N. Grigorov, unpublished results, 1996, referred to in .

22. S. Senoussi, C. Aguillon, Europhys. Lett. 12, 273 (1990).

23. X.K. Wang, R.P.H. Chang, A. Patachinski, J.B. Ketterson, J. Mater. Res. 9, 1578 (1994).

24. M. Szopper, E. Zipper, Int. J. Mod. Phys. B 9, 161 (1995).

25. I. Shlimak, V. Martchenkov, Solid State Commun. 107, 443 (1998).

26. L.N. Grigorov, T.V. Dorofeeva, A.V. Kraev, D.N. Rogachev, O.V. Demicheva, E.I. Shklyarova, Vysokomol. Soedin A 38, (1996) 2011. \[Polymer Science A 38, (1996) 1328\].

27. A.N. Ionov, V.A. Zakrevskiĭ, Pis’ma Zh. Tekh. Fiz. 26, No. 20, p. 36 (2000). \[Tech. Phys. Lett. 26, 910 (2000)\].

28. V.M. Kornilov, A.N. Lachinov, Synth. Met. 53, 71 (1992).

29. V.A. Zakrevskiĭ, A.N. Ionov, A.N. Lachinov, Pis’ma Zh. Tekh. Fiz. 24, No. 13, p. 89 (1998). \[Tech. Phys. Lett. 24, 539 (1998)\].

30. A.N. Ionov, V.A. Zakrevskiĭ, I.M. Lazebnik, Pis’ma Zh. Tekh. Fiz. 25, No. 17, p. 36 (1999). \[Tech. Phys. Lett. 25, 691 (1999)\].

31. V.I. Tsebro, O.E. Omel’yanovskiĭ, A.P. Moravskiĭ, Pis’ma Zh. Eksp. Teor. Fiz. 70, 457 (1999). \[JETP Lett. 70, 462 (1999)\].

32. V.I. Tsebro, O.E. Omel’yanovskiĭ, Phys. Usp. 43, 847 (2001).

33. G.-M. Zhao, Y.S. Wang, cond-mat 011268 (2001); Phil. Mag. B, to be published.

34. D. Djurek, Z. Medunić, A. Tonejc, M. Paljević, Physica C 351, 78 (2001).

Figure Caption

Fig. 1. Comparison of model calculations with diamagnetism associated with closed loops inferred from experiment for applied fields $`H_e`$ less than the fields at which the metamagnetic transition occurs, i.e. $`H_e<`$ 3360 Oe for sample 1, and $`H_e<`$ 1390 Oe for sample 2. The diamagnetic contribution to the net magnetisation has been inferred from the observed net magnetisation by correction for a probable positive contribution from conducting channels which do not form closed loops (see text).


text:
# Acknowledgements

## Acknowledgements

The authors thank A. Pilaftsis and C.E.M. Wagner for useful correspondence. This work was supported by the ‘Fonds zur Förderung der wissenschaftlichen Forschung’ of Austria, Project No. P13139-PHY, and by the European Community’s Human Potential Programme under contract HPRN-CT-2000-00149. W.P. is supported by ‘Fonds zur Förderung der wissenschaftlichen Forschung’ of Austria, Erwin Schrödinger fellowship Nr. J2095, and partly by the ‘Schweizer Nationalfonds’.

Figure Captions

Figure 1: Contours of the $`\stackrel{~}{\tau }_1`$ decay branching ratios $`B(\stackrel{~}{\tau }_1\tau \stackrel{~}{\chi }_1^0)`$ (a,b,c) and $`B(\stackrel{~}{\tau }_1\nu _\tau \stackrel{~}{\chi }_1^{})`$ (d,e,f) in the $`|A_\tau ||\mu |`$ plane for $`\mathrm{tan}\beta =3`$, $`M_2=200`$GeV, $`m_{\stackrel{~}{\tau }_1}=240`$GeV, $`m_{\stackrel{~}{\tau }_2}=250`$GeV, $`\phi _1=0`$, and ($`\phi _{A_\tau }`$, $`\phi _\mu `$) = (0,0) (a,d), ($`\pi /2`$,0) (b,e), and (0,$`\pi /2`$) (c,f) in the case $`m_{\stackrel{~}{\tau }_L}<m_{\stackrel{~}{\tau }_R}`$. The blank areas are excluded by the conditions (i) to (iii) given in the text and the inequality $`(m_{\stackrel{~}{\tau }_L}^2m_{\stackrel{~}{\tau }_R}^2)^2=(m_{\stackrel{~}{\tau }_1}^2m_{\stackrel{~}{\tau }_2}^2)^2(2|a_\tau m_\tau |)^20`$. The excluded region is different in each plot mainly due to the fact that the inequality depends on the phases ($`\phi _{A_\tau }`$, $`\phi _\mu `$).

Figure 2: Contours of the $`\stackrel{~}{\tau }_1`$ decay branching ratios $`B(\stackrel{~}{\tau }_1\tau \stackrel{~}{\chi }_1^0)`$ (a), $`B(\stackrel{~}{\tau }_1\tau \stackrel{~}{\chi }_2^0)`$ (b), and $`B(\stackrel{~}{\tau }_1\nu _\tau \stackrel{~}{\chi }_1^{})`$ (c) in the $`\phi _{A_\tau }\phi _\mu `$ plane for $`\mathrm{tan}\beta =3`$, $`M_2=200`$GeV, $`m_{\stackrel{~}{\tau }_1}=240`$GeV, $`m_{\stackrel{~}{\tau }_2}=255`$GeV, $`|A_\tau |=600`$GeV, $`|\mu |=350`$GeV, and $`\phi _1=0`$ in the case $`m_{\stackrel{~}{\tau }_L}<m_{\stackrel{~}{\tau }_R}`$.

Figure 3: $`\phi _\mu `$ dependence of the $`\stackrel{~}{\tau }_1`$ decay branching ratios for $`\mathrm{tan}\beta =3,30`$ (with $`\phi _1=0`$) (a) and $`\phi _1=0,\pi /2`$ (with $`\mathrm{tan}\beta =3`$) (b) with $`M_2=300`$GeV, $`m_{\stackrel{~}{\tau }_1}=240`$GeV, $`m_{\stackrel{~}{\tau }_2}=500`$GeV, $`|A_\tau |=600`$GeV, $`|\mu |=200`$GeV, and $`\phi _{A_\tau }=0`$ in the case $`m_{\stackrel{~}{\tau }_L}<m_{\stackrel{~}{\tau }_R}`$. In Fig.a the solid and dashed lines are for $`\mathrm{tan}\beta `$ = 3 and 30, respectively. In Fig.b the solid and dashed lines are for $`\phi _1`$ = 0 and $`\pi /2`$, respectively.

Figure 4: $`\phi _{A_\tau }`$ dependence of the $`\stackrel{~}{\tau }_2`$ decay branching ratios for $`\phi _\mu =0`$ (solid lines) and $`\pi /2`$ (dashed lines) in two cases: (a) small $`\stackrel{~}{\tau }`$-mass-splitting case with $`m_{\stackrel{~}{\tau }_1}=240`$GeV, $`m_{\stackrel{~}{\tau }_2}=260`$GeV, $`\mathrm{tan}\beta =6`$, $`M_2=500`$GeV, $`|A_\tau |=600`$GeV, $`|\mu |=150`$GeV, $`\phi _1=0`$, and $`m_{\stackrel{~}{\tau }_L}m_{\stackrel{~}{\tau }_R}`$ (for which $`m_{\stackrel{~}{\nu }_\tau }m_{\stackrel{~}{\tau }_2}`$), and (b,c) large $`\stackrel{~}{\tau }`$-mass-splitting case with $`m_{\stackrel{~}{\tau }_1}=240`$GeV, $`m_{\stackrel{~}{\tau }_2}=500`$GeV, $`\mathrm{tan}\beta =30`$, $`M_2=400`$GeV, $`|A_\tau |=900`$GeV, $`|\mu |=800`$GeV, $`\phi _1=0`$, and $`m_{\stackrel{~}{\tau }_L}<m_{\stackrel{~}{\tau }_R}`$ (for which $`m_{\stackrel{~}{\nu }_\tau }m_{\stackrel{~}{\tau }_1}`$). In the latter case (Figs.b and c) we have $`(m_{H_1^0},m_{H_2^0},m_{H_3^0})=(117.4,159.2,159.5)`$(GeV) and $`(117.7,150.6,151.2)`$(GeV) for $`\phi _\mu =0\text{and}\pi /2`$, respectively. The branching ratios $`B(\stackrel{~}{\tau }_2\stackrel{~}{\tau }_1Z^0)(10\%)`$ and $`B(\stackrel{~}{\tau }_2\stackrel{~}{\nu }_\tau W^{})(20\%)`$, which are rather insensitive to $`\phi _{A_\tau }`$, are not shown in Figs.b and c.

Figure 5: $`\phi _\mu `$ dependence of the $`\stackrel{~}{\nu }_\tau `$ decay branching ratios for $`\phi _1=0`$ (solid lines) and $`\pi /2`$ (dashed lines) with $`\mathrm{tan}\beta =3`$, $`M_2=500`$GeV, $`m_{\stackrel{~}{\tau }_1}=240`$GeV, $`m_{\stackrel{~}{\tau }_2}=500`$GeV, $`|A_\tau |=600`$GeV, $`|\mu |=150`$GeV, $`\phi _{A_\tau }=0`$, and $`m_{\stackrel{~}{\tau }_L}<m_{\stackrel{~}{\tau }_R}`$ (for which $`m_{\stackrel{~}{\nu }_\tau }m_{\stackrel{~}{\tau }_1}`$).

Fig.2

Fig.3


text:
# Phase diagram of dilute polyelectrolytes: Collapse and redissolution by association of counterions and coions

## I Introduction

Dilute solutions of strongly charged polymer electrolytes (polyelectrolytes) such as single and double stranded DNA , polystyrene sulphonate , and other biological and synthetic polymeric systems  exhibit phase and/or structural transitions upon addition of multivalent salts to the solutions. First, the addition of multivalent salt causes a collapse of the extended structure of flexible polyelectrolytes or the bundling in mono- or multi- molecular structures of semiflexible and rigid electrolytes that segregate into a dense phase, a phenomenon which in the case of DNA is also know as DNA condensation. Upon further addition of salt, the electrolytes redissolve acquiring again extended or unbundled states. There are many different examples of this type of system and it is clear that while the strong charge of the electrolytes is required for the transitions, the precise geometry of the polymers is not a crucial factor. The fact that the transitions are driven only by a combination of electrostatic interactions with entropic and steric effects has been shown by many computer simulations.

The first transition to a collapsed state is explained by the onset of attractive interactions between almost neutral structures formed by electrolytes whose charge is almost compensated by associated multivalent counterions.  The redissolution transition indicates that as the concentration of multivalent salt increases the system finds more energetically favorable conformations without packaging its structure. The model presented in this paper emerges from two previously proposed ideas regarding the redissolution transition. In a previous work  it was argued that a determining factor for the redissolution is the nature of the salty environment into which the electrolyte redissolves. In particular, it was emphasized that at large multivalent salt concentrations the thermodynamics of the multivalent counterions and monovalent coions is affected by the Bjerrum association  where the ions of the salt within the solution become strongly correlated with each other. Thus, an scenario for the redissolution is that the multivalent ions associated with the electrolyte prefer, without leaving the polymer, to be exposed to such environment. On the other hand, it has been argued by the group of Shklovskii  that the main reason for the redissolution is the decrease in energy available to the system by overcharging the electrolytes . Since the collapsed structure has to be effectively neutral, the electrolytes redissolve so that more counter-ions can associate to them. Both these ideas can be combined and the result provides an interesting picture of different possible charged states for the redissolved electrolyte.

The central point of this paper to exhibit, by means of a simple model, that for large amounts of multivalent and monovalent salt added to the system, the redissolution transition is marked by the onset of reassociation of the multivalent salt ions on the bulk of the solvent. In its redissolved state the electrolyte acquires a large number of multivalent counter-ions that overcompensate its bare charge but, also, associated with them, a large number of coions that form a second associated layer. Due to the large amount of screening of the electrostatic interaction near redissolution conditions, the decrease in energy due to the reassociation of coions with the multivalent counterions at the electrolyte surface is of magnitude similar (but smaller) to the energy of reassociation in the bulk of the solution. Thus, this second layer of associated coions can only exist if Bjerrum association is already present in the bulk of the solution. This second layer enhances the decrease in energy of the electrolyte due to the overcharging of the first layer. The presence of the coions reduces the net charge of the region near the electrolyte and thus decreases the Coulombic barrier to multivalent ions overcharge. Furthermore, in such a state, there is a greater number of particles in strongly correlated states thus decreasing the energy per electrolyte and making its redissolution more favorable. The total charge of the near region of the electrolyte that contains all these associated charges can have either charge sign, but in all cases the charge from the multivalent ions overcompensates the bare charge of the electrolyte. This multiple-species association has been observed in recent computer simulations by Tanaka and Grosberg  where the electrolyte is a small colloidal sphere and by Holm and Kremer  where the electrolyte is a rigid rod . The possibility of this second associated layer has been known for long time  but its implications for these phase transitions has not been fully explored. Figure 1 presents an sketch of different possible charged states of the near region of an electrolyte. Consider, for example, a section of electrolyte with total bare charge $`12`$, and multivalent counter-ions of charge $`z=+3`$. At low amounts of multivalent-salt the total charge of the system is negative, but partially compensated by (a) positive monovalent or (b) multivalent ions. The collapsed state (c) is neutral with the electrolyte charge compensated by positive multivalent ions, where different electrolytes or sections of electrolytes share the positive charges. In the redissolved state the charge of the electrolyte is overcompensated by multivalent ions but the total associated charge can be (d) positive, (e) neutral, or (f) negative due to the association of a second layer of negative ions.

A second goal of this article is to present a simple construction of the phase diagram of the multi-species electrolyte system. While conditions for collapse or redissolution implicitly assume some form of phase transition its global structure has not been previously considered in detail. This article proposes, on the basis of simple phenomenological observations, a construction of the structure of the phase diagram, the location of the coexisting phases, and the slopes of tie-lines.

In the rest of the paper, all free energies are measured in units of $`kT`$, with $`k`$ the Boltzmann constant and $`T`$ the temperature. All concentrations are presented as number densities. Lengths are measured with respect to a microscopic length unit $`a`$ that is the average size of monomers of the polymer and ions, which are assumed to be of roughly the same size. The strength of the electrostatic interaction is quantified by the Bjerrum number $`B=e^2/\epsilon akT`$ which is energy of interaction of two elementary charges at its nearest possible distance in water (or the medium considered).

## II Structure of the Phase Diagram.

### A Boundaries of the phase diagram

The basic components of the system are: a solvent, a dilute polyelectrolyte with monomer concentration $`\varphi `$, monovalent salt $`s`$ and a multivalent salt $`\rho `$. Since each of these components dissociate into different ions, we actually start with six different ionic components in solution. To reduce this number, it is possible to consider the case in which the charges of the ions are $`1`$ for each monomer of the polymer, $`+1`$ for the original counterions of the polymer, $`1`$ and $`+1`$ for the ions of the monovalent salt, and $`+z`$ and $`1`$ for the ions from the multivalent salt. It is further assumed that all of the $`+1`$ ions are identical, and similarly for the $`1`$ ions. Thus there are only four components, besides the solvent: dissociated polymer at concentration $`p=\varphi `$, positive multivalent ions $`m=\rho `$, positive monovalent ions $`t=\varphi +s`$, and negative monovalent ions $`u=zm+s`$. A further constraint is that each phase formed has to be electroneutral so that the condition $`p+zm+tu=0`$ is always satisfied. These restrictions allow for a description of all the ionic species in terms of only three concentrations, either of salts $`(\varphi ,\rho ,s)`$, or ions ($`p`$,$`m`$,$`t`$).

The allowed values of the concentrations of ions can be described as follows. Each ionic particle density should be non-negative, and in the reduced representation $`(p,m,t)`$ the requirement $`u0`$ becomes the constraint $`tp+zm0`$. Carrying this constraint to the salt representation shows that $`(\varphi ,\rho ,s)`$ should satisfy

$`0`$ $``$ $`\varphi ,`$ (1)
$`0`$ $``$ $`\rho ,`$ (2)
$`\mathrm{min}(z\rho ,\varphi )`$ $``$ $`s.`$ (3)
This last inequality has the following interpretation. Consider the case of a state with a small amount of solvent and a dense mixture of only polymer ions with multivalent ions, so that $`(p,m,t)=(x,x/z,0)`$ (and $`u=0`$). This state is physically realizable but its description in terms of the salt variables is $`(\varphi ,\rho ,s)=(x,x/z,x)`$. To obtain this state one can start with a mixture of polyelectrolyte and multivalent salt and then extract monovalent salt, that is, extract both positive monovalent counterions (given by the polyelectrolyte), and negative monovalent ions (given by the salt). The boundaries of the possible states of the system in the salt variables representation are sketched in Figure 2.

### B Location of the dense phase

Experiments have shown the phase segregation of neutral complexes with little amount of solvent from the bulk of the dilute polyelectrolyte solution. The precise location of this dense phase in the salt or ionic components diagram has not been precisely determined experimentally. Here it is only necessary for our purposes to assume that this phase has zero or very small amount of monovalent ions present and that is rich in polyelectrolyte neutralized with multivalent ions with a very small amount of solvent. The dense phase is not a unique point in the phase diagram, but it can be described as a region in the neighborhood of the point$`(p,m,t)=(\mathrm{\Phi },\mathrm{\Phi }/z,0)`$,or $`(\varphi ,\rho ,s)=(\mathrm{\Phi },\mathrm{\Phi }/z,\mathrm{\Phi })`$,where $`\mathrm{\Phi }`$, is the (large) concentration of monomers in this phase and, as explained above, the negative value of $`s`$ reflects the depletion of monovalent ions from these states. To simplify the model, however, it can be assumed that the dense phase consists of only one point, namely $`(\varphi ,\rho ,s)=(\mathrm{\Phi },\mathrm{\Phi }/z,\mathrm{\Phi })`$.

### C The tie lines

The precipitation of the collapsed neutral structures from the dilute solution indicates that the system has entered into a coexisting phase. A point $`(\varphi ,\rho ,s)`$ inside the coexistence region decomposes into two coexisting points with compositions $`(\mathrm{\Phi },\mathrm{\Phi }/z,\mathrm{\Phi })`$, namely the dense state, and a second polymer dilute point of the form $`(\varphi x\mathrm{\Phi },\rho x\mathrm{\Phi }/z,s+x\mathrm{\Phi })`$. The assumption that most of the high polymer density points of the region of coexistence lie in a very small region and are modeled as one single point immediately gives rise to this specific decomposition, and it is clear that the slope of the tie-lines is always determined by this construction.

To summarize, the phenomenological observation that the dense polymer phase is neutral with the polyelectrolyte charge neutralized by multivalent salts leads to a phase diagram of the form shown in Figure 3 including the tie lines shown. (The choice of cross section for different values of the $`t`$ variable are meant to aid visualization of the diagram and do not reflect actual energetic considerations. An example of actual cross-sections in logarithmic variables is presented in the Figure 4).

## III Free Energy Functional

To be able to write down an effective free energy functional that describes the polyelectrolyte phase transition it is necessary to consider coherently but separately the contributions to the free energy from three different regions. There is a dense precipitated phase that is relatively poor in solvent and rich in polymer and multivalent ions. Within the dilute phase, there are two distinct regions: the one occupied by the polymer and its associated charges, and the region away from the polymer occupied only by solvent and ions.

There are many variants of the problem all of which present similar behavior. The same type of phase diagram occurs for flexible and rigid polyelectrolytes of different lengths. It is simplest to write a model for a collection of rigid rods (say short or medium size DNA). In this case the main barrier arising from the structure of the polyelectrolyte to the formation of a collapsed state is purely entropic. If one considers very long semiflexible molecules, the collapse is monomolecular and the barrier is of a elastic nature instead of entropic, and requires the introduction of further parameters. The results are nevertheless similar.

### A Two States modeling

While it is clear that the region of phase coexistence does not separate two topologically disconnected regions in the full phase diagram, it is nevertheless useful to write down two separate free energy functionals for the states at the ends of a tie line. It is very hard to write down one single free energy functional that interpolates between the dense and dilute phase, especially for flexible electrolytes. When using different functionals for each phase, a phase transition can be properly determined only if both functionals are good approximations of the absolute free energy and refer to a same (arbitrary) zero energy point. In this case this goal is achieved by using similar approximations in the calculation of the free energy in the two states.

### B The Dense Phase

In the dense phase, the charges are very close to each other although it is likely that they retain a large amount of water. This phase is electroneutral and effectively forms a bulk region where the all the charges are strongly correlated with their near neighbors. In previous work  it was argued that a suitable free energy per unit volume $`F_d`$ for this state, is given by:

$$F_d(\mathrm{\Phi },\mathrm{\Phi }/z,\mathrm{\Phi })=M\mathrm{\Phi }.$$
(4)
where $`M`$ in an effective Madelung constant that depends on the specific geometry of the polymer and ions used but that reflects the strongly correlated state of the charges in that phase.

### C The dilute region

The dilute phase is not homogeneous at a microscopic level. The charge distribution is very different near the polyelectrolyte compared to regions away from it. The free energy of this phase can be approximated by separating a region near the polymer from the rest of the solvent. Since the polymer concentration is very dilute, this region can be considered as a separate thermodynamic system in equilibrium with an ion reservoir. The amount of charge per monomer associated to the polymer from each of the ionic species is $`f_m`$, $`f_u`$, $`f_t`$. The average ionic concentrations in the bulk of the solution excluding regions occupied by polymer is then $`(m,t,u)_{eff}=(mf_mp/z,tf_tp,uf_up)`$.

The free energy density for the system polyelectrolyte plus solvent can be written as:

$$F_t=\frac{\varphi }{N}\mathrm{ln}\frac{\varphi }{N}+\varphi F_e+F_s,$$
(5)
where $`F_t`$ is the total free energy density in the dilute phase, $`F_e`$ is the energy per electrolyte monomer due to the interactions between its own and associated charges, and $`F_s`$ is the energy of the solvent an the ions in solution. The values of the associated fractions $`f_i`$ are determined by minimization of this expression with respect to them.

The salty environment can be modeled as an ionic solution obeying a Debye-Huckel free energy where it is necessary to consider not only the multivalent positive, monovalent positive and negative ions but also an extra species, the reassociated ions. At the typical concentration at which the redissolution of the electrolyte takes place, there are important corrections to the free energy of the solution due to this phenomenon. The simplest way to describe the association is to assume that a finite fraction $`f_b`$ of all the multivalent ions in solution are associated with the same number of positive monovalent ions thus forming a new species of reduced charge $`(z1)`$ with concentration $`m_a=f_bm_{eff}`$. While in general there can be further reassociation to states of charge $`z2`$ etc., only the first reassociation will be treated explicitly here. Upon reassociation, the magnitude of the energy of a pair of associated ions is $`g_b`$. The new concentrations of non associated ions are $`m_n=m_{eff}m_a`$, and $`u_n=u_{eff}m_a`$. The free energy for the solvent is then:

$$F_s=\psi _i\mathrm{ln}\psi _i\frac{1}{12\pi }\kappa ^3g_bm_a$$
(6)
where $`\kappa `$ is the screening length of the system and is given by

$$\kappa ^2=\underset{i}{}4\pi z_i^2\psi _i,$$
(7)
where the ionic densities $`\psi _i`$ with charges $`z_i`$ are $`m_a`$, $`m_n`$, $`t_{eff}`$ and $`u_n`$.

The free energy of the near region of the polyelectrolyte, per monomer of electrolyte can be expressed as an expansion up to quadratic order on the fractions of associated charges. A good choice for this functional is

$$F_e=\frac{G_2}{2}q^2g_mf_mg_tf_tg_uf_u+g_cf_mf_u.$$
(8)
To argue for such functional it is convenient to separate within the near region of the electrolyte, the contributions to the free energy coming from a small region near each positive multivalent or monovalent ion, and the interaction of between these ions with others further away. The long range screened interaction between charges gives rise to the quadratic contribution $`G_2q^2/2`$, where the effective charge per monomer is $`q=1f_m+f_uf_t`$. The coefficient of interaction depends on the geometry of the system and for simplicity is chosen to have the form $`G_2=G_1B\mathrm{ln}((1+\kappa )/\kappa )`$, where $`G_1`$ is a geometric factor. For large screening lengths $`\kappa ^11`$, this term reflects the interaction of each charge with the logarithmic field of a charged rod of length $`\kappa ^1`$. At very small screening lengths $`\kappa ^1<1`$ the function goes to $`0`$.

When the effective electrolyte charge is near zero the linear terms in the free energy reflect the correlation energy of each charge with its near neighbors. At zero effective charge, the amount of energy gained from the addition of a unit charge from a given ion is simply $`g_m`$, $`g_s`$, or $`g_t`$. The association of a positive monovalent charge has as prerequisite the presence of multivalent negative ions associated to the electrolyte, so that $`f_u<f_m`$. The reduction in energy for this association is required to be comparable (but smaller) to the energy of association for these pairs of ions in the bulk so that $`g_u<g_b`$ the positive charge is brought into contact with a positive associated ion but it is also in the neighborhood of the negatively charged polyelectrolyte.

The reassociation of the multivalent ions and negative coions influences the strong correlations between associated multivalent ions. As a larger number of monovalent ions enter the near region, the multivalent ions reduce its correlations since their repulsion is diminished by the presence of the negative ions associated with them. This element is reflected in the correlation reduction term $`g_c`$. More such terms describing corrections to the correlations can be written but it turns out that the inclusion of just this one is sufficient to recover a phase diagram that is similar to those observed in experiments. Omission of this term leads to much more extended regions of coexistence.

## IV Determination of Coexistence Regions

### A Equilibrium between phases

A state with concentrations $`(p,m,t)`$ decomposes into two states if the following free energy, for some value of $`x`$, $`0x1`$, is smaller than the total energy of that state.

$$F_{comp}=xF_d+F_t(\varphi x\mathrm{\Phi },\rho x\mathrm{\Phi }/z,s+\mathrm{\Phi })F_t(\varphi ,\rho ,s).$$
(9)
The boundary of the coexistence region is determined by points at which the condition above is an equality and the required value of $`x`$ is zero. Expansion of $`F_{comp}`$ with respect to $`x`$ leads to

$$M\mu <0$$
(10)
where the effective chemical potential pre monomer for the addition of a neutral electrolyte complex into the dilute region is

$$\mu =_\varphi F_t+(1/z)_\rho F_t_sF_t.$$
(11)
On the basis of the assumption made, the final condition for the determination of the coexistence boundary does not explicitly depend on the concentrations in the dense phase given by $`\mathrm{\Phi }`$, but only on the bulk energy per electrolyte charge $`M`$ and the slope of the tie-line. The specific value of $`\mathrm{\Phi }`$ is, nevertheless, important to determine the condensed electrolyte fraction $`x`$.

### B Collapse

In a typical experiment performed to study the phase transitions of this type of system, the amount of multivalent salt present is changed from $`0`$ to some maximum value in small increments. The first transition occurs (for low monovalent salt) when the amount of charge provided by the multivalent salt is approximately the same as the amount of polyelectrolyte. If there is a large amount of monovalent salt, the transition requires a larger amount of multivalent salt. Thus, it is possible to approximate the location of the transition as follows: the condition Eq.9 is satisfied almost as soon as full charge compensation of the electrolyte by multivalent ions occur. The chemical potential $`\mu _m`$ of the ions in the solution, disregarding the Debye Huckel contribution, is $`\mathrm{ln}(mf_m\varphi )`$ with $`f_m1`$. Thus, for low monovalent salt concentrations the condition is simply that there must be enough ions to compensate the bare charge of the electrolyte:

$$zm\varphi .$$
(12)
At larger monovalent salt concentrations, the multivalent ions must first displace associated positive monovalent ions from the electrolyte. The reduction in the energy from the association of the monovalent and multivalent ions is, per charge, approximately the same and thus the chemical potential in the solvent for the multivalent ions should be smaller than that of the monovalent ions. The condition for the charging of the electrolyte with multivalent ions, and a first approximation to the collapse transition condition is:

$$\frac{1}{z}\mathrm{ln}m\mathrm{ln}(s+\varphi ).$$
(13)
### C Redissolution

In solutions where even at large salt concentrations the free energy of the ions is dominated by their translational entropy and there is little correlation between charges, the redissolution transition occurs when the decrease in energy per electrolyte due to multivalent ions overcharge is larger than the bulk correlation energy of the dense phase. However, if many monovalent ions are present, the onset of association of the positive monovalent with negative multivalent ions marks a sharp decrease in the energy of the electrolyte. Not only there is a decrease in the energy of the positive charges but also it is possible to associate many more ions. The repulsive interaction is decreased as the near region can be near neutral while having acquired large number of positive ions.

To first approximation, the onset of the association of positive ions occurs when the magnitude of the translational entropy of the free monovalent ions is comparable to the energy of association $`g_b`$. At larger monovalent ion concentrations, the charging of the electrolyte by positive multivalent ions together with monovalent coions occurs when the energy of association $`g_b`$ is capable of compensating the translational entropy of the monovalent ions. Thus we have

$$\mathrm{ln}u=\mathrm{ln}(zm+s)g_b.$$
(14)
This expression provides a rough approximation for the location of the second transition. For fixed amounts of added monovalent salt, the transition location is approximately a horizontal line in the ($`p`$,$`m`$) or ($`\varphi `$, $`\rho `$) plane.

### D  A numerical example.

The results of a numerical minimization of the model free energy appear in Figure 3. There it is shown, for different values of the monovalent salt concentration, the predicted regions of coexistence in the $`(\varphi ,\rho )`$ plane. The parameters used for this example are $`B=3.0`$, $`g_m=9.0`$, $`g_u=9.0`$, $`g_t=9.6`$, $`g_b=1.5`$, $`g_c=1.8`$, $`M=9.05`$. The curve containing points $`a`$, and $`b`$ has the lowest amount of added monovalent salt $`s=\mathrm{1.8\hspace{0.17em}10}^5`$ (intended to represent a 0.01 M molar concentration), the other two curves are for concentrations $`s=\mathrm{9.0\hspace{0.17em}10}^5`$, and $`s=1.810^4`$ (0.05 M, and 0.1 M). The largest salt concentration leads to smaller coexistence regions. The marked pairs of points $`(a,b)`$, $`(c,d)`$, and $`(e,f)`$ correspond to the location of the pairs of collapse-redissolution transitions that occur when for fixed amounts of monovalent salts and fixed amounts of polymer, multivalent salt is added to the system.

In this example, at the boundary of the redissolved region for the different values of monovalent salt added, the states $`d`$, $`e`$ and $`f`$ are all overcharged with positive multivalent ions but the net charge of all these states is negative. For the choice of parameters made, all states of the system that are not decomposed (lying outside the coexistence curve for fixed amount of monovalent salt) have a total negative charge. Different choices of parameters can lead to different charged states and, for example, in systems where the energy of reassociation of the multivalent salt is very small, all the dissociated states should have a positive effective charge, and the phase diagram has the shape proposed in work of Nguyen an Shklovskii , where the lines of redissolution and collapse do not join to form a bounded coexistence region. The model presented here also predicts that there is a maximum value of monovalent salt for which there is a coexistence region. As in the scheme of Figure 2, the coexistence region has a limited extent in the $`s`$, $`t`$ or $`u`$ axes.

For the case of largest amount of monovalent salt added, (with coexistence curve containing the points $`(e,f)`$), the dotted line divides the region of states that, while effectively negative, have enough multivalent ions to overcompensate the electrolyte bare charge. Separating curves for other salt concentration show similar behavior, always originating near the ”tip” of the diagram. The coexistence region sections are closed regions but their ”left” boundaries lie at very small concentrations that are typically experimentally inaccessible, and are not shown in Figure 4.

In extremely dilute solutions there is no large difference between the original concentration of ions in the bulk of the solution and the effective concentrations obtained by subtracting the associated ions. As the amount of polymer increases, the depletion from the bulk reduces the amount of available ions and states with similar amounts of associated ions appear only at higher multivalent salt concentrations. This leads to an upward curving trend in the redissolution curves.

## V Conclusions and Future work

It is still an important but unfinished problem to refine the description of the layer of associated ions by both providing a concrete definition of the layer and by testing possible (mean field) models of the free energy contribution of the layer to the total energy of the system. In the context of integral equation approaches this goal can be achieved by locating specific properties of the ion distribution functions , but ideally such definition ought to be model independent and experimentally useful. Nevertheless, even with arbitrary definitions of the layer, a free energy expansion of the form Eq. 8 is a useful instrument for the analysis of this type of system.

One would like to reduce the large number of parameters used in the present description of the system as many properties of real systems cannot be precisely measured or have a poor experimental definition. To use the ideas presented here in the context of the real system requires the collapse of information about the interaction between molecules and monomers into a few simple constants. These quantities reflect the average radii of ions and are supposed to include the steric effects of water as well as the specific location of charges in the electrolytes. A careful parametrization of simulation results, on the basis of the charge of the near region of an electrolyte was, for example, presented in the work of Messina et al. . In model systems studied though simulations or integral equations is possible to reduce the number of these parameters and calculate them on the basis of simpler geometric parameters. For comparison with experimental data it is, however, preferable to retain these constants as adjustable parameters.

A couple experimental tasks are suggested by this work as important to improve our understanding of these systems: the presentation of the data available from the point of view of phase separation of ionic species (instead of salts), a more systematic mapping of the effective charge of the electrolytes in the redissolved region of the phase diagram, and determination of the ionic composition of the associated layer (in this respect there has been progress by means of osmotic pressure measurements ). Future analysis of experimental data using the model presented in this article should provide a map, even if only a rough one, of the properties of specific ionic species with the parameters used. Finally, the addition of further ingredients that can directly modify the chemical potential of the different charged species should also induce dramatic changes in the location a of the transitions. From the point of view of the approach presented here, these modifications would be reflected in changes of the phenomenological parameters used in the model.

## Acknowledgements

The author thanks many useful talks with Monica Olvera de la Cruz, Eric Raspaud, Boris Shklovskii, Toan Nguyen, Motohito Tanaka, Mark Stevens, and John Widom.