text:
# Generalizing the Debye-Hückel equation in terms of density functional integral

## Abstract

We discuss the validity of generalized Debye-Hückel (GDH) equation proposed by Fisher et al. from the functional integral point of view. The GDH theory considers fluctuations around prescribed densities of positive and negative charges. Hence we first formulate a density functional integral expression for the canonical system of Coulomb gas, and also demonstrate that this is a dual form to the Sine-Gordon theory. Our formalism reveals the following: (i) The induced charge distribution around supposed density favors not only the cancellation of additional electrostatic potential like the original DH theory, but also the countervailing of chemical potential difference between imposed and equilibrium value. (ii) As a consequence apparent charge, absent in the GDH equation, comes out in our generalized equation. (iii) That is, the GDH equation holds only in special cases.

Much attention in recent years has been paid to novel phenomena of ionic fluids and charged colloids, such as unexpected behaviors of critical exponents in ionic fluids , (micro)phase separation of charged colloids , attraction between similarly charged objects , and so on. Of the numerous theoretical attempts for explaining these, this paper will be particularly concerned with a relevant tool for studying the ionic criticality: generalized Debye-Hückel (GDH) theory proposed by Fisher et al. . Then let us first see what is generalized.

One will find two extensions from the GDH equations (1) and (2) given below. Before doing so, however, the setting is to be described: Consider here the restricted primitive model (RPM), consisting of two oppositely charged, but otherwise identical, sets of $`NN_+=N_{}`$ hard spheres of diameter $`a`$ and charge per particle $`\pm e`$, immersed in a medium of dielectric constant $`ϵ`$. Also denote by $`\mathrm{\Psi }_+(𝐫;𝐫_1)`$ the mean electrostatic potential in the $`k_BT/e`$ unit (the others being the same) at a general point $`𝐫`$ when the positive charge 1 is fixed at $`𝐫_1`$, and by $`\mathrm{\Psi }_{imp}(𝐫)`$ the imposed electrostatic potential determined from prescribed densities of positive and negative charges, i.e., $`n_+(𝐫)`$ and $`n_{}(𝐫)`$ as $`^2\mathrm{\Psi }_{imp}(𝐫)=4\pi l_BQ(𝐫)`$ with putting that $`l_Be^2/ϵk_BT`$ and $`Q=n_+n_{}`$.

Fisher et al. advocate  that the ”local induced potential” defined by $`\stackrel{~}{\mathrm{\Psi }}_+(𝐫;𝐫_1)\mathrm{\Psi }_+(𝐫;𝐫_1)\mathrm{\Psi }_{imp}(𝐫)`$ be to satisfy the GDH equations as follows:

$`^2\stackrel{~}{\mathrm{\Psi }}_+(𝐫;𝐫_1)`$ $`=`$ $`4\pi l_B[\delta (𝐫𝐫_1)Q(𝐫)],|𝐫𝐫_1|a,`$ (1)
$`^2\stackrel{~}{\mathrm{\Psi }}_+(𝐫;𝐫_1)`$ $`=`$ $`\kappa ^2(𝐫)\stackrel{~}{\mathrm{\Psi }}_+(𝐫;𝐫_1),|𝐫𝐫_1|a,`$ (2)
implying modifications that are in order:

* A remarkable difference from the original DH equation  is seen in Eq. (2) where the Debye-Hückel screening length $`\kappa ^1`$ is generalized to be spatially dependent on the supposed densities as $`\kappa ^24\pi l_B[n_+(𝐫)+n_{}(𝐫)]`$.
* There also appears, unlike the primary version, the second term on the right hand side of (1) representing an effective ”cavity source” term.

The validity of these two generalizations are what we would like to investigate using the language of density functional integral.

For convenience of later discussion, let us further detail the GDH theory  exploiting the above equations (1) and (2). According to this theory, the Helmholtz free energy $`F`$ (in the $`k_BT`$ unit) as sum of the imposed and induced free energy, $`F=F_{imp}+F_{ind}`$, is given in the following density functional form:

$`F_{imp}`$ $`=`$ $`{\displaystyle 𝑑𝐫𝑑𝐫^{}\frac{l_B}{2}Q(𝐫)\frac{1}{|𝐫𝐫^{}|}Q(𝐫^{})}+{\displaystyle 𝑑𝐫n_+\mathrm{ln}n_+}+n_{}\mathrm{ln}n_{}n_+n_{},`$ (3)
$`F_{ind}`$ $`=`$ $`{\displaystyle \underset{\sigma =+,}{}}{\displaystyle 𝑑𝐫_\mathrm{𝟏}n_\sigma (𝐫_1)_0^1𝑑\lambda \mathrm{\Phi }_\sigma (𝐫_1;\lambda e)},`$ (4)
where $`\mathrm{\Phi }_\sigma (𝐫_1;\lambda e)`$ is the mean induced electrostatic potential at the site $`𝐫_1`$ of a fixed ion defined as $`\mathrm{\Phi }_\sigma (𝐫_1;e)lim_{𝐫𝐫_\mathrm{𝟏}}\left[\stackrel{~}{\mathrm{\Psi }}_\sigma (𝐫;𝐫_\mathrm{𝟏})\sigma l_B/|𝐫𝐫_\mathrm{𝟏}|\right](\sigma =+,)`$. The GDH theory especially imposes $`n_\sigma `$ $`(\sigma =+,)`$ on such a simple undulation as

$$n_\sigma =\overline{n}[1+\mathrm{\Lambda }\mathrm{cos}(𝐤𝐫+\theta _\sigma )].$$
(5)
Consequently the quadratic terms of $`\mathrm{\Lambda }`$ in the free energy difference, $`F(\{n_\sigma \})F(\{\overline{n}\})`$, gives the Fourier transform of various correlation functions by taking either phase, $`\mathrm{\Delta }\theta |\theta _+\theta _{}|=0\text{or}\pi `$, on a case-by-case basis. Thus this approach, in contrast to the original DH analysis, has succeeded in yielding both charge-charge oscillatory correlations at high densities (where $`\mathrm{\Delta }\theta =\pi `$ is taken) and density-density correlations (where $`\mathrm{\Delta }\theta =0`$) that exhibit a divergent correlation length at criticality.

Indeed, the validity of the GDH theory is strongly suggested not only by the above usefulness of this method, but also from satisfying both the Stillinger-Lovett sum rule  and the exact low-density limiting-law for various correlation lengths . In terms of functional integral, however, the GDH equations (1) and (2) are never trivial, due to the irrelevance of the Sine-Gordon (SG) mapping . The details follow: With use of the chemical potentials $`\mu _\sigma `$ $`(\sigma =+,)`$ and the charge density operator $`\widehat{q}(𝐫)`$ defined as $`\widehat{q}\widehat{\rho }_+\widehat{\rho }_{}`$ by the number density operators $`\widehat{\rho }_\pm _{i=1}^N\delta [𝐫𝐫_i^\pm ]`$, the grand partition function $`\mathrm{\Xi }`$ reads

$$\mathrm{\Xi }=\underset{N}{\overset{\mathrm{}}{}}\frac{e^{\mu _+N}}{N!}\underset{N}{\overset{\mathrm{}}{}}\frac{e^{\mu _{}N}}{N!}Z,$$
(6)
where

$$Zd\{𝐫^{2N}\}\mathrm{exp}\left[\frac{l_B}{2}𝑑𝐫𝑑𝐫^{}\frac{\widehat{q}(𝐫)\widehat{q}(𝐫^{})}{|𝐫𝐫^{}|}\right]$$
(7)
with putting that $`d\{𝐫^{2N}\}_{i=1}^Nd𝐫_i^+_{j=1}^Nd𝐫_i^{}`$. In the conventional SG mapping, eqs. (6) and (7) are Hubbard-Stratonovich transformed  to the functional integral of potential field, not of density. Therefore, within the SG frame, it is impossible to evaluate fluctuations around imposed number densities (not around prescribed potential) as the GDH theory does.

To generalize the original Debye-Hückel equation as Fisher et al. propose, it is hence indispensable to transform expressions (6) and (7) to the density functional integral formalism of the RPM. The former part in the remainder is then devoted to both its formulation and comparison with the SG theory. In the final, this formalism will reveal that the GDH equations (1) and (2) hold only in special cases.

Recently we have shown  that a dual method to the above Hubbard-Stratonovich way enables one to transform the configurational integral expressions (6) and (7) to the following density functional integral representation:

$`\mathrm{\Xi }={\displaystyle \underset{\sigma =+,}{}}{\displaystyle D\rho _\sigma \mathrm{exp}(H_\mathrm{\Xi }\{\rho _\pm \})},`$ (8)
$`H_\mathrm{\Xi }\{\rho _\pm \}{\displaystyle \frac{l_B}{2}}{\displaystyle 𝑑𝐫𝑑𝐫^{}\frac{q(𝐫)q(𝐫^{})}{|𝐫𝐫^{}|}}+{\displaystyle \underset{\sigma =+,}{}}{\displaystyle 𝑑𝐫\rho _\sigma \mathrm{ln}\rho _\sigma }\rho _\sigma \rho _\sigma \mu _\sigma ,`$ (9)
where $`D\rho _\sigma _{\{𝐫\}}d\rho _\sigma `$, and $`q`$ and $`\rho _\sigma `$ are c-numbers of the corresponding operators.

What is required in eqs. (3) and (4), however, is the Helmholtz free energy. We then move to the canonical partition function $`Z`$, via performing the following contour integral,

$$Z=\underset{\sigma =+,}{}\frac{1}{2\pi i}𝑑\lambda _\sigma \frac{\mathrm{\Xi }}{\lambda _\sigma ^{N+1}},$$
(10)
where $`\lambda _\sigma =e^{\mu _\sigma }`$ is now a complex variable. Substituting eqs. (8) and (9) into the integrand, we have

$`Z`$ $`=`$ $`{\displaystyle \underset{\sigma =+,}{}}{\displaystyle D\rho _\sigma \mathrm{exp}\left[H_\mathrm{\Xi }\{\rho _\pm \}\underset{\sigma =+,}{}𝑑𝐫\mu _\sigma \rho _\sigma \right]\frac{1}{2\pi i}𝑑\lambda _\sigma \lambda _\sigma ^{1[N{\scriptscriptstyle 𝑑x\rho _\sigma (x)}]}}.`$ (11)
Since one finds from the Cauchy’s integral theorem

$`{\displaystyle \frac{1}{2\pi i}}{\displaystyle 𝑑\lambda _\sigma \lambda _\sigma ^{1[N{\scriptscriptstyle 𝑑x\rho _\sigma (x)}]}}`$ (12)
$`=\{\begin{array}{cc}1\hfill & \text{if }𝑑x\rho _\sigma (x)=N\hfill \\ 0\hfill & \text{otherwise},\hfill \end{array}`$ (15)
eq. (11) is formally rewritten as

$`Z`$ $`=`$ $`{\displaystyle \underset{\sigma =+,}{}}{\displaystyle D\rho _\sigma \mathrm{exp}\left[H_Z\{\rho _\pm \}\right]\delta \left[𝑑𝐫\rho _\sigma (𝐫)N\right]},`$ (16)
$`H_Z`$ $``$ $`{\displaystyle \frac{l_B}{2}}{\displaystyle 𝑑𝐫𝑑𝐫^{}\frac{q(𝐫)q(𝐫^{})}{|𝐫𝐫^{}|}}+{\displaystyle \underset{\sigma =+,}{}}{\displaystyle 𝑑𝐫\rho _\sigma \mathrm{ln}\rho _\sigma }\rho _\sigma .`$ (17)
We have thus established the framework of which a priori use  has been made so far.

Before entering into the main question, let us also see that the density functional integral formalism gives back the SG theory of the canonical version. To this end we first introduce electrostatic potential field $`\psi (𝐫)`$ via using the equivalence between eqs. (16) and (17) and the following form:

$`Z={\displaystyle D\psi \underset{\sigma =+,}{}D\rho _\sigma \mathrm{exp}\left[H_{SG}\{\rho _\pm ;\psi \}\right]\delta \left[𝑑𝐫\rho _\sigma (𝐫)N\right]},`$ (18)
$`H_{SG}\{\rho _\pm ;\psi \}{\displaystyle 𝑑𝐫}{\displaystyle \frac{1}{8\pi l_B}}(\psi )^2+(\rho _+\rho _{})\psi +{\displaystyle \underset{\sigma =+,}{}}{\displaystyle 𝑑𝐫\rho _\sigma \mathrm{ln}\rho _\sigma }\rho _\sigma .`$ (19)
In these expressions, quadratic fluctuations of the density field around the saddle-point path $`\{\rho _\pm ^{sp}\}`$ are negligible as has been shown elsewhere . Therefore Gaussian approximation to $`\{\rho _\pm \}`$ reduces to the substitution of the Boltzmann distribution, $`\rho _\sigma ^{sp}=\rho _\sigma ^0\mathrm{exp}(\sigma \psi )(\sigma =+,)`$, into eq. (19), yielding

$`Z={\displaystyle D\psi \mathrm{exp}\left[H_{SG}\{\rho _\pm ^{sp};\psi \}\right]},`$ (20)
$`H_{SG}\{\rho _\pm ^{sp};\psi \}={\displaystyle 𝑑𝐫\frac{1}{8\pi l_B}(\psi )^2}+{\displaystyle \underset{\sigma =+,}{}}\rho _\sigma ^0e^{\sigma \psi },`$ (21)
where the total number invariance leads to $`\rho _\sigma ^0=N/𝑑𝐫e^{\sigma \psi }`$. The above transformation demonstrates that our formalism, eqs. (16) and (17), is a dual expression to the canonical SG theory, though there also exists a merit that our Hamiltonian (17) has two arguments (number densities of positive and negative charges) unlike the SG theory.

We now return to the starting formulae, (16) and (17), and consider the effect of quadratic fluctuations around the prescribed number densities $`n_\pm `$ on the electrostatic potential. Setting that $`\stackrel{~}{\rho }_\sigma =\rho _\sigma n_\sigma (\sigma =+,)`$ and $`\stackrel{~}{q}=\stackrel{~}{\rho }_+\stackrel{~}{\rho }_{}`$, and separating the imposed free energy $`F_{imp}`$, we have

$`Z=e^{F_{imp}}{\displaystyle \underset{\sigma =+,}{}}{\displaystyle D^{}\stackrel{~}{\rho }_\sigma e^{\stackrel{~}{H}_Z\{\stackrel{~}{\rho }_\pm \}}\delta \left[𝑑𝐫\stackrel{~}{\rho }_\sigma (𝐫)\right]},`$ (22)
$`\stackrel{~}{H}_Z\{\stackrel{~}{\rho }_\pm \}{\displaystyle \frac{l_B}{2}}{\displaystyle 𝑑𝐫𝑑𝐫^{}\frac{\stackrel{~}{q}(𝐫)\stackrel{~}{q}(𝐫^{})}{|𝐫𝐫^{}|}}+{\displaystyle \underset{\sigma =+,}{}}{\displaystyle 𝑑𝐫\frac{\stackrel{~}{\rho }_\sigma ^2}{2n_\sigma }}+\stackrel{~}{\rho }_\sigma (\sigma \mathrm{\Psi }_{imp}+\mathrm{ln}n_\sigma ),`$ (23)
where fixing positive charge 1 at $`𝐫_1`$ is taken into account by changing the definition of the integral measure as $`D^{}\stackrel{~}{\rho }_\sigma _{|𝐫_1𝐫|a}d\stackrel{~}{\rho }_\sigma `$, following the treatment of Fisher et al.. With use of expressions (22) and (23), however, Gaussian integration over $`\stackrel{~}{\rho }_\sigma `$ is not straightforward due to the spatial dependence of the prescribed densities $`n_\sigma (\sigma =+,)`$ other than conventional cases .

We then take a detour to introduce the induced potential $`\stackrel{~}{\psi }\psi \mathrm{\Psi }_{imp}`$ similarly to eqs. (18) and (19):

$`Z=e^{F_{imp}}{\displaystyle D\stackrel{~}{\psi }\underset{\sigma =+,}{}𝑑\nu _\sigma D^{}\stackrel{~}{\rho }_\sigma e^{\stackrel{~}{H}_{SG}\{\stackrel{~}{\rho }_\pm ;\nu _\pm ;\stackrel{~}{\psi }\}}},`$ (24)
(25)
$`\stackrel{~}{H}_{SG}\{\stackrel{~}{\rho }_\pm ;\nu _\pm ;\stackrel{~}{\psi }\}{\displaystyle 𝑑𝐫}{\displaystyle \frac{1}{8\pi l_B}}(\stackrel{~}{\psi })^2+(\stackrel{~}{\rho }_+\stackrel{~}{\rho }_{})\stackrel{~}{\psi }`$ (26)
$`+{\displaystyle \underset{\sigma =+,}{}}{\displaystyle ^{}}𝑑𝐫{\displaystyle \frac{\stackrel{~}{\rho }_\sigma ^2}{2n_\sigma }}+\stackrel{~}{\rho }_\sigma (\sigma \mathrm{\Psi }_{imp}+\mathrm{ln}n_\sigma ){\displaystyle 𝑑𝐫i\stackrel{~}{\rho }_\sigma \nu _\sigma },`$ (27)
where use has been made of the identity that $`\delta [𝑑𝐫\stackrel{~}{\rho }_\sigma (𝐫)]=𝑑\nu _\sigma \mathrm{exp}(i\nu _\sigma 𝑑𝐫\stackrel{~}{\rho }_\sigma )`$, and $`^{}𝑑𝐫`$ denotes the integration with omitting the region $`\mathrm{\Delta }r_1|𝐫_1𝐫|a`$. In these expressions Gaussian integration over $`\stackrel{~}{\rho }_\sigma `$ and $`\nu _\sigma `$ becomes trivial, though the resulting form is somewhat complicated:

$`Z=e^{F_{imp}}{\displaystyle D\stackrel{~}{\psi }e^{\stackrel{~}{H}_{SG}\{\stackrel{~}{\rho }_\pm ^{sp};\nu ^{sp};\stackrel{~}{\psi }\}}},`$ (28)
$`\stackrel{~}{H}_{SG}\{\stackrel{~}{\rho }_\pm ^{sp};\nu _\pm ^{sp};\stackrel{~}{\psi }\}`$ (29)
$`={\displaystyle 𝑑𝐫}{\displaystyle \frac{(\stackrel{~}{\psi })^2}{8\pi l_B}}+[\delta (𝐫𝐫_\mathrm{𝟏})Q(𝐫)]\gamma (\mathrm{\Delta }r_1)\stackrel{~}{\psi }`$ (30)
$`{\displaystyle \underset{\sigma =+,}{}}{\displaystyle ^{}}𝑑𝐫{\displaystyle \frac{n_\sigma }{2}}(\sigma \stackrel{~}{\psi }+\sigma \mathrm{\Psi }_{imp}+\mathrm{ln}n_\sigma )^2`$ (31)
$`+{\displaystyle \underset{\sigma =+,}{}}{\displaystyle \frac{1}{2}}{\displaystyle \frac{\left\{\stackrel{~}{\rho }_\sigma ^{fix}+^{}𝑑𝐫n_\sigma (\sigma \stackrel{~}{\psi }+\sigma \mathrm{\Psi }_{imp}+\mathrm{ln}n_\sigma )\right\}^2}{^{}𝑑𝐫n_\sigma }},`$ (32)
where the step function $`\gamma (\mathrm{\Delta }r_1)`$ is 1 for $`\mathrm{\Delta }r_1a`$ and $`0`$ otherwise, and we define $`\stackrel{~}{\rho }_\sigma ^{fix}(\sigma =+,)`$ as $`\stackrel{~}{\rho }_+^{fix}[\delta (𝐫𝐫_\mathrm{𝟏})n_+(𝐫)]\gamma (\mathrm{\Delta }r_1)`$ and $`\stackrel{~}{\rho }_{}^{fix}n_{}(𝐫)\gamma (\mathrm{\Delta }r_1)`$.

Equations (28) and (29) imply that the saddle-point paths satisfy the following relations:

$$\stackrel{~}{q}^{sp}=\underset{\sigma =+,}{}\sigma n_\sigma (\sigma \stackrel{~}{\psi }+\sigma \mathrm{\Psi }_{imp}+\mathrm{ln}n_\sigma i\nu _\sigma ^{sp}),$$
(34)
and

$$i\nu _\sigma ^{sp}=\frac{\stackrel{~}{\rho }_\sigma ^{fix}+^{}𝑑𝐫n_\sigma (\sigma \stackrel{~}{\psi }+\sigma \mathrm{\Psi }_{imp}+\mathrm{ln}n_\sigma )}{^{}𝑑𝐫n_\sigma },$$
(35)
where $`\stackrel{~}{q}^{sp}=\stackrel{~}{\rho }_+^{sp}\stackrel{~}{\rho }_{}^{sp}`$. We can easily check that these reduce to $`\stackrel{~}{q}^{sp}=2\overline{n}\stackrel{~}{\psi }`$ and $`i\nu _\sigma ^{sp}=\mathrm{ln}\overline{n}`$ in the conventional case of $`n_\sigma =\overline{n}`$. Such reduction demonstrates that our formulation correctly includes the original DH theory , and that $`i\nu _\sigma ^{sp}(\sigma =+,)`$ merely correspond to the chemical potentials $`\mu _\sigma `$ in equilibrium.

Thus expression (34) with (35) provides the physical insight into fluctuations around prescribed densities as follows: in general the induced charge density distribution $`\stackrel{~}{q}^{sp}`$ around supposed density favors not only the cancellation of additional electrostatic potential $`\stackrel{~}{\psi }`$ like the original DH theory, but also the countervailing of the chemical potential difference between imposed and equilibrium value,

$$\mathrm{\Delta }\mu _\sigma \sigma \mathrm{\Psi }_{imp}+\mathrm{ln}n_\sigma i\nu _\sigma ^{sp}.$$
(36)
Finally let us write down the DH-like equations generalized via density functional integral formalism. Saddle-point approximation to eqs. (28) and (29) yields

$`^2\stackrel{~}{\mathrm{\Psi }}_+(𝐫;𝐫_1)=4\pi l_B[\delta (𝐫𝐫_1)Q(𝐫)],|𝐫_\mathrm{𝟏}𝐫|a,`$ (37)
$`^2\stackrel{~}{\mathrm{\Psi }}_+(𝐫;𝐫_1)=\kappa ^2(𝐫)\stackrel{~}{\mathrm{\Psi }}_+(𝐫;𝐫_1)4\pi l_Bq_{ap}(𝐫),|𝐫_\mathrm{𝟏}𝐫|a,`$ (38)
where the saddle-point path of $`\{\stackrel{~}{\psi }\}`$ is identified with $`\stackrel{~}{\mathrm{\Psi }}_+`$ used in the GDH equations (1) and (2), and $`q_{ap}_{\sigma =+,}\sigma n_\sigma \mathrm{\Delta }\mu _\sigma `$ is the apparent charge arising from the chemical potential difference $`\mathrm{\Delta }\mu _\sigma `$ given by (36). Note also that eq. (38) corresponds to the equation that $`^2\stackrel{~}{\mathrm{\Psi }}_+=4\pi l_B\stackrel{~}{q}^{sp}`$ with use of eq. (34). Comparing our equations (37) and (38) with the GDH ones (1) and (2), one immediately finds that the former equations are different from the latter, in that the apparent charge $`q_{ap}`$ comes out on the right hand side of (38), though both the spatially dependent screening length $`\kappa ^1(𝐫)`$ and the cavity source term \[$`4\pi l_BQ`$ in eq. (37)\] can be reproduced indeed.

Accordingly the remaining problem is to investigate when the apparent charge $`q_{ap}`$ in eq. (38) disappears. A trivial condition is the chemical equilibrium $`\mathrm{\Delta }\mu _\sigma =0`$ where the imposed densities satisfy the Boltzmann distribution $`n_\sigma =\mathrm{exp}(\mu _\sigma \sigma \mathrm{\Psi }_{imp})`$. In such cases, though, what is imposed is not densities but potential, and the associated equations, already verified , are directly obtainable from the conventional SG theory, not via our formalism; this is not the case with us. A relevant case is to impose $`n_\sigma `$ on relation (5) with $`\mathrm{\Delta }\theta =0`$ as the GDH theory does, where $`\mathrm{\Psi }_{imp}=0`$ and the apparent charge $`q_{ap}`$ is negligible under the constraint that $`\nu _+^{sp}=\nu _{}^{sp}`$.

We thus conclude from the functional integral point of view the following: The GDH equations (1) and (2) proposed by Fisher et al. are valid for the prescribed densities of $`n_\sigma =\overline{n}[1+\mathrm{\Lambda }\mathrm{cos}(𝐤𝐫)]`$, i.e., $`\mathrm{\Delta }\theta =0`$, where the density-density correlation length has been extracted . However, in the other case of $`n_\sigma =\overline{n}[1+\sigma \mathrm{\Lambda }\mathrm{cos}(𝐤𝐫)](\sigma =+,)`$, i.e., $`\mathrm{\Delta }\theta =\pi `$, where both the charge-charge correlation length and the Lebowitz length have been derived , we have

$$q_{ap}2\overline{n}\mathrm{\Psi }_{imp}+Q+𝒪[\mathrm{\Lambda }^2],$$
(39)
and hence apparent charge is to be considered explicitly in the calculation of $`\stackrel{~}{\mathrm{\Psi }}_+`$; evaluating to what extent the additional term changes the previous results  remains a future problem.

We acknowledge the financial support from the Ministry of Education, Science, Culture, and Sports of Japan.


format: text

text:
# Phase separation and the possibility of orbital liquid states in CMR manganites. A 139La NMR study.

## Abstract

139La NMR spin-lattice relaxation rate $`1/T_1`$ and rf enhancement experiments provide evidence that the low temperature regime of the ferromagnetic (FM) phase of La1-xCaxMnO3 segregates into highly-conductive and poorly-conductive FM regions, associated with differences in the orbital structure. Remarkably, phase separation is accompanied with the appearance of an extra NMR signal from FM regions with vanishingly small magnetic anisotropy. This feature has been attributed to the appearance of regions with strong orbital fluctuations, resembling droplets of an orbital liquid within the inhomogeneous FM matrix.

The origin of the electronic properties of mixed valence manganites with the general formula R1-xDxMnO3 (R=rare earth, D=Ca, Ba, Sr), poses one of the most exciting open problems in the physics of strongly correlated electron systems. Initially, these properties were determined in the framework of the double exchange (DE) model , which is based on the strong Hund’s coupling between hopping $`e_g`$ and underlying $`t_{2g}`$ electrons in successive Mn4+,3+ sites. However, in the last years experiments have shown that DE is inadequate for the full description of the complex properties of these systems. For example, in the low dopping regime an important question concerns the origin of the ferromagnetic insulating (FMI) to ferromagnetic conductive (FMC) phase transition, which at first sight appears to contradict the DE model. Recent $`x`$-ray resonant scattering  and NMR  experiments have shown that this transition is controled rather by orbital than by spin degrees of freedom. Another question concerns the nature of the mixed FMI-FMC phase, which has been observed close to the corresponding phase boundary . According to theoretical modells including orbital ordering (OO) and Jahn-Teller distortions, such phase coexistence might indicate the presence of electronic phase separation , . However, it is still experimentally unclear whether the observed mixed phases possess the same intrinsic electronic properties as the constituting pure phases, or whether they form microscopically a ”third” novel electronic state.

In this letter, we demonstrate experimentally for the first time that phase separation, occuring in the low temperature regime of FM La based manganites, is accompanied with a significant change of the $`e_g`$ electron mobility at a local level. Specifically, we show that by cooling the high temperature ”homogeneous” FMC phase breaks into FMC and FMI regions with enhanced - respectively reduced - electron mobility in respect to the high temperature phase. Most important, the formation of the mixed phase is accompanied with the detection of FM regions with extremely low magnetic anisotropy, which is explained as showing the formation of regions with strong orbital fluctuations.

In order to label the various magnetic structures according to their electronic transport properties and OO we have employed 139La spin lattice relaxation, and 139La radiofrequency (rf) enhancement experiments. Particularly, we have shown that the abrupt decrease of the 139La NMR spin-lattice relaxation rate ($`1/T_1`$) in La1-xCaxMnO3 by decreasing temperature or increasing $`x`$, provides a direct measure of delocalization of the $`e_g`$ electrons. On the other hand, rf enhancement experiments may be used in order to monitor subttle changes of the magnetic anisotropy at a local level. In such experiments the NMR signal intensity $`I`$ is recorded as a function of the indensity $`H_1`$ or the duration $`\tau `$ of the applied rf field . The obtained curves follow an assymetric bell-shaped law with maximum at $`n\gamma H_1\tau =2\pi /3`$, which allows the calculation of the rf enhancement factor $`n`$ . Considering that the rf enhancement reflects the coherent response of the electron magnetic moments to an external rf fileld, it is clear that any variation of the magnetic anisotropy field $`H_A`$ will be directly reflected on $`n`$ .

Experiments were performed on polycrystalline samples prepared by annealing stoichiometric amounts of the corresponding oxides in air at $`1300`$ to $`1400`$ $`{}_{}{}^{0}C`$. All samples were then characterized structurally at room temperature with a D$`500`$ Siemens $`x`$-ray diffractometer, and magnetically with a SQUID magnetometer. The obtained crystallographic and magnetic data were found to be in accordance with literature. 139La NMR spectra in zero external magnetic field were acquired by applying a two pulse spin-echo technique, with pulse widths $`t_{p1}=t_{p2}=0.6`$ $`\mu sec`$ at very low rf power level, due to the very strong rf enhancement that characterizes FM materials . The obtained 139La spectra consist of broad lines located at $`20`$ MHz, which unambiguously are attributed to FM regions with fully polarized Mn spin octants , , . $`T_1`$ was then measured at the peak of the spectra, by applying a saturation recovery technique and fit with a multiexponential recovery law as in previous works .

Figure 1, exhibits 139La $`1/T_1`$ vs. $`x`$ measurements of La1-xCaxMnO3 at three different temperatures. It is observed that by increasing $`x`$, $`1/T_1`$ decreases abruptly by almost two orders of magnitude, at $`x0.175`$, which according to the literature defines the phase boundary between the FMI and FMC regimes of the $`Tx`$ phase diagram. This is a remarkable result, which shows that $`e_g`$ electrons play a major role in 139La spin-lattice relaxation phenomena. In order to describe qualitatively this behaviour, we have used the general approach invented by Warren for relaxation studies of metal-insulator transitions in liquid semiconductors . In general, nuclear spin-lattice relaxation in CMR manganites is determined by fluctuations of the local magnetic field $`H_{loc}`$ at the nuclear sites, (i) due to $`e_g`$ electron hopping, and (ii) due to interaction with spin waves. However, in systems with spin-polarized conduction bands, like CMR manganites , spin wave scattering is determined by weak two-magnon processes, and 139La $`1/T_1`$ should be dominated by direct $`e_g`$ electron spin flips, if energetically allowed. The latter assumption is supported by recent tight binding calculations in FM CMR manganites , which show that in the presence of electron-phonon coupling the minority $`t_{2g}()`$ band attains a small finite DOS at $`E_F`$, and spin flip between $`e_g()`$ and $`t_{2g}()`$ energy states is in principle possible. Following Warren’s approximation, and by assuming an exponentially decaying autocorrelation function for the fluctuations of $`H_{loc}`$, we may easily arrive at the expression  $`1/T_1TN(E_F)^2\tau _e`$/$`\left\{\mathrm{}N(E_F)\right\}`$, where $`\tau _e`$ has the significance of the lifetime residence of $`e_g`$ electrons on the Mn sites.

In case of delocalised electron states, where the electron mean-free path $`\mathrm{}>>a`$ ($`a`$ is the distance to the nearest neighbour), $`\tau _e=a/v_F=\mathrm{}N(E_F)`$ , and the above relation transforms to the standard Korringa relation, $`1/T_1TN^2(E_F)`$. However, in case of charge localization, i.e. for $`\tau _e>>\mathrm{}N(E_F)`$, $`1/T_1`$ is sufficiently higher than for delocalised electron states. Evidently, the abrupt reduction of $`1/T_1`$ in Figure 1 provides a direct evidence of conductivity enhancement, due to delocalization of the $`e_g`$ electron states. Such a behaviour is reminisence of a sharp metal-insulator transition of the Anderson type, and not of a percolative transition as proposed by many authors. Most notably, the onset of metallicity for $`x0.175`$ is characterized by a strong amplification of the rf enhancement factor, as recently reported by Dho et al. . This is clearly observed in the inset of Figure 1, which exhibits 139La NMR $`I`$ vs. $`H_1`$ curves of La1-xCaxMnO3 for $`x=0.1`$, and $`0.33`$, at $`T=5`$K. For $`x=0.33`$ the maximum of the $`I`$ vs. $`H_1`$ curve shifts to extremely low $`H_1`$ values, which indicates a considerable reduction of the magnetic anisotropy . Considering that (i) $`H_A`$ is mainly determined by the $`LS`$ coupling, and (ii) according to previous works on Sr-dopped systems  the metal-insulator transition correlates with OO changes, we argue that the observed change in the rf enhancement reflects a substantial modification in the orbital structure.

Figure 2 demonstrates 139La $`1/T_1`$ vs. $`T`$ measurements for $`x=0.25`$. It is observed that by decreasing temperature $`1/T_1`$ decreases almost exponentially, in agreement with previous works in relevant materials . For $`T<30`$K $`1/T_1`$ exhibits an almost stepwise decrease, which may be assigned to enhancement of the $`e_g`$ electron hopping below that temperature. Remarkably, the steep descent of $`1/T_1`$ is accompanied with the appearance of 55Mn NMR signal from FMI regions, as shown in ref. . Such FMI regions can not be detected by 139La NMR for $`x0.175`$, because signals from FMI regions with extremely low $`T_1`$s are masked from the much stronger signals from FMC regions. Apparently, the overall resistivity enhancement below $`30`$K (inset of Figure 2) may be attributed to restriction of the conductivity channels below that temperature. These results demonstrate that for $`T30`$K the system splits into regions with enhanced - respectively reduced - conductivity in comparison to the high temperature FMC phase.

A similar $`1/T_1`$ behaviour is followed by $`x=0.33`$ and $`0.5`$, as observed in Figure 3. For reasons of comparison, the 55Mn $`1/T_1`$ vs. $`T`$ curve for $`x=0.33`$, measured at the peak of the FMC signal at $`\nu 375`$ MHz , is also shown in the inset of Figure 3. In case of $`x=0.33`$, the deviation from the high temperature $`1/T_1`$ vs. $`T`$ law is realized below $`80`$K, whereas for $`x=0.5`$ a similar behaviour is observed; however, the experimental error of the $`1/T_1`$ data points above $`100`$K does not allow precise determination of the pertinent temperature, which defines the inflection point. It is worth to notice that for $`x=0.5`$ the reduction of $`1/T_1`$ at low temperatures takes place deep into the antiferromagnetic phase, hence it concerns enhancement of the electron mobility within the FMC islands that persist below $`T_N`$ , .

In order to examine whether the temperature dependence of the local electron mobility is associated with modifications of the OO, detailed $`I`$ vs. $`H_1`$ measurements for $`x=0.25`$, $`0.33`$, and $`0.5`$ were performed as a function of temperature (Figure 4). In case of $`x=0.25`$ and $`0.33`$, by decreasing temperature the maximum of the $`I`$ vs. $`T`$ curves is clearly observed to shift to lower $`H_1`$ fields. Remarkably, the abrupt reduction of $`1/T_1`$ is accompanied with the appearance of a second maximum in the $`I`$ vs. $`H_1`$ curves at extremely low $`H_1`$ values. This is a direct evidence that the low temperature regime of the FMC phase contains at least two distinct FM electron spin configurations, differing in $`H_A`$ and therefore $`OO`$ . By increasing $`x`$, the low-$`H_1`$ maximum remains unshifted, whereas the initial maximum shifts to higher rf power levels. It is also notable that for $`x=0.25`$ and $`0.5`$, the second maximum appears simultaneously with the onset of the mixed phase, i.e. below $`30`$K for $`x=0.25`$ , and $`140`$K for $`x=0.5`$ . These results corroborate with (i) recent 119Sn Moessbauer measurements in $`1\%`$ Sn-doped La0.5Ca0.5MnO3, which indicate the appearance of a third FM signal component for $`T<120K`$ , and (ii) magnetoresistance experiments in thin films of La0.67Ca0.33MnO3 , which show that below $`100`$K magnetoresistance anisotropy changes rapidly from two-fold to four-fold.

In view of the up-to-now discussion, the following picture may be envisaged: For $`x0.175`$ in the low temperature regime of the $`Tx`$ magnetic phase diagram the system breaks into FMC regions with high conductivity and FMI regions with low conductivity, which differ mainly in the orbital structure. In the FMC regions more than one orbital configurations may be realized, as implied by the subsequent shift of the prime maximum in the $`I`$ vs. $`H_1`$ curves, by varying dopping or temperature. This is in conformity with recent theoretical calculations , which have shown that by increasing $`x`$ the FM state is realized, whereas OO may change continuously from orbital $`G`$ $`(x^2y^2)/(3z^2r^2)`$ near $`x=0`$, to orbital $`C`$ $`(x^2y^2)/(3z^2r^2)`$ for $`x0.3`$, and orbital A $`[(3z^2r^2)+(x^2y^2)]/[(3z^2r^2)(x^2y^2)]`$ for $`0.4x0.7`$. Remarkably, the calculated energy difference among these states is very small , suggesting the presence of strong orbital fluctuations. In such a case, the formation of islands of an orbital liquid with extremely low $`H_A`$ is possible. This might give an explanation to the appearance of the maximum at extremely low $`H_1`$ field in the $`I`$ vs. $`H_1`$ experiments. In this scenario, by increasing temperature the onset of Jahn-Teller distortions  might lift the degeneracy of the various states, and the orbital liquid would freeze in the high temperature FMC orbital configuration.

In conclusion, 139La $`1/T_1`$ and rf enhancement measurements in La-based manganese perovskites provide clear evidence that differences in OO underly phase segregation tendencies into FM states with high (respectively poor) $`e_g`$ electron mobility. At low temperatures where the effect of Jahn-Teller distortions is minimized, NMR experiments suggest that degeneracy in the energy of these states may lead to the appearance of orbital liquid regions, in agreement with recent theoretical predictions .

Acknowledgment This work has been partially supported by the INTAS project 97-30253.


text:
# Buckling of rods with spontaneous twist and curvature

## 1 Equilibrium equations

Consider an elastic rod of length $`L`$ and cross-section $`𝒮`$ which has two axes of symmetry. We assume that $`L`$ $`a`$ where $`a`$ is the largest dimension associated with the cross-section. Denote by $`s`$ the arc-length which parameterizes the centerline of the rod in the stress-free reference configuration. The radius vector $`𝐑_0(s)`$ describes the spatial position of the centerline and $`𝐭_0=d𝐑_0/ds`$ the unit tangent vector in the reference state. The unit vectors $`𝐮_0(s)`$ and $`𝐯_0(s)`$ are directed along the axes of symmetry of the cross-section. We would like to emphasize that $`𝐮_0`$ and $`𝐯_0`$ are associated with physical attributes of rods of small but non-vanishing thickness and, in general, do not coincide with the normal and the bi-normal vectors which appear in the description of three-dimensional configurations of infinitely thin geometrical lines. Differentiation of the vectors $`𝐭_0`$, $`𝐮_0`$ and $`𝐯_0`$ with respect to $`s`$ results in

$$\frac{d𝐭_0}{ds}=\omega _{10}𝐮_0+\omega _{20}𝐯_0,\frac{d𝐮_0}{ds}=\omega _{10}𝐭_0+\omega _{30}𝐯_0,\frac{d𝐯_0}{ds}=\omega _{20}𝐭_0\omega _{30}𝐮_0,$$
(1)
where $`\omega _{k0}(s)`$ are generalized intrinsic torsions. The position of an arbitrary point of the rod is given by the Lagrangian coordinates $`\{\xi _k\}`$, where $`\xi _1,\xi _2`$ are Cartesian coordinates in the cross-sectional plane with unit vectors $`𝐮_0`$ and $`𝐯_0`$, and $`\xi _3=s`$,

$$𝐫_0(\xi _1,\xi _2,s)=𝐑_0(s)+\xi _1𝐮_0(s)+\xi _2𝐯_0(s).$$

We assume that the centerline is inextensible, which means that $`s`$ remains the arc-length in the deformed state as well. The position of the centerline in the deformed configuration is determined by its radius vector $`𝐑(s)`$. The position of an arbitrary point of the cross-section with the arc-length $`s`$ is

$$𝐫(\xi _1,\xi _2,s)=𝐑(s)+\xi _1𝐮(s)+\xi _2𝐯(s),$$
(2)
where $`𝐮`$ and $`𝐯`$ are unit vectors directed along the axes of symmetry of the cross-section, and $`𝐭=d𝐑/ds`$ is the unit tangent vector. The functional form of Eq. (2) is based on the assumption that after deformation any cross-section remains planar and perpendicular to the centerline of the rod. By analogy with Eq. (1), we write

$$\frac{d𝐭}{ds}=\omega _1𝐮+\omega _2𝐯,\frac{d𝐮}{ds}=\omega _1𝐭+\omega _3𝐯,\frac{d𝐯}{ds}=\omega _2𝐭\omega _3𝐮,$$
(3)
where $`\omega _k(s)`$ are generalized torsions.

Using the methods of Ref. , we calculate the displacement gradient in the deformed state $`𝐫_0/𝐫`$ and the strain tensor. The mechanical energy is constructed as a quadratic form in this strain tensor. Expressions for the components of the stress tensor are derived by differentiating this energy with respect to the corresponding components of the strain tensor. The constitutive equations relate the components of the internal moment $`𝐌=M_u𝐮+M_v𝐯+M_t𝐭`$ to deviations of the generalized torsions from their spontaneous values,

$$M_u=A_2(\omega _2\omega _{20}),M_v=A_1(\omega _1\omega _{10}),M_t=A_3(\omega _3\omega _{30}),$$
(4)
where $`A_1=E_1I_1`$, $`A_2=E_1I_2`$, $`A_3=E_2(I_1+I_2)`$, $`I_1`$ and $`I_2`$ are moments of inertia with respect to the two symmetry axes of the cross-section, and $`E_1`$ and $`E_2`$ are elastic moduli. Denote the internal force by $`𝐅=F_u𝐮+F_v𝐯+F_t𝐭`$, the external distributed force per unit length by $`𝐪=q_u𝐮+q_v𝐯+q_t𝐭`$, and the external distributed moment per unit length by $`𝐦=m_u𝐮+m_v𝐯+m_t𝐭`$. Upon some algebra, we eliminate the components of the internal force in the cross-sectional plane and write down the conditions of mechanical equilibrium,

$`{\displaystyle \frac{dF_t}{ds}}+\omega _1{\displaystyle \frac{dM_v}{ds}}\omega _2{\displaystyle \frac{dM_u}{ds}}+\omega _3(\omega _1M_u+\omega _2M_v)+\omega _1m_v\omega _2m_u+q_t=0,`$
$`{\displaystyle \frac{dM_t}{ds}}\omega _1M_u\omega _2M_v+m_t=0,`$
$`{\displaystyle \frac{d}{ds}}({\displaystyle \frac{dM_v}{ds}}+\omega _3M_u+\omega _2M_t+m_v)+\omega _3({\displaystyle \frac{dM_u}{ds}}\omega _3M_v+\omega _1M_t+m_u)`$
$`\omega _1F_tq_u=0,`$
$`{\displaystyle \frac{d}{ds}}({\displaystyle \frac{dM_u}{ds}}\omega _3M_v+\omega _1M_t+m_u)\omega _3({\displaystyle \frac{dM_v}{ds}}+\omega _3M_u+\omega _2M_t+m_v)`$
$`+\omega _2F_t+q_v=0.`$ (5)
Given boundary conditions, the four nonlinear differential equations (5) together with the constitutive relation (4), uniquely determine the four unknown functions, $`F_t`$, $`\omega _1`$, $`\omega _2`$ and $`\omega _3`$. The generalized torsions $`\omega _k`$ are substituted into Eqs. (3) which are resolved to determine the spatial configuration of the centerline and the rotation of the cross-section (twist) of the deformed rod. The tangential force $`F_t`$ plays the role of a Lagrange multiplier which enforces the inextensibility condition.

## 2 Stability of a rectilinear rod

We begin with the analysis of stability of a spontaneously twisted rectilinear rod. In the stress-free reference state, the rod’s centerline is straight, and its cross-sections rotate around the centerline by some angle $`\mathrm{\Omega }=\mathrm{\Omega }(s)`$. We set $`A_1<A_2`$, which means that the vector $`𝐮_0`$ coincides with the smaller of the two principal axes of inertia. The vectors $`𝐮_0`$, $`𝐯_0`$ and $`𝐭_0`$ are given by $`𝐮_0=\mathrm{cos}\mathrm{\Omega }𝐞_1+\mathrm{sin}\mathrm{\Omega }𝐞_2`$, $`𝐯_0=\mathrm{sin}\mathrm{\Omega }𝐞_1+\mathrm{cos}\mathrm{\Omega }𝐞_2`$, and $`𝐭_0=𝐞_3`$, where $`𝐞_k`$ are unit vectors of some Cartesian coordinate frame associated with the stress-free rod. These equalities and Eq. (1) imply that

$$\omega _{10}=0,\omega _{20}=0,\omega _{30}=\frac{d\mathrm{\Omega }}{ds}.$$

We will assume in the following that $`\omega _{30}`$ is independent of $`s`$.

When some compressive force $`P`$ (smaller than the critical force for buckling) is applied to the ends of an inextensible rod, the rod does not deform, the moments $`M_u`$, $`M_v`$ and $`M_t`$ vanish, and only a longitudinal internal force $`P`$ is generated. We now assume that, in addition to the load $`P`$, small perturbations of external distributed forces and external distributed moments are present. These perturbations cause small displacements of the rod that consist of small displacements of the centerline, $`\mathrm{\Delta }𝐑=x_1𝐞_1+x_2𝐞_2`$, and small rotations of the cross-sections around the centerline $`\alpha `$. Under this perturbation, the vectors $`𝐮_0`$, $`𝐯_0`$ and $`𝐭_0`$ are replaced by $`𝐮`$, $`𝐯`$, $`𝐭`$, which are substituted into Eq. (3) to obtain expressions for the perturbed generalized torsions $`\omega _k`$ (linear in the perturbations $`x_1`$, $`x_2`$ and $`\alpha `$). Substitution of these expressions into Eq. (4) yields the perturbations of internal moments,

$$M_u=A_2\left(\frac{d^2x_1}{ds^2}\mathrm{sin}\mathrm{\Omega }\frac{d^2x_2}{ds^2}\mathrm{cos}\mathrm{\Omega }\right),M_v=A_1\left(\frac{d^2x_1}{ds^2}\mathrm{cos}\mathrm{\Omega }+\frac{d^2x_2}{ds^2}\mathrm{sin}\mathrm{\Omega }\right),M_t=A_3\frac{d\alpha }{ds}.$$

Substituting these expressions into Eqs. (5), we find that the equations for the perturbation of the longitudinal force and for the rotation $`\alpha `$ decouple from the equations for the curvatures $`\omega _1`$ and $`\omega _2`$, which govern bending of the rod and determine it stability with respect to buckling.

When the external load $`P`$ reaches its critical value, the rod becomes unstable against arbitrarily small external perturbations of moments and forces and, therefore, buckling can be interpreted as the point at which a non-vanishing solution to Eqs. (5), in which these forces and moments are set to zero, first appears. This yields

$`A_1{\displaystyle \frac{d^2\omega _1}{ds^2}}+(PA_1\omega _{30}^2)\omega _12A_2\omega _{30}{\displaystyle \frac{d\omega _2}{ds}}=0,`$
$`A_2{\displaystyle \frac{d^2\omega _2}{ds^2}}+(PA_2\omega _{30}^2)\omega _2+2A_1\omega _{30}{\displaystyle \frac{d\omega _1}{ds}}=0.`$ (6)
Equations (6) are supplemented by the geometric relations which follow from Eqs. (3):

$$\frac{d^2x_1}{ds^2}=\omega _1\mathrm{cos}\mathrm{\Omega }\omega _2\mathrm{sin}\mathrm{\Omega },\frac{d^2x_2}{ds^2}=\omega _1\mathrm{sin}\mathrm{\Omega }+\omega _2\mathrm{cos}\mathrm{\Omega }.$$

The resulting linear equations are solved subject to the boundary conditions corresponding to hinged ends:

$$x_1(0)=x_1(L)=x_2(0)=x_2(L)=0,\frac{d^2x_1}{ds^2}(0)=\frac{d^2x_1}{ds^2}(L)=\frac{d^2x_2}{ds^2}(0)=\frac{d^2x_2}{ds^2}(L)=0.$$

When the spontaneous twist $`\omega _{30}`$ vanishes, Eqs. (6) turn into conventional equations for stability of an elastic beam . In the general case, Eqs. (6) may be resolved with respect to $`\omega _1`$ and $`\omega _2`$. Each of these functions obeys a linear fourth-order differential equation which contains only even derivatives with respect to $`s`$ and can, therefore, be written as a linear combination of $`\mathrm{sin}`$($`\beta _ks`$) and $`\mathrm{cos}(\beta _ks)`$, where the “wavenumbers” $`\beta _1`$ and $`\beta _2`$ are the positive roots of a characteristic equation that guarantees the existence of non-vanishing solutions $`\omega _1(s)`$ and $`\omega _2(s)`$ to Eqs. (6). Using standard methods, we derive a transcendental equation which relates the dimensionless critical force $`p^{}=P_{\mathrm{cr}}/P_{\mathrm{cr}}`$, where $`P_{\mathrm{cr}}=\pi ^2A_1/L^2`$ is the Euler critical force for an untwisted rod , to the total spontaneous twist angle $`\mathrm{\Omega }(L)=\omega _{30}L`$. The numerical solutions of these equations for various values of the cross-section asymmetry $`\nu =I_2/I_1`$ are plotted in Figure 1. Note that at the point of the instability, all generalized torsions $`\omega _k`$ do not vanish, which implies that, unlike the planar Euler instability, buckling of a spontaneously twisted rod is intrinsically three-dimensional.

For a circularly symmetric cross-section ($`\nu =1`$), the Euler result holds for arbitrary values of the spontaneous twist, i.e., $`p^{}`$ is independent of $`\mathrm{\Omega }`$. For any $`\nu >1`$, the critical force increases initially with $`\mathrm{\Omega }(L)`$ and reaches a maximum at $`\mathrm{\Omega }(L)=2\pi `$. Upon further increase in the total angle of twist, the critical force oscillates (with additional maxima at integer multiples of $`\pi `$) and eventually tends to some limiting value which depends on $`\nu `$. To understand the origin of this non-monotonic behavior, we investigated the behavior of critical wavenumbers $`\beta _k`$. Numerical calculations show that the difference between the critical wavenumbers at buckling, $`\beta _2\beta _1`$, oscillates as a function of $`\mathrm{\Omega }(L)`$, with maxima occurring at integer multiples of $`\pi `$, and tends asymptotically to $`2\pi `$. In this asymptotic limit the critical wavenumbers approach the “resonance” conditions $`\beta _kL`$ $`\omega _{30}L\pm \pi `$ ($`k=1,2`$). We conclude that the oscillatory dependence of the critical force on the twist angle arises as a result of the interplay between the two modes. Whenever the total angle of twist coincides with an integer multiple of $`\pi `$, the corresponding axes of symmetry of the cross–sections at both ends of the rod (at $`s=0`$ and $`s=L`$) become aligned, just as in the twist-free case.

## 3 Stability of a ring

We now turn to study the stability of a circular ring whose centerline in the stress-free state is a planar circle with radius $`R_0`$. We introduce cylindrical coordinates $`\{r,\varphi ,z\}`$ with unit vectors $`𝐞_r`$, $`𝐞_\varphi `$, $`𝐞_z`$, where the axis $`z`$ is perpendicular to the plane of the ring. The initial position of the centerline in the stress-free reference state is determined by the radius vector $`𝐑_0=R_0𝐞_r`$. We assume that in its stress free reference state the rod is characterized by a twist angle $`\mathrm{\Omega }=\mathrm{\Omega }(s)`$ between the vectors $`𝐞_r`$ and $`𝐮_0`$. The unit vectors $`𝐭_0`$, $`𝐮_0`$ and $`𝐯_0`$ can be expressed in terms of $`𝐞_r`$, $`𝐞_\varphi `$ and $`𝐞_z`$ as

$$𝐭_0=𝐞_\varphi ,𝐮_0=\mathrm{cos}\mathrm{\Omega }𝐞_r+\mathrm{sin}\mathrm{\Omega }𝐞_z,𝐯_0=\mathrm{sin}\mathrm{\Omega }𝐞_r+\mathrm{cos}\mathrm{\Omega }𝐞_z.$$

Substitution of these expressions into Eqs. (1) results in

$$\omega _{10}=\frac{1}{R_0}\mathrm{cos}\mathrm{\Omega },\omega _{20}=\frac{1}{R_0}\mathrm{sin}\mathrm{\Omega },\omega _{30}=\frac{d\mathrm{\Omega }}{ds}.$$

We study the stability of the twisted ring under a radially inward-directed force per unit length $`𝐪=q𝐞_r`$ that acts on the outer perimeter of the ring. As long as the force is smaller than some critical value, it is balanced by a tangential internal force $`F_t=qR_0`$ and no displacement arises. We consider a small deviation of the centerline from its unperturbed position, $`𝐑_0𝐑=𝐑_0+y_r𝐞_r+y_\varphi 𝐞_\varphi +y_z𝐞_z`$, and a small perturbation of the rotation angle about the centerline, $`\mathrm{\Omega }\mathrm{\Omega }+\alpha `$. We express the perturbed unit vectors $`𝐭`$, $`𝐮`$ and $`𝐯`$ in terms of the small deviations $`y_\varphi `$, $`y_z`$ and $`\alpha `$ ($`y_r`$ is eliminated using the inextensibility condition $`|𝐭|=1`$) and use Eq. (3) to relate these deviations to the perturbed spontaneous torsions $`\omega _k`$. The subsequent analysis is qualitatively similar (but much more cumbersome) to that leading to Eqs. (6). We derive three coupled differential equations with variable coefficients for the functions $`y_\varphi (s)`$, $`y_z(s)`$ and $`\alpha (s)`$. An analytical solution of these equations can be obtained only in the special case of equal moments of inertia, $`I_1=I_2`$ (e.g., for a circular cross-section). In this case we find that buckling takes place when the critical force per unit length reaches a critical value

$$q_{\mathrm{cr}}=\frac{3E_1I_1}{R_0^3}.$$

As expected for a circular rod, the critical force does not depend on the spontaneous twist $`\mathrm{\Omega }`$. In this case, the onset of instability with respect to in-plane buckling of a circular rod coincides with that for out-of-plane buckling.

In the general case, $`I_1<I_2`$, we have to resort to approximate methods. For simplicity, we set $`\mathrm{\Omega }(s)=Ns/R_0`$, where $`N`$ is a positive integer. A reasonable estimate which is expected to capture the qualitative aspects of the effect of spontaneous twist on the critical force, can be obtained using the Galerkin method . The perturbations $`y_\varphi `$, $`y_z`$ and $`\alpha `$ are presented in the form $`X_n\mathrm{sin}(ns/R_0)`$ with a fixed $`n`$. We substitute these expressions into the equilibrium equations, multiply both sides by $`\mathrm{sin}(ns/R_0)`$ and require that the difference between the left and the right hand sides of the resulting equalities vanishes on the average (i.e., when integrated over the ring contour, from $`s=0`$ to $`2\pi R_0`$). Upon some algebra, we arrive at the following expression for the critical compressive force per unit length

$$q_{\mathrm{cr}}=\frac{E_2(I_1+I_2)}{2R_0^3}\left[3\frac{I_2I_1}{I_2+I_1}(N^2+1)\right].$$
(7)
Note that for $`N^23(I_1+I_2)/(I_2I_1)1`$ the ring becomes unstable with respect to buckling, even in the absence of external forces. Because the equations for the perturbations are coupled, the instability has a three dimensional character and buckling is not confined to the plane of the ring.

It is of interest to compare these conclusions with the results of Ref. , where the stability of the stress-free configuration of a circular ring was analyzed and it was found that, in general, there were no non-trivial solutions to the equilibrium equations with periodic boundary conditions. The apparent discrepancy stems from the fact that the problems under consideration are rather different: while we examine the stability of a ring with spontaneous twist under compression in the plane of the ring, Ref.  analyzes the stability of a ring without spontaneous twist under the action of an external torque applied to a fixed cross-section.

## 4 Conclusions

Buckling instabilities have been studied for thin rods with spontaneous twist and non-circular cross-sections. Two cases are analyzed in detail: (i) a rectilinear rod compressed by forces applied to its ends, and (ii) a circular ring compressed by uniformly distributed pressure acting in its plane. It is demonstrated that for a rectilinear rod, an increase in the spontaneous twist angle leads to the growth of the critical force and stabilizes the rod against buckling. The effect increases with the ratio of principal moments of inertia, $`\nu =I_2/I_1>1`$. For a circular ring, an increase in the spontaneous twist leads to a monotonic decrease in the critical force. The rate of the decrease grows with $`\nu `$. For sufficiently large values of the spontaneous twist, the ring may become unstable to buckling even in the absence of external forces.

### Acknowledgement

We would like to thank S. Panyukov for helpful comments. AD gratefully acknowledges financial support by the Israeli Ministry of Science through grant 1202–1–98. YR’s work is supported by a grant from the Israel Science Foundation.


text:
# Optimal Mutation Rates in Dynamic Environments

## I Introduction

In any given environment the vast majority of mutations that have any effect on the fitness of a biological organism are deleterious. One might expect the damaging effect of non-zero mutation rates to imply that when under evolutionary control the lowest mutation rate compatible with physiological constraints should be selected for. However, when examined experimentally bacteria and viruses (and indeed all organisms) have significant non-zero rate, the magnitude and diversity of which have failed to find satisfactory theoretical explanation. Some results from a number of experiments measuring the mutation rates of a selection of small DNA-based organisms are shown in Table I.

Despite the huge variation in genome length over four orders of magnitude the mutation rate per genome and replication $`\mu _G`$ remains constant roughly within a factor of roughly $`2`$ (which is at the same level as the estimated accuracy of the figures). As pointed out by Drake and others  this constancy in $`\mu _G`$ is surprising since DNA/RNA repair and transcription are primarily local processes that act on individual bases. Thus the data strongly suggest that point mutation rates for the different organisms have evolved towards individual optimal values that result in almost constant genomic copying fidelity.

In this paper we develop a model of the evolution of mutation rates based on changing environments. The evolved point mutation rate of this model scales so that the genomic copying fidelity is approximately independent of genome length and insensitive to other parameters in the model. The evolved mutation rates are also of the same magnitude as observed in Table I for biologically plausible parameter settings. We also suggest a possible explanation for the high mutation rates of RNA viruses. Simulations confirm the predictions of the model.

## II Evolving Mutation Rates

It is impossible to perfectly maintain and copy genetic information. All molecules, including DNA and RNA are thermodynamically unstable, and their physical structure and hence the information they encode changes over time. In addition the binding sites of enzymes such as DNA polymerase are not perfectly specific and errors will be introduced during replication. Lowering the error rate requires the use of increasingly complex proof-reading and repair mechanisms, all of which ultimately impose an energetic, and hence fitness, cost on the organism. We can expect a balance to develop between the pressure to lower mutation rates due to the fitness cost of deleterious mutants and the physiological cost of high copying accuracy . Such a balance certainly provides an ultimate lower limit to the mutation rate of all organisms but explaining the concstancy in genomic copying fidelity using such arguments causes unnatural assumptions on the relation between cost of local copying fidelity and genome length. There is also little experimental evidence that mutation rates are actually determined by such a balance.

When viewed as a whole the genome encodes not only proteins that directly influence its reproductive or survival ability, but also the copying fidelity with which the genome reproduces. For example some viroids contain genes that are translated into surface coat proteins while others genes code for the replicase enzymes that perform the copying of its genetic material. In more complex organisms additional genes may encode for modifiers of the accuracy of copy and repair enzymes, usually increasing mutation rates , but sometimes resulting in a decrease . These modifiers can have large or small effects on mutation rate and affect individual bases or the entire genome .

One consequence of this flexibility of mutation rates and their encoding is that if there are random changes (mutations) in genes determining the mutation rate then the copying fidelity will itself undergo Darwinian evolution.

## III Population Genetics in Changing Environments

When comparing two haploid genomes, the one with lower mutation frequency will produce offspring that are on average more closely related to itself. This means that for an asexual haploid replicator evolving on a static fitness landscape the optimal mutation rate for a sequence whose fitness is already globally maximal is zero. If the fitness peak moves, however, the situation changes: to avoid extinction a genome with an initially superior fitness is forced to accept a non-zero mutation rate to survive. This leads to a non-trivial optimal copying fidelity.

Kimura formalized the evolutionary effect of a changing environment by considering the genetic load of a population : the proportion by which the population fitness is decreased in comparison with an optimum genotype. Genetic load results from a number of competing factors; most notably the *mutational load* due to the deleterious effects of most mutations and the *segregational load* due to the temporary reduction in fitness that occurs whenever the selective environment changes. Assuming that a population minimizes the genetic load, the optimal mutation rate can be calculated. Using a descrete time model, i.e. a model where there is no overlap between generations, with one fitness determining locus the optimal mutation rate becomes:

$$\mu _{opt}=\frac{1}{\tau }$$
(1)
where $`\tau `$ is the number of generations between environmental changes. This model only considers the effect of mutations on the population and is therefore based on group selection.

Later population genetic models that examined competition between genetic modifiers of the mutation rate demonstrated that (for haploids with a single fitness determining locus) a non-zero mutation rate comes to dominate a population evolving in an oscillating environment . These models are not built on group selection. However a general and simple to interpret multi-locus modifier model does not exist.

## IV The Model

We will explore a more general model of the evolution of mutation rates in a dynamic environment. Consider a population of haploid genomes where a genome consists of two separated parts, one coding for the fitness and one coding for the probability per base $`\mu `$ of an error occurring during copying. There is complete linkage (no recombination) between the sections of the genome that encode the mutation rate and those that determine the fitness. We also assume that the fitness determining region is of fixed length $`\nu `$. In general we are interested in the fates of certain genomes $`g_i`$ which have a (possibly time-dependent) fitness advantage $`\sigma (t)`$ over all other sequences. We call these genomes master-sequences. The genomic copying fidelity of the fitness determining region of each strain $`g_i`$ is $`Q_i=(1\mu _i)^\nu `$, the index $`i`$ refers to the mutation rate of the strain, different strains have different mutation rates but identical fitness $`\sigma `$. We assume that mutations do not affect the copying fidelity, only the fitness. Changes to the mutation rates occur on a time-scale significantly slower than the time it takes for the population to reach equilibrium. During a period when a specific sequence has superior fitness compared to the background (i.e. between environmental shifts) the changes in the relative concentrations $`x_i`$ of the master-sequences are described by the replicator equation

$`\dot{x}_i(t)`$ $`=`$ $`Q_i\sigma (t)x_i(t)f(t)x_i(t)`$ (2)
where $`f(t)=\sigma (t)_jQ_jx_j(t)`$ normalizes the relative concentrations of the master-sequence strains. Mutations from background sequences onto the strains with optimal fitness are ignored. Since we are only interested in competition between master-sequences the background is not explicitly expressed in these equations.

The environment changes as follows: for a time $`t[0,\tau _1]`$ one genotype has superior fitness, followed by a new gene-sequence for time $`t[\tau _1,\tau _1+\tau _2]`$, etc. The notation is chosen so that $`\tau `$ denotes lengths of time intervalls. We assume that the initial concentration of the new master-sequences $`x_i`$ immediately after the shift (at time $`t_a=_{i=1}^m\tau _i+ϵ`$, where $`m`$ denotes shifts of the fitness-peak and $`ϵ`$ is am infinitely small time-period) are proportional to the concentrations of the old master-sequence before the shift (at $`t_b=_{i=1}^m\tau _iϵ`$)

$`x_i(t_a)`$ $`=`$ $`h(\mu _i)x_i(t_b)`$ (3)
It is reasonable to assume that $`h(\mu _i)`$ is a function with Taylor-expansion in the mutation rate $`\mu `$

$`h(\mu )`$ $`=`$ $`{\displaystyle \underset{j=k_m}{\overset{\mathrm{}}{}}}\alpha _j\mu ^j`$ (4)
where $`k_m`$ is a measure of the environmental change, i.e. the number of point mutations needed to transform the old superior sequence into the new. This basically means that $`k_m`$ is the Hamming distance from the old peak to the new at shift $`m`$. The constants $`\alpha _j`$ are combinatorial factors. It will turn out that the optimal mutation rate is independent of these factors.

To analyze the long term behavior of this system we make a change of variables $`y_i(t)=e^{_0^tf(s)𝑑s}x_i(t)`$. The new system of differential equations is linear and the equations are decoupled (due to the assumption that the selective dynamics is significantly faster than the changes in mutation rate), it is therefore easy to find the analytical solution:

$`y_i(t)`$ $`=`$ $`y_i(0)e^{Q_i_0^t\sigma _m(s)𝑑s}`$ (5)
Since $`x_i`$ is propotional to $`y_i`$, maximizing the growth of $`y_i`$ and $`x_i`$ are equivalent. After a suitably long time interval the population will be completely dominated by genomes that have a mutation rate closest to the optimal value $`\mu _{opt}`$ which maximizes the long term growth of the strain

$$\underset{\mu }{\mathrm{max}}\left(\mathrm{\Pi }_mh(\mu )e^{(1\mu )^\nu \sigma _m\tau _m}\right)$$
(6)
where $`_m`$ denotes a time average during time-period $`m`$. Setting the derivative of this expression to zero and using Eq. 4 we find the optimal copying fidelity to be approximately

$`\mu _{opt}`$ $`=`$ $`{\displaystyle \frac{k}{\nu \sigma \tau }}`$ (7)
where $``$ denotes a time average over all time periods. We also assume no correlation between $`\sigma _m`$ and $`\tau _m`$. Since the genome lengths is large $`\nu 1`$, the optimal copying fidelity and mutation rate per genome become:

$`Q_{opt}`$ $`=`$ $`e^{\frac{k}{\sigma \tau }}`$ (8)
$`\mu _G`$ $`=`$ $`{\displaystyle \frac{k}{\sigma \tau }}`$ (9)
Thus we find that the genomic optimal copying fidelity is independent of the genome length for fairly general types of environmental change in both the advantage of the fittest genotype $`\sigma (t)`$ and the size of environmental shifts $`h(\mu )`$.

## V Simulations

To confirm the theoretical derivations we simulated the evolution of replicators in continuous time on a moving single peaked landscape using a birth-death process. Each time unit in the continuous time replicator equation is the mean replacement time of the population and could therefore be identified as a generation. In the simulation each generation is devided into $`N`$ time-steps (where $`N`$ is the population size). At each of these time-steps a single individual is selected to copy and mutate. Individuals are selected wita h probability proportional to their relative fitness, which is given by $`\sigma `$ or $`1`$ on the single-peaked landscape. Thus a master-sequence of strain $`g_i`$ (with mutation rate $`\mu _i`$) is chosen with probability $`\frac{x_i\sigma }{f}`$. This copy replaces a randomly chosen individual in the existing population which is then discarded. Thus the population is replaced one by one in discrete birth-death events. In the limit of large population size the dynamics of this simulation approaches the continuous time replicator equation.

The fitness peak is changed every $`\tau `$ generations to one of its nearest neighbors. For the binary genomes used here it accomplished by flipping a randomly chosen bit in the definition of the fitness peak.

The population was first seeded with a diverse range of mutation rates and the population was allowed to evolve while these rates were kept fixed. This is a true test of $`\mu _{opt}`$, since the fastest growing sequence should come to dominate. In general the population converged to the strain with mutation rate closest to the theoretically predicted $`\mu _{opt}`$. Figure 1 shows the mean mutation rate of the population $`\overline{\mu }`$ evolving down towards the theoretically predicted optimum $`\mu _{opt}\frac{1}{\nu \sigma \tau }=0.00445`$. From about generation 800 the variance in mutation rates in the population is larger than the fluctuations in the mean and the evolution of rates has effectively ended.

Simulations were also made to study the effects of more rapidly changing mutator dynamics. In these simulations errors in the copying process not only introduce changes in the fitness determining genotype, but also result in offspring with slightly different mutation rates than their parents, i.e. the mutation rate is allowed to evolve. The mutation rate was treated as a continuous variable which had Gaussian noise introduced during the copying process.

Fig. 2 shows the evolution of mutation rates in detail in a population with a reasonably fast rate of change of mutation rates. This simulation has the same landscape parameters as Fig. 1. The mean mutation rate fluctuates around the optimum. For mutation rates close to the optimum fluctuations in selection are significantly larger than the selective advantages of one mutation rate over another. In this region the evolution of mutation rates is effectively neutral and thus the mean mutation rate conducts a random walk about the optimum. We also note that the population typically spends more time with mutation rates above the optimum than below. This is mainly a finite population size effect. When the peak moves and the population size is limited there is a relatively large probability that there will be no individuals representing a master-sequence with very low mutation rate on the new peak. This leads to a temporary increase in mutation rate in the population after an environmental shift.

## VI Biological Implications

In nature the existence, and value, of an optimum mutation rate that results from a changing environment depends on many different parameters: the time between shifts in the selective environment, the complex structure of the fitness-landscape, the genome length, co-evolutionary effects, the strength of selection, neutrality in the fitness landscape and fluctuations due to finite population sizes etc. One must therefore be careful when comparing the results of a simple model, such as the one we have presented in this paper, and biological measurements. Nonetheless it is this range of possible differences between organisms and the complexity of their evolutionary environments that leads us to consider the possibility that simple laws of biology — such as the scaling of point mutation rates with genome length — are likely to have quite simple explanations that do not depend on the details of the particular organism. It is therefore worth comparing the results of the model presented in this paper with the biological data.

For low mutation rates Eq. 9 is relatively insensitive to changes in the average fitness or size and frequency of environmental changes, as shownin Fig. 3. This insensitivity of the optimal genomic mutation rates to evolutionary parameters is important, since the bacteria and phages illustrated in table I are most unlikely to live in environments with the same types of time-dynamics and time-scales. In Fig. 3 we see that the sensitivity to one of the parameters in the model, $`\sigma `$ or $`\tau `$, depends strongly on in which region the other parameter is. For most realistiv populations we may expect the selective advantage $`\sigma `$ to be weak, maybe on average less than $`2`$. The predicted mutation rate is then highly insensitive to the average time between shifts in the fitness landscape, e.g. $`\sigma =2`$ gives $`\tau [110,200]`$ for the organisms listed in Table I. It is also reasonable to assume the fitness landscapes of the organisms listed in Table I to be more similar to each other than to higher eukaryotes and since our predictions as to $`Q_{opt}`$ are rather insensitive to the details of $`\sigma (t)`$, $`\tau `$ and $`h(\mu )`$ we would expect many organisms to have approximately the same mutation rate per genome (within an order of magnitude). This is what we observe for simple DNA-based organisms.

## VII RNA viruses

The lytic RNA viruses consistently show an extremely high mutation rate — orders of magnitude larger than that of any DNA viruses of similar size. This rate of around one substitution per genome per generation is inconsistent with the analysis conducted above for mutation rates evolving in a changing selective environment. Such high rates imply implausible values for the dynamic environment parameters.

As an explanation for the high mutation rates observed in many RNA viruses and the mutation rate scaling with genome length it has been suggested that these viruses have evolved the highest mutation rate possible to be able to adapt to a rapidly changing environment. The maximal mutation rate is then given by the error-threshold, which was first discussed in a model by Eigen et al. . It basically states that on a singled peaked fitness landscape an organism must have high enough copying fidelity so that its relative superiority in reproduction rate multiplied by the probability of reproducing onto a perfect copy of itself must be larger than one, otherwise there will be no effective selection for the genotype. It has later been shown that the error-threshold can rather easily be generalized to include effects of a dynamic environment . From this argument it is however not clear why RNA viruses should evolve towards the error-threshold while DNA based organism tend to have much lower mutation rates (by orders of magnitude). In this section we will combine the error-threshold with the model presented in this paper to suggest a possible explanation to the difference in observed mutation rates between DNA and RNA based organisms.

The dynamic environment model presented in this paper applies to organisms where the copying fidelity is encoded in a part of the genome that has little or no effect on fitness. In many viruses this may not be appropriate, partly because the proteins involved in mutagenesis may have a multitude of functions but also because the relatively high selective pressure towards short genome lengths will result in the overlap and multiple use of genetic material where possible. This give rise to a different possibility for the evolution of optimal mutation rates and might help explain the large differences between the observations for RNA and DNA based organisms.

We suggest that for organisms which have strong overlaps between genes coding for the mutation rate and genes coding more directly for reproductive advantage there is no effective selection for lower mutation rates, as long as the mutation rate is below the error threshold. This argument is based on the assumption that most mutations are deleterious in terms of fitness, and that the relative fitness advantage on the local peak results in stronger selection pressure than the pressure towards lower mutation rates. We also assume that evolution of mutation rates usually affect regions of the genome where the organism need mutations to be able to adapt ot changes in the environment. If these assumptions apply we expect a population to have mutation rates close to the error-threshold. Changes to mutation rate is transient, assuming that the organism is not pushed beyond the error-threshold.

For this hypotheses to apply, viruses with high mutation rate (mainly RNA viruses) should have overlapping genes regulating mutation frequency as well as reproduction rate, whereas organisms with low mutation rates (such as those listed in Table I) should not have overlapping reading frames in their genomes. There are observations that support this, but it is unclear whether the correlation is strong enough for this hypothesis to be valid.

## VIII Conclusions

In this paper we have studied the evolution of mutation rates in a population of multi locus genomes. The genomic mutation rate $`\mu _G`$ leading to the greatest long term growth of a strain (the optimal rate) was analytically determined for reasonably general peak shifts and time-dependent replication rates $`\sigma (t)`$

$$\mu _G\frac{k}{\sigma \tau }$$

where $`k`$ is the mean Hamming distance between successive fitness optima and $`\tau `$ is the mean time between shifts. These optimal rates were quantitatively confirmed by computational simulations of populations whose mutation rates were allowed to evolve.

These continuous time multi-locus replicator models predict the kind of scaling of point-mutation rate with genome length that has been observed in some bacteria and viruses/phages and puzzled over for years. When combined with the consequences of the multiple use/pleiotropic encoding of copying machinery these models of the evolution of mutation rate in dynamic environments also suggest why lytic RNA viruses may have rates at or about the error-threshold.

We would like to thank Claes Andersson and Erik van Nimwegen for useful discussions. Thanks are also due to Mats Nordahl who has given valuable comments on the manuscript. Nigel Snoad and Martin Nilsson were supported by SFI core funding grants. N.S. would also like to acknowledge the support of Marc Feldman and the Center for Computational Genetics and Biological Modeling at Stanford University while preparing this manuscript.


text:
# On the criterion for Bose-Einstein condensation for particles in traps

\[

## Abstract

We consider the criterion for Bose condensation for particles in a harmonic trap. For a fixed angular momentum, the lowest energy state for a cloud of bosons with attractive interactions is the ground state of the cloud with all the angular momentum in the center-of-mass motion, and the one-particle reduced density matrix generally does not have a single large eigenvalue, but a number of them, suggesting that the state is an example of a fragmented condensate (Wilkin, Gunn, and Smith, Phys. Rev. Lett. 80, 2265 (1998)). We show that a convenient way to describe correlations in the system is by defining an internal one-particle reduced density matrix, in which the center-of-mass motion is eliminated, and that this has a single eigenvalue equal to the number of particles for the problem considered here. Our considerations indicate that care is necessary in formulating a criterion for Bose-Einstein condensation.

PACS numbers: 03.75.Fi, 05.30.Jp, 67.40.Db, 67.40.Vs

\]

In investigations of the properties of the helium liquids, rotation has proved to be a valuable probe, and this has led to the discovery of vortex lines in liquid 4He and in the superfluid phases of liquid 3He. In the discussion of superfluidity and Bose-Einstein condensation it is usually assumed that particles are confined in a box, with a uniform potential throughout its interior, and with infinitely high walls. This situation is very different from that encountered in studies of dilute atomic vapors in traps, where the potential is generally closer to a harmonic oscillator one. We therefore consider the properties of rotating atomic clouds of bosons at zero temperature. Among possible ways of putting angular momentum into the system are the creation of vortex lines  and excitation of surface modes and other internal excitations . For traps with an axis of symmetry, another possibility is to excite the center-of-mass motion, which for modes having rotational quantum number $`m_z=\pm 1`$ about the symmetry axis is equivalent to excitation of the simplest surface mode. The purpose of this paper is to consider the criterion for Bose condensation in a cloud of particles in a trap.

A criterion for Bose-Einstein condensation was proposed by Penrose , and subsequently elaborated by Penrose and Onsager , and by Yang . Consider the one-particle reduced density matrix. This is defined in terms of the many-body wave function $`\psi (𝐫_1,\mathrm{}𝐫_N)`$, where the $`𝐫_i`$ are the coordinates of the particles, by the equation

$$\rho ^{(1)}(𝐫,𝐫^{})=𝑑𝐫_2\mathrm{}𝑑𝐫_N\psi (𝐫,𝐫_2\mathrm{}𝐫_N)\psi ^{}(𝐫^{},𝐫_2\mathrm{}𝐫_N).$$
(1)
The density matrix may be expanded in terms of its eigenfunctions $`\chi _j(𝐫)`$ with eigenvalues $`\lambda _j`$:

$$\rho ^{(1)}(𝐫,𝐫^{})=\underset{j}{}\lambda _j\chi _j^{}(𝐫^{})\chi _j(𝐫).$$
(2)
The condition generally adopted for the existence of Bose condensation is that there should be one eigenvalue that is of order the number of particles, in the limit when the number of particles tends to infinity. However, as we shall indicate below, this can be misleading.

In an instructive paper Wilkin et al.  considered bosons in a harmonic trap, and showed that for weak attractive forces the lowest state with a given angular momentum is one in which all the angular momentum is associated with the center-of-mass motion. This corresponds to excitation of a dipolar (magnetic quantum number $`m_z=\pm 1`$) surface mode of oscillation . An analysis of the eigenvalues of the one-particle density matrix revealed that there is not a single large eigenvalue, but rather there can be many eigenvalues of a comparable size. This led the authors to suggest that the state of the rotating system corresponds to that of a fragmented condensate of the type proposed by Nozières and Saint James . However, this conclusion seems surprising, because the state is just the ground state of the cloud, which is Bose-Einstein condensed, executing center-of-mass motion. The problem is to identify a quantity that can reveal the Bose condensation.

Let us consider N spinless bosons in a harmonic trap. This system has the important property that the center-of-mass and relative motions of the particles are separable, and therefore the wave function may be written in the form

$$\psi (𝐫_1,\mathrm{}𝐫_N)=\psi _{\mathrm{cm}}(𝐑)\psi _{\mathrm{rel}}(𝐪_1,𝐪_2,\mathrm{},𝐪_N),$$
(3)
where the center-of-mass coordinate is given by $`𝐑=(1/N)\mathrm{\Sigma }_{i=1}^N𝐫_i`$ and $`𝐪_i=𝐫_i𝐑`$ is the coordinate of particle $`i`$ relative to the center of mass. Note that the definition of the center-of-mass coordinate implies that only $`N1`$ of the $`𝐪_i`$ are independent. Let us now introduce a quantity which we call the internal density matrix

$`\rho _{\mathrm{int}}(𝐪_i,𝐪_i^{})`$ $`=`$ $`{\displaystyle 𝑑𝐑\psi (\{𝐪_i+𝐑\})\psi ^{}(\{𝐪_i^{}+𝐑\})}`$ (4)
$`=`$ $`\psi _{\mathrm{rel}}(\{𝐪_i\})\psi _{\mathrm{rel}}^{}(\{𝐪_i^{}\}).`$ (5)
This has the property that the center-of-mass motion has been eliminated, and consequently it depends only on coordinates relative to the center of mass. One may also construct reduced density matrices in terms of the relative coordinates in the standard way. For example, the internal one-particle reduced density matrix is

$$\rho _{\mathrm{int}}^{(1)}(𝐪_1,𝐪_1^{})=𝑑𝐪_2\mathrm{}𝑑𝐪_{N1}\rho _{\mathrm{int}}(𝐪_i,𝐪_{}^{}{}_{i}{}^{}).$$
(6)
Consequently the internal density matrix for the lowest energy state of a cloud of bosons with weak attractive interactions and nonzero angular momentum in a harmonic trap is the same as for the ground state in the absence of rotation. It has one large eigenvalue which is equal to $`N`$, to within terms of unity, the deviation from $`N`$ being due to the constraint that the sum of the coordinates relative to the center of mass must vanish.

Our discussion indicates that caution is required in formulating a criterion for Bose-Einstein condensation in finite systems, and that the one generally used can be misleading. For systems in harmonic traps the internal density matrix introduced above can be a useful concept, because it removes the effects of the center-of-mass motion which, for a rotating cloud, tends to spread out the weight of the usual one-particle reduced density matrix over a number of states.

We are grateful to Ben Mottelson for helpful conversations.


text:
# REFERENCES

\[ Comment on “The Transition Temperature of the Dilute Interacting Bose Gas” and on “Transition Temperature of a Uniform Bose Gas”
\] The transition temperature of a uniform, weakly interacting Bose gas has received considerable interest in the past months, with at least three contributions appearing in these pages alone. In  Baym et al. argued that because of infrared divergences, the shift in temperature due to the weak repulsive interaction relative to a free Bose gas cannot be obtained in perturbation theory, no matter how weak the interaction. To regularize the infrared divergences, they derived a self-consistent equation for the energy spectrum of the elementary excitations and estimated the temperature shift as

$$\frac{T_\mathrm{c}T_0}{T_0}=c_0\left(na^3\right)^{1/3},$$
(1)
with $`c_02.9`$. Here, $`T_0=(2\pi \mathrm{}^2/mk_\mathrm{B})[n/\zeta (3/2)]^{2/3}`$, with $`\zeta `$ the Riemann zeta function, is the transition temperature of an ideal Bose gas with particle number density $`n`$, and $`a`$ is the s-wave scattering length characterizing the repulsive interaction between the particles.

In  Huang addressed the same issue using a virial expansion of the equation of state. He arrived at the conclusion that the shift in $`T_\mathrm{c}`$ increases with the square root of the scattering length, rather than linearly with $`a`$.

In this Comment we wish to point out that first, contrary to the claim by Baym et al. , the shift in temperature can be calculated perturbatively, as was done in . And second, that the regime considered by Huang  does not describe the Bose-Einstein condensed phase.

Specifically, he confines himself to the regime where the fugacity $`z=\mathrm{exp}(\mu /k_\mathrm{B}T)`$ is smaller than unity, corresponding to negative values of the chemical potential $`\mu `$. This is because the used expression for the pressure as function of $`z`$ becomes complex for $`z>1`$. However, the broken-symmetry phase of a weakly interacting Bose gas is, unlike He-II, characterized by a positive value of $`\mu `$. In particular, at $`T=T_\mathrm{c}`$ perturbation theory gives  $`\mu =8\pi \mathrm{}^2an/m`$ to linear order in $`a`$. Note that in a free Bose gas $`\mu `$ is either negative (in the gas phase) or zero (in the condensed phase).

In  the shift in temperature, among other things, was obtained perturbatively by calculating the pressure to the one-loop order in a grand canonical ensemble as in , and expanding it in a high-temperature series. The justification of the high-temperature expansion was given a posteriori by the observation that the leading term in the expression for the transition temperature is of order $`(\mu /a)^{2/3}`$, where we recall that $`\mu `$ is an independent variable in a grand canonical ensemble. Since this is large for $`a`$ small, the high-temperature expansion is consistent with the weak-coupling assumption of perturbation theory.

The emerging infrared divergences were regularized by analytic continuation, frequently used to regularize ultraviolet divergences. Finally, using the equation of state to swap the chemical potential for the particle number density as independent variable—which is more appropriate from an experimental view point—we found that the shift in temperature is given by Eq. (1), with

$$c_0=\frac{8}{3}\frac{\zeta (\frac{1}{2})}{\zeta ^{1/3}(\frac{3}{2})}2.83.$$
(2)
The estimate of Baym et al. is remarkably close to this perturbative result.

Earlier path-integral Monte Carlo simulations  restricted to small (up to 216) particle numbers gave a value $`c_0=0.34\pm 0.06`$ almost an order of magnitude smaller that the perturbative result (2). More recent simulations  on larger systems containing of the order of $`10^4`$ particles, gave a value $`c_0=2.3\pm 0.25`$ in reasonable agreement with Eq. (2).

In closing, we remark that in  also a similar high-temperature expansion and regularization of infrared divergences as discussed here was applied to the BCS theory of superconductivity, and shown to reproduce, among other things, the known result for the BCS transition temperature.

I thank G. Baym, K. Huang and F. Laloë for useful correspondence, and M. Krusius for helpful advice. This work was funded in part by the EU sponsored programme Transfer and Mobility of Researchers under contract No. ERBFMGECT980122.

Adriaan M. J. Schakel
Low Temperature Laboratory
Helsinki University of Technology
P.O. Box 2200, FIN-02015 HUT, Finland
and
Institut für Theoretische Physik
Freie Universität Berlin
Arnimallee 14, 14195 Berlin, Germany
e-mail: schakel@physik.fu-berlin.de

PACS numbers: 03.75.Fi, 05.30.Jp