text:
# Complexation between oppositely charged polyelectrolytes: beyond the Random Phase Approximation

## 1 Introduction

Most polymer blends are incompatible and phase separate in a wide range of their phase diagram into two phases each containing essentially one of the components of the mixture. The physical origin of the segregation of the polymers is their low mixing entropy of order $`1/N`$, where $`N`$ is the number of monomers per chain, which cannot overcome any weak enthalpic interaction. A slight dilution of the two polymers in a solvent does not change the qualitative picture of the phase separation, except that the two phases in coexistence are diluted. The chemical mismatch between two polymers is measured by a Flory-Huggins parameter $`\chi _{AB}`$ that gives the contribution to the free energy $`F_{enthalpic}=kT\chi _{AB}\varphi _A\varphi _B`$ where the monomer concentrations of polymers $`A`$ and $`B`$ are denoted by $`\varphi _A`$ and $`\varphi _B`$. For most polymer pairs, $`\chi _{AB}`$ is positive leading to incompatibility. There exists however some cases where $`\chi _{AB}`$ is negative and where the polymers are compatible and form a single homogeneous phase. A way to enhance the compatibility in solution is to add some electrical charges on the polymers. The simplest case is a mixture between a polyelectrolyte and a neutral polymer. In addition to the ionic charges along the polyelectrolyte backbone, the solution contains small counterions that insure the electroneutrality. Their entropy is sufficient to overcome the enthalpic interaction between the backbones since it is proportional to the monomer concentration independent of the polymer length: a phase separation is unfavorable as the counterions loose too much entropy by being confined in one of the phases. One can also directly take advantage of the electrostatic interaction by charging both polymers with ions of opposite charges. Such experiments have been performed by Djadoun et al. by mixing acidic and basic copolymers. If the number of charged groups along the copolymers is large enough, a single homogeneous phase is observed. The corresponding ternary phase diagram is sketched in Figure 1.a. It is qualitatively well described by the theory of Khokhlov et al. for low charge densities along the polymer. For higher polymeric charges, the shape of the phase diagram is completely different, and is sketched in Figure 1.b: the two polymers form complexes leading to a polymer-solvent type of phase separation. In this case the Flory parameter is no longer relevant since the phase behavior is governed by electrostatics.
Polyelectrolyte complexes have been studied extensively for thirty years . Nevertheless a clear theoretical description is not yet available. A so-called ”host” model has been proposed by Nordmeier et al.; the complexation is described as an ion binding between the polymers. This model is also invoked in a number of experiments  and should be realistic for highly charged semi-flexible polyelectrolytes. The first attempt to describe complexes between fully flexible weakly charged polyelectrolytes is due to Borue et al.. In this approach, the complex is stabilised by attractive interactions induced by charge fluctuations. We have used this model recently to describe the complexation between adjacent layers in polyelectrolyte multilayers .
Our purpose in this paper is to describe with the same model the transition from an incompatibility phase diagram (Fig.1a) to a complexation phase diagram (Fig.1b) when the charge density along the polymer chains is increased. The general formalism is introduced in section 2. It relies on a field theoretical formulation of the Random Phase Approximation. In section 3, we apply the method to a simple electrolyte as an example of the new results brought by the theory. Then it is applied to a charged polymer mixture where the two types of chains show chemical mismatch. We discuss the calculated structure factors and the phase diagram of the solution in section 4.2. In the last section, we present some concluding remarks and discuss some possible issues.

## 2 General theoretical framework

Our aim in this section is to introduce the formalism that we use throughout the paper. We introduce a Hamiltonian $`H[\varphi (𝐫)]`$ that can be written as a functional of the local field $`\varphi `$ (below $`\varphi `$ is the local concentration or the composition of the polymer mixture). The thermodynamic properties of this system are given by the partition function :

$$Z[h(𝐫)]=D\varphi (𝐫)\mathrm{exp}[\frac{H[\varphi (𝐫)]𝑑𝐫h(𝐫)\varphi (𝐫)}{kT}]$$
(1)
where we introduce an external field $`h(𝐫)`$. For the rest of this paper we use $`kT`$ as our energy unit ($`kT=1`$) to simplify the notations. The $`n`$-point correlation functions of the field are:

$$G(𝐫_1,𝐫_2,\mathrm{}𝐫_n)\varphi (𝐫_1)\varphi (𝐫_2)\mathrm{}\varphi (𝐫_n)=\frac{1}{Z[0]}\frac{\delta ^nZ[h(𝐫)]}{\delta h(𝐫_1)\delta h(𝐫_2)\mathrm{}\delta h(𝐫_n)}|_{h(𝐫)=0}$$
(2)
In this paper we focus on the two-point correlation function (the structure factor). It can be calculated as a derivative of the thermodynamic potential $`\mathrm{\Gamma }[\overline{\varphi }(𝐫)]`$ which is the Legendre transform of the free energy $`F[h(𝐫)]=\mathrm{log}Z[h(𝐫)]`$:

$$\mathrm{\Gamma }[\overline{\varphi }(𝐫)]F[h(𝐫)]+𝑑𝐫h(𝐫)\overline{\varphi }(𝐫)$$
(3)
We call $`\overline{\varphi }(𝐫)`$ the thermodynamic average of the field $`\varphi (𝐫)`$. The correlation function is obtained as:

$$G^1(𝐫_1,𝐫_2)=\frac{\delta ^2\mathrm{\Gamma }[\overline{\varphi }(𝐫)]}{\delta \overline{\varphi }(𝐫_1)\delta \overline{\varphi }(𝐫_2)}|_{\overline{\varphi }(𝐫)}$$
(4)
For any realistic Hamiltonian, the calculation of the path integral involved in Eq.1 cannot be performed exactly. The simplest approximation is to compute the partition function at the saddle point of its integrand: this is the mean field approximation. It neglects the contribution of all other paths to the integral in the partition function. This approximation leads to rather poor results when correlation effects are important. One can improve the approximation by taking into account the contribution of the fluctuations around the saddle point path at the gaussian level. This is the so-called Random Phase Approximation (RPA). We present here a standard route to derive the RPA thermodynamic potential.
The field can be decomposed into two parts $`\varphi (𝐫)=\varphi _{MF}(𝐫)+\delta \varphi (𝐫)`$ where the mean field part stems from the minimisation equation:

$$h(𝐫)=\frac{\delta H[\varphi (𝐫)]}{\delta \varphi (𝐫)}|_{\varphi _{MF}}$$
(5)
Upon expansion of the Hamiltonian around the saddle point at second order in the fluctuations of the field, the partition function is rewritten as:

$$Z[h(𝐫)]=Z_{MF}[h(𝐫)]D\delta \varphi (𝐫)\mathrm{exp}[\frac{1}{2}𝑑𝐫_1𝑑𝐫_2\delta \varphi (𝐫_1)G_{RPA}^1(𝐫_1,𝐫_2)\delta \varphi (𝐫_2)]$$
(6)
where we introduce the mean field partition function and the mean field two-point correlation function:

$`Z_{MF}[h(𝐫)]`$ $`=`$ $`\mathrm{const}.\mathrm{exp}[H[\varphi _{MF}(𝐫)]+{\displaystyle 𝑑𝐫h(𝐫)\varphi _{MF}(𝐫)}]`$
$`G_{RPA}^1(𝐫_1,𝐫_2)`$ $`=`$ $`{\displaystyle \frac{\delta ^2H[\varphi (𝐫)]}{\delta \varphi (𝐫_1)\delta \varphi (𝐫_2)}}|_{\varphi _{MF}(𝐫)}`$ (7)
Note that $`G_{RPA}^1`$ depends a priori on the mean field path $`\varphi _{MF}`$. One can then calculate the RPA thermodynamic potential:

$$\mathrm{\Gamma }[\overline{\varphi }(𝐫)]=H[\overline{\varphi }(𝐫)]\mathrm{log}\left[D\delta \varphi (𝐫)\mathrm{exp}[\frac{1}{2}𝑑𝐫_1𝑑𝐫_2\delta \varphi (𝐫_1)G_{RPA}^1(𝐫_1,𝐫_2)\delta \varphi (𝐫_2)]\right]$$
(8)
with $`\overline{\varphi }(𝐫)=\frac{\delta F[h(𝐫)]}{\delta h(𝐫)}`$. Using Eq.4, the inverse two-point correlation function, or the inverse structure factor reads:

$`G^1(𝐫_1,𝐫_2)=G_{RPA}^1(𝐫_1,𝐫_2)`$
$`{\displaystyle \frac{\delta ^2}{\delta \overline{\varphi }(𝐫_1)\delta \overline{\varphi }(𝐫_2)}}\left\{\mathrm{log}\left[{\displaystyle D\delta \varphi (𝐫)\mathrm{exp}[\frac{1}{2}𝑑𝐫^{}𝑑𝐫^{\prime \prime }\delta \varphi (𝐫^{})G_{RPA}^1(𝐫^{},𝐫^{\prime \prime })\delta \varphi (𝐫^{\prime \prime })]}\right]\right\}`$
The RPA formalism is of common use in polymer physics when one has to deal with rather dense systems where concentration fluctuations are known to be weak,. One can calculate both the thermodynamic potential and the structure factors. The free energy is useful to discuss macroscopic phase separation, but the knowledge of the structure factor is necessary if one wants to study properties at a finite lenth scale such as a micro-phase separation. For this purpose, several theoretical approaches only consider the first term of Eq.2, which corresponds to a mean field approximation for the inverse structure factor: fluctuation effects are neglected. In most cases, this is sufficient to discuss the phase separations of interest.
However for some coulombic systems, correlation effects might be relevant. Consider for example the case of a simple monovalent electrolyte. We show in the next section that neglecting fluctuations in the RPA calculation of the structure factor, which is the spirit of the Debye-Hückel theory of electrolytes (DH), leads to a constant density-density structure factor that does not satisfy all the sum rules imposed by the electroneutrality. This observation is at the origin of the Generalized Debye-Hückel (GDH) theory of Lee and Fisher which produces a $`q`$-dependent structure factor .
One can qualitatively understand why a mean field approximation misses some relevant properties in coulombic system: because of the electroneutrality condition, no pure coulombic term can appear in a mean field thermodynamic potential in the absence of external field, since fluctuations around the electroneutral state are not included.
The RPA calculation of the thermodynamic potential is also called the one loop approximation, since it corresponds to a resumation of all the one loop diagrams of the associated field theory for this potential. In the following, we call the result of Eq.2 the one loop approximation (OLA) of the inverse structure factor, to distinguish it from the standard RPA approximation, that includes only the zeroth order term in a loop expansion.

## 3 Structure factor of a simple electrolyte

We consider in this section a simple monovalent electrolyte solution in water composed of discrete point-like ions embedded in a continuous dielectric background of permittivity $`ϵ`$. For a sake of simplicity, we neglect the effect of the hard core potential since we only want to show the effect of the concentration fluctuations and not to discuss quantitatively the phase diagram. Simple electrolytes have been investigated quite recently by Netz et al.  and Frusawa et al.  independently using field theoretical methods. The N-particle problem has been mapped exactly onto a field theory with two fields which correspond to collective variables: the first one is a density field and the second one is an auxiliary field which acts as a complex potential on the ions. At the saddle point of the auxiliary field, the partition function is written $`Z=Dn_+(𝐫)Dn_{}(𝐫)\mathrm{exp}[H[n_+(𝐫),n_{}(𝐫)]]`$ with the Hamiltonian:

$`H[n_+(𝐫),n_{}(𝐫)]`$ $`=`$ $`{\displaystyle }d𝐫\{n_+(𝐫)(\mathrm{log}n_+(𝐫)1)+n_{}(𝐫)(\mathrm{log}n_{}(𝐫)1)`$
$`+{\displaystyle \frac{l_B}{2}}{\displaystyle }d𝐫^{}{\displaystyle \frac{(n_+(𝐫)n_{}(𝐫))(n_+(𝐫^{})n_{}(𝐫^{}))}{|𝐫𝐫^{}|}}\}`$
The Bjerrum length is defined here by $`l_B=\frac{e^2}{4\pi ϵ}`$; it measures the strength of electrostatic interactions. This hamiltonian can be used as the starting point of a density functional theory where the equilibrium density fields result from a balance between entropy and electrostatics. At this level, we neglect the fluctuations of the auxiliary field and investigate the effect of composition fluctuations. We show below that we recover this way the results of the Generalized Debye-Hückel theory if we include the composition fluctuations at the level of a one loop approximation, and the classical results of the Debye-Hückel theory if we neglect composition fluctuations. Let $`n_\sigma (𝐫)=\overline{n}_\sigma (𝐫)+\delta n(𝐫)`$ with $`\sigma =+,`$. Then using the definition of section 2, the partition function is:

$$Z=Z_{MF}D\delta n_+(𝐫)D\delta n_{}(𝐫)\mathrm{exp}[\mathrm{\Delta }H[\overline{n}_+,\overline{n}_{},\delta n_+,\delta n_{}]]$$

where the effective hamiltonian is defined as

$$\mathrm{\Delta }H[\overline{n}_+,\overline{n}_{},\delta n_+,\delta n_{}]=\frac{1}{2}\underset{\sigma ,\sigma ^{}}{}d𝐫d𝐫^{}\delta n_\sigma (𝐫)G_{RPA\sigma ,\sigma ^{}}^1(𝐫,𝐫^{})\delta n_\sigma ^{}(𝐫^{})]$$

The RPA inverse structure factor matrix is given by:

$$G_{RPA\sigma ,\sigma ^{}}^1(𝐫,𝐫^{})=\left(\begin{array}{cc}\frac{\delta (𝐫𝐫^{})}{\overline{n}_+(𝐫)}+\frac{4\pi l_B}{|𝐫𝐫^{}|}& \frac{4\pi l_B}{|𝐫𝐫^{}|}\\ \frac{4\pi l_B}{|𝐫𝐫^{}|}& \frac{\delta (𝐫𝐫^{})}{\overline{n}_{}(𝐫)}+\frac{4\pi l_B}{|𝐫𝐫^{}|}\end{array}\right)$$
(11)
One can easily generalize Eq.2 to a multi-component system to obtain the inverse structure factor matrix in the one loop approximation:

$`G_{\sigma ,\sigma ^{}}^1(𝐫,𝐫^{})`$ $`=`$ $`G_{RPA\sigma ,\sigma ^{}}^1(𝐫,𝐫^{})`$ (12)
$`+{\displaystyle \frac{\delta ^2\mathrm{\Delta }H}{\delta \overline{n}_\sigma (𝐫)\delta \overline{n}_\sigma ^{}(𝐫^{})}}+{\displaystyle \frac{\delta \mathrm{\Delta }H}{\delta \overline{n}_\sigma (𝐫)}}{\displaystyle \frac{\delta \mathrm{\Delta }H}{\delta \overline{n}_\sigma ^{}(𝐫^{})}}{\displaystyle \frac{\delta \mathrm{\Delta }H}{\delta \overline{n}_\sigma (𝐫)}}{\displaystyle \frac{\delta \mathrm{\Delta }H}{\delta \overline{n}_\sigma ^{}(𝐫^{})}}`$
where the average is performed with the gaussian weight $`\mathrm{exp}[\mathrm{\Delta }H[\overline{n}_+,\overline{n}_{},\delta n_+,\delta n_{}]]`$ of the fluctuating variables $`\delta n_+,\delta n_{}`$. The average densities $`\overline{n}_+,\overline{n}_{}`$ are constant in a homogeneous phase. Note that the densities must be considered as constant only at the step of averages, after the formal derivations with respect to inhomogeneous densities; otherwise the contributions of the fluctuations vanishes.
It is convenient to study the Fourier transform of the structure factors. Each element of the matrix is then rewritten as:

$`G_{\sigma ,\sigma ^{}}^1(𝐪)=G_{RPA\sigma ,\sigma ^{}}^1(𝐪)`$
$`{\displaystyle \frac{1}{2\overline{n}_\sigma ^2\overline{n}_\sigma ^{}^2}}{\displaystyle \frac{d^3𝐪_1}{(2\pi )^3}G_{RPA\sigma ,\sigma ^{}}(𝐪_1)G_{RPA\sigma ,\sigma ^{}}(𝐪𝐪_1)}+{\displaystyle \frac{\delta _{\sigma ,\sigma ^{}}}{\overline{n}_\sigma }}{\displaystyle \frac{d^3𝐪_1}{(2\pi )^3}G_{RPA\sigma ,\sigma ^{}}(𝐪_1)}`$
where the gaussian average has been already performed. All the integrals involved in Eq.3 can be calculated analytically by using Feynman integral techniques. The diagonal terms are diverging at high $`q`$. This divergence can be removed by substracting the contribution of the neutral system (perfect gas). The structure factor matrix of the electrolyte solution reads:

$$G_{\sigma ,\sigma ^{}}^1(𝐪)=\left(\begin{array}{cc}\frac{1}{\overline{n}}+\frac{4\pi l_B}{𝐪^2}\mathrm{\Delta }G^1(𝐪)& \frac{4\pi l_B}{𝐪^2}\mathrm{\Delta }G^1(𝐪)\\ \frac{4\pi l_B}{𝐪^2}\mathrm{\Delta }G^1(𝐪)& \frac{1}{\overline{n}}+\frac{4\pi l_B}{𝐪^2}\mathrm{\Delta }G^1(𝐪)\end{array}\right)$$
(14)
The electroneutrality condition imposes the equality of homogeneous densities $`\overline{n}_+=\overline{n}_{}=\overline{n}`$. The correction term is then given by:

$$\mathrm{\Delta }G^1(𝐪)=\frac{2\pi l_B^2}{q}\mathrm{arctan}\left(\frac{q}{2\kappa }\right)$$
(15)
with the Debye-Hückel length defined by $`\kappa ^2=8\pi l_B\overline{n}`$. The effect of the fluctuation contribution can be drawn from the density-density structure factor that is calculated from the partial structure factors $`G_{\overline{n}_{tot},\overline{n}_{tot}}=G_{\overline{n}_+,\overline{n}_+}+G_{\overline{n}_{},\overline{n}_{}}+2G_{\overline{n}_+,\overline{n}_{}}`$.

$$G_{\overline{n}_{tot},\overline{n}_{tot}}(𝐪)=\frac{2\overline{n}}{12\overline{n}\mathrm{\Delta }G^1(q)}$$
(16)
If we neglect the fluctuations at all orders,the density-density structure factor is constant. The fluctuations induce a decrease of the structure factor for increasing wave vectors. Notice that the value value at $`q=0`$ is consistent with the Debye-Hückel polarization energy $`F_{pol}=\frac{\kappa ^3}{12\pi }`$. This is expected since the results of the OLA and the RPA are the same for the free energy.
This natural generalization of the Debye-Hückel theory to structure factors must be compared to the Generalised Debye-Hückel theory proposed by Lee and Fisher (. Starting from the same constatation on the behaviour of the density-density structure factor in the Debye-Hückel theory, they solve the Debye-Hückel equation for the electrostatic potential with modulated densities. When the amplitude of the modulation vanishes at the end of the calculation, they obtain a $`q`$-dependent structure factor. Notice that they take into account the finite radius of the ions that we neglect here. Nevertheless their procedure is very similar to ours. In fact, the correlation length that we can define from the density-density structure factor is exactly the same in the limit of vanishing ion radius. More precisely the $`q=0`$ susceptibility and the correlation length are in the limit of vanishing ion radius:

$`\chi ^1(0)=1{\displaystyle \frac{\kappa l_B}{4}},\xi ^2={\displaystyle \frac{\chi (0)l_B}{48\kappa }}`$ (17)
in agreement with the results of Lee and Fisher. It is clear from Eqs.16 and 17 that the ”One loop approximation” is valid as long as $`\frac{\kappa l_B}{4}<1`$. If it is not the case, one has to consider higher orders in the loop expansion. The perturbation parameter of the theory is $`\zeta =\frac{\kappa l_B}{4}`$.
We do not discuss further the properties of the structure factor because it goes beyond the scope of this paper. We note however that simple tests do not reveal any inconsistencies, in particular with the Stillinger-Lovett sum rules .

## 4 Complexation and incompatibility in mixtures of oppositely charged polyelectrolytes

### 4.1 Inverse structure factor matrix

In this section we generalise the results obtained in the previous section for simple electrolytes to polyelectrolyte complexes which are mixtures of polyelectrolytes of opposite charges. We use the structure factors to study the stability of the polyelectrolyte mixture and compare our results to the experimental work of Djadoun et al. mentioned in the introduction. In this work, the charge of acidic and basic copolymers is progressivley increased. The neutral backbones of those copolymers exhibit a strong incompatibility. Starting from an incompatibility phase diagram at zero charge, the compatibility range (solubilisation) increases with the charge and finally turns into a complexation phase diagram where the two copolymers are in the same neutral dense phase in equilibrium with almost pure water because of the strong electrostatic attraction between them.
In order to describe this system, we consider a mixture of two polyelectrolytes A and B having the same linear charge density $`f`$ but of opposite charges . The charge density is smeared out along the backbone and the solvent is modelled as a continuous dielectric background of permittivity $`ϵ`$. A concentration $`n_++n_{}=2n`$ of small point-like ions of opposite charges is added to the solution. We study only a symmetric system and the monomer concentrations are $`c_+=c_{}=c`$. Each chain is composed of $`N`$ monomers of size $`a`$. If the electrostatic interactions are switched off, the short range interactions between polymers are characterised by three Flory parameters $`\chi _{AS},\chi _{BS}`$ and $`\chi _{AB}`$. For a sake of simplicity, we will only consider both polymers are in a $`\theta `$ solvent and fix the values $`\chi _{AS}=\chi _{BS}=\frac{1}{2}`$. The last parameter $`\chi _{AB}\chi `$ characterises the chemical mismatch between the two polymers and is varied independently.

The total Hamiltonian of the solution can then be written as:

$`H[c_+(𝐫),c_{}(𝐫),n_+(𝐫),n_{}(𝐫)]`$ $`=`$ $`{\displaystyle }d𝐫\{n_+(𝐫)(\mathrm{log}n_+(𝐫)1)+n_{}(𝐫)(\mathrm{log}n_{}(𝐫)1)`$ (18)
$`+{\displaystyle \frac{a^2}{24}}\left({\displaystyle \frac{|c_+(𝐫)|^2}{c_+(𝐫)}}+{\displaystyle \frac{|c_{}(𝐫)|^2}{c_{}(𝐫)}}\right)`$
$`+{\displaystyle \frac{l_B}{2}}{\displaystyle 𝑑𝐫^{}\frac{\rho _z(𝐫)\rho _z(𝐫^{})}{|𝐫𝐫^{}|}}`$
$`+\chi c_+(𝐫)c_{}(𝐫)\}`$
The charge density is defined here as $`\rho _z(𝐫)=f(c_+(𝐫)c_{}(𝐫))+n_+(𝐫)n_{}(𝐫)`$. The first term of Eq.18 is the translational entropy of the small ions. The second term is the conformational entropy of the polymer chains in the limit of infinite chain length. In this limit the conformational entropy can be calculated from a ground state dominance approximation and it is given by the Lifshitz formula . In the limit of infinite chain length, the translational entropy of the polymer chains is negligibly small. The third term is the electrostatic interaction between the charged components of the solution. Finally the last term models the enthalpic interactions that give rise to the incompatibility between the polymers.

From this model Hamiltonian we want to calculate various structure factors in order to discuss the stability limits of the solution and its phase diagram. As for simple electrolytes, the calculation is performed at the level of the one loop approximation presented in section 2.
For a sake of simplicity, we include the enthalpic interaction between the polymers only at the mean field level, or equivalently at the RPA level in the calculation of the structure factor. We believe that this approach gives the correct physical result because we expect fluctuations to be dominated by the connectivity of the chains. This calculation can also be used as a starting point to study various systems involving gaussian polymer chains in external potentials (in order to take into account excluded volume or solvent effects for example). Note that the contribution of the Flory parameter to the fluctuation can also be computed but the calculations are lengthy and much more tedious.
We first ignore the short range interactions, $`\chi =0`$. In order to find the fluctuation contributions to the structure factor, let $`c_\sigma (𝐫)=\overline{c}_\sigma (𝐫)+\delta c_\sigma (𝐫)`$ and $`n_\sigma (𝐫)=\overline{n}_\sigma (𝐫)+\delta n_\sigma (𝐫)`$. By expanding the Hamiltonian up to second order in the composition fluctuations, we calculate the thermodynamic potential and the resulting inverse structure factor. Within the one loop approximation, the inverse structure factor reads:

$`G_{\rho ,\sigma ,\sigma ^{}}^1(𝐫,𝐫^{})=G_{RPA\rho ,\sigma ,\sigma ^{}}^1(𝐫,𝐫^{})`$ (19)
$`+{\displaystyle \frac{\delta ^2\mathrm{\Delta }H}{\delta \overline{\rho }_\sigma (𝐫)\delta \overline{\rho }_\sigma ^{}(𝐫^{})}}+{\displaystyle \frac{\delta \mathrm{\Delta }H}{\delta \overline{\rho }_\sigma (𝐫)}}{\displaystyle \frac{\delta \mathrm{\Delta }H}{\delta \overline{\rho }_\sigma ^{}(𝐫^{})}}{\displaystyle \frac{\delta \mathrm{\Delta }H}{\delta \overline{\rho }_\sigma (𝐫)}}{\displaystyle \frac{\delta \mathrm{\Delta }H}{\delta \overline{\rho }_\sigma ^{}(𝐫^{})}}`$
with the short-hand notation $`\overline{\rho }_\sigma =(\overline{c}_+,\overline{c}_{},\overline{n}_+,\overline{n}_{})`$. The quantity $`\mathrm{\Delta }H`$ is the second order term in the expansion of the Hamiltonian $`H`$. As in the previous section, the average is performed with the gaussian weight $`\mathrm{exp}[\mathrm{\Delta }H]`$ of the fluctuating densities. The average densities $`\overline{\rho }`$ can in the calculation of the averages be set equal to their uniform values in the solution, $`c`$ and $`n`$ respectively for the monomers and the small ions. The RPA inverse structure factor has been already calculated by several authors , it is given in reciprocal space by:

$$G_{RPA\rho ,\sigma ,\sigma ^{}}^1=\left(\begin{array}{cccc}\frac{q^2a^2}{12\overline{c}}+\frac{4\pi l_Bf^2}{q^2}& \frac{4\pi l_Bf^2}{q^2}& \frac{4\pi l_Bf}{q^2}& \frac{4\pi l_Bf}{q^2}\\ \frac{4\pi l_Bf^2}{q^2}& \frac{q^2a^2}{12\overline{c}}+\frac{4\pi l_Bf^2}{q^2}& \frac{4\pi l_Bf}{q^2}& \frac{4\pi l_Bf}{q^2}\\ \frac{4\pi l_Bf}{q^2}& \frac{4\pi l_Bf}{q^2}& \frac{1}{\overline{n}}+\frac{4\pi l_B}{q^2}& \frac{4\pi l_B}{q^2}\\ \frac{4\pi l_Bf}{q^2}& \frac{4\pi l_Bf}{q^2}& \frac{4\pi l_B}{q^2}& \frac{1}{\overline{n}}+\frac{4\pi l_B}{q^2}\end{array}\right)$$
(20)
We give below the values of the corrections to the RPA in Eq.19 for the various components of the inverse structure factor matrix.
The monomer-monomer components read:

$`\mathrm{\Delta }G_{c\sigma ,\sigma ^{}}^1=`$ (21)
$`{\displaystyle \frac{a^4}{576}}\{{\displaystyle \frac{2}{c_\sigma ^2c_\sigma ^{}^2}}{\displaystyle }{\displaystyle \frac{d^3𝐪_1}{(2\pi )^3}}(𝐪_1.(𝐪𝐪_1)𝐪^2)^2G_{RPAc\sigma ,\sigma ^{}}(𝐪_1)G_{RPAc\sigma ,\sigma ^{}}(𝐪𝐪_1)`$
$`+{\displaystyle \frac{48\delta _{\sigma ,\sigma ^{}}}{c_\sigma ^3a^2}}{\displaystyle }{\displaystyle \frac{d^3𝐪_1}{(2\pi )^3}}(𝐪^2+𝐪_1^2)G_{RPAc\sigma ,\sigma ^{}}(𝐪_1)\}`$
The correction term to the RPA for the salt-salt components of the inverse structure factor matrix are given by the same formula Eq.3 as for a simple electrolyte:

$`\mathrm{\Delta }G_{n\sigma ,\sigma ^{}}^1=`$
$`{\displaystyle \frac{1}{2\overline{n}_\sigma ^2\overline{n}_\sigma ^{}^2}}{\displaystyle \frac{d^3𝐪_1}{(2\pi )^3}G_{RPAn\sigma ,\sigma ^{}}(𝐪_1)G_{RPAn\sigma ,\sigma ^{}}(𝐪𝐪_1)}+{\displaystyle \frac{\delta _{\sigma ,\sigma ^{}}}{\overline{n}_\sigma ^3}}{\displaystyle \frac{d^3𝐪_1}{(2\pi )^3}G_{RPAn\sigma ,\sigma ^{}}(𝐪_1)}`$
Note however that the RPA structure factor is different from the simple electrolyte case.
Finally the corrections to the crossed salt-monomer terms of the inverse structure factor matrix are given by:

$$\mathrm{\Delta }G_{c_\sigma ,n_\sigma ^{}}^1=\frac{a^2}{24c_\sigma ^2n_\sigma ^{}^2}\frac{d^3𝐪_1}{(2\pi )^3}(𝐪_1.(𝐪𝐪_1)𝐪^2)G_{RPAc_\sigma ,n_\sigma ^{}}(𝐪_1)G_{RPAc_\sigma ,n_\sigma ^{}}(𝐪𝐪_1)$$
(23)
All the integrals involved in Eqs.21, 4.1 and 23 can be calculated using Feynman integrals techniques. The analytical expressions of these integrals are reported in appendix A.

We can therefrom calculate the inverse structure factor at the level of the one loop approximation. As the inverse structure factor matrix is symmetric, we only give the independent terms:

$`G_{OLAc_+,c_+}^1=G_{OLAc_{},c_{}}^1`$ $`=`$ $`{\displaystyle \frac{q^2a^2}{12c}}+{\displaystyle \frac{4\pi l_Bf^2}{q^2}}+{\displaystyle \frac{f^2}{c}}\alpha (q)f^4\beta (q)`$ (24)
$`G_{OLAc_+,c_{}}^1`$ $`=`$ $`{\displaystyle \frac{4\pi l_Bf^2}{q^2}}f^4\beta (q)`$ (25)
$`G_{OLAn_+,n_+}^1=G_{OLAn_{},n_{}}^1`$ $`=`$ $`{\displaystyle \frac{1}{n}}+{\displaystyle \frac{4\pi l_B}{q^2}}\delta (q)`$ (26)
$`G_{OLAn_+,n_{}}^1`$ $`=`$ $`{\displaystyle \frac{4\pi l_B}{q^2}}\delta (q)`$ (27)
$`G_{OLAc_+,n_+}^1=G_{OLAc_{},n_{}}^1`$ $`=`$ $`{\displaystyle \frac{4\pi l_Bf}{q^2}}f^2\gamma (q)`$ (28)
$`G_{OLAc_+,n_{}}^1=G_{OLAc_{},n_+}^1`$ $`=`$ $`{\displaystyle \frac{4\pi l_Bf}{q^2}}f^2\gamma (q)`$ (29)
The quantities $`\alpha ,\beta ,\gamma ,\delta `$ are related to the Feynman integrals of appendix A by:

$`\alpha (q)`$ $`=`$ $`{\displaystyle \frac{24\pi l_B}{a^2}}I_2(q)`$
$`\beta (q)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{48\pi l_B}{a^2}}\right)^2I_1(q)`$ (30)
$`\gamma (q)`$ $`=`$ $`{\displaystyle \frac{96(\pi l_B)^2}{a^2}}I_4(q)`$
$`\delta (q)`$ $`=`$ $`8(\pi l_B)^2I_3(q)`$
### 4.2 Phase diagram and discussion of the experimental results

We have obtained in the preceding section the inverse structure factor matrix in the particular basis $`\{c+,c_{},n_+,n_{}\}`$. In order to discuss the phase diagram of the system, it is more convenient to calculate the structure factor matrix for the variables $`\{c_{tot},c_z,n_{tot},n_z\}`$ denoting respectively the total monomer density ($`c_{tot}=c_++c_{}`$), the total monomer charge density ($`c_z=f(c_+c_{})`$), the total small ion density ($`n_{tot}=n_++n_{}`$) and the total small ion charge density ($`n_z=n_+n_{}`$). Indeed the density-density structure factor is essential to predict a phase separation such a complexation which is a polymer-solvent phase separation; it is associated to a singular behavior of this structure factor. On the other hand, the incompatibility between the chains, is a polymer-polymer phase separation, and is associated to a singularity in the charge-charge structure factor, since the polymer charge is the only way to distinguish the two types of chains within our model. Previous studies have focused on this last structure factor because it exhibited a singularity while the density-density structure factor was regular. We show below that the effect of the composition fluctuations is to induce an attraction between the chains that counterbalances in some parameter range the chemical mismatch between the chains. This effect was ignored in all the previous studies that stay at the RPA or mean field level for the structure factor.

The change of variables and the inversion of the matrix $`G_{OLA}^1(q)`$ can be performed analytically. At this point we add at the RPA level the Flory parameter $`\chi `$. The structure factor bloc matrix for the polymer variables then reads:

$`G_{OLAc_\sigma ,c_\sigma ^{}}(q)=`$ (33)
$`\left(\begin{array}{cc}\frac{1}{\frac{1}{2c}\left(\frac{q^2a^2}{12}+f^2\alpha (q)2f^4c\beta (q)\right)+\frac{\chi }{2}}& 0\\ 0& \frac{1}{\frac{1}{2f^2c}\left(\frac{q^2a^2}{12}+f^2\alpha (q)\right)+\frac{4\pi l_B}{q^2+\kappa ^2}\frac{\chi }{2f^2}}\end{array}\right)`$
where $`\sigma =(tot,z)`$. In this last equation, we have neglected higher order contributions obtained in the matrix inversion involving the square of correction terms since these orders in the loop expansion have been neglected from the very beginning.
If we neglect $`\alpha (q)`$ and $`\beta (q)`$ in the structure factor, we recover the results of Khokhlov et al: their discussion is limited to the charge-charge structure factor. We present here a similar analysis by omitting first the effect of the density-density correlations on the phase diagram. The phase diagram can be drawn in a very transparent and elegant way by using the dimensionless variables introduced first by Borue et al. So far, we have introduced two characteristic lengths: the Debye-Hückel screening length $`\kappa ^1`$ and the correlation length of a dense polyelectrolyte mixture $`q_{}^1`$ ($`q_{}^4=\frac{96\pi l_Bf^2c}{a^2}`$). The first length is associated to the screening of electrostatic interactions by the small ions and the second length is associated to the screening of the same interactions by the polymers. The enthalpic interactions between the neutral backbones of the chains are measured by the Flory parameter $`\chi `$ to which we can associate the mean field length $`\xi _\chi `$ defined by :

$$\xi _\chi ^2=\frac{12\chi c}{a^2}$$
(34)
The physical behaviour of the system is then governed by the following dimensionless ratios:

$`s`$ $`=`$ $`{\displaystyle \frac{\kappa ^2}{q_{}^2}}`$
$`t`$ $`=`$ $`{\displaystyle \frac{\xi _\chi ^2}{q_{}^2}}`$ (35)
The choice of $`q_{}^1`$ as the unit length is motivated by the fact that we are expecting the connectivities of the chains to play a dominant role.
The stability of the system against the polymer-polymer phase separation is ensured if $`1/G_{c_zc_z}(q)`$ is positive for all values of $`q`$. The spinodal line giving the stability limit is obtained when this condition is violated. If the first mode for which $`1/G_{c_zc_z}(q)`$ becomes negative is $`q_00`$, the solution shows a microphase separation. On the contrary, for $`q_0=0`$ there is macrophase separation. The partial phase diagram is shown in Figure 2. It does not take into account yet instabilities of the density-density structure factor. An analysis based on the stability of the homogeneous phase can only give the equation of the spinodal lines but it cannot be used to predict the precise morphology of the different phases after the transition. A complete phase diagram can only be obtained by calculating the free energy of the various possible phases and by looking for the phases that minimise it. In this paper we study only the stability limits of the homogeneous phase.
We checked by numerical evaluation of the structure factor with the exact formula that the correction term $`\alpha (q)`$ can be neglected. If $`s>1`$ the system is homogeneous for $`st<1`$, otherwise it is macroscopically phase separated: chains of opposite charge are segregated into two dilute phases. This is due to the incompatibility of the polymers. If $`s<1`$ it is possible to have microphase separation. The spinodal line for this transition is given by $`t=2s`$. The wave vector at the transition is $`q_0^2=q_{}^2\frac{ts}{2}`$. It is related to the wavelength of modulation of the mesophase and it tends to infinity as we approach the macroscopic spinodal line $`st=1`$. Notice that in the limit of high ionic strength ($`s1`$) which is the relevant case for m ost experiments, no microphase separation is expected.

The full phase diagram includes possible instabilities associated to the density-density structure factor. The numerical analysis of this structure factor shows that the correction term $`\alpha (q)`$ can again be neglected and that the term $`\beta (q)`$ can be approximated by a constant $`\beta (q)\beta (q=0)`$. The result of this analysis is that the composition fluctuations induce a macroscopic complexation between the oppositely charged polyelectrolytes when:

$$t<\frac{u}{(s+2)^{3/2}}$$
(36)
where we introduce the following dimensionless ratios:

$`u`$ $`=`$ $`{\displaystyle \frac{\xi _{mesh}}{q_{}^1}}`$ (37)
$`\xi _{mesh}`$ $`=`$ $`{\displaystyle \frac{3}{8\pi a^2c}}`$ (38)
The length $`\xi _{mesh}`$ is the mesh size of a gaussian transient network (gaussian chains in semi-dilute regime). This length is relevant to discuss the validity limits of the OLA or the RPA calculations. The main assumption of those approximations is that the chains are gaussian at all length scales. For simple polyelectrolyte solutions, this remains true as long as the electrostatic blob size <sup>1</sup><sup>1</sup>1The electrostatic blob size is the length scale at which the typical electrostatic energy of a subunit composed of $`\left(\frac{\xi _{el}}{a}\right)^2`$ becomes larger than the thermal energy. More precisely $`kT\frac{kTl_Bf^2}{\xi _{el}}\left(\frac{\xi _{el}}{a}\right)^4`$ is larger than the mesh size of the semi-dilute solution because the electrostatic interactions are screened by the transient network at large length scales. Up to irrelevant prefactors this condition is equivalent to $`u<1`$. In terms of polymer concentration this threshold is written $`ca^3\frac{l_B^{4/3}f^{2/3}}{a^{4/3}}`$. Therefore the approximation used in this paper is valid at moderately high concentration.
The phase diagram of the system is presented in Figure 3. Note that there is no microphase separation associated to the polyelectrolyte complexation. Although the variables used to describe the phase diagram look very physical, they are not well suited to describe the different phases in terms of concentration. In Figure 4 we translate the $`ts`$ diagram into a $`nc`$ diagram (all the other parameters being kept constant). The equations of the various spinodal lines are reported in Appendix B.
If we now come back to the experiments of Djadoun et al. discussed above, our model is able to describe qualitatively the transition from incompatibility to compatibility and finally complexation for increasing polymeric charges. Indeed the parameters $`t`$ and $`s`$ scale like:

$`t`$ $``$ $`{\displaystyle \frac{\chi c^{1/2}}{l_B^{1/2}af}}`$ (39)
$`s`$ $``$ $`{\displaystyle \frac{l_B^{1/2}na}{c^{1/2}f}}`$ (40)
so that the two parameters have the same dependence in the polymer charge density, except for prefactors. If we increase $`f`$, we move in the phase diagram on a line represented in Figure 3. Starting at low $`f`$, the solution is in the incompatibility range; increasing $`f`$ pushes it into the compatibility range and finally complexation takes place. A very rough model of the complex can be made by describing it as a collapsed polymer . This argument also shows that incompatibility becomes an irrelevant parameter once the complex is formed.

In real systems, the polymers are not symmetric in length and charge, but we believe that this will not affect significantly the general shape of our proposed phase diagram: the neutral dense phases should be similar in composition but we expect the supernatant to be composed of charged soluble species. The description of the supernatant will be the topic of a forthcoming paper.

## 5 Concluding remarks

We have proposed in this article description of polymer mixtures in a solvent that goes beyond the classical Random Phase Approximation and that allows to describe simultaneously in a certain range of parameters both the incompatibility between the two polymers and the formation of polyelectrolyte complexes when the polymers are charged. This one loop approximation has been used to calculate the structure factors. The Random Phase Approximation is commonly used to discuss phase stability in polymer systems but in the calculation of the structure factors, it misses important correlation effects. As pointed out by Haronska et al., the RPA is able to predict only repulsive electrostatic contributions to the mean field effective potential between polymers. Including fluctuations at the order of one loop allows for attractive effective interactions. Using a One Loop Approximation at the level of the structure factors, we were able to include the effect of the fluctuations for two systems: simple electrolytes and mixtures of incompatible oppositely charged polyelectrolytes. For simple electrolytes, we obtain a natural generalization of the Debye-Hückel theory which leads to results similar to the one proposed by Lee and Fisher. For polyelectrolytes blends, we were able to describe within a single model the experiments showing a transition from an incompatibility phase diagram to a complexation phase diagram. The transition towards polyelectrolyte complexation is predicted to take place in the range of validity of the OLA, namely in the limit of weakly charged chains
The main limitation of this work is that for high enough charge densities the one loop approximation breaks down: the system becomes strongly fluctuating and one should use another formalism. A model based on the formation of ion pairs seems appropriate in this case; whe hope to be able to set up such a model in a future publication.

Acknowledgments: We are grateful for discussions to I. Erukhimovitch (Moscow state University). This research was supported by the Deutsche Forschung Gemeinschaft through the Schwerpunkt program ”Polyelektrolyte”.

## Appendix A

### A.1 Exact expression for the monomer-monomer correction terms

By replacing the expression of the RPA structure factor into Eq.21, we obtain for the off-diagonal term:

$$\mathrm{\Delta }G_{c+,}^1=\frac{1}{2}\left(\frac{48\pi l_Bf^2}{a^2}\right)^2I_1(𝐪)$$
(41)
$`I_1(𝐪)`$ is a Feynamn-like integral defined by:

$`I_1(𝐪)=`$ (42)
$`{\displaystyle \frac{d^3𝐪_1}{(2\pi )^3}\frac{(𝐪_1.(𝐪𝐪_1)𝐪^2)^2}{𝐪_1^2(𝐪^2𝐪_1^2)}\frac{1}{[𝐪_1^2(𝐪_1^2+\kappa ^2)+q_{}^4][(𝐪𝐪_1)^2((𝐪𝐪_1)^2+\kappa ^2)+q_{}^4]}}`$
with the characteristic length $`q_{}^1`$ defined by $`q_{}^4=\frac{96\pi l_Bf^2c}{a^2}`$. It can be interpreted as the correlation length of a dense mixture of oppositely charged polyelectrolytes in the absence of small ions . By reducing the integrand of Eq.42 to simple fractions we can expressed the whole integral as a sum of Feynman diagrams. We give here only the expansion at low $`q`$ because the discussion on the possibility of macro- or microphase separation is relevant only in this range, while the exact formula is rather lengthy:

$$I_1(q)=\frac{1}{8\pi q_+q_{}(q_++q_{})^3}+q^2\frac{15q_+^4+75q_+^3q_{}+116q_+^2q_{}^2+75q_+q_{}^3+15q_{}^4}{96\pi q_+^3q_{}^3(q_++q_{})^5}$$
(43)
where the two wave vectors $`q_+,q_{}`$ verify $`q_+^2q_{}^2=q_{}^4`$ and $`(q_++q_{})^2=\kappa ^2+2q_{}^2`$.

The diagonal terms exhibit a divergence at high $`q`$ which can be removed by substracting the contribution of the neutral system as in section 3. The convergent term is given by:

$$\mathrm{\Delta }G_{c\sigma ,\sigma }^1=\frac{24\pi l_Bf^2}{a^2c}I_2(𝐪)\frac{1}{2}\left(\frac{48\pi l_Bf^2}{a^2}\right)^2I_1(𝐪)$$
(44)
with the new integral:

$`I_2(𝐪)=`$
$`{\displaystyle \frac{d^3𝐪_1}{(2\pi )^3}\frac{1}{𝐪_1^2(𝐪^2𝐪_1^2)}\left\{\frac{(𝐪_1.(𝐪𝐪_1)𝐪^2)^2}{(𝐪𝐪_1)^2((𝐪𝐪_1)^2+\kappa ^2)+q_{}^4}\frac{(𝐪^2+𝐪_1^2)(𝐪𝐪_1)^2(𝐪.𝐪_1)^2}{𝐪_1^2(𝐪_1^2+\kappa ^2)+q_{}^4}\right\}}`$
The low $`q`$ expansion of this last integral is given by:

$$I_2(q)=\frac{q^2}{6\pi q_+q_{}(q_++q_{})}$$
(46)
### A.2 Exact expression for the salt-salt correction terms

The diagonal and off-diagonal correction terms in the salt-salt matrix are equal after substraction of the unphysical divergence and the common value reads:

$$\mathrm{\Delta }G_{n\sigma ,\sigma ^{}}^1=8\pi ^2l_B^2I_3(q)$$
(47)
with the integral:

$$I_3(q)=\frac{d^3𝐪_1}{(2\pi )^3}\frac{𝐪_1^2(𝐪𝐪_1)^2}{[𝐪_1^2(𝐪_1^2+\kappa ^2)+q_{}^4][(𝐪𝐪_1)^2((𝐪𝐪_1)^2+\kappa ^2)+q_{}^4]}$$
(48)
Although we calculated exactly this integral we do not give here the result because it is lengthy and we do not need it for the discussion of the structure factor as it will appear clearly in the section 4.2.

### A.3 Exact expression for the salt-monomer correction terms

As in the salt-salt matrix, the diagonal and off-diagonal correction terms of the salt-monomer matrix are equal to:

$$\mathrm{\Delta }G_{cn\sigma ,\sigma ^{}}^1=\frac{96(\pi l_Bf)^2}{a^2}I_4(q)$$
(49)
with the following integral

$$I_4(q)=\frac{d^3𝐪_1}{(2\pi )^3}\frac{\left(𝐪^2𝐪_1.(𝐪𝐪_1)\right)}{[𝐪_1^2(𝐪_1^2+\kappa ^2)+q_{}^4][(𝐪𝐪_1)^2((𝐪𝐪_1)^2+\kappa ^2)+q_{}^4]}$$
(50)
As in the preceding subsection, we will not give the precise evaluation of the integral for the sake of clarity.

## Appendix B

Figure 4 represents the $`nc`$ translation of the $`ts`$ phase diagram found in the main text. We have redefined the small ion concentration by substracting the polymer counterions; it is thus such that at $`n=0`$ there are still some small ions, the respective counterions of the polycations and polyanions. With this convention, the macroscopic spinodal line given by $`st=1`$ is written as:

$$\stackrel{~}{n}=\frac{l_Bf^2}{\stackrel{~}{\chi }a}f\stackrel{~}{c}$$
(51)
with the following rescaled parameters:

$`\stackrel{~}{n}`$ $`=`$ $`na^2l_B`$ (52)
$`\stackrel{~}{c}`$ $`=`$ $`ca^2l_B`$ (53)
$`\stackrel{~}{\chi }`$ $`=`$ $`{\displaystyle \frac{\chi }{a^3}}`$ (54)
The spinodal corresponding to microphase separation ($`t=2s`$) is given by:

$$\stackrel{~}{n}=\frac{6^{1/2}f\stackrel{~}{c}^{1/2}}{\pi ^{1/2}}\left(f+\frac{3a\stackrel{~}{\chi }}{2\pi l_B}\right)\stackrel{~}{c}$$
(55)
Finally the spinodal of complexation ($`t=\frac{u}{(s+2)^{3/2}}`$) is given by the formula:

$$\stackrel{~}{n}=\frac{3l_B^{1/3}f^2}{2^{4/3}\pi ^{2/3}\stackrel{~}{\chi }^{2/3}\stackrel{~}{c}^{1/3}}\frac{\mathrm{2\hspace{0.17em}3}^{1/2}f\stackrel{~}{c}^{1/2}}{(2\pi )^{1/2}}f\stackrel{~}{c}$$
(56)
At fixed $`c,n`$ and $`\chi `$, the critical value of the charge density for polyelectrolyte complexation is given by the solution of the quadratic Eq.56 in $`f`$.

## Figure Caption

1a.Typical phase diagram for a solution of oppositely charged polyelectrolytes at very low charge densities. In a wide range of this phase diagram the polymers undergo a polymer-polymer phase separation, the neutrality being ensured by the small ions. The thick lines are the tie lines indicating coexisting phases.
1b. Typical phase diagram for a solution of oppositely charged polyelectrolytes with higher charge densities. In a wide range of this phase diagram there is coexistence of dense phases containing both polymers and dilute phases with almost no polymers: there is thus complexation between the polyelectrolytes. The thick lines are the tie lines.

Partial $`ts`$ phase diagram: the different phases are indicated with their relative boundaries. This phase diagram can be associated to the representation used in Figure 1a: the polymers are either incompatible or compatible (soluble).

Complete $`ts`$ phase diagram: the effects of composition fluctuations are taken into account, resulting in the possibility of complexation for low $`\chi `$. When the charge density $`f`$ is increased, the other parameters being imposed, the system follows the double arrow sketched on the phase diagram, from the incompatibility region, through the compatibility region and finally through the complexation region.

Complexation phase diagram in the concentration variables. The value of $`\stackrel{~}{\chi }`$ is fixed ($`\stackrel{~}{\chi }=0.001`$) and the complexation line is shown for various values of $`f`$ (from bottom to top $`f=1/50,1/20,1/10`$).


format: text

text:
# The Riemannium

## Abstract

The properties of a fictitious, fermionic, many-body system based on the complex zeros of the Riemann zeta function are studied. The imaginary part of the zeros are interpreted as mean-field single-particle energies, and one fills them up to a Fermi energy $`E_F`$. The distribution of the total energy is shown to be non-Gaussian, asymmetric, and independent of $`E_F`$ in the limit $`E_F\mathrm{}`$. The moments of the limit distribution are computed analytically. The autocorrelation function, the finite energy corrections, and a comparison with random matrix theory are also discussed.

According to the Riemann hypothesis, the complex zeros of the function $`\zeta (s)=n^s`$ (defined for $`\mathrm{Re}s>1`$, and by analytic continuation to the rest of the complex plane) are located on the critical line $`s_\mu =1/2+iE_\mu `$, with $`E_\mu `$ real. There is a natural connexion between this hypothesis and quantum mechanics that originates on a spectral interpretation of the complex zeros $`s_\mu `$. Indeed, a general strategy to prove the hypothesis was suggested by Hilbert and Pólya, who proposed to find an hermitian operator whose eigenvalues are precisely the $`E_\mu `$’s. This operator may be viewed as the quantization of an hypothetical classical dynamical system, which determines the “Riemann dynamics”. Although the classical Hamiltonian corresponding to the Riemann dynamics is not known, there are evidences for its existence. When the $`E_\mu `$ are interpreted as quantum eigenvalues, both their statistical properties  as well as a semiclassical interpretation of their density  indicate that the Riemann dynamics is fully chaotic and has no time-reversal symmetry. The semiclassical interpretation provides in fact much more detailed information (cf Eq.(2)), in particular concerning the classical periodic orbits. This leads to a somewhat paradoxical situation: although the “Riemann Hamiltonian” is not known, the detailed properties of its dynamics are. From the perspective of a dynamicist, the Riemann zeta thus offers the rare opportunity of a chaotic motion for which the relevant dynamical information is simple and explicitly known. Not to mention the impressive amount of existing numerical data on the complex zeros , as well as the number-theoretic background of $`\zeta (s)`$, also useful in the present context.

Many aspects of the Riemann dynamics have been investigated (see  for recent review articles). Here we explore new facets of the problem related to a different use of the complex zeros as suggested by physical analogies. We consider the location $`E_\mu >0`$ of the complex zeros of $`\zeta (s)`$ as the single-particle levels of a fermionic many-body system. In the mean field approximation, the ground-state total energy is obtained by filling the single-particle levels from the lowest Riemann zero up to a “Fermi energy” $`E_F`$. We are interested in the properties of such a Fermi gas. Following nuclear physics terminology, we refer to this “element” as the Riemannium.

At zero temperature, the properties of the Riemannium are described by the grand potential

$$\mathrm{\Omega }(E_F)=_0^{E_F}(EE_F)\rho (E)𝑑E=_0^{E_F}𝒩(E)𝑑E.$$
(1)
Here $`\rho (E)=_\mu \delta (EE_\mu )`$ is the spectral density of the complex Riemann zeros interpreted as quantum eigenvalues and $`𝒩(E)=\rho (E)𝑑E`$ its counting function. $`\mathrm{\Omega }(E)`$ corresponds to the sum of the imaginary part of the Riemann zeros satisfying $`0<E_\mu E_F`$, using $`E_F`$ as reference energy . It is therefore the total energy of the system. To calculate the grand potential we make use of the decomposition of the spectral density in smooth plus oscillatory parts, $`\rho =\overline{\rho }+\stackrel{~}{\rho }`$. The former, as a Weyl series, has an explicit asymptotic expansion for large $`E`$. The oscillating term is an interference sum over the prime numbers $`p=2,3,\mathrm{}`$

$$\stackrel{~}{\rho }(E)=\frac{1}{\pi }\mathrm{Re}\underset{p}{}\underset{r=1}{\overset{\mathrm{}}{}}\frac{\mathrm{log}p}{p^{r/2}}\mathrm{exp}(irE\mathrm{log}p).$$
(2)
The comparison of this equation with the semiclassical Gutzwiller trace formula for the spectral density of a dynamical system , $`\stackrel{~}{\rho }(E)=_{po}A_{po}\mathrm{cos}(S_{po}(E)/\mathrm{})`$ where $`po`$ are the classical periodic orbits, shows that each prime number labels an unstable periodic orbit of a fully chaotic system of action $`S_p=ET_p`$, period $`T_p=r\mathrm{log}p`$, Lyapounov stability $`\lambda _p=1`$, repetitions labeled by $`r`$, and $`\mathrm{}=1`$ . An unusual fact of the Riemann dynamics that plays an important role in what follows is the independence with respect to the energy $`E`$ of the periods and Lyapounov exponents (notice also the minus sign in front of Eq.(2)).

The spectral density can be integrated twice to obtain the smooth and oscillatory contributions to the grand potential, $`\overline{\mathrm{\Omega }}(E_F)=E_F^2\mathrm{log}(E_F/2\pi )/4\pi +3E_F^2/8\pi 7E_F/8\mathrm{log}E_F/48\pi +cte`$, and

$$\stackrel{~}{\mathrm{\Omega }}(E_F)=\frac{1}{\pi }\mathrm{Re}\underset{p}{}\underset{r=1}{\overset{\mathrm{}}{}}\frac{\mathrm{exp}(iE_Fr\mathrm{log}p)}{r^2p^{r/2}\mathrm{log}p}.$$
(3)
We are interested in the statistical properties of $`\stackrel{~}{\mathrm{\Omega }}`$ as a function of the Fermi energy $`E_F`$. From Eq.(3) we have $`\stackrel{~}{\mathrm{\Omega }}=0`$. The average is done over an energy window which is small compared to $`E_F`$ but contains several oscillations of $`\stackrel{~}{\mathrm{\Omega }}`$ (the typical scale of oscillation will be given below).

As we will now see, the Riemannium has very peculiar properties. The most important ones demonstrated here are: i) the distribution of $`\stackrel{~}{\mathrm{\Omega }}`$, denoted $`P(\stackrel{~}{\mathrm{\Omega }})`$, is independent of energy as $`E_F\mathrm{}`$; ii) $`P(\stackrel{~}{\mathrm{\Omega }})`$ is non-Gaussian and asymmetric; iii) all the moments of $`P(\stackrel{~}{\mathrm{\Omega }})`$ may be computed with very good accuracy from Eq.(3), their value being dominated by the contribution of the lowest prime numbers; iv) the asymmetry of the distribution is due to an interference effect between repetitions of periodic orbits; v) the autocorrelation $`\stackrel{~}{\mathrm{\Omega }}(E_F)\stackrel{~}{\mathrm{\Omega }}(E_F+ϵ)`$ is a non-decaying irregular oscillatory function; vi) the finite energy corrections to $`P(\stackrel{~}{\mathrm{\Omega }})`$ are universal and well described by a circular unitary ensemble of random matrices ($`CUE`$).

Some of the statistical properties of $`\stackrel{~}{\mathrm{\Omega }}`$ were explored in the past. Selberg  computed the even moments of the distribution. Up to an additive constant and a global change of sign, $`\stackrel{~}{\mathrm{\Omega }}(E)`$ coincides with the function denoted $`S_1(t)`$ by him. More recently, Odlyzko  has numerically calculated the first four moments, but no comparison has been made with the analytical results of Selberg.

On Fig. 1 is displayed the distribution $`P(\stackrel{~}{\mathrm{\Omega }})`$ whose properties are now discussed. We begin by computing the second moment. Equation (3) allows to express it as an integral over the period of the orbits

$`\stackrel{~}{\mathrm{\Omega }}^2=(1/2\pi ^2){\displaystyle _0^{\mathrm{}}}𝑑TK(T)/T^4,`$ (4)
$`K(T)=`$ $`{\displaystyle \underset{i,j}{}}A_iA_j\mathrm{cos}[E_F(r_i\mathrm{log}p_ir_j\mathrm{log}p_j)]\delta (T\overline{T}).`$ (5)
Use has been made here of the definition of the form factor, $`K(T)=4\pi _0^{\mathrm{}}𝑑ϵ\mathrm{cos}(ϵT)\stackrel{~}{\rho }(E_F)\stackrel{~}{\rho }(E_F+ϵ)`$, with the amplitudes defined as $`A_i=\mathrm{log}p_i/p_i^{r_i/2}`$, and $`\overline{T}=(T_{p_i}+T_{p_j})/2`$. Rigorous arguments valid for times $`T_{min}TT_H`$, where $`T_{min}=\mathrm{log}2`$ is the shortest period of the system and $`T_H=2\pi \overline{\rho }=\mathrm{log}(E_F/2\pi )`$ is the Heisenberg time at $`E_F`$, and heuristic for longer ones, show  that the form factor of the complex zeros of $`\zeta (s)`$ tends, as $`E_F`$ goes to infinity, to the corresponding function $`K_{GUE}`$ of Gaussian random matrices with unitary symmetry (GUE). The latter behaves as $`K_{GUE}(T)=T`$ for $`T<T_H`$, and $`K_{GUE}(T)=T_H`$ for $`T>T_H`$. However, the replacement $`K=K_{GUE}`$ in Eq.(4) and its extrapolation to short times is of no use to understand the behavior of $`\stackrel{~}{\mathrm{\Omega }}`$, because the integral diverges. In other words, the second moment is dominated by the (non-universal) short periodic orbits whose contribution has to be considered explicitly. We therefore keep as leading-order approximation the diagonal part $`(i=j)`$ in the double sum (4). This gives the convergent sum

$$\stackrel{~}{\mathrm{\Omega }}_0^2=\frac{1}{2\pi ^2}\underset{p}{}\underset{r=1}{\overset{\mathrm{}}{}}\frac{1}{r^4p^r\mathrm{log}^2p}7.9\times 10^2.$$
(6)
This equation exhibits two of the main properties of $`P(\stackrel{~}{\mathrm{\Omega }})`$. The first one is its asymptotic independence with respect to energy. This fact is not generic for dynamical systems, since for example for a Fermi gas in a chaotic cavity $`\stackrel{~}{\mathrm{\Omega }}^2`$ grows linearly with the Fermi energy . It is due to the independence of the periods and Lyapounov exponents with respect to energy in the Riemann dynamics. The second property is its non-universality (the sum (6) depends on the particular properties of the short orbits). The longer, statistically universal, periodic orbits provide next-to-leading order corrections to Eq.(6). For times $`TT_{min}`$ the form factor $`K_{GUE}`$ can be used in Eq.(4), leading to the correction (independent of the prime numbers)

$$\stackrel{~}{\mathrm{\Omega }}^2=\stackrel{~}{\mathrm{\Omega }}_0^2\frac{1}{12\pi ^2\mathrm{log}^2(E_F/2\pi )}.$$
(7)
We have numerically checked the validity of Eq.(7) down to values of $`E_F`$ of the order of a few thousands.

The third basic feature of the distribution is its asymmetry, as revealed by the third moment $`\stackrel{~}{\mathrm{\Omega }}^3`$, whose computation is now sketched. The third power of $`\stackrel{~}{\mathrm{\Omega }}`$ computed from Eq.(3) involves products of three cosines containing as argument the action of three different prime numbers. This product may be expressed as a sum of cosines involving the sum and differences of the actions. The term involving the sum of actions has zero average. As before, the remaining sum is dominated by the smallest primes. Due to the non-commensurability of the periods of the different primes we restrict moreover to the approximation $`p_i=p_j=p_k=p`$. The oscillating factors in the sum now have the typical form $`\mathrm{cos}[E(r_i+r_jr_k)\mathrm{log}p]`$. The only terms of this type having a non-zero average are those satisfying $`r_k=r_i+r_j`$. Since there are three different possibilities for choosing the backward repetition $`r_k`$, we finally obtain the convergent sum

$`\stackrel{~}{\mathrm{\Omega }}_0^3`$ $`={\displaystyle \frac{3}{4\pi ^3}}{\displaystyle \underset{p}{}}{\displaystyle \underset{r_i,r_j=1}{\overset{\mathrm{}}{}}}[r_i^2r_j^2(r_i+r_j)^2p^{r_i+r_j}\mathrm{log}^3p]^1`$ (9)
$`5.78\times 10^3.`$
The asymmetry of $`P(\stackrel{~}{\mathrm{\Omega }})`$ is therefore related to a simple but non-trivial interference phenomenon occurring between two (different or not) repetitions of a given orbit compensated by a backward one. The extreme case where only $`p=2`$ and its repetitions are taken into account in Eq.(9) yields already $`90\%`$ of its value.

The higher moments of the distribution are computed likewise. In the limit $`E_F\mathrm{}`$ to leading order they are given by expressions similar to, though more complicated than, Eq.(9). Table I shows a comparison of the analytical results with numerical data up to the sixth moment. The agreement between both calculations is astonishing. Our numerical results also coincide for the lowest moments with those obtained in .

An additional characterization of the Riemannium comes from the autocorrelation of the total energy $`𝒞_\mathrm{\Omega }(ϵ)=\stackrel{~}{\mathrm{\Omega }}(E_F)\stackrel{~}{\mathrm{\Omega }}(E_F+ϵ)/\stackrel{~}{\mathrm{\Omega }}_0^2`$. We evaluate this quantity to leading order by using the expansion (3) and making a diagonal approximation. This leads to

$$𝒞_\mathrm{\Omega }(ϵ)=\frac{1}{2\pi ^2\stackrel{~}{\mathrm{\Omega }}_0^2}\underset{p}{}\underset{r=1}{\overset{\mathrm{}}{}}\frac{\mathrm{cos}(ϵr\mathrm{log}p)}{r^4p^r\mathrm{log}^2p}.$$
(10)
Due to the dominance of the shortest orbits in the sum (10) $`𝒞_\mathrm{\Omega }`$ is an irregularly fluctuating non-decaying function, as illustrated on Fig. 2. The fundamental period $`2\pi /\mathrm{log}29.06`$ associated to the shortest orbit can be clearly observed.

It is interesting to make a comparison of the previous results with random matrix theory. As already emphasized, the $`GUE`$ is not really appropriate because no short-time scale equivalent to $`T_{min}`$ is built-in in the theory (and this produces a divergence of the moments). We therefore concentrate on the $`CUE`$ which, contrary to the Gaussian one, has an inherent short-time scale (although more appropriate measures, not relevant for the present discussion, have been recently proposed ). The following analysis reveals also some striking similarities existing between the Riemann zeros and eigenvalues of circular ensembles. Consider an $`N\times N`$ unitary matrix $`U`$ describing the time-periodic dynamical evolution of a quantum system. We fix for simplicity the periodicity to one. The Floquet spectrum of $`U`$ is given by the eigenvalue equation $`U\psi _\alpha =\mathrm{exp}(i\theta _\alpha )\psi _\alpha `$, $`\alpha =1,\mathrm{},N`$. $`\psi _\alpha `$ are the (stroboscopic) eigenstates and $`\theta _\alpha `$ the eigenphases. The spectral density on the unit circle is $`\rho _{CUE}(\theta )=_\alpha \delta (\theta \theta _\alpha )`$. The $`2\pi `$–periodicity of this function leads to a decomposition in smooth plus oscillatory parts, $`\rho _{CUE}=N/2\pi +\stackrel{~}{\rho }_{CUE}`$, with

$$\stackrel{~}{\rho }_{CUE}(\theta )=\frac{1}{\pi }\mathrm{Re}\underset{k=1}{\overset{\mathrm{}}{}}\mathrm{Tr}U^k\mathrm{exp}(ik\theta ).$$
(11)
As was done for Eq.(2), a direct comparison of Eq.(11) with Gutzwiller trace formula allows to make a “naïve” semiclassical interpretation of this equation. We look at it as a sum over the periodic orbits of a classical map labelled by the index $`k=1,2,\mathrm{}`$, having period $`T_k=k`$, action $`S_k=k\theta `$, (local) energy $`\theta `$, stability amplitude $`A_k=\mathrm{Tr}U^k`$, and $`\mathrm{}=1`$ (repetitions are degenerated with the fundamental periods). The interpretation of $`\theta `$ as the energy satisfies, as it should, the classical relation $`T_k=S_k/\theta `$. We recover in this analogy one of the basic properties of the Riemann dynamics: the independence of the periods with respect to energy. Short orbits correspond to $`T_kT_{min}=1`$, while long ones contributing to the density at the scale of the mean level spacing have period $`T_kT_H=2\pi \overline{\rho }=N`$. If now an ensemble of unitary matrices is considered, Eq.(11) acquires a statistical meaning since the prefactors $`\mathrm{Tr}U^k`$ have now a distribution. It is well known through, e.g., the prime number theorem and the Hardy-Littlewood conjecture, that asymptotically the statistical properties of large prime numbers are those required to mimic random matrix statistics . So for long times Eqs.(2) and (11) are very similar from a statistical point of view. This is not the case for times $`TT_{min}`$. To have a quantitative description we look at the $`CUE`$ distribution of the grand potential, obtained by integrating $`\stackrel{~}{\rho }_{CUE}`$ twice with respect to $`\theta `$,

$$\stackrel{~}{\mathrm{\Omega }}_{CUE}=(1/\pi )\mathrm{Re}\underset{k}{}(\mathrm{Tr}U^k/k^2)exp(ik\theta ).$$
(12)
The second moment is

$$\stackrel{~}{\mathrm{\Omega }}_{CUE}^2=(1/2\pi ^2)\underset{k}{}|\mathrm{Tr}U^k|^2/k^4.$$
(13)
The average is done over the ensemble of matrices. This expression is the $`CUE`$ analog of Eq.(4) ($`K_{CUE}=|\mathrm{Tr}U^k|^2`$). Using that $`|\mathrm{Tr}U^k|^2=k`$ if $`kN`$ and $`|\mathrm{Tr}U^k|^2=N`$ if $`k>N`$ , we get

$$\stackrel{~}{\mathrm{\Omega }}_{CUE}^2=\frac{\zeta (3)}{2\pi ^2}\frac{1}{12\pi ^2N^2}+𝒪(1/N^4),$$
(14)
to be compared to Eq.(7). As for the Riemannium, the leading term $`\zeta (3)/2\pi ^26.1\times 10^2`$ is controlled by short times (small values of $`k`$) in the sum (13). Because the short-time structure differ in both cases, the constants do not coincide. But the independence of the leading terms with respect to $`N`$ (or energy) is not a coincidence and comes from the structural similarities between Eqs.(2) and (11). Moreover, the next-to-leading order corrections are the same if one identifies  the corresponding Heisenberg times $`N=\mathrm{log}(E_F/2\pi )`$. Both corrections come from times of order $`T_H`$, where the statistical properties agree. This situation is in contrast with the results obtained, for instance, for the distribution of $`\mathrm{log}\zeta (1/2+iE)`$, which has been shown to agree to leading order with the (finite $`N`$) $`CUE`$ distribution, followed by non-universal corrections .

The qualitative agreement between $`\stackrel{~}{\mathrm{\Omega }}^2`$ and $`\stackrel{~}{\mathrm{\Omega }}_{CUE}^2`$ is lost when considering odd powers of $`\stackrel{~}{\mathrm{\Omega }}_{CUE}`$. The latter can be shown to vanish in the limit $`N\mathrm{}`$. More generally, and contrary to the Riemann case, in that limit the distribution of $`\stackrel{~}{\mathrm{\Omega }}_{CUE}`$ tends to a Gaussian. This can be seen from Eq.(12). On the one hand we have demonstrated that $`\stackrel{~}{\mathrm{\Omega }}_{CUE}`$ is asymptotically dominated by the lowest contributions in the sum. But $`\mathrm{Tr}U^k`$ is known to be, for $`k`$ finite and $`N\mathrm{}`$, Gaussian independent distributed . Therefore $`\stackrel{~}{\mathrm{\Omega }}`$ has also a limiting normal distribution.

The techniques employed here to analyze the distribution of the grand potential are not specific to the Riemann dynamics and can be applied in a wider context to the theory of Fermi gases. Our results indicate that in general the grand potential of a Fermi gas, as well as other physical quantities derived from it, are controlled by the shortest non-universal classical orbits. This dominance leads, via the interference mechanisms illustrated here for the grand potential of the Riemannium, to non-Gaussian asymmetric distributions for the associated quantity. Some related results illustrating this general phenomenon were obtained by Tsang , who demonstrated non-Gaussian asymmetric distributions for the error term in the mean square formula of $`|\zeta (1/2+iE)|`$, and for the fluctuations of the number of lattice points inside a curve, namely a circle (see also ) and an hyperbola (known as the circle problem and Dirichlet’s divisor problem, respectively). The last two problems have a dynamical interpretation in terms of the fluctuations of the spectral counting function of an associated integrable system. The mechanism leading to the asymmetry in  is very similar to the one found here. In contrast, the fluctuations of the counting function in the Riemann case are known to be asymptotically Gaussian  and are not dominated by the short orbits. This difference in the behavior of the counting function can be traced back to the different short-time behavior of the form factor for integrable and chaotic systems, that produces a dominance of the short orbits in the former. We find however that integrals of the counting function (the grand potential being the first one, cf Eq.(1)) are always dominated by the short orbits, irrespective of the (chaotic or integrable) nature of the underlying classical dynamics. On the mathematical side, the present results can be extended in a straightforward manner to general Dirichlet’s $`L`$-functions.

In summary we have considered, guided by physical analogies, new properties of the Riemann zeros that reveal new aspects of their dynamical interpretation as quantum eigenvalues of a classically chaotic system. Our results are relevant in the theory of Fermi systems, as well as in the general context of quantum chaotic motion. Concerning the Riemann hypothesis and the search of the Hilbert–Pólya Hamiltonian, the present study suggests that time–periodic dynamical evolutions have to be considered as serious candidates.

We are specially grateful to D. Hejhal, J. Keating and A. Odlyzko for fruitful discussions and suggestions. The Laboratoire de Physique Théorique et Modèles Statistiques is a Unité Mixte de Recherche de l’Université Paris XI and CNRS.


text:
# Polarization dynamics of femtosecond pulses propagating in air

## I Introduction

There has been a great interest in long distance femtosecond pulse propagation in air in recent years motivated in part by potential applications in remote sensing and laser-induced lightning. The first experimental observations of highly localized, high-intensity filaments propagating over distances that exceed their corresponding Rayleigh lengths by orders of magnitude motivated the efforts to understand the phenomenon . Due to the violent formation process, and the extreme time and spatial scales of the filaments, the computer simulations and analytic approaches turned out to be vital tools to grasp the underlying physics . Several models have been proposed. The first explanation suggested a stationary waveguide formed by competing effects of nonlinear self-focusing and linear defocusing by underdense plasma generated in the most intense part of the pulse. An alternative model employed a notion of the “moving focus” to explain how the geometrical focus of a beam is transformed into a long filament. More recently, a dynamic spatial replenishment model emerged from numerical investigation by Mlejnek et al.. The essential feature of the femtosecond propagation of infrared (IR) pulses is its dynamics that makes it possible that the localized filaments propagate over long distances and do not suffer significant energy loses. The basic mechanism involves a dynamic balance between the nonlinear self-focusing and defocusing by free electrons generated via multi-photon absorption by the high-intensity filaments. When the self-focusing leading edge of the pulse starts to generate plasma, it thus creates a defocusing “lens” for the trailing portion of the pulse. This has two effects; First, it limits the losses due to absorption in the plasma and, second, it prevents the major portion of the pulse from experiencing self-focusing collapse. After the plasma generating leading portion of the pulse exhausts its energy, the strength of the defocusing lens decreases, and the self-focusing starts over again, this time in the “next temporal slice” of the pulse. The whole process can repeat several times, depending on the total energy of the pulse. The qualitative features of the dynamical spatial replenishment remain valid also in transversely wide pulses, that break-up into multiple filaments .

Up to the present, most of the work on the experimental side and all numerical studies have been concerned with the case of linearly polarized input pulses. Recently, Petit et al.  studied the effects of the polarization on the propagation of femtosecond IR pulses. They have measured luminescence from the plasma generated in the filaments to show that the polarization of the pulse plays an important role in the plasma generation. Due to the highly dynamic nature of the filament formation and propagation, it is natural to expect a rich polarization dynamics in femtosecond pulses. In this initial study, we restrict ourselves to femtosecond pulses with modest peak powers that retain their initial radial transverse symmetry and are just sufficient to produce several refocusing events within a single pulse. The initial polarization of the pulse is varied and the polarization state is recorded along the propagation path. Our results indicate a tight correlation between the evolution of the pulse waveform along the propagation distance and certain global polarization parameters. Thus, the measurements of the polarization state could provide yet another tool to extract information on dynamics of pulses as well as an opportunity to correlate experiment and theory.

The remainder of the paper is organized as follows. Our basic model is described in Sec. II, and details of our simulated experiment are given in Sec. III. Plasma production, polarization dynamics, and the associated supercontinuum generation are then discussed in Secs. IV-VI. Finally, our summary and conclusions are given in Sec. VII.

## II Model Equations

Since the typical transverse dimension of the self-focusing collapsing filaments are a few tens to hundred microns in the IR wavelength region, it is a good approximation to consider the optical field as transverse. We describe the complex optical amplitude in terms of two circularly polarized components, $`^\pm `$. The choice of the circularly polarized basis is motivated by the fact that in this basis the nonlinear birefringence is “diagonal” and preserves the power in each component, which in turn makes its implementation easier. Our model is a straightforward extension of the scalar version we used in our previous studies . It takes into account the effects of diffraction, normal group velocity dispersion (GVD), multi-photon ionization (MPI) and avalanche plasma generation, defocusing of light by the generated free electrons, and instantaneous and delayed cubic nonlinearity:

$`{\displaystyle \frac{^\pm }{z}}=`$ $`{\displaystyle \frac{i}{2k}}_{}^2^\pm {\displaystyle \frac{ik^{\prime \prime }}{2}}{\displaystyle \frac{^2^\pm }{t^2}}`$ (1)
$``$ $`{\displaystyle \frac{\sigma }{2}}(1+i\omega \tau )\rho ^\pm {\displaystyle \frac{\beta ^{(K)}}{2}}||^{2K2}^\pm `$ (2)
$`+`$ $`i{\displaystyle \frac{2\omega }{3c}}(1f)n_2(|^\pm |^2+2|^{}|^2)^\pm `$ (3)
$`+i{\displaystyle \frac{2\omega }{3c}}fn_2`$ $`\left[{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑t^{}R(tt^{})(|(t^{})^\pm |^2+2|(t^{})^{}|^2)\right]^\pm .`$ (4)
Here $`\omega `$ is the optical frequency, $`||^2=|^+|^2+|^{}|^2`$ the combined intensity of the two circularly polarized components, $`k=\omega /c`$, $`k^{\prime \prime }=^2k/\omega ^2`$, $`\rho `$ is the electron density, $`\sigma `$ is the cross-section for inverse bremsstrahlung, $`\tau `$ is the electron collision time, $`\beta ^{(K)}`$ is the $`K`$-photon absorption coefficient, and the nonlinear change in refractive-index for a linearly polarized continuous wave (cw) field is $`n_2||^2`$. The corresponding critical power for self-focusing collapse for linearly polarized fields is $`P_{cr}=\lambda _0^2/2\pi n_2`$, or $`P_{cr}=`$1.7 GW for our parameters. In contrast, the critical power for self-focusing for circularly polarized fields $`P_{cr}^\pm =1.5P_{cr}`$ is 1.5 times that for linearly polarized fields . The normalized response function (characterized by the resonance frequency $`\mathrm{\Omega }`$ and the decay $`\mathrm{\Gamma }`$)

$`R(t)`$ $`=`$ $`\theta (t)\mathrm{\Omega }^2e^{\mathrm{\Gamma }t/2}{\displaystyle \frac{\mathrm{sin}(\mathrm{\Lambda }t)}{\mathrm{\Lambda }}},\mathrm{\Lambda }=\sqrt{\mathrm{\Omega }^2\mathrm{\Gamma }^2/4},`$ (5)
accounts for delayed nonlinear effects, and $`f`$ is the fraction of the cw nonlinear optical response which has its origin in the delayed component, and we denoted the Heaviside step function by $`\theta (t)`$. In the present model, we chose the relative weight of the “self” and “cross” non-linear birefringent terms the same as for the instantaneous Kerr effect and nonlinear Raman effect, namely, the cross effect has a weight twice that of the self effect as is appropriate for an isotropic medium .

The optical field amplitude equations are completed by a simple equation describing the evolution of the plasma density:

$`{\displaystyle \frac{\rho }{t}}`$ $`=`$ $`{\displaystyle \frac{\sigma }{E_g}}\rho ||^2+{\displaystyle \frac{\beta ^{(K)}||^{2K}}{K\mathrm{}\omega }}a\rho ^2,`$ (6)
where the first and second terms represent the avalanche and multi-photon effects, while the last one models the plasma recombination. Note, that at the time-scales relevant for the present study, the only practically important contribution comes from the the multi-photon term. For the pulse powers we use, it is sufficient to include oxygen alone as a source for MPI, since its corresponding multiphoton order is lower than that for nitrogen. We use the Keldysh theory formula to calculate the MPI rate .

Explanations of symbols that appear in our model equations are listed in Table 1 together with the values used in our simulation.

## III Numerical simulation

To study the role of the initial polarization state of the pulse on its subsequent propagation, we performed a series of simulations. In all runs, the initial pulse was chosen to be a Gaussian wave-packet (in time and space) characterized through its central wavelength $`\lambda =775`$ nm, pulse waist $`w=0.7`$ mm and temporal duration $`\tau _{\mathrm{FWHM}}=200`$ fs. We fixed the initial peak intensity as $`1.0\times 10^{16}`$ W m<sup>-2</sup> which is a relatively modest value: for a linearly polarized input this corresponds to an input peak power of $`P=7.7`$ GW$`=4.5P_{cr}`$, whereas for a pure circular polarization $`P=7.7`$ GW$`=3.0P_{cr}^\pm `$. At these powers, there are typically two to three refocusing events in the pulse propagation, and the numerics can be reliably controlled. At higher powers under a perfect axial symmetry of the pulse, it may be necessary to extend the model beyond the nonlinear Schrödinger equation (NLSE), and to include correction terms that allow to handle pulses with a very broad spectrum. To check our numerical procedures and model implementation, we performed most of the runs in two or three different resolutions. The data we present were obtained with the time-domain resolution of 0.24 fs. This resolution is sufficient to capture spectra several hundredth of nanometers wide. We provide more details of our testing procedures concerning the spectral resolution and algorithm validity in Sec. VI devoted to supercontinuum generation.

Below, we present our results for a series of runs that differ in the initial polarization state of the pulse. We change the polarization from the linear, through elliptic to circular, to see how it affects the dynamics of the filaments. In previous work on nonlinear propagation in fibers, the Stokes parameter formalism has been employed to classify the polarization dynamics for plane wave fields . Here we employ space and time averaged Stokes parameters $`(s_0,s_1,s_2,s_3)`$ as a description of the polarization dynamics of the propagating pulses. In particular, our Stokes parameters are calculated numerically according to the prescription

$`s_0`$ $`=`$ $`(_{\pi /2}+_0)`$ (7)
$`s_1`$ $`=`$ $`(_{\pi /4}_{\pi /4})/s_0`$ (8)
$`s_2`$ $`=`$ $`(_{\pi /2}_0)/s_0`$ (9)
$`s_3`$ $`=`$ $`(_{\mathrm{CR}}_{\mathrm{CL}})/s_0`$ (10)
where $`_\beta `$ is the total energy detected after passing the pulse through a polarizer of state $`\beta `$

$$_\beta =2\pi _{\mathrm{}}^+\mathrm{}_0^R|_\beta |^2(r,t)r𝑑r𝑑t$$
(11)
with $`_\beta `$ the field resolved along the polarizer direction, and $`R`$ the radius of a detection aperture chosen to select the the most intense region of the pulse around the filament. For a linear polarizer $`\beta `$ is the angle of the polarizer, whereas for a circular polarizer $`\beta =CR,CL`$ corresponding to right and left circular polarization settings. Here we chose $`R=0.1`$ mm, and perform the measurement of the Stokes parameters “in the near field”. We note that $`s_0`$ is the total energy of the pulse detected over the aperture, and the remaining Stokes parameters are calculated as differences between the detected energy with different polarizer settings normalized to the total energy. Thus our prescription has a direct experimental interpretation and should therefore be of utility. The above prescription for determining the Stokes parameters also reduces to the usual definitions in the limit of long pulses of broad transverse extent.

In the numerical simulations to be presented we fix the initial Stokes parameter $`s_2`$ equal to zero, and vary $`s_1`$ and $`s_3`$ between zero and one to vary the initial polarization from linear through elliptic to pure circular polarization. Besides the polarization state, we also recorded the data pertaining to the plasma generation inside the pulse, and generation of the supercontinuum light. We start our discussion with plasma generation.

## IV Plasma production

As the femtosecond pulse undergoes multiple self-focusings, the amount of plasma generated by its high-intensity portions reflects the spatio-temporal shape of the pulse. The total number of generated electrons as well as the maximal plasma densities exhibit peaks along the propagation distance. These peaks coincide with the locations of self-focusing collapses within the pulse, each peak being produced by a different temporal portion of the pulse. Figure 1 shows the plasma generation for three different initial polarizations of the pulse. The trend that one can see is quite in line with what is expected based on the functional form of the nonlinear birefringence. Namely, as we change the initial polarization from linear through elliptic to circular, the amount of the generated plasma decreases. Also, the onset of filament formation is delayed for the circularly polarized pulse because the critical power for self-focusing is higher for circularly polarized pulse as noted earlier. In other words, keeping the input peak power the same for different polarizations, we effectively decrease the self-focusing power of circularly polarized pulses. This is also the reason why the number of refocusing events can be higher in a close-to-linear or linear polarization than in a circularly polarized pulse. While the overall plasma production depends on the polarization state, the typical dimensions of the filaments are not very sensitive to it. That can be seen from the Figure 1 which shows the longitudinal extent of the plasma columns. The transverse dimensions of the plasma channels can be estimated from the ratio of the two curves shown in the figure as the square root of the ratio between the total number of electrons and the maximal plasma density. This characteristic dimension of the plasma channel is shown in Fig. 2 for three different polarizations. Though there are small variation between different initial polarizations, the thickness of the plasma channels is always roughly 60 microns in the most dense parts. The plasma channel generated by the circularly polarized pulse seems to be more “homogeneous”, exhibiting less thickness variation along the propagation distance.

## V Polarization dynamics

While the dynamics of the plasma generation and its dependence on the polarization described in the previous Section is straightforwardly linked to the structure of the equations governing the optical field evolution, the polarization dynamics seems to be more difficult to interpret. Figures 3,4,5 show the Stokes parameters and the polarization degree as functions of the propagation distance for the three different initial polarizations we discussed in the previous Section. An interesting feature is the difference between the “stability” of initial linear and circular polarization. Figure 3 shows that a small perturbation to the linear polarization in the initial pulse leads to an increasing deviation of the polarization state from the initial one. In this sense, the linear polarization appears to be unstable, as the polarization measured after the filament formation can significantly differ from the initial one. Naturally, the rate of divergence for two close but not identical initial conditions decreases with the decreasing input power. On the other hand, in the case of almost circular polarization shown in Fig. 4, the pulse polarization state doesn’t change that dramatically. Though there is a small decrease of the polarization degree, one can say that final polarization stays close to the initial one even after two refocusings of the pulse. Thus, the circular polarization seems to be more stable against small perturbations than the linear polarization. Figure 5 shows an interesting case of an initially elliptic polarization. Note that the Stokes parameter $`s_3`$, which measures the degree of circular polarization, only exhibits small variations, while the other two parameters decrease significantly after their initial increase in the first collapse. That means that the light focused in the second collapse is predominantly circularly polarized. This observation is confirmed by examining the spatio-temporal polarization pattern within the pulse. This is an observation that may not be expected based on the previous results concerning cw self-focusing of polarized pulses . Namely, in a situation close to a continuous wave regime, one can argue that the weaker circular component experiences a stronger focusing “lens” because of the factor two in the birefringence cross-term, and that eventually leads to equal intensities of both circular components and, therefore, linear polarization of the central filament. However, the important point here is that the femtosecond light filaments under consideration are extremely dynamic objects. The resulting polarization distributions strongly depend on the spatial and temporal location within the pulse, and any interpretation based on steady-state-like considerations becomes invalid. Namely, there is a delay between when the light encounters the focusing “lens” and when it actually reaches the focus. This delay interferes with the temporal profile of the pulse, which typically exhibits multiple peaks that may be just a few femtosecond long. As a consequence, the above simple argument is not sufficient to capture all essential features of the phenomenon. The tendency of the predominantly circular polarization of the most intense portion of a filament was also observed in our simulations that were not restricted to axial symmetry . A wide beam with a random perturbation breaks up into multiple filaments that exhibit polarization properties similar to those we discuss here. However, one has to keep in mind, that in both cases, axisymmetric as well as fully spatially resolved, our simulation modeled pulses with relatively small energy fluence when compared to some current experiments. Therefore, it would be extremely interesting to see what happens to the polarization of the central filament in a pulse that has enough energy for many self-focusing events and also has the transverse profile clean enough to preserve its axial symmetry.

We conclude this Section with yet another presentation of the polarization dynamics data we have shown above. Namely, we want to show that the polarization changes closely reflect the self-focusing events within the pulse and, consequently, the locations where most of the plasma is generated. Figure 6 shows the root-mean-square rate of the change of the Stokes vector along the propagation distance

$$\frac{ds}{dz}=\sqrt{\left(\frac{ds_1}{dz}\right)^2+\left(\frac{ds_2}{dz}\right)^2+\left(\frac{ds_3}{dz}\right)^2}$$
(12)
for the case of elliptic initial polarization of the pulse. The curve shown corresponds to the data depicted in Fig. 5 and in Fig. 1b). Note, that the maxima of the rate of the polarization change closely follow those in the plasma production curve. We thus see that the multiple self-focusing events in the single pulse leave their signature on the polarization. This could provide another way, besides the indirect plasma density observations, to visualize the dynamics of the spatial replenishment.

## VI Supercontinuum generation

After contrasting the behavior of pulses polarized close to linear and circular from the points of view of plasma generation and of their polarization dynamics, we want to discuss the effects of polarization on the supercontinuum generation. However, before presenting our results, we feel a note concerning some technical questions is in order. The explosive spectral broadening in the supercontinuum generation in femtosecond pulses is a rather subtle phenomenon from the point of view of numerical simulation. Clearly, one needs a sufficient resolution in the time (spectral) domain to capture the broad spectrum, but the resolution may not be the only issue here. One has to check how broad is the spectral region within which the model and its numerical implementation describes the wave propagation correctly. It is expected, that at extreme powers, correction terms beyond the basic NLSE (see e.g. Ref. ) need to be included in the field equation. To ensure that we work in the regime where the correction terms may be neglected, and to assess the spectral band over which our numerics works well, we performed some comparative simulations. The choice of the reference frequency (wavelength) around which the NLSE is built is in principle arbitrary, though it is obviously most appropriate to choose it close to the central frequency of the modeled pulse. This means, that simulations that only differ in the choice of the reference frequency should give the same results. We have compared simulation with the reference frequency shift of 150 nanometers off the central wavelength of the pulse, and obtained a very close match of the spectra in the region from 500 to 1200 nm and over four decades in spectral intensity. Thus, in this interval we can trust the spectra extracted from our simulations. We would like to point out that this is in fact a rather strong test for the overall numerical implementation of the solver. But most important, the close agreement between the spectra shows that the correction terms beyond NLSE do not play an important role in our regime of modest input powers; This is because one can interpret them as corrections that partially restore the invariance of the original wave equation with respect to the choice of the (physically meaningless) reference frequency – that is exactly what our test shows is not needed for the conditions in our simulations.

Figure 7 shows comparison of the spectral broadening of two pulses that differ only in their initial polarization. The way we extracted the spectra from the pulse waveforms corresponds to measurement in a near field with the same aperture we used for polarization characterization. The figure shows spectra “measured” after the last self-focusing collapse, after which the pulse will eventually diffract and there will not be more supercontinuum generation. One can see that the pulse which was initially polarized close to linear exhibits a much stronger supercontinuum generation. However, note that we compare pulses with the same peak power, and what we see here is an effect similar to the plasma generation. Supercontinuum generation strongly depends on the available power and the natural measure of that power is in units of critical power for self-focusing. From that point of view the circularly polarized pulse is weaker and that is the main reason that it exhibits less spectral broadening. However, this situation represents a reasonable experimental setup in which only the polarization is changed.

To get a feeling regarding the role of the group velocity dispersion in the supercontinuum generation, we performed most of the simulation runs also with a higher group velocity dispersion parameter. It turns out that increasing GVD by an order of magnitude leads to a strong suppression of the continuum production. We speculate that it may be one of the reasons that, at least in some experiments, there is only little spectral broadening in the ultraviolet (UV) femtosecond pulses because the GVD value of air is significantly higher in the UV region.

The findings from our numerical simulations should be accessible to experimental testing. However, extreme caution should be exercised when trying to compare experimental and simulational spectra. The spectra we present are taken from “one shot”. They exhibit modulations typical for supercontinuum generation in gases . While the characteristic “frequency” of the modulation is rather reproducible, the exact spectral shape is not. Thus, even small fluctuation in the parameters of the pulses will result in a suppression of fine features in the multiple-shot experimental spectra.

## VII Conclusions

We have performed a computer simulation study of the effects of the initial pulse polarization on its propagation and filamentation dynamics. In agreement with the experiment , we have found that the filamentation onset is reached earlier for a linearly polarized pulse than in a circularly polarized pulse of the same peak power. However, in some cases, the experiment indicates that the circularly polarized pulses create higher plasma densities in comparison with linearly polarized pulses. Our simulations suggest the opposite, but one has to keep in mind that the experimental measurements and our simulation pertain to rather different conditions, including much higher power and focusing in the experiment. In our simulations, we also see more self-focusing events with linearly polarized pulses than with circularly polarized ones.

We have observed a rich spatio-temporal polarization dynamics. Naturally, the limited range of the full parameter space explored in our simulations prevents us from drawing general conclusions, but we believe some tendencies are already discernible. First, the initially circular polarization seems to be stable in a sense that a small perturbation of the polarization state doesn’t grow significantly. On the other hand, when a small polarization perturbation is applied to the linearly polarized pulse, it grows and the polarization degree of the central filament decreases significantly. The growth rate of the deviation is expected to increase for higher powers. An interesting case is the one of an initially elliptic polarization. We have observed that the center of the filament is almost purely circularly polarized after subsequent self-focusing collapses within the pulse. Apparently, we have here a rather different situation than in the self-focusing in nanosecond pulses, which tend to create linearly polarized filaments independently of initial polarization state. The usual argument, that the the weaker circular component experiences stronger self-focusing which finally leads to balancing the power of the circularly polarized components and thus to linear polarization, can’t be applied in the femtosecond pulses. The evolution of the polarization state along the propagation distance is extremely dynamic and hardly possible to describe in simple “static” terms. The multiple refocusing within a pulse and the defocusing effect of the generated plasma play a major role. The important feature of the whole process is that apart from the relatively small energy losses due to plasma generation and radiation, the energy in each circular component remains conserved. Thus, the main mechanism that results in changing polarization pattern within the pulse is the spatio-temporal energy redistribution within each circular component. As a consequence, the polarization state of the whole pulse is very complicated and, therefore, any projection onto global quantities like Stokes parameters has to be interpreted in relation to details of the measurement (aperture, near vs. far field, collecting angle, …).

We have also looked at spectra “measured” after the last self-focusing collapse for different initial polarizations. In accordance with our observation about the plasma generation, we see much stronger supercontinuum generation in linearly polarized than in circularly polarized pulses.

Finally, we have seen that rate of change of the polarization state in the center of the filament is closely correlated with self-focusing and plasma generation. Thus, the polarization offers, in principle, an alternative way to investigate the dynamics of the spatial replenishment in femtosecond pulses.

In this work, we have concentrated on investigating “global” quantities to characterize the femtosecond pulse propagation that should be experimentally accessible, at least in principle. Naturally, the question is how much can be done practically. To measure the Stokes parameters evolution along the propagation distance, for example, a very good reproducibility of the initial pulse would be required. However, the “final output” polarization state and spectra, which should be easier to measure, also carry a lot of signatures about the inner dynamics of the femtosecond pulse propagation.

## VIII Acknowledgments

Work supported by AFOSR grant no. F4962-00-1-0312, AFOSR DURIP grant no. F4962-00-1-0190, and in part by Boeing. M.K. was partly supported by the GASR grant VEGA 2/7174/20.