text:
# Soft Listeria : actin-based propulsion of liquid drops

## Abstract

We study the motion of oil drops propelled by actin polymerization in cell extracts. Drops deform and acquire a pear-like shape under the action of the elastic stresses exerted by the actin comet. We solve this free boundary problem and calculate the drop shape taking into account the elasticity of the actin gel and the variation of the polymerization velocity with normal stress. The pressure balance on the liquid drop imposes a zero propulsive force if gradients in surface tension or internal pressure are not taken into account. Quantitative parameters of actin polymerization are obtained by fitting theory to experiment.

Actin polymerization is a key element in the motility of most cells and bacteria. The bacteria Listeria monocytogenes are propelled inside cells by the growth of a soft elastic comet made of a filamentous actin network. Actin polymerizes at the back of the advancing bacterium. The biochemistry of the comet formation is now well-understood pollard . Theoretical approaches have been proposed to explain the physical mechanism of force production. They differ by the scale at which they describe the mechanism. Molecular models by Mogilner et al Mogilner  consider the Brownian flexibility of growing actin filaments, whereas Carlsson Carlsson  concentrates on the effect of branching and growth of the actin network. Gerbal et al Gerbal  analyze at a mesoscopic scale the elastic stresses exerted by the deformed comet gel on the bacterium resulting in the propulsive force prost .

A further step in understanding the propulsion mechanism is provided by the study of biomimetic experimental systems where Listeria are replaced by solid spherical beads on which actin polymerization promoters theriotpnas ; anne ; Vincent ; welch  are attached. These beads mimic closely the natural propulsion mechanism of Listeria with comet tail formation, after the breaking of the initial spherical symmetry theriotncb ; anne . Biomimetic systems allow for a systematic variation of the parameters and thus for a quantitative comparison to theory.

The aim of this letter is to demonstrate both theoretically and experimentally the mechanism of force production due to the elastic stresses exerted by the actin comet gel. A newly designed experimental system is made of oil drops with actin polymerization promoters attached on their surface. Once placed in cell extracts, such an oil drop moves by actin polymerization and deforms under the action of the elastic stresses exerted by the gel. The same squeezing effect is observed on endosomes Tauton  driven by actin comets and synthetic vesicles covered with the bacterial protein ActA VanOudenaarden ; giardini . However, for liquid drops, the knowledge of the surface tension and the constant volume allow for a quantitative analysis of the observed shape and thus the determination of the elastic stress distribution on the drop surface.

An emulsion of oil drops is obtained by sonicating ($`5`$ seconds, $`90`$ W) a mixture of edible oil, Isio4, and buffer ($`2.3\%`$ oil in borate buffer $`100`$mM, pH $`8.5`$). The actin polymerization promoter VCA is derived from the Wiskott-Aldrich syndrom protein (WASP) and purified as described in Fradelizi et al Fradelizi . VCA is adsorbed onto the oil drops by incubating 10 $`\mu `$l of the emulsion with a $`0.2`$mg/ml VCA solution in borate buffer. A volume of $`0.2`$ $`\mu `$l of the emulsion coated with VCA is added to $`15`$ $`\mu `$l of HeLa extracts prepared as explained in Noireaux et al Vincent  with a final protein concentration of $`19`$mg/ml. As usual the extracts are supplemented with 0.019mg/ml G-actin, 0.06mg/ml rhodamin actin, 1mM ATP, 27mM creatine phosphate, 1mM DTT.

The sample is observed by bright field or fluorescence microscopy using an Olympus BX51(Germany)and Metamorph software (Princeton Instruments, USA). The drop radii range from $`1.5`$ to $`5.5`$ $`\mu `$m. As with Listeria, actin polymerizes only at the interface between the drop and the comet and depolymerizes at the back of the comet. When placed in extracts, actin is first polymerized on the drops with a spherical symmetry. After roughly one hour, symmetry is broken for c.a. $`70\%`$ of the drops. They develop an actin comet, deform into a pear-like shape and move. Smaller drops are less deformed than larger drops. An example of a moving drop is shown on Fig.1.

At the beginning of the experiment, the actin polymerization factor VCA is uniformly distributed around the spherical drop. After deformation of the drop and formation of the comet, fluorescence intensity measurements using FITC-labeled VCA fitc  show that $`90\%`$ of the VCA is found on the comet side of the drop, as seen on Fig.2.

This means either that VCA has been displaced from the interface to the bulk of the extract except where the gel is present, or that all VCA has been collected by the gel during the symmetry breaking process. In any event this also means that the VCA surface density is comparable to the density of filament extremities at the surface. The average distance $`d`$ between filament extremities is larger than the close-packing distance. A lower bound of $`d`$ is $`10`$nm. The surface tension change due to the presence of VCA at the interface is then of order $`kT/d^2<4\text{x}10^2`$mN/m, more than $`100`$ times smaller than the oil-extract tension that we measure by the pendant drop method $`\gamma _0=4\pm 0.6`$ mN/m.

The experimental shape of the drop is sketched on Fig.3.

The front part is a spherical cap of radius $`R`$ not covered by actin. The radius $`R`$ is different from the radius of the undeformed spherical drop $`R_0`$ that fixes the volume $`4\pi R_0^3/3`$. The back part of the drop, surrounded by the comet, has a blunted cone-like shape (rotationally symmetric around the direction of motion). In order to give a theoretical description of the drop shape, we parameterize it by the liquid thickness $`h`$ and the local angle $`\theta `$ between the tangent to the shape and the direction of motion. The spherical cap and the cone match at the triple line between the drop, the comet and the surrounding solvent; we call the tangent angle at the contact point $`\theta _m`$ where the corresponding oil thickness is $`h_m=R\mathrm{cos}\theta _m`$. The pressure inside the drop varies from point to point since the drop motion induces an internal flow. At the interface with the extract, in the spherical part, it is given by Laplace’s law $`P_{in}=2\gamma _0/R+P_0`$, where $`P_0`$ is the pressure in the surrounding liquid. At any point inside the drop, the pressure differs from $`P_{in}`$ and reads $`P=P_{in}+\delta P`$. At a point of thickness $`h`$ the local stress balance along the normal of the drop is given by the local Laplace’s law

$$\frac{2\gamma _0}{R}+\delta P(h)=\gamma \left(\frac{\mathrm{cos}\theta (h)}{h}+\frac{d\mathrm{cos}\theta (h)}{dh}\right)\sigma _{nn}(h)$$
(1)
where $`\sigma _{nn}`$ is the normal stress at the surface of the comet, exerted by the drop ($`\sigma _{nn}`$ is positive, dilative stress, if the comet pulls on the drop and negative, compressive stress, if the comet pushes the drop). Although we argued that the surface tension gradients are small, we consider here that the surface tension $`\gamma =\gamma _0+\delta \gamma (h)`$ varies along the interface; it is constant in the spherical part with a value $`\gamma _0`$ and it is continuous at the contact line. The total elastic force $`F_e`$ exerted by the comet on the drop is obtained by integrating the projection of the normal elastic stress on the direction of motion.

$$F_e/2\pi =_0^{h_m}𝑑hh\left(\mathrm{cos}\theta \frac{d\delta \gamma }{dh}+\delta P\right)$$
(2)
If we ignore both the pressure gradient inside the drop and the surface tension gradient, the propulsive force equals zero. This result holds for any axisymmetric drop shape whatever the elastic stress distribution along the surface, independent of any model for the comet elasticity and the actin polymerization velocity. This hypothesis, used in references VanOudenaarden ; giardini , is unable to produce an estimate of the propulsive force. The important consequence of Eq.2 is that an experimental measurement of the propulsive force must take into account the surface tension gradient and the flow inside the drop as discussed later.

We now proceed with the determination of the drop shape and the stress distribution. Given that the surface tension variation is small along the drop contour and neglecting hydrodynamic effects, we consider, for local equations, that the surface tension and the internal pressure are constant. The elastic stresses in the gel influence the polymerization kinetics; polymerization is normal to the surface of the drop and it is accelerated by a dilative stress and slowed down by a compressive stress. Classical rate theories kramers  predict a polymerization velocity $`v_p`$ varying as a Boltzmann law

$$v_p(h)=v_p^0\mathrm{exp}[\sigma _{nn}(h)/\sigma _0],\sigma _0kT/a^2\delta $$
(3)
where $`a`$ is the distance between actin polymerization promoters on the drop surface, $`\delta `$ is of the order of the size of an actin monomer and $`v_p^0`$ is the polymerization velocity in the absence of stress.

The last equation determining the shape of the drop is the conservation of the gel volume upon polymerization. In a first approximation, we assume both that the gel density is constant and that the comet is a perfect cylinder. In a steady state, the drop advances at a constant velocity $`V`$. The local gel thickness $`e`$ shown on Fig.3 is then such that $`de/ds=\mathrm{tan}\theta `$ where $`s`$ is the length along the drop contour. With these simplifying approximations, the local polymerization velocity is related to the advancing velocity by $`v_p(h)=V\mathrm{sin}\theta (h)`$.

This approximation does not allow the determination of the drop velocity $`V`$ since the propulsive force vanishes. We thus find a family of solutions for the drop shape parametrized by the advancing velocity. We determine the drop shape using the measured advancing velocity $`V=0.15\pm 0.03\mu m`$/min.

The shape of the comet-drop interface departs from a pure cone when the Boltzmann factor is significantly larger than 1. This defines the size of the blunted region as $`\mathrm{}\gamma a^2\delta /kT`$ which leads to $`\sigma _0=\gamma /\mathrm{}`$. In the following we consider $`ϵ\mathrm{}/R`$ as a small number.

For $`h`$ of order $`R`$, the elastic stress $`\sigma _{nn}(h)`$ is small and can be neglected. Then Eq.3 describes a perfect cone with angle $`\theta =\theta _0`$ such that $`V\mathrm{sin}\theta _0=v_p^0`$. This result is independent of the polymerization law giving $`v_p`$ as a function of $`\sigma _{nn}(h)`$. The measure of $`\theta _0`$ gives thus access to the polymerization velocity in the absence of elastic stress.

In the blunted region, $`h`$ is of order $`\mathrm{}`$ and we neglect the Laplace pressure term on the left hand side of equation 1. In the vicinity of the rear point, the drop profile is given by $`h^2=4\mathrm{}z/\mathrm{log}(1/\mathrm{sin}\theta _0)`$.

The normal stress in the rear region of size $`\mathrm{}`$ is positive and of order $`\gamma /\mathrm{}`$. At the rear point, $`\sigma _{nn}(h=0)=\gamma \mathrm{log}(1/\mathrm{sin}\theta _0)/\mathrm{}`$. In the conical region, the stress reads at lowest order in $`ϵ`$

$$\sigma _{nn}(h)=(\gamma /R)\left(2+R\mathrm{cos}\theta _0/h\right)$$
(4)
It is positive at the back of the drop (pulling the drop backwards) and negative in the front part of the cone (pushing the drop forwards) as qualitatively observed in reference VanOudenaarden . It vanishes for a thickness $`h=(R\mathrm{cos}\theta _0)/2=h_m/2`$. This is in accordance with the prediction made in ref. Gerbal  that the actin gel could pull at the rear of Listeria and explains the pear-like deformation that we now describe quantitatively.

A more detailed description of the drop profile is obtained by solving numerically equations 1 and 3. We use the radius $`R`$ as a unit length; the drop profile depends on the two dimensionless parameters, $`ϵ`$ and $`\mathrm{sin}\theta _0=V/v_p^0`$. On Fig.4, we show a comparison of the calculated and experimental profiles for the drop of Fig.1. The experimental profile has been digitized and adjusted by a continuous curve. The best fitting parameters are $`ϵ=0.049`$ and $`\mathrm{sin}\theta _0=0.58`$ ($`\theta _0=35.6^{}`$); this gives a determination of the length $`\mathrm{}=0.125\mu m`$, of the stress $`\sigma _0=32nN/\mu m^2`$ and of the polymerization speed in the absence of stress $`v_p^0=1.4`$nm/s footnote . The local normal stress is obtained from the experimental drop shape by using equation  with a constant interfacial tension $`\gamma _0`$ and neglecting the variation of the internal pressure $`\delta P`$. The length $`\mathrm{}`$ being constant, smaller drops corresponding to larger $`ϵ`$ are less deformed, in agreement with experiments.

In the description proposed so far, the curvature of the interface between the comet and the drop has a discontinuity at the triple line. In the spherical region, the curvature of the interface is $`2/R`$; in the conical region, the interface is curved only in one direction and the curvature is $`1/R`$. The local forces normal to the interface are still balanced and the elastic normal stress in the comet at the triple line is $`\sigma _{nn}(h_m)=\gamma /R`$. (see insert of figure 4). At the triple line, the gel thickness $`e`$ vanishes and cannot create a finite stress. This explains why on figure 4, the theoretical and experimental normal stresses are in disagreement in the very vicinity of the triple line. We now assume that the gel density remains constant but that the comet shape is not a perfect cylinder in the region close to the triple line. Volume conservation then imposes a polymerization velocity $`v_p=V\mathrm{cos}\theta de/ds`$ where $`s`$ is the length along the interface. When a gel element is created in a time $`dt`$, it is stretched by an amount $`\delta u=\delta t(V\mathrm{sin}\theta +v_p)`$; the tensile stresses $`\sigma _{ii}`$ in the azimuthal direction ($`i=\varphi `$) and in the tangential direction along the interface ($`i=t`$) are then increased. At the level of scaling laws, we write this increase as $`\frac{d\sigma _{ii}}{ds}=\frac{E}{R_iV\mathrm{cos}\theta }(V\mathrm{sin}\theta v_p)`$ where $`R_i`$ is the radius of curvature of the interface in the direction $`i`$ and E an elastic modulus of the comet. In the vicinity of the triple line, we use the thin shell approximation and relate the tensile and normal stresses in the gel $`\sigma _{nn}=e(\frac{\sigma _{tt}}{R_t}+\frac{\sigma _{\varphi \varphi }}{R_\varphi })`$. In this boundary layer, the two tensile stresses can be considered as constant; the matching to the pear-like shape imposes that the azimuthal tensile stress $`\sigma _{\varphi \varphi }`$ vanishes. Defining the dimensionless tensile stress $`\stackrel{~}{\sigma }_t=\mathrm{}\sigma _{tt}/\gamma `$, the normal stress in the boundary layer is then calculated as

$$\sigma _{nn}=\frac{\gamma }{R}\frac{e\stackrel{~}{\sigma }_t}{\mathrm{}+e\stackrel{~}{\sigma }_t}$$
(5)
As the thickness of the gel vanishes, the normal stress vanishes as expected. The boundary layer, where the comet is deformed, has a thickness $`\mathrm{}/\stackrel{~}{\sigma }_t`$ of order $`\mathrm{}`$ and is thus small in the limit where $`ϵ`$ is small. When $`e`$ is large, further away from the triple line, the comet reaches a cylindrical shape and the normal stress is $`\gamma /R`$; one can consider the comet as a perfect cylinder as done above.

Our experimental observations on the shape of liquid drops propelled by actin polymerization are well described by the theoretical model based on a local normal force balance and on a Boltzmann variation of the polymerization velocity with normal stress. The results are robust if we use, for the polymerization velocity, the mathematical forms suggested by simulations on flat surfaces Carlsson . We demonstrate that the elastic propulsive force cannot be calculated from an elastic stress distribution that ignores both pressure variations inside the drop and surface tension gradients. We have estimated the surface tension gradient by assuming that the actin polymerization promoter density profile along the interface follows the gel elastic deformation. This leads to a propulsive force of order $`F_e2\pi \gamma _0ϵ_G/E`$ where $`ϵ_G`$ is the Gibbs elastic modulus of the interface. With reasonable values of the parameters, we find a propulsive force of order $`100`$pN. This is however a lower bound since there could exist other contributions to the surface tension gradient. The advancing velocity of the drop results from a balance between the propulsive force and the friction force between the comet and the drop but a precise study of the velocity selection will require a more refined analysis.

In our experiments, as well as in reference VanOudenaarden , the motion stops when the drop (or the liposome) becomes spherical. The large stress at the back of the drop could lead to “cavitation” or rupture of the links between the drop and the comet. The elastic stresses exerted on the drop then relax and provoke the experimentally observed arrest. We expect the stress distribution to be similar for other types of actin propelled objects such as the bacteria Listeria or solid beads. This would explain the observation of hollow comets of Listeria Tauton . A final output of our analysis is an estimate of the polymerization velocity of an actin gel in the absence of stress and its variation with stress.

###### Acknowledgements.

J.F.J. is grateful to A. van Oudenaarden for sending a copy of reference VanOudenaarden . O.C. thanks the European Network PHYNECS and Ministerio de Educación, Cultura y Deporte for financial support. We thank J.Casademunt and K.Zeldovitch (Institut Curie) for help in the numerical work and V.Noireaux for initiating the experiments.


format: text

text:
# Constraining the neutrino magnetic moment with anti-neutrinos from the Sun

On leave from \]Institute of Terrestrial Magnetism, Ionosphere and Radio Wave Propagation of the Russian Academy of Sciences (IZMIRAN), 142190, Troitsk, Moscow, Russia

On leave from \]IZMIRAN

## Abstract

We discuss the impact of different solar neutrino data on the spin-flavor-precession (SFP) mechanism of neutrino conversion. We find that, although detailed solar rates and spectra allow the SFP solution as a sub-leading effect, the recent KamLAND constraint on the solar antineutrino flux places stronger constraints to this mechanism. Moreover, we show that for the case of random magnetic fields inside the Sun, one obtains a more stringent constraint on the neutrino magnetic moment down to the level of $`\mu _\nu <\mathrm{few}\times 10^{12}\mu _B`$, similar to bounds obtained from star cooling.

preprint: IFIC/03-49

The latest KamLAND result :2003gg  greatly improves the expected sensitivity limit on a possible antineutrino component in the solar flux from 0.1 % Kamland:proposal  of the solar boron $`\nu _e`$ flux to $`2.8\times 10^2`$ % at the 90 % C.L., about 30 times better than the recent Super-K limit Gando:2002ub .

The presence of anti-neutrinos in the solar flux may signal the existence of SFP conversions induced by non-vanishing neutrino transition magnetic moments Schechter:1981hw ; Akhmedov:uk  or, alternatively, neutrino decays in models with spontaneous violation of lepton number Schechter:1981cv . Here we focus on the more likely case of anti-neutrinos produced by SFP conversions. Several solar magnetic field models are possible, characterized by different assumptions pertaining to their magnitude, location and typical scales Kutvitskii ; others ; Friedland .

Previous to the latest KamLAND limit on the solar anti-neutrino flux, the first KamLAND evidence of reactor anti-neutrino disappearance Eguchi:2002dm  had already excluded SFP scenarios as solutions to the solar neutrino problem Barranco:2002te . However, even with the recent SNO salt results Ahmed:2003kj  which confirm the simplest three-neutrino oscillation picture Maltoni:2003da , a neutrino magnetic moment could still play a role as a sub-leading effect. In order to illustrate this, we have performed a $`\chi ^2`$ analysis taking into account all the solar experimental data plus the KamLAND disappearance data for the same self-consistent profile for the magnetic field in the convective zone used in previous analyses Miranda:2000bi . The results are shown in fig. 1, were we have taken different values for the product $`\mu _\nu B`$ (neutrino magnetic moment in units of $`10^{11}\mu _B`$ and maximum magnetic field value in convective zone in units of MG) and plotted the allowed regions at 90 %, 95 % and 99 % C. L.

One can see how the allowed LMA-MSW region of oscillation parameters (upper left panel) is modified in the presence of SFP conversions (upper right panels) in Fig. 1. The lower panels in this figure give the $`\mathrm{\Delta }\chi ^2`$ profiles with respect to $`\mathrm{\Delta }m_{\text{sol}}^2`$ and $`\mathrm{sin}^2\theta _{\text{sol}}`$, from where one can determine the corresponding allowed ranges. Given the current laboratory best limit on the neutrino magnetic moment $`\mu _{\overline{\nu }_e}<1.0\times 10^{10}\mu _B`$ Daraktchieva:2003dr ; Grimus:2002vb , and the fact that magnetic field amplitude inside the convective zone should not exceed 0.1-0.3 MG Moreno  one can see that, although limited, there is room left for sub-dominant SFP conversions.

We now turn to the recent KamLAND result on the $`\overline{\nu _e}`$ ’s from the Sun. The Collaboration has reported that the antineutrino flux $`\mathrm{\Phi }_{\overline{\nu _e}}`$ is less than $`3.7\times 10^2`$ cm -2 s -1 at 90 % C. L. which corresponds to a 0.028 % of the solar 8 B $`\nu _e`$ flux (in the energy window between 8.3 MeV and 14.8 MeV). It is this last result that will substantially limit the possibility of subleading SFP component in the neutrino conversion mechanism, establishing the robustness of the oscillation hypothesis.

Within our generalized picture (LMA-MSW+SFP), after the MSW $`\nu _e\nu _\mu `$ conversion takes place in the central region of the Sun, $`\overline{\nu _e}`$’s are produced due to the magnetic moment conversion $`\nu _\mu \overline{\nu _e}`$. To analyse quantitatively the restrictions imposed by this result we consider three models for the solar magnetic field:

1. regular magnetic fields, both in the convective (CZ) Kutvitskii  and radiative (RZ) Friedland  zones of the Sun
2. convective zone random magnetic fields Bykov:1998gv

For definiteness we consider the simplest aproximate two-neutrino picture, which is justified in view of the stringent limits that follow mainly from reactor neutrino experiments Maltoni:2003da . From the resulting $`4\times 4`$ form of the neutrino evolution equation Schechter:1981hw ; Akhmedov:uk , we compute the expected $`\overline{\nu _e}`$ yield for each magnetic field model. For the case of CZ magnetic fields, this $`4\times 4`$ form of the neutrino evolution further decouples into LMA-MSW conversions deep in the Sun followed by (approximate) vacuum SFP conversions Schechter:1981hw  inside the CZ. In contrast, for the RZ case there is no decoupling of the neutrino evolution, due to the large strength of the magnetic field 111Although possible, this model is not favored on physical grounds..

In order to obtain a conservative bound on $`\mu _\nu B`$ we have calculated the minimum $`\overline{\nu _e}`$ yield as the oscillation parameters vary within the acceptable range (at the 90 % C.L.) in the absence of magnetic field (pure LMA-MSW case). Such minimum predicted anti-neutrino fluxes for the case of regular field profiles are shown in the curves depicted in Fig. 2, while the recent KamLAND limit is indicated by the lower horizontal lines. For comparison, we present the previous Super-K bound, indicated by the upper horizontal lines. Note that in both cases of regular field profile we can only constrain the product $`\mu _\nu B`$, as opposed to the intrinsic neutrino magnetic moment 222For constraints on the intrinsic neutrino magnetic moment see Ref. Grimus:2002vb . . Grimus:2002vb .

A more interesting picture for SFP scenario is obtained if we consider the case of turbulent random magnetic fields inside the Sun. As we will see, the reason for this is twofold: (i) we will be able to fix, to some extent, the dependence on the magnetic field, in contrast to previous analysis of the random magnetic field presented in Ref. Bykov:1998gv ; Torrente-Lujan:1998sy , where an extra parameter $`L_0`$ appeared, characterizing the scale of the random magnetic field cells, and (ii) this model leads to more stringent limits on the neutrino magnetic moment by itself.

In the following we give a brief discussion of our new approach to the problem (a more detailed analysis will be published elsewhere Inprep ). The mean magnetic field value over the solar disc is of the order of 1 G and magnetic field strength in the solar spots reaches 1 kG. It is commonly accepted that fields measured at the solar surface are substantially weaker than those near the bottom of the convective zone (CZ) where they are supposed to be generated by a dynamo mechanism Yoshimura ; Sokoloff . In dynamo theory, the mean magnetic field is accompanied by a small-scale random magnetic field, which is not directly traced by sunspots or any other tracers of solar activity. By small scales we denote a typical scale of about the solar granule size ($``$ 1000 km) Stix  where the rms magnetic field amplitudes (generated either as a by-product of a large-scale dynamo mechanism or directly by a small-scale dynamo) in the range of 50-100 kG are reasonable.

For LMA oscillations the MSW conversion occurs well below the CZ resulting in a coherent mixture of two neutrino flavours at the bottom of the convective zone. When neutrinos cross the randomly fluctuating magnetic field of the CZ it is expected that $`\nu _{eL}`$ will convert to $`\overline{\nu }_{\mu R}`$ and $`\nu _{\mu L}`$ will convert to $`\overline{\nu }_{eR}`$ as a result of the CZ magnetic field. For the case of random CZ magnetic fields these populations will tend to equilibrate. However, given the laboratory upper limits on $`\mu _\nu `$ and the finite depth of the CZ, this equilibration is not achieved. The relaxation can be viewed as a small-amplitude random walk of the neutrino polarization vector leading to a small appearance of electron and muon antineutrinos at the solar surface (cf. Loeb:1989dr  for the vacuum case, and also discussion in Bykov:1998gv ). In what follows we describe the main features of this calculation and give the main result.

Because of the randomness of the underlying magnetic fields, the spin flavour evolution looses coherence, i.e., instead of evolving neutrino wave functions over the whole CZ and then obtaining the final probabilities, one has to compute the probabilities in each correlation cell and afterwards add them. As a result, independently of the random magnetic field model, the appearance of antineutrinos is proportional to the relevant Fourier harmonic of the transverse two-point magnetic correlation function, with space period equal to the LMA-MSW neutrino oscillation length (cf. Loeb:1989dr  and discussion in Bamert:1997jj  for the MSW matter noise).

We will assume for definiteness that random magnetic field evolution is due to the highly developed steady-state MHD turbulence treated within the Kolmogorov scaling theory Kolmogorov ; Landau . The magnetic Reynolds number within the CZ is very large, $`R_m10^8`$ Sokoloff , and the corresponding inertial range is very wide, $`l_{\mathrm{diss}}<l<L_{\mathrm{max}}`$ where the outer scale is supposed to be $`L_{\mathrm{max}}1000`$ km, about the solar granule size, and $`l_{\mathrm{diss}}=L_{\mathrm{max}}R_m^{3/4}1`$ m. The rms magnetic field $`b_l`$ is assumed to scale as $`b_ll^{1/3}`$. This implies that, after fixing the maximum field amplitude at the outer scale $`L_{\mathrm{max}}`$ it is straightforward to obtain the rms field at the neutrino oscillation scale, which is about several hundreds kilometers for LMA-MSW case, well within the inertial range where scaling arguments are valid Inprep . Here we show only the final result for neutrino mass eigenstate probabilities at the surface of the Sun:

$$|\nu _{1L}|_R_{}^2=P_1(1\eta ),|\nu _{2R}|_R_{}^2=P_1\eta ,$$
(1)
$$|\nu _{2L}|_R_{}^2=P_2(1\eta ),|\nu _{1R}|_R_{}^2=P_2\eta ,$$
(2)
(note that these obey unitarity) where the common factor $`\eta `$ is given by

$$\eta 0.3\frac{\mu _{\nu }^{}{}_{}{}^{2}\overline{b}_{\mathrm{max}}^2}{\delta ^2}\frac{\delta L}{(\delta L_{\mathrm{max}})^{2/3}}S^2.$$
(3)
Here $`P_i=P_{iL}=|\nu _{iL}(r=0.7R_{})|^2`$ are the probabilities that solar neutrinos reach the bottom of the CZ in a given mass state, $`L`$ is the CZ width, $`\delta =\mathrm{\Delta }m^2/4E`$, and $`S^2`$ is a rms magnetic field profile shape factor

$$S^2=\frac{1}{L}_0^L𝑑z\frac{\overline{b}^2(z)}{\overline{b}_{max}^{}{}_{}{}^{2}}.$$
(4)
Note that $`S=1`$ for constant shape of rms field and of the order of unity for other profiles, e.g. $`S=0.579`$ for the “smooth” profile of Bykov:1998gv , and $`S=0.577`$ for a triangle profile.

The different factors in Eq. 3 can be explained if we keep in mind that $`\delta `$ is inversely proportional to the neutrino oscillation length, $`\delta =\pi /\lambda _{osc}`$ . The ratio $`\overline{b}_{\mathrm{max}}^2/(\delta L_{\mathrm{max}})^{2/3}=\overline{b}_{\lambda _{\mathrm{osc}}}^2`$ is the squared rms field at the scale $`l=1/\delta \lambda _{\mathrm{osc}}`$, and the ratio $`\mu _\nu ^2\overline{b}_{\lambda _{\mathrm{osc}}}^2/\delta ^2`$ determines the fraction of neutrinos which experience (on average) spin-flavour conversion to the antineutrino states in a correlation cell of length $`l\lambda _{\mathrm{osc}}`$. Finally $`\delta LL/\lambda _{\mathrm{osc}}`$ is the number of correlation cells along the neutrino trajectory and the factor $`S^2`$, as already mentioned, accounts for the specific shape of the rms field profile. (The numerical factor 0.3 comes from the integration along the trajectory).

We can rewrite Eq. 3 in normalized units as

$$\eta 3\times 10^3\mu _{11}^2\epsilon ^2S^2\left(\frac{7\times 10^5eV^2}{\mathrm{\Delta }m^2}\right)^{5/3}\left(\frac{E}{10\mathrm{M}\mathrm{e}\mathrm{V}}\right)^{5/3}$$
(5)
where $`\mu _{11}`$ is the magnetic moment in units of $`10^{11}\mu _B`$, and the ratio $`\epsilon =(b/100\mathrm{kG})(L_{\mathrm{max}}/1000\mathrm{km})^{1/3}`$.

Although the ratio $`\epsilon `$ is not known precisely, we can obtain a good estimate assuming the equipartition between kinetic energy of hydrodynamic pulsations and the rms magnetic energy at the upper (most energetic) scale $`L_{max}`$ Sokoloff

$$\frac{\rho \overline{v^2}}{2}\frac{\overline{b}_{\mathrm{max}}^2}{8\pi }$$
(6)
Taking $`v3\times 10^4\mathrm{cm}/\mathrm{s}`$, $`\rho 1g/\mathrm{cm}^3`$ we obtain $`b_{\mathrm{max}}100`$ kG. In what follows we will consider the range $`b_{\mathrm{max}}50100`$ kG.

This specific behaviour implies, on physical grounds, that $`\epsilon `$ is not allowed to vary substantially, $`0.5<\epsilon <1`$, with our assumptions. As mentioned, for different profiles we found that the parameter $`S^2`$ lies in the range between 0.5 and 1. With this in mind, we can see that the product $`k=\epsilon S`$ lies in the interval $`0.25<k<1`$. Since the overall $`\nu _e\overline{\nu _e}`$ conversion probability depends linearly on $`\eta `$, it follows that, by comparing the resulting $`\overline{\nu _e}`$ yield to the recent KamLAND bound one can constrain the value of $`\mu _\nu `$ to within a factor 4. More precisely, we have computed the solar anti-neutrino yield for this random magnetic field model for all allowed values of neutrino oscillation parameters in the LMA-MSW region. Our results are shown in fig. 3 as a function of $`\mu _\nu `$, where the width of the band corresponds to the freedom in choosing the parameter $`k`$. It is easy to see that, indeed, the constraint we obtain is better than those that hold for the case of regular magnetic fields. Moreover, from our estimate for the minimum value of $`k`$, we can obtain an upper bound for $`\mu _\nu `$ which lies in the range $`\mu _\nu 1.5\times 10^{12}\mu _B`$ to $`\mu _\nu 5\times 10^{12}\mu _B`$.

We have obtained this result for the Kolmogorov theory. However, it can be easily generalized for arbitrary power-law exponent $`p`$ (for Kolmogorov theory $`p=5/3`$). We have checked that our bound on $`\mu _\nu `$ is rather robust with respect to possible changes of the power-law exponent $`p`$ within the interval from 1 to 2 Inprep .

In short, we have discussed the impact of recent solar neutrino data on the spin-flavor-precession (SFP) mechanism of neutrino conversion. The recent KamLAND bound on the solar anti-neutrino yield leads to a constraint on the product of Majorana neutrino transition magnetic moment and some magnetic field strength in different models for the solar magnetic field. For the case of random magnetic fields inside the Sun, one can obtain a direct constraint on the intrinsic neutrino magnetic moment of $`\mu _\nu <\mathrm{few}\times 10^{12}\mu _B`$, similar to bounds obtained from star cooling. Comparing these constraints with Fig. 1 one sees that the robustness of the oscillation interpretation of current solar neutrino data against possible magnetic-field-induced transitions is firmly established.

We thank V. B. Semikoz and D. D. Sokoloff for useful discussions. This work was supported by Spanish grant BFM2002-00345, by European RTN network HPRN-CT-2000-00148, by European Science Foundation network grant N. 86, by INTAS YSF grant 2001/2-148 and MECD grant SB2000-0464 (TIR). TIR and AIR were partially supported by the Presidium RAS and CSIC-RAS grants. OGM was supported by CONACyT-Mexico and SNI.


text:
# Extended tachyon field, Chaplygin gas and solvable 𝑘-essence cosmologies

## I Introduction

Cosmological inflation has become an integral part of the standard model of the universe and perhaps the only known mechanism which can dynamically solve the flatness and the horizon problem of the universe . It has gained certain support from the recent observations of the cosmic microwave background anisotropies -. Although particle physics, in particular M/String theory, provides several very weakly coupled scalar fields which are natural inflaton candidates, there exists no clearly preferred inflationary model. A link between string theory and inflation was investigated in Ref.  where the authors introduce the $`k`$-inflation and show that the $`k`$-field may drive an inflationary evolution starting from rather generic initial conditions. Also, the process of rolling tachyon field has been extensively used to implement power-law accelerated expansion scenarios -. In Ref.  it was pointed out that rolling tachyons can contribute a mass density to the universe that resembles classical dust. This has brought a new understanding of the role of the tachyon in string theory. The tachyonic matter, then, might provide an explanation for inflation at the early epochs and could contribute to some new form of cosmological dark energy  at late times . These facts suggest to enlarge the theoretical background exploring some new possibilities, in this sense, we will introduce an extended tachyon field (ETF) which supply accelerated and superaccelerated expanding scenarios.

In Ref.  it was argued that the coincidence problem may be solved assuming a universe filled with a viscous fluid and dark energy modelled with a tachyonic field or a Chaplygin gas. The conjecture that dark matter and dark energy can be unified by using a generalized Chaplygin gas obeying an exotic equation of state has been investigated in several works in view of the cosmological observations -. The present situation is somewhat controversial, with some tests indicating good agreement between observational data and the theoretical predictions of the model and others ruling out the model as an actual possibility of description for our Universe . There were found some differences with the observed CMB and mass power spectrum data . Due to these discrepancies, it may be useful to consider other candidates for dark energy, as for instance, the extended Chaplygin gas (ECG) generated from the ETF. Another interesting possibility appears selecting a $`k`$-field model whose equation of state is like that of two fluids, one obeying a baryotropic equation of state with constant baryotropic index $`\gamma `$, and the other is an ECG. This ”modified Chaplygin gas” interpolates between a power-law phase and a de Sitter phase.

From theoretical and experimental points of view it is important to find an exact shape of the potential, for instance, tachyonic inflation has been studied using phenomenological potentials that have not been derived from string theory and can be related to the so called “k-Inflation” , . Such $`k`$-fields, described by a non-standard kinetic term, has been one of the recent suggested candidates to play the role of same unknown component of the universe. One of the purpose of introducing $`k`$-essence is to provide a dynamical explanation of cosmological observations which does not require the fine-tuning of initial conditions -. In this sense, it may be considered as an alternative to quintessence, which requires a scalar field that slowly rolls down a potential to depict the observed acceleration of the present universe. Also, it was argued that in certain dynamical regimes the $`k`$-essence is equivalent to quintessence and it may prove difficult to distinguish between the two fields. In the light of these results, it seems clear that it is worth to search the links between scalar field and $`k`$-essence cosmologies, in especial, whether they are kinematically or dynamically equivalent.

In section II, we introduce the ETF and find the power-law solution generated by an inverse square potential. At the end, we define the ECG and show the general behaviour of the background geometry. In section III, we find the first integral of the $`k`$-field equation for an inverse square potential and the general solution of coupled Einstein-$`k`$-field equations for a linear function $`F`$. For power-law solutions we demonstrate that the linear $`k`$-field model, with constant function $`F`$, is isomorphic to the model with divergent sound speed. Finally, we introduce the modified Chaplygin gas. In section IV, we show the kinematical equivalence of the $`k`$-field driven by a tachyonic potential and a linear function $`F`$, with the scalar field driven by an exponential potential. The conclusions are stated in section V.

## II Extend tachyon field and Chaplygin gas

The spatially flat homogeneous and isotropic space-time is described by the Friedmann-Robertson-Walker (FRW) line element

$$ds^2=dt^2+a^2(t)\left[dx_1^2+dx_2^2+dx_3^2\right],$$
(1)
where $`a(t)`$ is the scale factor. This metric allows a particular calculational simplicity, on account of both the high degree of symmetry and the single metric degree of freedom. We assume that $`8\pi G=1`$.

Let us present the Lagrangian of the $`k`$-essence field

$$=V(\varphi )F(x),x=g^{ik}\varphi _i\varphi _k,$$
(2)
where $`V(\varphi )`$ is a positive defined potential, $`F(x)`$ is an arbitrary function of $`x`$, $`\varphi `$ is the $`k`$-field and $`\varphi _i=V(\varphi )/x^i`$. Associating the energy-momentum tensor of the $`k`$-field,

$$T_{ik}=V(\varphi )\left[2F_x\varphi _i\varphi _kg_{ik}F\right],F_x=\frac{dF}{dx},$$
(3)
with the energy-momentum tensor of a perfect fluid, we compute the energy density $`\rho `$, the pressure $`p`$ and the baryotropic index $`\gamma =1+p/\rho `$ of this equivalent fluid

$$\rho =V(\varphi )[F2xF_x],p==V(\varphi )F,$$
(4)
$$\gamma =\frac{2\dot{H}}{3H^2}=\frac{2xF_x}{F2xF_x}.$$
(5)
From Eq. (4), the Einstein field equations are

$$3H^2=V[F2xF_x],$$
(6)
$$\dot{H}=VxF_x,$$
(7)
and the conservation equation reads

$$\dot{\rho }+3H(\rho +p)=0,$$
(8)
where $`H=\dot{a}/a`$ is the expansion rate. Substituting the Eq. (4) into the conservation equation (8), we find the field equation for the $`k`$-field

$$[F_x+2xF_{xx}]\ddot{\varphi }+3HF_x\dot{\varphi }+\frac{V^{}}{2V}[F2xF_x]=0.$$
(9)
The stability of the $`k`$-essence with respect to small wavelength perturbations requires that the effective sound speed

$$c_s^2=\frac{p_x}{\rho _x}=\frac{F_x}{F_x+2xF_{xx}},$$
(10)
be positive. However, in Ref.  it was shown that a positive sound speed is not a sufficient condition for the theory to be stable. For the tachyon field, which is obtained from the $`k`$-field selecting $`F=(1\dot{T}^2)^{1/2}`$, with $`\dot{T}^2=x=\dot{\varphi }^2`$ the sound speed is $`c_s^2=1\gamma >0`$, where $`\gamma =\dot{T}^2`$. Let us look for the set of functions $`F`$, such that, $`c_s^2`$ is proportional to the sound speed of the tachyon field. We express this proportionality as

$$c_s^2=\frac{1\gamma }{2r1},$$
(11)
where $`r`$ is a real constant, with $`r<1/2`$ for $`\gamma >1`$ and $`r>1/2`$ for $`\gamma <1`$. From Eqs. (5), (10) and (11), we obtain a differential equation for the function $`F`$

$$(1r)FF_x+(2r1)xF_x^2+xFF_{xx}=0.$$
(12)
Integrating for $`r=0`$, we obtain the set of functions $`F_0(x)=x^{\gamma _0/2(\gamma _01)}`$, with $`\gamma _0`$ constant, investigated in ; while for $`r=1`$ we find three types of solution $`F_1=(1\dot{T}_1^2)^{1/2}`$, $`F_1=(1+\dot{T}_1^2)^{1/2}`$, and $`F_1=(\dot{T}_1^21)^{1/2}`$. The first function corresponds to the ordinary tachyon field $`T=T_1`$. Its energy density $`\rho _1=VF_1^1`$ and pressure $`p_1=VF_1`$ leads to the relation $`p_1=(\gamma _11)\rho _1`$, where $`\gamma _1=\dot{T}_1^2`$. The other two functions generate two ”new tachyon field” which will be investigated below. For the remaining values of $`r`$ the general solution of Eq. (12) is given by $`F_r^{2r}=c_1(x)^r+c_2`$, where $`c_1`$ and $`c_2`$ are arbitrary integration constants, the baryotropic index (5) is $`\gamma _r=c_1(x)^r/c_2`$ and we have $`F_r^{2r}=c_2(1\gamma _r)`$. From Eqs. (4) and (5), the energy density of the $`k`$-field can be rewritten as $`\rho =VF/(1\gamma )`$. So, without loss of generality, we can split the solutions as follow: for $`\gamma _r<1`$ and $`r>1/2`$, it is necessary that $`F_r>0`$ and $`c_2>0`$, with two options, $`0<\gamma _r<1`$ and $`c_1<0`$ or $`\gamma _r<0`$ and $`c_1>0`$; for $`\gamma _r>1`$ and $`r<1/2`$, we have $`F_r<0`$, $`c_2`$ and $`c_1>0`$. Finally

$$F_r=(1\gamma _r)^{1/2r},\gamma _r=\pm \dot{T}_r^{2r}<1,$$
(13)
$$F_r=(\gamma _r1)^{1/2r},\gamma _r=\dot{T}_r^{2r}>1,$$
(14)
where $`\gamma _r`$ is the extended baryotropic index and the $`T_r`$-field linked with $`F_r`$ will be called ETF. Inserting Eqs. (13) and (14) into (4), we get the energy density $`\rho _r`$ and the pressure $`p_r`$ of each ETF

$$\rho _r=V(1\dot{T}_r^{2r})^{(12r)/2r},p_r=V(1\dot{T}_r^{2r})^{1/2r},0<\gamma _r<1,$$
(15)
$$\rho _r=V(1+\dot{T}_r^{2r})^{(12r)/2r},p_r=V(1+\dot{T}_r^{2r})^{1/2r},\gamma _r<0,$$
(16)
$$\rho _r=V(\dot{T}_r^{2r}1)^{(12r)/2r},p_r=V(\dot{T}_r^{2r}1)^{1/2r},1<\gamma _r,$$
(17)
with

$$p_r=(\gamma _r1)\rho _r.$$
(18)
The Eq. (15) comprise perfect fluids which generalize the normal tachyon field, which is obtained for $`r=1`$. The fluids represented by the Eq. (16) with negative pressure and negative baryotropic index describe phantom cosmologies, while those described by the Eq. (17) give rise to nonaccelerated expanding evolutions. In the limit of large $`r`$, the fluids (15) and (16) satisfy the equation of state $`p=\rho =V`$ acting like a variable cosmological constant depending on the k-field. Besides, the fluids (17) fulfill the equation of state $`p=\rho =V`$ behaving like stiff matter in the same limit. These exotic fluids satisfy the relations

$$p_r=\frac{V^{2r/(2r1)}}{\rho _r^{1/(2r1)}},\gamma _r<1,$$
(19)
$$p_r=\frac{V^{2r/(2r1)}}{\rho _r^{1/(2r1)}},\gamma _r>1,$$
(20)
which, for a constant potential $`V=V_0`$, become exotic equations of state which extend that of generalized Chaplygin gas.

### II.1 Power-law expansion for the $`T_r`$-field

We begin investigating the atypical behaviour of the $`T_{1/2}`$-field. Its sound speed diverges and the $`k`$-field equation (9) becomes a first order equation. The equation (5) can be integrated and its solution gives a link between the expansion rate and $`T_{1/2}`$

$$H=\frac{2}{3T_{1/2}},$$
(21)
where the integration constant has been chosen to set the singularity at $`T_{1/2}=0`$. For $`r=1/2`$, the three equations (15), (16) and (17) lead to $`3H^2=V`$. So, combining them with Eq. (21) we obtain the inverse square potential

$$V_{1/2}=\frac{4}{3T_{1/2}^2},$$
(22)
which is the uniquely allowed potential for $`T_{1/2}`$. As $`T_{1/2}`$ is not controlled by the $`k`$-field equation (9), it can be chosen freely and this choice determines, after integrating Eq. (21), the form of the scale factor. Although the potential (22) diverges at $`T_{1/2}=0`$, it reasonably mimics the behaviour of a typical potential in the condensate of bosonic string theory. One expects the potential to have a unique local maximum at the origin and a unique global minimum away from the origin at which $`V`$ vanishes. In the most interesting case the global minimum is taken to lie at infinity. Obviously more complicated potentials may be contemplated but this is the simplest case to begin with.

It will be useful to investigate the existence of accelerated and superaccelerated expanding solutions when the $`T_r`$-field is driven by the potential $`V_r=V_0/T_r^2`$. To do that, let us consider an evolution of the form $`a=t^n`$ and a $`T_r`$-field, such that $`T_rt`$. Power law solutions are very important because they can be always obtained from any function $`F`$ with an inverse square potential or from polynomial functions $`F=(x)^N`$ with any potential. As these solutions are usual ingredients in the quintessense and $`k`$-essence models they would allow to recognize and compare the differences between the two cosmological models. The complete solution for the ETF is given by

$$T_r=\left(\frac{2}{3n}\right)^{1/2r}t,n=\frac{1}{3}\left[1+\sqrt{1+9\beta ^2}\right],0<\gamma _r<1,$$
(23)
$$T_r=\left(\frac{2}{3n}\right)^{1/2r}t,n=\frac{1}{3}\left[1\sqrt{1+9\beta ^2}\right],\gamma _r<0,$$
(24)
$$T_r=\left(\frac{2}{3n}\right)^{1/2r}t,n=\frac{1}{3}\left[1\pm \sqrt{19\beta ^2}\right],1<\gamma _r,$$
(25)
where

$$\beta ^2=\left[\frac{3}{V_0}\left(\frac{2}{3}\right)^{1/r}\right]^{2r/(12r)}.$$
(26)
The sound speed associated with the above solutions is given by

$$c_s^2=\frac{12/3n}{2r1}.$$
(27)
Assuming a positive sound speed bounded by the light speed $`0c_s^21`$ we can see that $`r>11/3n`$ is required for $`n>2/3`$ or $`n<0`$ and $`r<11/3n`$ for $`0<n<2/3`$.

An accelerated power-law expanding universe ($`n>1`$), can be described by the set of ETF having $`r>2/3`$ and it is represented by the solutions (23). All these models are cinematically and dynamically different because the scale factor and the $`k`$-field are strongly depending of $`r`$ as can be seen from Eqs. (23) and (26). In particular, for $`r=1`$ the Eq. (23) reduces to the solution of the tachyon field found in . The solutions (24) with $`n<0`$ describe phantom cosmologies and we have to select ETF with $`r>1`$ to reach this scenario

### II.2 Extended Chaplygin gas

Recently it was proposed a set of simple cosmological models based on the use of peculiar perfect fluids . In this simple model the universe is filled with the so called Chaplygin gas, which is a perfect fluid characterized by the following equation of state $`p=A/\rho `$, where $`A`$ is a positive constant. It describes a transition from a decelerated cosmological expansion to a cosmic accelerated de Sitter stage. Other possibility is the inhomogeneous Chaplygin gas, which is able to combine the roles of dark energy and dark matter , and the generalized Chaplygin gas model discussed in Ref. , having two free parameters $`p=A/\rho ^\alpha `$ with $`0<\alpha 1`$.

These cosmological models may give a unified macroscopic phenomenological description of dark energy and dark matter and generalize the usual $`\mathrm{\Lambda }`$CDM models. On the other hand, the Chaplygin gas can be considered as the simplest tachyon cosmological models where the tachyon field is a purely kinetic $`k`$-essence model with a constant potential. In the same way, we will show that the generalized Chaplygin gas can be conceived as the simplest ETF model driven by a constant potential. This identification has the advantage of producing a variety of new Chaplygin gases, same of which leads to superaccelerated scenarios as we will see in this subsection.

Coming back to our Eqs. (19) and (20), we see that in the case of a constant potential $`V(T_r)=V_0`$, with $`r=1`$ and $`0<\gamma _r<1`$, the exotic equation of state (19) represents a Chaplygin gas . Such equation may be the consequence of a scalar field with a non-standard kinetic term, e.g., the string theory motivated tachyon field . For $`1<r`$ and $`0<\gamma _r<1`$, the equation of state (19) represents a perfect fluid called generalized Chaplygin gas . Finally an ECG will be characterized by the equation of state (19) with $`1/2<r<1`$, or by the Eq. (20) with $`r<1/2`$. Their properties will be investigated below. Using Eqs. (19)-(20) and the relativistic energy conservation equation (8), we obtain the energy density of the ECG:

$$\rho =V_0\left[1+\left(\frac{a_0^3}{a^3}\right)^{2r/(2r1)}\right]^{(2r1)/2r},0<\gamma _r<1,$$
(28)
$$\rho =V_0\left[1\left(\frac{a_0^3}{a^3}\right)^{2r/(2r1)}\right]^{(2r1)/2r},\gamma _r<0,$$
(29)
$$\rho =V_0\left[1+\left(\frac{a_0^3}{a^3}\right)^{2r/(2r1)}\right]^{(2r1)/2r},1<\gamma _r.$$
(30)
The Eq. (28) with $`r>1`$ gives the energy density of the generalized Chaplygin gas in terms of the scale factor interpolating between a dust dominated phase where the energy density is $`\rho V_0(a_0/a)^3`$, and a de Sitter phase where $`\rho V_0`$. While during the intermediate stage it could be interpreted as a mixture of two-fluids, one of which is the cosmological constant and the other is a perfect fluid with equation of state $`p\rho `$. The additional free parameter $`r`$ of the generalized Chaplygin gas can be used to compare it with observational data.

In any other case, Eqs. (28) with $`1/2<r<1`$, (29) and (30) represent new perfect fluids which are interesting from the cosmological point of view. The scale factor generated by the source (29) is nonsingular and has a bounce at the minimum $`a=a_0`$. The universe begins from a contracting era and ends in a superaccelerated stage. Such cosmologies may be interpreted as universes filled with baryotropic fluids having a negative constant baryotropic index and violating the weak energy condition $`\rho >0`$, $`\rho +p>0`$. The models will be dubbed phantom cosmologies following the standard terminology. Phantom matter can apparently be accommodated by current observations , and it can be based on the motivation provided by string theory . It looks interesting to admit that the origin of dark energy should be searched within a fundamental theory, as string theory. However, at present there is no consensus whether a universe violating the weak energy condition should generically possess a future singularity or big rip .

The source (30) leads to two types of solutions according to the values of the parameter $`r`$. For $`r<0`$ the universe has a finite time span, interpolating between two dust dominated phases and has a maximum at $`a=a_0`$. For $`0r1/2`$ the universe begins at a singularity with a finite scale factor, $`a=a_0`$, and ends in a dust dominated phase. From the cosmological point of view these solutions are not relevant because they do not describe the present observed accelerated expansion stage.

## III Solvable k-essence cosmologies

In this section, we will show some cases where the coupled Einstein-$`k`$-field equations (6,9) admit a first integral or can be solved exactly. Expressing the energy density of the $`k`$-field as $`\rho =VF/(1\gamma )`$ and using the conservation equation (8), we get the $`k`$-field equation in terms of the baryotropic index $`\gamma `$

$$\left(\frac{\gamma }{\dot{\varphi }}\right)^.+3H\left(\frac{\gamma }{\dot{\varphi }}\right)(1\gamma )+\frac{V^{}}{V}(1\gamma )=0.$$
(31)
where $`V^{}=dV/d\varphi `$. This form of writing the field equation allows to show that the first integral of the $`k`$-field equation (9) for any function $`F`$ is giving by

$$\frac{\gamma }{\dot{\varphi }}=\varphi ^1\left(\frac{2}{3H}+\frac{c}{a^3H^2}\right),$$
(32)
where we have ssumed an inverse square potential

$$V=\frac{V_0}{\varphi ^2},$$
(33)
and $`c`$ is an arbitrary integration constant. The usual linear field solution $`\varphi =\varphi _0t`$ along with the power law scale factor $`a=t^{2/3\gamma }`$ (with constant $`\gamma `$) is obtained from the last equation for a vanishing integration constant. Combining Eqs. (5), (7) and (33) with Eq. (32), it can be rewritten as

$$\dot{\varphi }F_x\left[H+\frac{3c}{2a^3}\right]\frac{\varphi }{V_0}=0,$$
(34)
or as

$$V_0F_x\dot{H}=\left[H+\frac{3c}{2a^3}\right]^2.$$
(35)
Then, for the tachyonic potential (33), Eqs. (32), (34), (35) are three different forms of writing the first integral of the $`k`$-field equation (9) or (31).

At this point, using Eqs. (5), (33) and (32), it is interesting to see the coupled Einstein-$`k`$-field equations (6,9) as a system of differential equations for the function $`F`$. Hence, integrating we obtain

$$F=\frac{3h_0^2}{V_0}+\sqrt{x}\left[b+\frac{3c}{2V_0}\frac{(2h_0+ch)h}{(x)^{3/2}}𝑑x\right],$$
(36)
$$H\varphi =\frac{3c}{2}\frac{d\varphi }{a^3}=h_0+ch,$$
(37)
where $`h_0`$ and $`b`$ are arbitrary constants and $`\dot{h}=3\dot{\varphi }/2a^3`$. In the particular case $`c=0`$, the Eq. (36) becomes

$$F^{\mathrm{}}=\frac{3h_0^2}{V_0}+b\sqrt{x},$$
(38)
which, after a redefinition of the constants, turns into the function $`F_{1/2}`$ that generates the ETF, $`T_{1/2}`$, (see Eqs. (13,14) for $`r=1/2`$). For this ”divergent” $`k`$-essence theory the sound speed (10) diverges and the $`k`$-field equation (9) becomes a first order equation that is only consistent with an inverse square potential. This divergent model is related to the linear $`k`$-field model, $`\varphi =\varphi _0t`$, obtained by evaluating $`F`$, $`F_x`$ at $`x=x_0=\varphi _0^2`$ into the Einstein-$`k`$-field equations (6,9). The linear $`k`$-field model driven by the potential (33) leads to the power-law solutions $`a=t^\lambda `$ with

$$\lambda =\frac{1}{3}\frac{f+2\varphi _0^2f^{}}{\varphi _0^2f^{}},V_0=\frac{\lambda }{f^{}},$$
(39)
where we have defined $`f=F(\varphi _0^2)`$, $`f^{}=F_x(\varphi _0^2)`$.

From Eqs. (38,39) it is easy to show that the linear model is isomorphic to the divergent one. This can be done by constructing a one-to-one mapping between these two models. In fact, choosing in the divergent model the constants $`h_0^2=\lambda ^2\varphi _0^2`$ and $`b=2\lambda \varphi _0/V_0`$ we find the same power-law solutions obtained from the linear model. In Ref.  it was suggested that this might be the reason as to why a model with a diverging sound speed leads to serious problems as discussed in a recent paper . In addition, assuming a series expansion of the function $`F(x)`$ around $`x=x_0`$, the background cosmology is completely determined by the first two coefficients $`(f,f^{})`$ of the expansion of $`F`$ and the value of $`\varphi _0`$. Hence, the model is insensitive to the remaining coefficients in the expansion of the function $`F`$ and both the linear and divergent models should be considered as equivalent. This means that the power-law solutions generated by an inverse square potential possesses a degeneracy. Possibly this degeneracy may be removed by perturbating the solutions.

### III.1 Inverse square potential and linear function $`F`$

From the $`k`$-field equation (35), it can be seen that a linear function

$$F=1+mx,$$
(40)
where $`m`$ is a constant, decouples the dynamics of the background geometry from the dynamics of the $`k`$-field . Therefore, this choice clearly introduces a rather small degree of nonlinearity into the dynamical equations allowing us to obtain the general solution of the coupled Einstein-$`k`$-field equations. On the other hand, the function (40) mimics the behaviour of other models. For instance, when $`x1`$ the tachyonic function $`F=(1+x)^{1/2}`$ can be approximated by $`F1+x/2`$ and it has the form (40) . In Ref.  it was introduced a set of models where $`F`$ admits a power series expansion similar to (40). This form is reminiscent of a Born-Infeld action with higher order corrections in $`x`$, and particular cases were investigated in . Hence, the knowledge of the general solution for the linear function (40) should be of interest, at least for understanding the asymptotic behaviour of many others models generated by an analytical function $`F(x)=F(0)+F_x(0)x+\mathrm{}\mathrm{}`$. In fact, keeping the first order term in the expansion of $`F`$, it adopts the linear form (40), after redefining the potential $`V`$ to set the constant $`F(0)=1`$.

Combining Eqs. (35) and (40), we obtain the following nonlinear second-order differential equation for the scale factor

$$\frac{d^2s}{d\tau ^2}+s^\sigma \frac{ds}{d\tau }+\frac{1}{4}s^{2\sigma +1}=0,\sigma =3mV_0,$$
(41)
where we have used the new variables $`s`$ and $`\tau `$, defined by

$$s=a^{3/\sigma },\tau =\frac{3c}{mV_0}t.$$
(42)
The general solution of Eq. (41) can be found changing to the nonlocal variables $`z`$ and $`\eta `$, defined by

$$z=\frac{s^{\sigma +1}}{\sigma +1},\eta =s^\sigma 𝑑\tau ,\sigma 1,$$
(43)
$$z=\mathrm{ln}s,\eta =\frac{d\tau }{s},\sigma =1.$$
(44)
Then, in these new variables, the Eq. (41) becomes a linear homogeneous differential equation with constant coefficients

$$\frac{d^2z}{d\eta ^2}+\frac{dz}{d\eta }+\frac{\sigma +1}{4}z=0,\sigma 1,$$
(45)
equivalent to a damped harmonic oscillator equation, and

$$\frac{d^2z}{d\eta ^2}+\frac{dz}{d\eta }+\frac{1}{4}=0,\sigma =1.$$
(46)
On the other hand, expressing the Eq. (34) for the $`k`$-field in terms of the independent variable $`\eta `$, we get

$$\frac{d\varphi }{\varphi d\eta }=\frac{3da}{\sigma ad\eta }+\frac{1}{2}.$$
(47)
Once $`z(\eta )`$ is known from Eq. (45), one can compute $`s(\tau )`$ from Eq. (43), the scale factor $`a(\tau )`$ from (42), and the $`k`$-field $`\varphi (\tau )`$ from equation (34). Following this procedure and inserting the general solutions of Eqs. (45) and (47) into the Einstein equation (6), we find the scale factor, the $`k`$-field and the relationship amongst the integration constants

$$a=\left[\sqrt{B}\text{e}^{\frac{\eta }{2}}\mathrm{sinh}\left(\frac{\sqrt{\sigma }}{2}\eta +\eta _0\right)\right]^{\sigma /3(\sigma +1)},\sigma <1,$$
(48)
$$a=a_0\text{e}^{\eta /12+V_0\text{e}^\eta /81c^2\varphi _0^2},\sigma =1,$$
(49)
$$a=\left[\sqrt{B}\text{e}^{\frac{\eta }{2}}\mathrm{cosh}\left(\frac{\sqrt{\sigma }}{2}\eta +\eta _0\right)\right]^{\sigma /3(\sigma +1)},1<\sigma <0,$$
(50)
$$a=\left[\sqrt{B}\text{e}^{\frac{\eta }{2}}\mathrm{sin}\left(\frac{\sqrt{\sigma }}{2}\eta +\eta _0\right)\right]^{\sigma /3(\sigma +1)},0<\sigma ,$$
(51)
the $`k`$-field is given by

$$\varphi =\varphi _0a^{3/\sigma }\text{e}^{\eta /2},$$
(52)
where

$$B=\frac{4(\sigma +1)V_0}{27c^2\varphi _0^2},$$
(53)
$`\eta _0`$ and $`\varphi _0`$ are arbitrary integration constant.

For $`\sigma <1`$, the solution expands from a singularity as $`t^{1/3}`$ and ends as $`t^{\sigma /3}`$. When $`\sigma >3`$ the scale factor displays a power-law inflationary scenario. For $`1<\sigma <0`$, the universe expands from a singularity as $`t^{1/3}`$ and its final behaviour is given by $`t^{1/3}`$. For $`\sigma >0`$ the solution represents a contracting universe which begins at a finite time, reaches a minimum where it bounces, exhibiting a final superaccelerated expansion $`a(t)(t_0t)^{\sigma /3}`$. The universe has a finite time span and bounces when the $`k`$-field satisfies the condition $`\dot{\varphi }^2=1/m`$.

As this model displays an accelerated expanding stage at late times it may be an interesting alternative to describe the epoch where dark energy dominates.

### III.2 The explicit solution

For $`\sigma =4`$ or $`mV_0=4/3`$, we can solve the Eq. (41) by means of the substitution

$$s^4=\frac{1}{2}\frac{v^4}{v^4𝑑\tau },$$
(54)
so that Eq. (41) reads as

$$\ddot{v}=0.$$
(55)
Inserting its solution into (54), using (42) and integrating (47), we get the general solution for the evolution and the $`k`$-field satisfying the Eq. (6):

$$a(\tau )=\left(\frac{3}{2}\right)^{1/3}\left[\frac{9c^2V_0}{16\varphi _0^2}t^4ct\right]^{1/3},$$
(56)
$$\varphi (\tau )=\frac{2}{3t^{1/2}}\left[\frac{\varphi _0^2}{c}+\frac{9V_0}{16}t^3\right]^{1/2},$$
(57)
where we have chosen $`c<0`$ to set the singularity at $`t=0`$. The scale factor (56) exhibits a transition from $`at^{1/3}`$ to an accelerated expansion where $`at^{4/3}`$. Curiously, it coincides with the solutions found in  where it was considered a FRW spacetime filled with a scalar field driven by an exponential potential. In the next section we will investigate this relation between the scalar field and the $`k`$-essence field. The $`k`$-field diverges at the singularity as $`\varphi t^{1/2}`$ and behaves as $`\varphi t`$ for large time. It has a minimum and a turning point where the kinetic energy vanishes.

Another set of solutions can be found when $`c=0`$. Here, the Eq. (35) reduces to $`\dot{H}=H^2/mV_0`$ and we get the power-law expansion $`a=t^{mV_0}`$ and a linear $`k`$-field $`\varphi =\varphi _0t`$.

### III.3 The polynomial function $`F_\gamma (x)=(x)^{\frac{\gamma }{2(\gamma 1)}}`$

This polynomial function yields a constant baryotropic index $`\gamma `$ and the power-law expansion $`a=t^{2/3\gamma }`$ . For $`F_\gamma `$ the general solution of Eq. (34) is

$$\varphi =\left[b+4\varphi _0t^{(2\gamma )/\gamma }\right]^{\gamma /(2\gamma )},\gamma 2,$$
(58)
$$\varphi =\varphi _0t^b,\gamma =2,$$
(59)
where $`b`$ and $`\varphi _0`$ are integration constants. Consistency between these solutions and the Eq. (6) gives the following relation between the integration constants,

$$V_0=\frac{4(1\gamma )}{3\gamma ^2}(4\varphi _0)^{\gamma /(1\gamma )},\gamma 2$$
(60)
$$b=\pm \frac{1}{\sqrt{3V_0}},\gamma =2.$$
(61)
For large cosmological time the $`k`$-field behaves as $`\varphi t`$. More details about the cosmological model generated by the set $`F_\gamma `$ can be found in Ref. .

### III.4 Constant potential case

We have seen that the generalized Chaplygin gas model was proposed as unified dark matter. It is derived from a Lagrangian containing non-standard kinetic-energy terms (i.e., non-quadratic) and can be considered as driven by a constant potential. Below we show that, even in absence of any potential energy term, a general class of models exists, including the ECG, connecting a dust dominated era at early times with an accelerated expansion stage at late times.

For a constant potential $`V=V_0`$, the $`k`$-field equation (31) has the first integral

$$\frac{\gamma }{\dot{\varphi }}=\frac{c}{a^3H^2},$$
(62)
or, after using Eqs. (5) and (7), it turns into

$$a^3F_x\dot{\varphi }=\frac{3c}{2V_0},$$
(63)
where $`c`$ is an arbitrary integration constant. Also, combining them with the Friedmann equation (6), the baryotropic index associated with this kind of $`k`$-essence can be written in a more convenient form

$$\gamma =\frac{1}{1+2V_0^2a^6FF_x/9c^2}.$$
(64)
From the last equation we see that for a large set of model, i.e., the class of model generated by the set of functions $`F`$ satisfying the condition $`a^6FF_x1`$ at early time, the universe is dust dominated in the beginning. At intermediate times it behaves as it were filled with a perfect fluid with equation of state $`p\rho `$. Finally, the universe ends in an accelerated expansion scenario. So, these alternative models play the same role that the generalized Chaplygin gas i.e., interpolating between dark matter at early time and dark energy at late time.

Now, we investigate a simple kinetic $`k`$-essence model generated by a functon $`F`$ satisfying the more general condition $`a^6FF_xconstant`$ at early time. This model is generated by the following function

$$F=\frac{1}{(2\alpha 1)V_0}\left[2\alpha \alpha _0\sqrt{x}(x)^\alpha \right],$$
(65)
where $`\alpha `$ and $`\alpha _0`$ are two real constants. The energy density and pressure of the $`k`$-field are calculated from Eqs. (4) and (65)

$$\rho =(x)^\alpha ,p=\frac{1}{(2\alpha 1)}\left[2\alpha \alpha _0\sqrt{x}(x)^\alpha \right],$$
(66)
the equation of state is

$$p=\frac{1}{2\alpha 1}\left[\rho \frac{2\alpha \alpha _0}{\rho ^{1/2\alpha }}\right],$$
(67)
and the sound speed becomes

$$c_s^2=\frac{1}{2\alpha 1}\left[1\frac{\alpha _0}{\rho ^{(2\alpha 1)/2\alpha }}\right].$$
(68)
Solving the conservation equation (8), we obtain the energy density

$$\rho =\left[\alpha _0+\frac{c_0}{a^3}\right]^{2\alpha /(2\alpha 1)},$$
(69)
where $`c_0`$ is a redefinition of the integration constant $`c`$. Restricting to positive sound speed, the source (69) and the ECG induce a similar evolution of the scale factor. However, there is a significant change when both constant $`\alpha _0`$, $`c_0`$ are positive and $`\alpha >1/2`$, because in this case the baryotropic index becomes

$$\gamma =\frac{2\alpha }{2\alpha 1}\left[1\frac{\alpha _0}{\rho ^{(2\alpha 1)/2\alpha }}\right].$$
(70)
So, near the singularity, the energy density diverges and the baryotropic index behaves as $`\gamma 2\alpha /(2\alpha 1)`$, indicating that the model begins to evolve with a power-law dominated phase, where the scale factor is $`at^{(2\alpha 1)/3\alpha }`$. In this limit the sound speed (68) behaves as $`c_s^2(2\alpha 1)^1`$. For large $`\alpha `$, the model is initially dust dominated, with approximated vanishing sound speed, approaching to the ECG generated by the source (28) in the limit of large $`r`$. At late times the model ends in a de Sitter stage. Such ”modified Chaplygin gas” may be considered as an alternative model to the generalized Chaplygin gas investigated in . It allows the evolution of the initial perturbations in the energy density into a nonlinear regime to form a gravitational condensate of particles that could play the role of cold dark matter. The cool dark matter condensates gravitationally into the regions where the pressure $`p0`$ and the $`k`$-field is close to $`\varphi _c=(2\alpha \alpha _0)^{1/(2\alpha 1)}`$. In this case, the model yields an energy density which scales like the sum of a non-relativistic dust component at $`\varphi =\varphi _c`$ with equation of state $`p=0`$ and a cosmological-constant-like component $`p=\rho `$.

## IV Linking scalar and $`k`$-essence cosmologies

The observed acceleration of the present Universe has been investigated assuming that the dark energy can be described by quintessence and more recently by $`k`$-essence . The last one involves an effective scalar field theory generated by a Lagrangian with a non-canonical kinetic term. Particular cases of $`k`$-essence are generalized Chaplygin gas and tachyon dark energy models. Quintessence and $`k`$-essence frameworks are usually based in a homogeneous scalar fields driven by an exponential potential in the case of quintessence or an inverse square potential in the case of k-essence. Both encounter the so–called coincidence problem, namely, why are the energy densities of dark energy and dark matter of the same order today. Standard quintessence model appear promising in this point, as it can solve this problem for flat FRW universes provided the dark matter component is assumed to be dissipative . The system is attracted to a stationary and stable solution characterized by the constancy of both density parameters at late times. In addition, a class of k-essence models has been claimed to solve the coincidence problem by linking the onset of dark energy domination to the epoch of matter domination . From these satisfactory results and taking into account that both models involve evolving scalar fields, we believe it is reasonable to explore whether quintessence and $`k`$-essence frameworks have same kind of similitudes. We will deal with this question by investigating in which cases quintessence resembes k-essence.

Another interesting aspect to be considered is related with current observations which would indicate that the universe is superaccelerated and it could be filled with a non-standard fluid that violates the weak energy condition. In this sense, it was recently proposed a kind of matter described by an homogeneous scalar field with negative kinetic energy term. The fluid with negative pressure obeys and equation of state of the form $`p=(\gamma 1)\rho `$, where $`\gamma `$ is taken negative and the models are known as phantom or ghost cosmologies -. The phantom and quintessence cosmologies can be investigated simultaneously by taking a scalar field with both signs of the kinetic term driven by an exponential potential. The dynamical equations of these cosmological models are

$$3H^2=\frac{1}{2}q\dot{\phi }^2+𝒱$$
(71)
$$\ddot{\phi }+3H\dot{\phi }+\frac{1}{q}\frac{d𝒱}{d\phi }=0,$$
(72)
$$𝒱(\phi )=𝒱_0\text{e}^{qA\phi },$$
(73)
where $`q`$, $`A`$ are real numbers and $`𝒱_0`$ is a positive constant. The exponential potential is interesting because it may be considered just as a limit of a more complex potential . For negative values of $`q`$ the above equations describe a phantom cosmology .

It will be demonstrated that the scale factor obtained from Eqs. (71)-(73) or from the $`k`$-essence model generated by the linear function $`F=1+mx`$, and driven by an inverse square potential (see Eqs. (6), (33), (35) and (40)) is the same. To this end, we sketch the procedure followed in the Ref. . From Eqs. (71)-(73), we find

$$\dot{H}=\frac{1}{2}q\dot{\phi }^2,$$
(74)
and the first integral of the Klein-Gordon equation (72)

$$\dot{\phi }=AH+\frac{c_1}{a^3},$$
(75)
where $`c_1`$ is an arbitrary integration constant. Now, inserting Eq. (75) into (74), we obtain the second-order differential equation for the scale factor

$$\frac{d^2S}{d\zeta ^2}+S^\nu \frac{dS}{d\zeta }+\frac{1}{4}S^{2\nu +1}=0,\nu =6/qA^2,$$
(76)
where we have used the new variables $`S`$ and $`\zeta `$, defined by

$$S=a^{3/\nu },\zeta =c_1qAt.$$
(77)
With the following identification of the parameters

$$mV_0=\frac{2}{qA^2},\frac{3c}{2}=\frac{c_1}{A},$$
(78)
the Eqs. (41) and (76) coincide. Therefore, both models are described by the same scale factor and they are geometrically equivalent. Also, from Eqs. (6,7), (40) and Eqs. (71,72), we get a relationship between both potentials

$$3H^2+\dot{H}=V(t)=𝒱(t),$$
(79)
showing they are the same function of the cosmological time $`t`$, so that

$$\frac{V_0}{\varphi ^2}=𝒱_0\text{e}^{qA\phi }.$$
(80)
After inserting the $`k`$-field (52) into the last equation we find the scalar field

$$\phi =\frac{1}{qA}\mathrm{ln}\phi _0+\frac{\eta }{qA}+A\mathrm{ln}a,$$
(81)
where $`\phi _0=𝒱_0\varphi _0^2/V_0`$. The Eq. (80) supplies the link between the scalar field and the $`k`$-essence field

$$\varphi =\frac{\varphi _0}{\phi _0^{1/2}}\text{e}^{qA\phi /2},$$
(82)
displaying that these models are dynamically no equivalent. Thus, the homogeneous quintessence and phantom fields are different than k-essence field.

Resuming, from the cinematical point of view, that is the background geometry, we have obtained exact equivalences and it is impossible to differentiate between quintessence, k-essence and phantom cosmologies because they share the same scale factor. To distinguish between them, it appears necessary to focus our attention on the scalar field.

## V Conclusions

We have investigated the set of $`T_r`$-fields whose effective sound speed is proportional to the sound speed of the tachyon field. They can be grouped into three types according whether they yield phantom expansion or accelerated expansion with or without inflation. We have shown these behaviours by finding exact power-law solutions for an inverse square potential and proved that the $`T_{1/2}`$-field is compatible only with this potential. Each $`T_r`$-field produces an ECG and the set of all these gases can be divided into three kinds, one of which contains the generalized Chaplygin gas and the others give rise to ”perfect fluids”, leading to new evolutions. There exists basically nonsingular bouncing solutions which interpolate between two superaccelerated stages or singular ones with a finite time span, as well as, peculiar singular solutions that begin with a finite scale factor. In this manner, both, the ETF and the ECG may be considered fair candidates to implement phantom cosmologies.

For an inverse square potential, we have found the first integral of the $`k`$-field equation for any function $`F`$ and shown that the coupled Einstein-$`k`$-field equations can be solved in some cases. In particular, the divergent $`k`$-essence theory generated by the $`T_{1/2}`$-field becomes an intrinsic component of all $`k`$-essence models. Therefore, for power-law expansions the linear $`k`$-field model driven by an inverse square potential and the divergent model are isomorphic. We have obtained the general solution of Einstein-$`k`$-field equations for a linear function $`F`$. From the kinematical point of view this model and the quintessence scalar field one driven by an exponential potential are the same. However, they are dynamically nonequivalent, because the $`k`$-field and the scalar field are linked by the Einstein equation, e.g., both potentials are the same function of the cosmological time.

For a constant potential, we have studied a $`k`$-essence field associated with a perfect fluid whose equation of state contains a term proportional to the energy density and the other has the form of an ECG. This model, essentially different than the ECG model, smoothly interpolates between a power-law dominated phase and a de Sitter phase. In this ”modified Chaplygin gas” scenario, it may be chosen the value of the $`k`$-field where initial perturbations condensate gravitationally in cold dark matter.

## VI acknowledgments

This work was supported by the University of Buenos Aires under Project X223 and the Consejo Nacional de Investigaciones Científicas y Técnicas. LPC thanks Alejandro S. Jakubi for useful discussions of this work.


text:
# Implicit Purification for Temperature-Dependent Density Matrices

## Abstract

An implicit purification scheme is proposed for calculation of the temperature-dependent, grand canonical single-particle density matrix, given as a Fermi-Dirac operator expansion in terms of the Hamiltonian. The computational complexity is shown to scale with the logarithm of the polynomial order of the expansion, or equivalently, with the logarithm of the inverse temperature. The system of linear equations that arise in each implicit purification iteration is solved efficiently by a conjugate gradient solver. The scheme is particularly useful in connection with linear scaling electronic structure theory based on sparse matrix algebra. The efficiency of the implicit temperature expansion technique is analyzed and compared to some explicit purification methods for the zero temperature density matrix.

Linear scaling electronic structure theory in combination with tight-binding, self-consistent Hartree-Fock or density functional theory has become a very powerful tool for studying complex large material systems . There are several ways to achieve a computational cost that scales linearly with system size. Here we focus on methods based on the single-particle density matrix for band-gap materials, the matrix elements of which decay exponentially with overlap distance. For large systems the number of matrix elements above some numerical threshold therefore scales as $`𝒪(N)`$. In these schemes the two major steps are the construction of the tight-binding, Fockian or Kohn-Sham Hamiltonian $`H(𝐫,𝐫^{})`$ and the calculation of the density matrix $`\rho (𝐫,𝐫^{})`$. The present article concerns aspects of the second problem, the construction of the density matrix.

The relation between the density matrix at $`T=0`$ and the Hamiltonian is given by the Heaviside step function

$$\rho =\theta (\mu IH).$$
(1)
In density matrix schemes this relation is approximated by constructing $`\rho `$ from $`H`$ using sparse matrix algebra, where each major operation computationally scales linearly with the system size, thanks to $`𝒪(N)`$ matrix sparsity. This can be achieved through constrained minimization schemes , spectral projections or purification methods , or by an expansion of the temperature dependent Fermi-Dirac function or similar step function approximations . Contour integral representations of the Fermi distribution with complex Padé polynomials as resolvents, that can be calculated $`𝒪(N)`$ implicitly , as well as combinations of various approaches have also been explored .

The quadratically convergent purification techniques have turned out to be some of the most efficient approaches for the construction of the density matrix, both in memory and speed , with a computational complexity, in terms of number of matrix-matrix multiplications necessary to reach convergence, that scales linearly with the logarithm of the inverse band gap and the degree of expansion, and with a numerical error that scales linearly with the threshold . The Fermi-Dirac operator expansion-based methods based on Chebychev expansion techniques are generally much slower, with the computational cost scaling at best with the square root of the degree of expansion . However, these methods have an important advantage; they can account for a finite temperature distribution of the density matrix. Here we propose an expansion scheme that combines the low logarithmic complexity and quadratic convergence of the purification schemes with the finite temperature Fermi-Dirac distribution. We show how this can be accomplished by an implicit purification scheme, based on a Padé approximation of the rescaled Fermi-Dirac function, with a computational complexity that scales logarithmically with the expansion order, or equivalently, the inverse temperature.

The Fermi-Dirac distribution ,

$$\mathrm{\Phi }_{\mathrm{FD}}(\epsilon ,\mu ,\beta )=\frac{1}{e^{\beta (\epsilon \mu )}+1},$$
(2)
occurs in statistical mechanics as the occupational distribution of fermions at finite temperatures. It converges to a step function with the step formed at the chemical potential $`\mu `$ when $`T0`$. The temperature dependent grand canonical density matrix is formally given by the operator relation

$$\rho (\beta )=\mathrm{\Phi }_{\mathrm{FD}}(H,\mu ,\beta ).$$
(3)
The single-particle energy of a fermion system at a finite temperature is given by

$$E_s=Tr[H\rho (\beta )]=\underset{i,j}{}\varphi _i|H|\varphi _j\varphi _j|\rho (\beta )|\varphi _i,$$
(4)
in some set of basis functions $`\varphi _i`$. In this formulation the expression for $`\rho (\beta )`$ does not have to be calculated explicitly; instead the Fermi-Dirac function can be expanded in Chebychev polynomials. The Chebychev expansion technique is one of the most efficient ways of approximating a function and the Chebychev functions $`T_n(x)`$ obey a simple two-step recurrence formula, ($`T_0=1`$, $`T_1=x`$),

$$T_{n+1}(x)=2xT_n(x)T_{n1}(x).$$
(5)
The products $`\mathrm{\Phi }_{\mathrm{FD}}(H,\mu ,\beta )|\varphi _i`$ can be calculated efficiently using the two-step recurrence formula using only matrix-vector multiplications . The Chebychev expansion technique has many advantages: for example, the costly matrix-matrix multiplications are avoided, and error accumulation is small. However, compared to linear scaling purification techniques it is fairly inefficient . The problem is the slow linear increase in polynomial order as a function of iterations in the two-step recurrence formula. To improve the efficiency, Liang et al.  recently suggested an alternative approach, where the expansion polynomials are not calculated by the two-step formula, but by a direct expansion. This can be achieved with a computational complexity that scales with the square root of the polynomial order of the expansion $`𝒪(\sqrt{n})`$. This limit is optimal for a general polynomial . However, by choosing the expansion with a particular set of Padé polynomials, we will show how the computational complexity can be reduced even further, scaling only linearly with the logarithm of the polynomial order $`𝒪(\mathrm{log}n)`$ or the inverse temperature $`𝒪(\mathrm{log}\beta )`$.

There are at least 19 different ways to calculate matrix exponentials . Here we use one particular technique based on a Padé approximation. Consider the exponential function

$$e^x=\left(e^{x/n}\right)^n=\left(\frac{e^{x/(2n)}}{e^{x/(2n)}}\right)^n.$$
(6)
A Taylor expansion to first order gives

$$e^x=\underset{n\mathrm{}}{lim}\left(\frac{2n+x}{2nx}\right)^n.$$
(7)
This Padé approximation can be used in the Fermi-Dirac function and for the rescaled chemical potential and inverse temperature , $`\mu ^{}=1/2`$ and $`\beta ^{}=4n`$,

$$\mathrm{\Phi }_{\mathrm{FD}}(x,\frac{1}{2},4n)\frac{(1x)^n}{x^n+(1x)^n}.$$
(8)
At higher values of $`n`$ the approximation becomes increasingly better. The choice $`\mu ^{}=1/2`$ is made to center the step of the Fermi-Dirac function at $`x=1/2`$. In the interval $`[0,1]`$ the approximation is a continuously decreasing function with a maximum of $`1`$ at $`x=0`$ and a minimum of $`0`$ at $`x=1`$. This is the interval in which the temperature dependent density matrix has its eigenvalues and it is the interval where the Fermi-Dirac distribution is well approximated already at fairly high temperatures. It is therefore the interval around which we chose to perform the expansion. This choice requires an initial rescaling of the Hamiltonian spectra around the interval $`[0,1]`$. Let

$$G_n(x)=\frac{x^n}{x^n+(1x)^n}.$$
(9)
then the Fermi-Dirac function for $`\mu ^{}=1/2`$ and $`\beta ^{}=4n`$ in the interval $`[0,1]`$ is approximated by

$$\mathrm{\Phi }_{\mathrm{FD}}(x,\frac{1}{2},4n)=\left[e^{4n(x1/2)}+1\right]^11G_n(x).$$
(10)
The polynomial order $`n`$ of the Padé approximation is proportional to the inverse temperature since $`n=\beta ^{}/4`$. This means that the lower the temperature the better the approximation. In practice, however, the approximation at the normalized energy interval $`[0,1]`$ is already very good at orders as low as $`n5`$. An example given in Fig. 1, which shows the Padé approximation $`1G_5(x)`$, is virtually identical to the corresponding Fermi-Dirac function with $`\beta ^{}=20`$ and $`\mu ^{}=1/2`$. The inset shows the error. With the interval $`[0,1]`$ equal to 1 Ry this example corresponds to a temperature of 7894 K. At lower temperatures the approximation becomes increasingly better. The Padé approximation in Eq. (10) is only one alternative, but, as will be shown below, it turns out to be particularly simple and efficient.

A major advantage with the Padé approximation in Eq. (10) is how efficiently we can calculate high orders of $`G_n`$. The computational complexity is very low thanks to the iterative relation

$$G_{k\times l}(x)=G_k(G_l(x)).$$
(11)
In an operator expansion this corresponds to purifications, projecting the eigenvalues towards 0 and 1. In contrast to a more general polynomial expansion such as the Chebychev expansion, which computationally scales at best with the square root of the polynomial order, $`𝒪(\sqrt{n})`$ , or as $`𝒪(n)`$ if the two-step recurrence formula is used, the iterative relation above makes it possible to reach the same order of expansion in only $`𝒪(\mathrm{log}n)`$ steps. The same low logarithmic complexity is found generally in purification expansion schemes that are based on iterative spectral projections. The Padé approximation of the Fermi-Dirac distribution can thus be used in a highly efficient purification scheme which calculates the finite temperature density matrix $`\rho (\beta )`$ at a set of normalized inverse temperatures $`\beta ^{}=4n`$. The purification algorithm can be described by

$$\begin{array}{cc}X_1\hfill & =F(H,\mu )\hfill \\ X_{i+1}\hfill & =G_m(X_i),i=1,2,\mathrm{},\mathrm{log}_m(n)\hfill \\ \rho (\beta )\hfill & =IX_{i+1}.\hfill \end{array}$$
(12)
The expansion order $`n`$ and thus the normalized inverse temperature $`\beta ^{}`$ must be chosen so that the number of iterations $`\mathrm{log}_m(n)`$, is an integer. The function

$$F(H,\mu )=\alpha (H\mu I)+0.5I$$
(13)
is a normalization function that rescales all the eigenvalues of $`H`$ to the interval $`[0,1]`$, with the chemical potential $`\mu `$ shifted to $`\mu ^{}=1/2`$. The chemical potential and spectral bound must thus be known in advance. The normalization factor

$$\alpha \frac{1}{2}\mathrm{min}[(\mu H_{\mathrm{min}})^1,(H_{\mathrm{max}}\mu )^1],$$
(14)
rescales the spectra and sets the temperature scale. The implicit purification scheme converges to the zero temperature density matrix for any value of $`\alpha >0`$, but the convergence is faster and the approximation is more accurate at higher temperatures if $`\alpha `$ is chosen to normalize the spectra around $`[0,1]`$. The spectral bounds $`H_{\mathrm{max}}`$ and $`H_{\mathrm{min}}`$ can be estimated by for example Lanczos’ algorithm or Gersgorin circles. Generally we have that the temperature  $`T=1/(\alpha k_B4n)`$, where $`n`$ is the accumulated expansion order in Eq. (12).

Because of the rational form of $`G_m(X_i)`$ the purification scheme is implicit. Assuming a finite orthogonal basis representation, a set of linear equations in Eq. (12) has to be solved in each step for $`i=1,2\mathrm{},\mathrm{log}_m(n)`$:

$$\left[X_i^m+(IX_i)^m\right]X_{i+1}=X_i^m,$$
(15)
which is given from the second step in Eq. (12) and from the definition of $`G`$ in Eq. (9) along with its nested iterative expansion property given in Eq. (11). Here we find another major advantage with our particular choice of Padé approximation. The left side system matrix $`A_i=[X_i^m+(IX_i)^m]`$ is symmetric and positive definite for symmetric $`X_i`$’s with their spectra belonging to $`[0,1]`$. In fact, with increasing $`i`$, the system matrix $`A_i`$ converges to the identity matrix $`I`$. The implicit equations are therefore very well suited for solutions with the linear conjugate gradient method , that in turn, can efficiently exploit the close approximation of $`X_i`$ to the unknown columns of $`X_{i+1}`$, which becomes increasingly more accurate and efficient towards the last iterations. Another possibly efficient alternative is the application of the sparse approximate inverse (AINV)  that can be expected to work well for this particular problem. However, this approach has not been explored in the present study.

Alternative implicit purification schemes can also be derived from various sign matrix expansions . Sign matrix expansions are equivalent to purification. The only difference is that spectral projections are performed in the interval $`[1,1]`$ instead of $`[0,1]`$, as in the case of purification.

To analyze the efficiency of the algorithm compared with explicit purification schemes, we have chosen an $`N\times N`$ model Hamiltonian with $`N`$ random diagonal elements. The overlap elements decay exponentially as a function of site separation on a randomly distorted simple cubic lattice. This model represents a Hamiltonian of an insulator that might occur, for example, with a Gaussian basis set in density functional theory or in various tight-binding schemes. The convergence is mainly determined by the occupation and the band gap. The test Hamiltonian was therefore modified such that all $`N`$ eigenvalues were uniformly distributed on $`[0,1]`$. In this case the band gap $`\mathrm{\Delta }g=1/N`$, independently of the fractional occupation $`f_{\mathrm{occ}}=N_e/N`$. This simplifies the analysis and comparison of the different methods, which otherwise are hard to perform for a less idealized set of material systems. After each iteration a numerical threshold $`\tau =1.0\times 10^7`$ was applied and convergence was determined when the error in energy $`|E_{\mathrm{approx}}E_{\mathrm{exact}}|<1.0\times 10^5`$ . The convergence criterion corresponds in practice to $`T=0`$. A comparison at $`T0`$ is necessary since the explicit purification schemes used in the comparison only give the zero temperature density matrix. At room temperature the computational effort with the implicit purification scheme is only slightly reduced because of the rapid convergence. The computational complexity was measured in number of matrix-matrix multiplications, where $`N`$ conjugate gradient steps, i.e. $`N`$ matrix-vector multiplications, are counted as one matrix-matrix multiplication.

Figure 2 shows the computational cost for various occupation factors. The implicit purification scheme of order two (IP), i.e. with $`m=2`$ in Eq. (15), using $`X_i`$ as initial approximations to $`X_{i+1}`$, is compared to the trace correcting scheme with second order polynomials (TC2) by Niklasson , the trace resetting asymmetric fourth order method (TRS4) by Niklasson et al. , the grand canonical scheme with fourth order projections (GC4) by Niklasson , the grand canonical McWeeny purification scheme (McW) , and finally the canonical scheme (PM) by Palser and Manolopolous . The grand canonical schemes that require prior knowledge of the chemical potential are indicated in the figure by bold italics.

For small band gaps, i.e. high values of $`N`$, and with prior knowledge of $`\mu `$, the asymmetric GC4 method is the most efficient technique. The best performing schemes that require no prior knowledge of $`\mu `$ are the TC2 and TRS4 schemes. The TC2 scheme is more memory efficient since it only needs to calculate a second order polynomial in each iteration and intermediate storage needed in higher order expansions is avoided. However, it can not deal with degeneracy and fractional occupancy, which are addressed with the TRS4 scheme . At low occupation the PM scheme becomes very inefficient. This sensitivity is not seen for any of the other schemes. The proposed implicit purification scheme is slower than the alternative explicit purification schemes except for the PM scheme at low occupancies. However, it is the only method that, for only a slightly increased computational cost, correctly gives the temperature dependent Fermi-Dirac distribution of the single-particle eigenstates. The implicit purification scheme scales with the logarithm of the expansion order. This is an important improvement over previous Fermi-Dirac operator expansion methods. Whereas in general, a polynomial can be calculated with a computational cost scaling at best with the square root of the order of the polynomial , we restrict the polynomial approximation to a nested form $`f(f(\mathrm{}f(x)\mathrm{}))`$. In this case a high order can be reached much more efficiently than for the general form. This is the key idea behind purification expansions.

In summary, we have proposed an implicit purification scheme for the calculation of the temperature dependent single-particle density matrix given as a Fermi-Dirac operator expansion in terms of the Hamiltonian. The method is useful in connection with linear scaling electronic structure theory and it has a computational complexity that scales with the logarithm of the inverse temperature $`𝒪(\mathrm{log}\beta )`$ or as the logarithm of the polynomial expansion order $`𝒪(\mathrm{log}n)`$.

Discussions with Matt Challacombe, Eric Chisolm, Siobhan Corish, S. Goedecker, C. J. Tymczak, and John Wills are gratefully acknowledged.