text:
# Star formation and chemical evolution of damped Lyman 𝛼 systems

## 1 Introduction

Damped Lyman-$`\alpha `$ systems (DLAs) are high column density (the HI column density is higher than $`10^{20.3}`$/cm2) absorbers detected in the optical spectra of quasars up to relatively high redshift (up to $`z5`$), and can be used to probe the cosmic chemical evolution, i.e. the globle evolution of gas, metallicity and star formation rates in the Universe. DLAs are now generally believed to be the progenitors of present-day galactic discs (Wolfe et al. 1986; Wolfe 1988). This view is supported by the arguments based on mass estimates, which are from the column densities and absorber sizes at redshifts $`z23`$, as well as by the kinematic features consistent with models of fast-rotating (at velocities of order 200 $`\mathrm{km}\mathrm{s}^1`$), thick discs (Prochaska & Wolfe 1997a,b, 1998; Wolfe & Prochaska 1998). At the same time, on the basis of \[$`\alpha `$/Fe\] and \[N/O\] abundance ratios and their large scatter, an origin of DLAs in dwarf and low surface brightness galaxies is also presented (Matteucci, Molaro & Vladilo 1997; Vladilo 1998; Jimenez, Bowen & Matteucci 1999). Because DLAs contain most of the cool, neutral gas in the universe at high redshift, they play an important role in studies of the galaxy formation, gas process, star formation and metal enrichment as well as the large scale structure. The pioneering work in this field was done by Meyer, Welty & York (1989) and Pettini, Boksenberg & Hunstead (1990). In the past few years, Pettini et al. (1994, 1997a) have devoted great efforts in observation of DLAs’ column densities and metal abundances. They found that, with the increase of observed data, DLAs are generally metal poor with a mean value only of 1/13 $`Z_{}`$ and are very young galaxies in the early stage of chemical enrichment, and the metallicity distribution of DLAs shows a large spread of nearly two magnitudes. They also pointed out that the abundance of Cr is relatively poorer than that of Zn due to the depletion of dust.

Lu et al. (1996) used the powerful Keck 10 m telescope to study the metal abundances of many heavy elements, and presented a very detailed analysis of DLAs including comparing their distribution of metallicities, electron density and kinematics with those of our Milky Way galaxy, etc. Boissé et al. (1998) presented HST spectra for a sample of six DLAs with intermediate redshift ($`z_a1`$), and concluded that the available observations may be biased against dust-rich absorbers and dust extinction causes a preferential selection of QSOs with intervening gas relatively poor in metals, dust and molecules. They also strongly suspect that surveys of DLAs down to fainter QSOs would reveal more systems with large HI column densities and high metallicities. M$`ø`$ller $`\&`$ Warren (1998) used HST images to estimate the cross-section-weighted mean radius of DLA absorbers at high redshift ($`z>2`$), that is they are smaller than 23.6 $`h^1`$ kpc and 12.7 $`h^1`$ kpc for $`q_0=0.0`$ and 0.5, respectively. They also found that for $`q_0=0.5`$ the space density of DLAs at high redshift is more than five times that of the spiral galaxies locally, but the problem will be ruled out if $`q_0`$ decreases. The observations of H$`\alpha `$ emission done by Bunker et al. (1999) and near-$`IR`$ and $`IR`$ done by Bechtold et al. (1998) also supported the conclusions of M$`\mathrm{ø}`$ller & Warren (1998) above.

Untill now, most of the researches in DLAs are for observational constrains. There has, however, been relatively little work on its simple prescription. It is worthy of emphasizing that the observation of HST (M$`\mathrm{ø}`$ller & Warren 1998) strongly suggests DLAs and LBGs may commonly associated with each other. Based on their successful disc galaxy formation model, Mo, Mao & White (1999) have investigated the global properties of LBGs, such as the sizes, luminosities, kinematics and star formation rates. They presented that the number of density and clustering properties of LBGs are consistent with them being the central galaxies of the most massive dark haloes present at $`z3`$.

In this paper, under an assumption that the DLAs are the central galaxies of less massive dark haloes present at $`z3`$, and they should inhabit haloes of moderately low circular velocity, we study how star formation and chemical enrichment may have proceeded in them by using the simplest models possible with the smallest number of parameters.

The outlines of this paper are as follows. In Section 2, we give a simple description of our models; The predicted distributions of metallicity and neutral hydrogen column density for DLAs are presented in Section 3, and some correlations will also be shown in this section; The summary and discussion are presented in Section 4.

## 2 Modelling damped Lyman-$`\alpha `$ systems

We model the assembly of DLAs in the context of the standard hierarchical picture (White & Rees 1978; White & Frenk 1991). Structure growth in this model is specified by the parameters of the background cosmology and by the power spectrum of initial density fluctuations. As an illustration, we show theoretical results for a CDM model with cosmological density parameter $`\mathrm{\Omega }_0=0.3`$, cosmological constant $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$. The power spectrum is assumed to be that given in Bardeen et al. (1986), with shape parameter $`\mathrm{\Gamma }=0.2`$ and with normalization $`\sigma _8=1.0`$. We denote the mass fraction in baryons by $`f_\mathrm{B}=\mathrm{\Omega }_\mathrm{B}/\mathrm{\Omega }_0`$, where $`\mathrm{\Omega }_\mathrm{B}`$ is the cosmic baryonic density parameter. According to the cosmic nucleosynthesis, the currently favoured value of $`\mathrm{\Omega }_\mathrm{B}`$ is $`\mathrm{\Omega }_\mathrm{B}0.019h^2`$ (Burles & Tytler 1998), where $`h`$ is the present Hubble constant in units of 100 $`\mathrm{kms}^1\mathrm{Mpc}^1`$, and so $`f_\mathrm{B}0.063h^2`$. Whenever a numerical value of $`h`$ is needed, we take $`h=0.7`$.

### 2.1 Galaxy formation

There are some traditional methods to investigate the galaxy formation and evolution. One can, for instance, assume the redshift of galaxy formation and star formation history, based on the stellar IMF and cosmological model, to trace the observations (Tinsley 1980; Pozzetti, Bruzual $`\&`$ Zamorani 1996; Jimenez, Padoan & Matteucci 1998). White and his collaborators (White $`\&`$ Frenk 1991; Kauffmann $`\&`$ White 1993; Kauffmann, White $`\&`$ Guiderdoni 1993; Kauffmann et al. 1999) construct hierarchical models of galaxy formation and evolution, and are successful in many observations such as the statistical properties of local and distant galaxies, etc. Madau, Pozzetti $`\&`$ Dickinson (1998) developed a method which focuses on the integrated light of galaxy populations as a whole to match the properties of the star formation history for the field galaxies.

In the present paper, we use the Mo, Mao $`\&`$ White (1998) (hereafter, MMW) galaxy formation model in assuming that central galaxies form when collapse of the protogalactic gas is halted either by its angular momentum or by fragmentation as it becomes self-gravitating. This model is in hierarchical cosmogonies, in which a disc is a thin and centrifugally supported structure with an exponential surface density profile. LBGs can be well understood in MMW with large circular velocity $`V_c`$ and small spin parameter $`\lambda `$ which correspond to compact objects with massive haloes (Mo et al. 1999). In this paper, we assume that DLAs are the central galaxies of less massive dark haloes present at redshifts $`z3`$, and they should inhabit haloes of moderately low circular velocity.

According to MMW, if the mass profile of the disc is taken to be exponential

$$\mathrm{\Sigma }(r)=\mathrm{\Sigma }_0\mathrm{exp}(r/R_\mathrm{d}),$$
(1)
then $`\mathrm{\Sigma }_0`$ and $`R_\mathrm{d}`$ that are the central surface density and the scale length of the exponential disc, respectively, and the rotation curve of the galaxy, are determined uniquely. Specifically,

$$R_\mathrm{d}8.8h^1\mathrm{kpc}(\frac{\lambda }{0.05})(\frac{V_\mathrm{c}}{250\mathrm{km}\mathrm{s}^1})[\frac{H(z)}{H_0}]^1$$
(2)
and

$$\mathrm{\Sigma }_0380h\mathrm{M}_{}\mathrm{pc}^2(\frac{m_\mathrm{d}}{0.05})(\frac{\lambda }{0.05})^2(\frac{V_\mathrm{c}}{250\mathrm{km}\mathrm{s}^1})[\frac{H(z)}{H_0}],$$
(3)
where $`m_\mathrm{d}`$ is the fraction of disk to halo mass, $`V_\mathrm{c}`$ the circular velocity of the halo, $`\lambda `$ the dimensionless spin parameter, $`H(z)`$ the Hubble constant at redshift $`z`$, $`H_0`$ the present-day Hubble constant (See detailed in MMW). Because the Hubble constant $`H(z)`$ increases with redshift, it is expected from equations (2) and (3) that the galaxy disc should be less massive and smaller but have a higher surface density at higher redshift. One can find if $`V_\mathrm{c}`$, $`\lambda `$ and cosmogony are given, $`R_\mathrm{d}`$, which dominates the disk profile, will not change.

For a given cosmogony, the mass function of dark matter haloes at redshift $`z`$ can be estimated from the Press-Schechter formalism (Press $`\&`$ Schechter 1974):

$$dN_\mathrm{h}(M_\mathrm{h},z)=\sqrt{\frac{2}{\pi }}\frac{\overline{\rho }_0}{M_\mathrm{h}}\frac{\delta _\mathrm{c}(z)}{\mathrm{\Delta }(M_\mathrm{h})}\frac{\mathrm{d}\mathrm{ln}\mathrm{\Delta }(M_\mathrm{h})}{\mathrm{d}\mathrm{ln}M_\mathrm{h}}$$

$$\times \mathrm{exp}[\frac{\delta _\mathrm{c}^2(z)}{2\mathrm{\Delta }^2(M_\mathrm{h})}]\frac{dM_\mathrm{h}}{M_\mathrm{h}},$$
(4)
where $`\delta _\mathrm{c}(z)=\delta _\mathrm{c}(0)(1+z)g(0)/g(z)`$ is the linear overdensity corresponding to collapse at redshift $`z`$, and $`g(z)`$ is the linear growth factor relative to an Einstein-de Sitter universe at $`z`$, $`\delta _\mathrm{c}(0)1.686`$, $`\mathrm{\Delta }(M_\mathrm{h})`$ is the $`rms`$ linear mass fluctuation at $`z=0`$ in a sphere which on the average contains mass $`M_\mathrm{h}`$ by $`M_\mathrm{h}=(4\pi /3)\overline{\rho }_0R^3`$, $`\overline{\rho }_0`$ is the mean mass density of the universe at $`z=0`$. The relation between the mass $`M_\mathrm{h}`$ and circular velocity $`V_\mathrm{c}`$ of a halo is in equation (2) of MMW, and a more detailed description of the PS formalism and the related issues can be found in the Appendix of MMW.

Based on the PS formalism and the $`\lambda `$ distribution, with the latter being a log-normal function with mean $`\overline{\mathrm{ln}\lambda }=\mathrm{ln}0.05`$ and dispersion $`\sigma _{\mathrm{ln}\lambda }=0.5`$ (MMW), we can generate Monte Carlo samples of the halo distribution in the $`M_\mathrm{h}`$-$`\lambda `$ plane at redshifts $`z=3`$. We can then use the galaxy formation model of Section 2.1 to transform the halo population into a DLA population. Finally, suppose that each halo at $`z3`$ has a central galaxy with a star formation rate (SFR) that is a monotonic function of halo mass, and only a negligible fraction of haloes hosts a second galaxy bright enough to be seen. The observed DLAs then correspond to haloes at $`z3`$. Here a lower-limit for $`V_\mathrm{c}`$ of 50 km/s of haloes is chosen, and a random distribution of galactic discs is taken into account.

It should be pointed out that the interaction between discs and bulges is not considered in the present study. This interaction will become more and more important for very compact objects which correspond to $`\lambda <0.025`$ (MMW). Here, we treat galaxies with $`\lambda <0.025`$ as $`\lambda =0.025`$. It can be easily estimated from Warren et al. (1992), Cole $`\&`$ Lacey (1996) and Steinmetz $`\&`$ Bartelmann (1995) that, there are only 10 percent galaxies with $`\lambda <0.025`$. This will not influence our main results.

### 2.2 Star formation rates

For star formation rates (SFR), we take the empirical Schmidt law (Kennicutt 1998)

$$\mathrm{\Sigma }_{\mathrm{SFR}}=a(\frac{\mathrm{\Sigma }_\mathrm{g}(r)}{1\mathrm{M}_{}\mathrm{pc}^2})^b\mathrm{M}_{}\mathrm{yr}^1\mathrm{pc}^2,$$
(5)
where

$$a=2.5\times 10^{10},b=1.4,$$
(6)
respectively. Here $`\mathrm{\Sigma }_{\mathrm{SFR}}`$ is the SFR per unit area, and $`\mathrm{\Sigma }_\mathrm{g}(r)`$ the gas surface density. This star formation law was derived by averaging the cool gas density over large areas on spiral discs and over starburst regions.

### 2.3 Closed box model

For a closed box model, the metallicity $`Z`$ of gas is (Binney & Merrifield 1997; Tinsley 1980)

$$ZZ_i=y\mathrm{ln}(\mathrm{\Sigma }_\mathrm{g}(r)/\mathrm{\Sigma }_{\mathrm{g0}}(r))=y\mathrm{ln}\mu ^1,$$
(7)
where $`Z_i`$ is the initial metallicity of cool gas that is taken to be 0.01$`Z_{}`$ in the present paper, $`y`$ the stellar yield that is obtained by normalizing the predicted distribution of metallicity to the mean value of 1/13 $`Z_{}`$ as presented by Pettini et al. (1997a) in this study (see next section), $`\mathrm{\Sigma }_{\mathrm{g0}}(r)`$ the initial gas surface density and $`\mu `$ the gas fraction, and

$$\frac{\mathrm{d}\mathrm{\Sigma }_\mathrm{g}(r)}{\mathrm{d}t}=(1R_\mathrm{r})\mathrm{\Sigma }_{\mathrm{SFR}},$$
(8)
where $`R_\mathrm{r}`$ is the returned fraction of stellar mass into the ISM, and we take $`R_\mathrm{r}=0.3`$ for a Salpeter IMF (see Madau et al. 1998). We assume that at initial time ($`t_{}=0`$), $`\mathrm{\Sigma }_{\mathrm{SFR}}=0`$ and $`\mathrm{\Sigma }_{\mathrm{g0}}(r)=\mathrm{\Sigma }_0(r)\mathrm{exp}(r/R_\mathrm{d})`$, with the latter being given by equation (1). After a period of $`T`$, the surface density of the survived cool gas can be obtained by equations (5) and (8), that is

$$\mathrm{\Sigma }_\mathrm{g}(r)=\mathrm{\Sigma }_{\mathrm{g0}}(r)(1+T/\tau )^{\frac{1}{b1}},$$
(9)
and the current SFR is

$$\mathrm{\Sigma }_{\mathrm{SFR}}=a\mathrm{\Sigma }_{\mathrm{g0}}^b(r)(1+T/\tau )^{\frac{b}{1b}},$$
(10)
where

$$\tau =\frac{\mathrm{\Sigma }_{g0}^{1b}(r)}{a(b1)(1R_\mathrm{r})}.$$
(11)
As we know, in the simple closed-box model of chemical evolution, all metals produced by stars are retained within the galaxy. The metallicity of the galaxy is simply determined by the stellar yield $`y`$, defined as the mass of metals produced per solar mass of long-lived stars that are formed (Tinsley 1980). It is worthy to note that, the instantaneous recycling approximation is not valid as regards the contributions to metal enrichment of either type Ia supernovae or of stars of a few solar masses (Binney & Merrifield, 1997). Pagel (1997) has refined to take into account the time delay between star formation and the ejection into the ISM of heavy elements.

In this paper, we use \[X/H\]111We use the normal notation where $`[X/Y]=\mathrm{log}(X/Y)\mathrm{log}(X/Y)_{}`$. for the logarithm metallicity, and all the histograms are normalized to unity.

### 2.4 Criterion for star formation

According to Toomre (1964), in the case of a gaseous disc, instability for star formation is expected if the surface density of gas exceeds the following critical value

$$\mathrm{\Sigma }_\mathrm{c}(r)=\frac{\alpha \kappa \sigma _\mathrm{g}}{3.36G},$$
(12)
where $`\alpha `$ is taken to be 2 (van der Kruit & de Grijs, 1999); and

$$\kappa =1.41\frac{V}{r}(1+\frac{r}{V}\frac{\mathrm{d}V}{\mathrm{d}r})^{1/2};$$
(13)
and $`\sigma _\mathrm{g}`$ is a gas velocity dispersion, we take it to be $`6\mathrm{km}\mathrm{s}^1`$ (Kauffmann 1996), which is consistent with the observational results at any redshift till now (see Larson 1998).

When this criterion is considered, the star formation occurs only in the region where the surface density of gas satisfies the critical value, not everywhere of the gas disc.

## 3 Results

In this section, we present the predicted normalized distributions of metallicity and HI column density for the DLA population. As explained by Pettini, King & Smith (1997b), \[Zn/H\] is a straight forward measure of the degree of metal enrichment analogous to the stellar \[Fe/H\]. So, we use the mean value of \[Zn/H\] as presented by Pettini et al. (1997a) to normalize the predicted distribution of metallicity for DLAs in order to obtain the steller yield. A random distribution of discs is taken into account.

### 3.1 Distribution of metallicity without considering the Toomre’s criterion

At first, we do not consider the criterion for star formation which is developed by Toomre (1964). Based on the Monte Carlo samples at $`z=3`$, we can calculate the metallicity $`Z`$ of DLA population by the formulas from (5) to (11) by assuming that, at initial time ($`t_{}=0`$), $`\mathrm{\Sigma }_{\mathrm{SFR}}=0`$ and $`\mathrm{\Sigma }_{\mathrm{g0}}(r)=\mathrm{\Sigma }_0\mathrm{exp}(r/R_\mathrm{d})`$; and the star formation time $`T`$ is 1 Gyr. $`\mathrm{\Sigma }_0`$ and $`R_\mathrm{d}`$ are decided by equations (2) and (3). The stellar yield is obtained by normalizing the predicted distribution of metallicity to the mean value of 1/13 $`Z_{}`$ as presented by Pettini (1997a). It is 0.22 $`Z_{}`$. The results are plotted in Fig. 1. For comparison, the metallicities of 18 DLAs at redshifts $`z3`$ $`(2.3<z<3.4)`$ that were obtained by Meyer, Welty & York (1989) and Pettini et al. (1994, 1997a) 222The solar system abundance of Zn is from the compilation by Anders & Grevesse (1989), i.e. $`\mathrm{log}(\mathrm{Zn}/\mathrm{H})_{}=7.35`$. are also presented. This figure shows that, the observational data that have a spread of nearly two orders of magnitude agree with our model, but the predicted metallicity distribution cannot spread as widely as the observations. The observational data are still quite sparse, and larger and more complete samples are needed. In order to test whether our model fits the abundances of DLAs, we provide the Kolmogorov-Smirnov Test. The results are that the maximum value of the absolute difference is 0.33, and the significance level probability is 69.94 percent.

### 3.2 Distribution of metallicity with considering the Toomre’s criterion

When we take into account the Toomre instability criterion of star formation, the star formation occurs only in the region where the surface density of gas must be larger than $`\mathrm{\Sigma }_c(r)`$. So, we first obtain the regions where the star formation can occur based on the formulas (1), (12) and (13), assuming that the disc has a flat rotation curve and rotational velocity is equal to the circular velocity of the surrounding dark matter halo. If the observation is not in these regions, the predicted distribution of metallicity for DLA population is initial one of the cool gas, i.e. 0.01 $`Z_{}`$. When the observation is just in these regions, the metallicity of the DLA population can be obtained as section 3.1. In this case, the stellar yield is 0.70 $`Z_{}`$. Fig. 2 presents the normalized model distribution of metallicity. For comparison, the metallicities of 18 DLAs at redshifts $`z3`$ $`(2.3<z<3.4)`$ that were obtained by Meyer, Welty & York (1989) and Pettini et al. (1994, 1997a) are also plotted. This figure shows that, the model data which spreads more than two magnitudes can match the observed ones (Meyer et al. 1989, Pettini et al. 1997a), although the model is a little metal rich. As in Fig. 1, we also provide the Kolmogorov-Smirnov Test. The results are that the maximum value of the absolute difference is 0.25, and the significance level probability is 84.75 percent.

### 3.3 Distribution of HI column density

HI column density is the major characteristics of DLAs, because it is the first physical parameter that we can obtain from observations of QSO absorption lines. If we assume that the survived gas from the star formation in a disc is dominated by the neutral hydrogen, the model column density of neutral hydrogen can be obtained by the formulas from (5) to (11). The predicted distribution of HI column density for DLAs without and with considering the Toomre’s criterion are shown in Figs. 3 and 4, respectively. For comparison, the HI column densities of 18 DLAs at redshifts $`z3`$ $`(2.3<z<3.4)`$ that were obtained by Meyer, Welty & York (1989) and Pettini et al. (1994, 1997a) are also plotted. These two figures present that, the observational neutral hydrogen column densities spread narrower than the predicted ones. The reason is, perhaps, that the extinction may results in a preferential selection of QSOs with systems displaying a low HI and H2 content.

In fact, Boissé (1995) has already noticed that, the absence of very large N(HI) values (e.g. $`10^{22}`$ cm-2) which would be quite easily recognized in the numerous low resolution optical spectra obtained so far, supports the reality of cut-off caused by extinction.

Surveys of DLAs samples drawn from the observation of fainter QSOs may reveal more systems with large N(HI), and spectroscopic surveys of fainter QSOs can lead to a more representative view of the whole diversity of the DLAs.

### 3.4 Relation between metallicity and HI column density

In order to present the predicted relation between metallicity and HI column density for DLAs, we randomly select 100 samples from the DLA population by their cross-sections because the observations of DLAs are based on the QSO absorption lines.

In Figs. 5 and 6, we plot the predicted distribution of \[Zn/H\] as a function of HI column density for DLAs without and with considering the Toomre’s criterion, respectively. For comparison, the observational data of 18 DLAs at redshifts $`z3`$ $`(2.3<z<3.4)`$ that were obtained by Meyer, Welty & York (1989) and Pettini et al. (1994, 1997a) are also plotted. The predicted correlations (especially, without considering the Toomre’s criterion) are that, the higher the column density, richer the metallicity. In fact, Boissé et al. (1998) have already emphasized that, in an unbiased sample, the large N(HI) values should correspond to the innermost parts of galaxies where the metallicity is presumable higher. The observational data show a different trend, i.e. the metallicity is richer with the decreasing of column density because of the extinction bias as pointed out by Boissé et al. (1998).

## 4 Summary and Discussion

In the present paper, we establish two simple scenarios to investigate the nature of DLAs at redshifts $`z3`$ based on the disc galaxy formation model developed by Mo, Mao & White (1998). One is without considering the Toomre’s criterion, another is with considering the Toomre’s criterion. The star formation rate is chosen to be in the form of the empirical Schmidt law (Kennicutt 1998). The approximation of instantaneous recycling is assumed. The predicted distribution of metallicity for DLAs spreads larger than one order of magnitude, which is consistent with the observations by Meyer, Welty & York (1989) and Pettini et al. (1994, 1997a).

Although the star formation time of 1 Gyr, which we adopted, is arbitrary, its order is consistent with that of observations (Bechtold et al. 1998). It is no problem for us to choose it as a typical time of star formation to investigate the global properties of DLAs.

In our models, the star formation time scale in the disc is implicitly assumed to be the same everywhere. This will lead to the slow consumption of gas in the central parts of the disc. In fact, as we know well, the star formation in the central region is very rapid, and the central region should correspond to higher metallicities and lower HI column densities. Furthermore, we treat all the survived gases as HI. There actually exist $`\mathrm{H}_2`$, ionized gases and other elements. This is also leads to a overestimation of HI column density, especially in the central region of the disc. In addition, the DLAs that correspond to have small $`V_\mathrm{c}`$, should form earlier and the merging could take place more frequently. So, the DLAs are to a certain extent more complicated.

There is only one output parameter in our models, i.e. the stellar yield, which relates to the star formation time scale and is obtained from the normalization of the predicted distribution of metallicity to the mean value of 1/13 $`Z_{}`$ as presented by Pettini et al. (1997a). The metallicity of the initial cool gas is assumed to be of 0.01 $`Z_{}`$. According to Pagel (1987), the stellar yield is 0.3 $`Z_{}`$ for the disc clusters and 0.4 $`Z_{}`$ for the solar neighborhood. The stellar yields in our models are also quite similar to that of LMC of some 0.2 $`Z_{}`$ (Binney $`\&`$ Tremaine 1987).

Finally, we must point out, although our models are very simple without many physical processes being considered, such as the interaction between discs and bulges, merging effects, gas reheating and supernovae feedback, etc., it can predict some properties of DLAs.

## Acknowledgments

We are indebted to the anonymous referee for many critical comments and helpful suggestions that have greatly improved our paper. We are also grateful to H. J. Mo for much help in finishing this paper. This project is partly supported by Chinese National Natural Foundation and the “Sonderforschungsbereich 375-95 für Astro-Teilchenphysik” der Deutschen Forschungsgemeinschaft.


format: text

text:
# Random sequential adsorption and diffusion of dimers and 𝑘-mers on a square lattice

## I Introduction

Random Sequential Adsorption (RSA) of hard particles on lattice is employed to describe a wide class of irreversible phenomena in physics, chemistry and biology.  The prototype of this process is the dimerization of monomers on a chain.  Each monomer occupies initially a site of the chain and can react irreversibly with one of its neighbors to form a dimer: irreversible dimerization takes place sequentially until no more pairs of neighbor monomers are available. According to Flory  this process can be viewed as the random sequential adsorption of non-overlapping dimers onto a one dimensional lattice: a dimer can be irreversibly deposited if two neighbor free sites are available. At a given stage only single isolated sites will be available and further deposition will be no more allowed. At that stage, the density of dimers is finite and takes a non random value in the limit of infinite system, the so called jamming coverage $`\theta _J`$. This study was followed by the introduction and investigation of a variety of RSA processes involving different objects, lattices and rules.  The quantities of interest are in general the time behavior of the coverage $`\theta (t)`$ and its jamming value. Exact solutions are available in few cases, typically for dimers in one dimension, whereas only numerical or approximate results exist for more complex objects and in higher dimension.

In recent years, the class of systems for which RSA yields effective description has widened by the observation that in many experimental processes particles also diffuse on the lattice after deposition. Diffusional Relaxation (DR) affects the spatial distribution of free sites, allowing for further deposition of particles whenever neighbor empty sites are created. This produces a different time behavior of $`\theta (t)`$ and different values of the asymptotic coverage $`\theta (\mathrm{})`$. The effect of DR on RSA models has been investigated by many authors (see  and Refs. therein), often displaying the emergency of non trivial behaviors both for $`\theta (t)`$ and $`\theta (\mathrm{})`$. In the case of dimers on a line, for instance, the introduction of DR produces a very slow increase of coverage with time, terminating eventually with the complete filling: $`\theta (\mathrm{})=1`$. Exact results  and numerical simulations  show that, after a transient depending on the initial conditions, $`1\theta (t)`$ decreases asymptotically as $`t^{1/2}`$, and thus an infinite time is needed to reach the completely filled state. The exponent $`1/2`$ is due to the relevance of statistical fluctuations in the formation of neighbors free sites, so that any computation based on mean field treatments is not effective. It is also for this reason that solutions of RSA problems with Diffusion (RSAD) are very rare and usually limited to one dimensional systems (see  and Refs. therein) and some particular cases.  Some two dimensional models have been recently investigated , and features sometimes different from the one dimensional case have been observed. Eisenberg and Baram have found a stable asymptotic coverage $`\theta (\mathrm{})<1`$ for the RSAD of crosses on a square lattice  whereas such a behavior is not observed in the case of squares , in which the complete coverage is eventually reached for an infinite system, although through coexistence of different phases. The pure RSA of dimers in two dimensions was first investigated by Nord and Evans  and was later extended also considering anisotropic deposition rates.  Grigera et al.  computed the asymptotic coverage $`\theta (\mathrm{})`$ in presence of diffusion and rotation of dimers and found that it is less than unity. They also observed indications for an exponential approach to the limit coverage $`\theta (\mathrm{})`$.

In this paper we present extensive Monte Carlo (MC) simulations of the RSAD of linear particles on a square lattice. In the systems under consideration hard particles, the $`k`$mers, are dropped on an initially empty square lattice and are allowed to diffuse. Each $`k`$mer occupies $`k`$ consecutive sites and cannot overlap to the others. We consider different initial conditions as well as a set of different deposition and diffusion rates and probabilities. Whereas the former only have influence on initial and intermediate stages of the coverage, the latter influence the late stage of dynamics and different values may lead to different final states. In fact, we show that the value of the final coverage depends on these parameter and, as expected, decreases with increasing $`k`$. In this paper the focus will be mainly on the case $`k=2`$, since simulations show that its main features are also shared by systems with higher values of $`k`$, while the time required for simulations rapidly increases with $`k`$. For this reason Secs. II-IV are dedicated mainly to the RSAD of dimers.

In Sec. II we define the model and the procedure for the Montecarlo simulation. We also test possible algorithm to save computational time when simulating such slow relaxation processes, and describe the one we adopted. Section III is devoted to investigate the dynamical features of the coverage process as well as of the final state and of their dependence on the deposition and diffusion rates. In Sec. IV the results are put in relation with the behaviour of the vacancies in the lattice, whose statistical behaviour is also studied. The $`k`$mers RSAD model with $`2<k5`$ is addressed in Sec. V, where similarities and differences with the case $`k=2`$ are pointed out. Finally sec. VI is devoted to a summary of our main findings and to some general conclusions.

## II Model and Computational Details

### A RSAD dynamics

The process that we have more extensively investigated in this paper consists in the deposition and diffusion of dimers on a square lattice containing $`N=L\times L`$ sites with periodic boundary conditions. Each dimer occupies two nearest neighbor sites, or equivalently a bond, and cannot overlap with other dimers. For $`t>0`$ it is free to diffuse along the direction of its axis towards empty sites. We have assumed deposition to occur with equal probability on horizontal and vertical bonds and we have considered two possible different initial conditions: $`a)`$ jamming coverage, i.e. $`\theta (0)=\theta _J`$ and $`b)`$ empty lattice, i.e. $`\theta (0)=0`$. In the case $`a)`$ the initial stage consists of pure RSA. Dimers are deposited on the lattice until no more free pair of sites is available and the jamming state is reached. The jamming coverage for this kind of model is known from previous work  to be $`\theta _J=0.90654`$. It is reached exponentially and is higher than the corresponding one in one dimension. At this stage deposited dimers are allowed to diffuse and more dimers can be possibly added when suitable space is made available. In case $`b)`$ the deposited dimers can move from the beginning, so that diffusion in the initial stages largely dominates the deposition process, giving rise to a large amount of vacancies.

We have computed the time behavior of the coverage and of the abundance of different kinds of vacancies in the above processes by means of MC simulations. At each time step deposition of dimers is attempted with probability $`p`$ and diffusion with probability $`1p`$. Both directions of diffusion are equiprobable. Deposition happens at rate $`R`$, i.e. $`pRN`$ dimers are tried to be deposited on the average in a MC step and a unitary jumping rate $`d=1`$ follows for diffusion by assuming unitary lattice spacing and simulation (micro) time discretization $`\delta t=1/N`$.

### B Saving time algorithms

The above dynamics is extremely time consuming to simulate, mainly in the late stage of the process, since due to the exclusion rules most of diffusion and deposition attempts fail. A large amount of computational time can however be saved and dynamics sped up by suitable algorithms that maximize the number of successful events. We have adopted the strategy of storing in memory the dimers that can successfully diffuse along a given direction, $`n_{difH}`$ and $`n_{difV}`$ and the bonds available for deposition $`n_{depH}`$ and $`n_{depV}`$ respectively for the horizontal and vertical directions.

By accounting for the different probabilities and rates at which these actions may take place, it turns out that on the average a successful event occurs with probability

$$p_{succ}=p\left(\frac{n_{depH}+n_{depV}}{2N}\right)+\frac{1p}{2}\left(\frac{n_{difH}+n_{difV}}{2N}\right).$$
(1)
Since the law of large numbers is expected to hold (at least for system size large enough), at each step of the simulation the MC time is increased by the amount

$$\mathrm{\Delta }t=\frac{1}{p_{succ}}.$$
(2)
We adopted the above method that produces a time saving of about one order of magnitude with respect to the usual MC simulation with equivalent results.

As a further check we have also tested the method adopted by Gan and Wang , which is based on the probability distribution of time lags $`\mathrm{\Delta }t`$ between two consecutive successful events

$$P(\mathrm{\Delta }t=m)=p_{succ}(1p_{succ})^{m1}$$
(3)
and also yields results equivalent to those of usual MC. The identity of results confirms the reliability of both of these methods; however time saving is much more efficient via Eq. (1), since in fact this latter method corresponds to an average of the former over a large number of successful events. Both of them require additional memory occupation for the storage of dimers and bonds that may change their state, but this does not seriously hinder simulations of large systems if compared with required computational time.

## III Time coverage

### A Asymptotic coverage

Mean field approximations are known to work very badly for RSAD of small particles in one dimension. Diffusional dynamics of dimers along the lattice can be described in terms of mobile vacancies, that annihilate through the deposition of a dimer when they meet on neighbor sites. As already mentioned, both exact solutions and numerical simulations show that the process of vacancies annihilation is dominated by random fluctuations of density rather than by average values, thus leading asymptotically to a coverage described by a random walk process, $`\theta (\mathrm{})\theta (t)t^{1/2}`$, opposed to the $`t^1`$ coverage estimated by mean field. Random walk process with regular exponent $`0.5`$ concerns the dynamics of vacancies, since particles suffer an anomalous random walk with exponent $`0.5`$, as discussed in Ref. .

Mean field is generally speaking expected to work better in higher dimensions. However studies on RSAD of different objects in $`d=2`$ found that fluctuations of isolated vacancies still dominate relaxational dynamics.  Investigations of Eisenberg and Baram  on RSAD of monomers particles with nearest neighbors and next nearest neighbors exclusion have shown that $`\theta (\mathrm{})\theta (t)t^{1/2}`$ to the leading order in the asymptotic regime. Systems made by superposition of several one dimensional lattices to form multilayer structures, that allow for both intra-layer and inter-layer diffusion, have been found to behave in a similar fashion.

Our simulations show that also the RSAD of dimers in $`d=2`$ displays a regime in which $`t^{1/2}`$ dominates, thus implying the recombination of one-site vacancies to be the main process, in analogy with the $`d=1`$ case (we will turn to this point in Sec. IV). This result has shown to be independent from the parameter values and initial conditions chosen for the simulation. Nevertheless in $`d=2`$ there is a stage of the covering process at which the dynamics definitely changes. At that stage we observe that the annihilation rate slows down and correspondingly the filling slowly stops. The system is prevented to reach the complete coverage by the formation of isolated vacancies that are trapped at one site or along closed loops. As firstly pointed out by Grigera et al.  this seems to be a genuine property of the system, since independent of the lattice size. The results displayed in Fig. 1, where systems of different sizes approach the same asymptotic coverage $`\theta (\mathrm{})0.998`$, support the idea of a dynamically jammed state related to local configurations of the system. We find, in agreement with the authors in , that finite-size effects are reasonably negligible for $`L40`$.

In Fig. 2 we show the coverage behavior for different choices of $`p`$. Different values of $`\theta (\mathrm{})`$ are reached, implying that different choices of $`p`$ and $`R`$ can lead to different final coverages and can therefore be regarded as different quench rates of an initially annealed system.

### B System evolution

Fig. 2 puts into evidence a somewhat counterintuitive feature, i.e. that faster coverages are obtained at lower deposition probabilities, showing that to this aim the process of vacancies generation is more important than deposition attempts themselves. In order to understand this behavior it is worth to consider the different processes that contribute to increase the system coverage. This can be done by describing the dynamics of the system on a wide time scale by a series of the form

$$1\theta (t)=\underset{i=0}{\overset{\mathrm{}}{}}r_it^{i/2}$$
(4)
A fit of this curve to the results of numerical simulations shows that the terms beyond the fourth in the series can be neglected (cfr. Table I) and is displayed in Fig. 2. As expected, the effect of an increased diffusion is to decorrelate the system, and thus to produce a faster coverage, as shown by higher weights for the $`t^1`$ term. The coverage depends on the initial condition only at early times.

Fig. 3 displays the time coverage when $`\theta (0)=0`$. After an initial stage ($`1<t<10^2`$) in which deposition attempts dominate and the density progressively increases, the system approaches the same regime as for $`\theta (0)=\theta _J`$.

A comparison of the two cases is made in Fig. 4. This time behavior can also be well described by curves of the form of Eq. (4), obtaining the coefficients of Table II. Coefficients relevant to the late stage process exhibit a similar trend with respect to those of Table I.

Another interesting situation, also often investigated in literature, is the Diffusion Limited Deposition (DLD). In this case $`p=1`$ and $`R=\mathrm{}`$, corresponding to the filling of two-site vacancies with dimers as soon as they are created by diffusion. In this case only one-site vacancies exist, making the coverage process much more effective with respect to the case of finite $`R`$, as evidenced in Fig. 5.

On the other hand the asymptotic coverage is higher for finite deposition rate, showing that the larger diffusion allowed becomes important at very high coverage values, since in that case the long time frozen configurations are more easily removable. Thus, in analogy with Figs. 2 and 3, different asymptotic coverages may be reached. Also for DLD a superposition of non analytic powers of time fits the curve quite well in all the dynamical range:

$$1\theta (t)b_0+b_1t^{1/2}+b_2t^1+b_3t^{3/2},$$
(5)
where $`b_00.0016`$, $`b_10.060`$, $`b_20.112`$ and $`b_30.137`$ as shown in Fig. 5.

## IV Vacancy dynamics

As mentioned above, RSAD can be described in terms of reaction and diffusion of vacancies on the lattice. In the general case of $`k`$mers with diffusion coefficient $`D`$ and deposition rate $`R`$, one deposition corresponds to the annihilation of $`k`$ vacancies on contiguous sites.

If one adopts this point of view the system behavior appears to be ruled by the dynamics of vacancies. In order to give a full account of this dynamics, knowledge of vacancies correlation functions is necessary. Nevertheless a qualitative understanding of some of the main dynamical features can be obtained also through the investigation of the sole vacancy population.

### A Initial jammed state

Fig. 6 compares the coverage behavior with the change in time of vacancy population for different values of $`p`$. Solid and dotted lines represent respectively the total fraction of empty sites, $`1\theta `$, and of one-site vacancies, $`\rho _1`$. At early and intermediate times these two quantities display some differences, more visible at lower values of $`p`$ where formation of many-site vacancies is enhanced by diffusion. For a jammed initial state, an increase in two-site vacancies takes place in the early stage of the process and, for low values of $`p`$, can last for several decades both in density ($`\rho _2`$) and time. In fact, two dimensional geometry allows for the formation of many different kinds of voids; besides vertical and horizontal two-site vacancies, there can be $`L`$ many-site vacancies with different orientations, square vacancies and so on.

A comparison of $`\rho _2`$ with the overall non-one-site vacancy density, $`1\rho _1\theta `$, is shown in Fig. 6. We observe that the abundance of the latter does not vanish for large times, reflecting the configurational richness of the two-dimensional process with respect to the unidimensional one where $`\rho _2`$ and $`1\rho _1\theta `$ coincide asymptotically.

If one expands $`\rho _2`$ and $`1\rho _1\theta `$ in powers of $`t^{1/2}`$ one finds, as for the coverage, that larger contributions from mean field terms appear at lower $`p`$ values in the asymptotic expansions:

$$\rho _2(t)=\underset{i=0}{\overset{\mathrm{}}{}}x_it^{i/2},$$
(6)
$$1\rho _1(t)\theta (t)=\underset{i=0}{\overset{\mathrm{}}{}}y_it^{i/2}.$$
(7)
The coefficients of the best fit of Eqs. (6) and (7) to the numerically determined curves for $`\rho _2`$ and $`1\rho _1\theta `$ are reported respectively in Table III and Table IV.

### B Initial empty lattice

An initially empty lattice produces early and intermediate stages in vacancy dynamics different from the jammed case. Whereas the late stage is still dominated by recombination of one-site vacancies, an initial stage exists in which the abundance of this kind of defects increases monotonically up to reach about the same density of many-site vacancies, as displayed in Fig. 7. On the other hand, right at that stage the latter’s population starts to be less than that of one-site vacancies. Interestingly single empty sites rule the dynamics for $`t>10^2`$, that is, as already explained, the same time range in which density no longer depends on initial conditions (cfr. Fig. 4). It can be noted that two-site vacancies grow at pace with one-site ones up to a maximum, where they start to decay at the same rate of larger ones. The location of this maximum as well as the intersection between $`2\rho _2`$ and $`1\rho _1\theta `$ tend to move to early times when $`p`$ increases, suggesting that enhanced deposition rate destroys many-site vacancies more efficiently. The functional forms that fit well $`\rho _2`$ and $`1\rho _1\theta `$ are the same as those for the initial jammed lattice case (cfr. Eqs. (6) and (7)).

## V RSAD of $`k`$mers

We now address the problem of RSAD of $`k`$mers on a square lattice, which is a generalization of the model described in Sec. II that corresponds to $`k=2`$. Since for increasing $`k`$ the simulations of the dynamics becomes more and more time consuming and finite-size effects more important, we concentrate on low values of $`k`$, namely $`k=3,4,5`$.

We have simulated the RSAD dynamics always starting with an initial jammed state and then switching on the diffusive part of the process. We checked that the jamming densities obtained as initial conditions after the pure RSA process were in agreement with the values found earlier by Bonnier et al It turns out that there were no discrepancies with their results within the numerical precision. The jamming densities obtained are summarized in the Table V.

In Fig. 8 we compare the temporal evolution of the $`k`$mers density for two values of $`p`$.

The major aspect to underline is the very sluggish dynamical behavior for increasing $`k`$, which is a common feature with the one dimensional case: the probability for $`k`$ vacancies to meet and annihilate decreases when $`k`$ gets larger. It can be noted that, as for dimers, a dynamical saturation at long times exists, although our data are less clear for $`k=5`$, since the time at which this plateau sets in is an increasing function of $`k`$. Such a saturation is absent in $`d=1`$, where the density follows a single power law for long times also for larger values of $`k`$. It appears therefore to be a feature typical of the two dimensional geometry and independent of the particle size, related to the complexity of the configurational space accessible to the system. A proof of that is that the curves in Fig. 8 can be fitted by the same functional form adopted for dimers (Eq. (4)) with the parameters of Table VI. It is seen that the $`t^{1/2}`$ term is more important at larger $`k`$, in agreement with a slow down in the coverage process.

Arguments based on diffusive reactions  have shown in the past that a critical dimension $`d_c=2/(k1)`$ should exist above which fluctuations would be negligible and mean field behavior ($`1\theta t^{1/(k1)}`$) would dominate. This would imply that for $`d=2`$ the mean field law should work for $`k3`$; nevertheless we have found no evidence of that, at least for $`k=3`$. In order to draw such a conclusion for larger $`k`$, longer runs would be required, since a change in the dynamical behavior could occur at very large times.

Finally we have compared the data for the density of one-site vacancies and of many-site vacancies, for $`k=3,4,5`$ and $`p=0.1`$ (Fig. 9). At variance with the case $`k=2`$, many-site vacancies tend now to prevail for a long time and disappear more slowly for larger $`k`$; the time behavior of their density can be well approximated by a single power law:

$$1\rho _1(t)\theta (t)t^{\beta _k},$$
(8)
with $`\beta _3=0.94`$, $`\beta _4=0.75`$ and $`\beta _5=0.61`$. Moreover, the crossover time $`t_c`$ after which one-site vacancies dominate increases with $`k`$, as shown in table VII of $`t_c`$ values, where we note that for higher $`k`$ values diffusive regime sets in at later time.

## VI Summary and conclusions

We studied irreversible deposition and diffusion of dimers and $`k`$-mers on a square lattice, using a driven algorithm which allowed to decrease the CPU time needed for Monte Carlo simulations. Both time coverage behavior and vacancy dynamics have been analyzed in detail as a function of initial conditions and deposition rate and probability. A wider variety of dynamical regimes with respect to the one dimensional case has been observed, reflecting the complexity of the configurational space accessible to dimers in $`d=2`$ with respect to $`d=1`$. In particular, the existence of dynamically frozen configurations, that only depend on the local environment, slows down the asymptotic dynamics and does not allow the system to reach the complete coverage. Different choices of deposition and diffusion rates and probabilities can lead to different final values of the asymptotic coverage. Simulation shows that the behavior of coverage $`\theta `$ strongly deviates from the mean field prediction. Deviations from mean field are less pronounced when the diffusion probability is higher. Furthermore $`\theta `$ is observed to depend on the initial conditions only in the early stage of the dynamics. Starting from the jammed state correlations are higher and the dynamics evolves more slowly. Vacancy dynamics is dominated not only by single-site vacancies, as in $`d=1`$, but also by larger vacancies.

It was also shown that the dynamical behavior is severely slowed down with the increase of the length $`k`$ of the particles. We do not have evidences of a mean field behavior, in spite of the fact that it would be expected to hold for $`k3`$. We have also pointed out that the dimensionality plays an important role in reaction-diffusion models with the result of a rich dynamical behavior in the investigated quantities.

As an open possibility for the future, we think that in order to better understand the dynamics of these models it would be useful to vary the kind of relaxation mechanisms, investigating the possibility of occurrence of collective effects in $`d=2`$ and $`d=3`$.

## Acknowledgments

C. F. would like to thank C.K. Gan and J.S. Wang for having received the code of their algorithm. A. P. enjoyed many useful discussions with M.J. de Oliveira.


text:
# On Some Exponential Potentials for a Cosmological Scalar Field as Quintessence

## I Introduction

In the last few years, new models of the universe have been built taking dark energy into account . Together with baryons, cold dark matter, photons and neutrinos, a fifth component has been added, the socalled quintessence field $`Q`$  (or, in general, the x-field  ). Generalizing ideas like that of a cosmic equation of state variable with a $`\mathrm{\Lambda }`$-term , with respect to a more usual cosmological constant $`\mathrm{\Lambda }`$, such a $`Q`$-field, even if still implying a negative pressure contribution to the total pressure of the cosmic fluid, is characterized by the fact that its equation of state is given by $`1<w_Qp_Q/\rho _Q<0`$, $`p_Q`$ and $`\rho _Q`$ being, respectively, the pressure and energy density of the $`Q`$-field. Actually, the interval $`1<w_Q0.6`$ is usually considered . As a matter of fact, when $`w_Q=1`$ we recover a constant $`\mathrm{\Lambda }`$-term , which can be regarded as a measure of vacuum energy density, leading to the well known discrepancy between theory and observations , based on the question of why $`\rho _Q`$ is so small with respect to typical particle physics scales. (But there are also mechanisms of relaxation of the cosmological constant during the initial inflationary stage, which could explain such a discrepancy; see , for instance.)

As well known, an interesting possibility to handle the presence of quintessence in the universe is to see it as given by a scalar field $`\phi `$ slowly rolling down its potential $`V(\phi )`$. If we define

$$\rho _\phi \frac{1}{2}\dot{\phi }^2+V(\phi );p_\phi \frac{1}{2}\dot{\phi }^2V(\phi ),$$
(1)
(dot indicating time derivative, and $`V(\phi )`$ being the potential for $`\phi `$), the slow rolling condition immediately gives $`p_\phi <0`$ and $`w_\phi p_\phi /\rho _\phi 1`$. With such a negative pressure, the universe evolves like in a sort of present day soft inflationary scenario, so allowing to explain observations on supernovae  and why vacuum and matter densities are today comparable (‘cosmic coincidence’ problem ). (Alternatively to quintessence, a negative pressure and an explanation of current observations can also be obtained in a Chaplygin cosmology .)

Many cosmological models with a dynamical scalar field have been proposed, showing scaling solutions, i.e., such that at some time $`\rho _m`$ and $`\rho _\phi `$ simultaneously depend on some powers $`n_1`$ and $`n_2`$ of the scale factor $`a`$, acting as attractors in the phase space. When $`n_1=n_2`$, we have the socalled self-tuning solutions , which are typically driven by exponential potentials. This kind of potential has been studied extensively , especially from a qualitative point of view (see  and references therein, for instance).

The simplest possibility is of course $`V(\phi )=\alpha e^{\lambda \phi }`$, which is often discarded (see discussion in ). More promising seems to be a combination of two terms $`V(\phi )=\alpha e^{\lambda \phi }+\beta e^{\lambda \phi }`$ .

In this paper, we consider a particular class of both these types from a different point of view, obtaining general exact solutions. This allows a very stringent comparison with experimental data on supernovae, so that also the first type seems to deserve further investigations; for the second type, we obtain a solution which can mimic very well the presence of a cosmological constant in the late evolution of the universe. Both of them are not, strictly speaking, scaling solutions, although this concept may be recovered in a more general sense.

Another experimental fact which we use is the strong evidence of a spatially flat universe . Thus, we set the scalar curvature $`k=0`$ in all our equations from the very beginning. However, as we shall see, the values of $`\mathrm{\Omega }_m`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }`$ ($`\mathrm{\Omega }\phi `$ in our case) derived from the experiments strongly depend on the model, so that some discussion is needed.

Mostly, the scalar field $`\phi `$ has been considered as minimally coupled to gravity, even if (more recently) a nonminimal coupling has also been introduced . Here, we will consider a very simple model consisting of a two-component cosmological fluid: matter and scalar field. ‘Matter’ means baryonic + cold dark matter, with no pressure, and the scalar field is minimally coupled and noninteracting with matter. Clearly, this model cannot be used from the very beginning of the universe, but only since decoupling of radiation and dust. Thus, we do not take into account inflation, creation of matter, nucleosynthesis, etc. The main shortcut is that we cannot really check for the tracker feature  of the $`\phi `$-solution.

In Sec. II we take a particular exponential potential into account, deriving the general exact solution of the cosmological equations and, thereby, cosmological parameters as functions of time, so allowing the comparison with observational data. Sec. III is devoted to the same kind of considerations for a potential given by a linear combination of two such exponentials. Conclusions are drawn in Sec. IV.

## II An exponential potential

### A Mathematical treatment

Let us consider a spatially flat, homogeneous and isotropic universe, filled with two noninteracting components only, i.e., pressureless matter (or dust) and a scalar field $`\phi `$, minimally coupled to gravity. The cosmological equations are then

$$3H^2=\frac{8\pi G}{c^2}(\rho _m+\rho _\phi ),$$
(2)
$$\dot{H}+H^2=\frac{4\pi G}{3c^2}(\rho _m+\rho _\phi +3(p_m+p_\phi )),$$
(3)
$$\ddot{\phi }+3H\dot{\phi }+V^{}(\phi )=0,$$
(4)
where prime indicates derivative with respect to $`\phi `$, $`H\dot{a}/a`$ is the Hubble parameter, $`p_m=w_m\rho _m`$ and $`p_\phi =w_\phi \rho _\phi `$ are the equations of state for matter and scalar field. Let us stress that $`w_\phi `$ is  not constant, and that we want to describe some features of cosmology after the decoupling.

We set $`w_m=0`$, so that $`\rho _m=Da^3`$. The parameter $`D\rho _{m}^{}{}_{0}{}^{}a_{0}^{}{}_{}{}^{3}`$ (the lower index ‘$`0`$’ indicating present day values) is the amount of matter. The equations can also be rewritten as

$$\left(\frac{\dot{a}}{a}\right)^2=\frac{8\pi G}{3c^2}(Da^3+\frac{1}{2}\dot{\phi }^2+V(\phi )),$$
(5)
$$2\frac{\ddot{a}}{a}+\left(\frac{\dot{a}}{a}\right)^2=\frac{8\pi G}{3c^2}\left(\frac{1}{2}\dot{\phi }^2V(\phi )\right),$$
(6)
$$\ddot{\phi }+3\left(\frac{\dot{a}}{a}\right)\dot{\phi }+V^{}(\phi )=0.$$
(7)
In this Section, we consider the potential

$$V(\phi )=B^2e^{\sigma \phi },$$
(8)
where $`B^2`$ is a generic positive constant and

$$\sigma ^2\frac{12\pi G}{c^2}.$$
(9)
(The minus sign in the exponential is irrelevant, since there is symmetry with respect to a change $`\phi \phi `$.)

This type of potential leads to a late time attractor in a scalar-field dominated situation ($`\mathrm{\Omega }_\phi =1,w_\phi =0.5`$) . Being aware of such a behaviour, anyway, we stress that we are especially interested in the contemporary or, at most, the recent past regimes, where the situation is different. Usually, associated with an exponential potential, a scalar field is considered such that $`\mathrm{\Omega }_\phi 8\pi G\rho _\phi /(3c^2H^2)`$ is practically constant during part of the matter-dominated era. This implies that assuming $`w_\phi constant`$ leads to a constant ratio of quintessence to matter energy density, so that $`\mathrm{\Omega }_\phi `$ (being $`0.15`$ at the beginning of matter-dominated era, due to nucleosynthesis ) must remain small forever . Mainly for such reasons, this kind of potential is not considered as suitable for a quintessence field.

The particular choice of Eq. (9) for $`\sigma `$ allows for general exact integration of equations. Such a choice was in fact used in the context of inflationary theory by us  and others , with a scalar field only.

Let us concentrate on the second order equations (6) and (7), while Eq. (5), which is a first integral, is considered as a constraint on the integration constants. Let us introduce the new variables $`u`$ and $`v`$, defined by the tranformation

$$a^3=uv;\phi =\frac{1}{\sigma }\mathrm{log}\frac{u}{v},$$
(10)
which is always invertible (the Jacobian being $`J=2/\sigma `$). We get for the potential

$$V(u,v)=B^2\frac{u}{v},$$
(11)
and Eqs. (6) and ( 7) become

$$\ddot{u}=0;\ddot{v}=\omega u,$$
(12)
where $`\omega =\sigma ^2B^2=12\pi GB^2/c^2>0`$. They are immediately integrated to

$`u(t)`$ $`=`$ $`u_1t+u_2,`$ (13)
$`v(t)`$ $`=`$ $`{\displaystyle \frac{1}{6}}u_1\omega t^3+{\displaystyle \frac{1}{2}}u_2\omega t^2+v_1t+v_2,`$ (14)
being $`u_1`$, $`u_2`$, $`v_1`$, and $`v_2`$ arbitrary integration constants. Taking into account Eq. (5), we find

$$\frac{2}{3}u_1v_1\frac{1}{3}\omega u_2\frac{4\pi G}{c^2}D=0.$$
(15)
Since $`D`$ is a physical parameter, it would be natural to use it as given and derive one of the other constants. But this complicates calculations without substantial advantages, as its value depends on the normalization of the scale factor $`a`$, which can be fixed arbitrarily. Thus, we will determine $`D`$ from Eq. (15)

$$D=\frac{c^2}{12\pi G}\left(2u_1v_1\omega u_2\right).$$
(16)
Being of course $`D>0`$, this gives limitations on the choice of $`u_1`$, $`u_2`$, $`v_1`$.

The well known  solution $`\phi =2/\sqrt{3(1+\alpha )}\mathrm{log}t`$, coming from the potential

$$V(\phi )=\frac{2(1\alpha )}{3(1+\alpha )^2}\mathrm{exp}\left(\sqrt{3(1+\alpha )}\phi \right),$$
(17)
is a very particular case of what we find in Eqs. (13) and (14 ). It can be obtained by setting

$$B^2=\frac{2(1\alpha )}{3(1+\alpha )^2};\sigma =\sqrt{3(1+\alpha )};u_2=v_2=v_1=0.$$
(18)
Eq. (16) then gives $`D=0`$, so that we get a model without matter, not really interesting in our context.

A more interesting possibility is given by the choice $`u_1=v_1=v_2=0`$, involving, from Eq. (16), $`u_2<0`$. This implies, in fact,

$$D=\frac{c^2\omega u_2}{12\pi G},a^3=\frac{1}{2}u_{2}^{}{}_{}{}^{2}\omega t^2,\rho _\phi =\frac{4}{\sigma ^2t^2}.$$
(19)
In this particular case, $`\rho _\phi `$ ($`a^3`$) scales as $`\rho _m`$. (But see below for a discussion on the scaling properties.)

Without any special assumptions on constants, we can get many important quantities as functions of $`u`$ and $`v`$ (we do not write them explicitly in terms of $`t`$, for sake of brevity)

$`\rho _\phi (u,v)`$ $`=`$ $`{\displaystyle \frac{1}{\sigma ^2}}\left({\displaystyle \frac{(\dot{u}vu\dot{v})^2}{2u^2v^2}}+\omega {\displaystyle \frac{u}{v}}\right),`$ (20)
$`p_\phi (u,v)`$ $`=`$ $`{\displaystyle \frac{1}{\sigma ^2}}\left({\displaystyle \frac{(\dot{u}vu\dot{v})^2}{2u^2v^2}}\omega {\displaystyle \frac{u}{v}}\right),`$ (21)
$`w_\phi (u,v)`$ $`=`$ $`{\displaystyle \frac{(\dot{u}vu\dot{v})^22\omega u^3v}{(\dot{u}vu\dot{v})^2+2\omega u^3v}},`$ (22)
$`H(u,v)`$ $`=`$ $`{\displaystyle \frac{\dot{u}v+u\dot{v}}{3uv}},`$ (23)
$`\mathrm{\Omega }_\phi (u,v)`$ $`=`$ $`{\displaystyle \frac{(\dot{u}vu\dot{v})^2+2\omega u^3v}{(\dot{u}v+u\dot{v})^2}},`$ (24)
$`\mathrm{\Omega }_m(u,v)`$ $`=`$ $`{\displaystyle \frac{24\pi GDuv}{c^2(\dot{u}v+u\dot{v})^2}}={\displaystyle \frac{4u\dot{u}v\dot{v}2\omega u^3v}{(\dot{u}v+u\dot{v})^2}},`$ (25)
the last equality coming from Eq. (16). It can be easily checked that $`\mathrm{\Omega }_m+\mathrm{\Omega }_\phi =1`$.

### B Physical considerations

We now pass to simplify the situation with a suitable choice of initial conditions. Of course, the following choice is not the only one possible. We shall see, anyhow, that it is able to reproduce observational data. Other more general choices may improve the situation, but we do not treat them here.

If $`t_{in}u_2/u_1`$ (which is always possible, being $`u_10`$), the scale factor $`a`$ is zero, and we can show that there is no other time $`t>t_{in}`$ when this occurs again. We can thus fix the time origin in such a way that $`a(0)=0`$. This condition has to be interpreted just as an arbitrary choice of the time origin. The real beginning (of physical meaning) for the model starts a little bit afterwards, at a time $`t_1`$. This delay is otherwise arbitrary, so that this setting does not seem to exclude important cases, as said before, and leads to a great simplification in the formulae. Now, $`a(0)=0`$ implies $`u_2=0`$ or $`v_2=0`$, or both. If we set only one of them to zero, we obtain $`\phi (0)=\mathrm{}`$ (which could be accepted, but is rather disturbing), and, most of all, $`\mathrm{\Omega }_\phi (0)=1`$, which would mean an initial scalar-field dominated universe, with a neglectable content of other types of matter. Of course, if we consider the situation in general, the scalar field does dominate. But, as already mentioned, in our case we start after decoupling time, when a matter-dominated behaviour seems to be more natural. So, if we set $`u_2=v_2=0`$, we get instead (in a matter-dominated situation, then)

$$\phi (0)=\frac{1}{\sigma }\mathrm{log}\frac{u_1}{v_1};\mathrm{\Omega }_\phi (0)=0;D=\frac{c^2}{6\pi G}u_1v_1.$$
(26)
We prefer to stick to this choice, so that we have

$$u(t)=u_1t;v(t)=\frac{1}{6}u_1\omega t^3+v_1t.$$
(27)
Let us now define a time scale $`t_s`$ such that $`H(t_s)=1/t_s`$, which is of the order of the age of the universe. This leads to consider a dimensionless time $`\tau t/t_s`$. From Eqs. (23) and (27) we get

$$t_s^2=\frac{6v_1}{\omega u_1}.$$
(28)
By means of these choices the formulae found above for the relevant cosmological parameters reduce to

$`\rho _\phi `$ $`=`$ $`{\displaystyle \frac{2(3+4\tau ^2)}{\sigma ^2t_s^2(1+\tau ^2)^2}},`$ (29)
$`p_\phi `$ $`=`$ $`{\displaystyle \frac{2(3+2\tau ^2)}{\sigma ^2t_s^2(1+\tau ^2)^2}},`$ (30)
$`w_\phi `$ $`=`$ $`{\displaystyle \frac{3+2\tau ^2}{3+4\tau ^2}},`$ (31)
$`a`$ $`=`$ $`(u_1v_1t_s^2(1+\tau ^2)\tau ^2)^{1/3},`$ (32)
$`(1+z)^3`$ $`=`$ $`{\displaystyle \frac{\tau _0^2(1+\tau _0^2)}{\tau ^2(1+\tau ^2)}},`$ (33)
$`H`$ $`=`$ $`{\displaystyle \frac{2(1+2\tau ^2)}{3t_s\tau (1+\tau ^2)}},`$ (34)
$`\mathrm{\Omega }_m`$ $`=`$ $`{\displaystyle \frac{1+\tau ^2}{(1+2\tau ^2)^2}},`$ (35)
where $`za(\tau _0)/a(\tau )1`$ is the redshift, and $`\tau _0`$ indicates the present time.

If we define dimensionless pressure and energy density

$$\stackrel{~}{p}_\phi \frac{\sigma ^2t_s^2}{2}p_\phi ;\stackrel{~}{\rho }_\phi \frac{\sigma ^2t_s^2}{2}\rho _\phi ,$$
(36)
we find the equation of state for the scalar field

$$\stackrel{~}{p}_\phi =\stackrel{~}{\rho }_\phi 12+6\sqrt{4\stackrel{~}{\rho }_\phi },$$
(37)
which is well approximated by

$$\stackrel{~}{p}_\phi =0.382\stackrel{~}{\rho }_\phi 0.196\stackrel{~}{\rho }_\phi ^2.$$
(38)
In Fig. 1, we compare the plots of the two functions in Eqs. (37) and (38), and show that the approximation is quite good; a comparison is also made with a straight line $`\stackrel{~}{p}_\phi =0.86\stackrel{~}{\rho }_\phi `$, where the coefficient has been obtained through a numerical approximation, as well those in Eq. (38).

From Eqs. (29) and (32) it is possible to derive

$$\delta _1\frac{d\mathrm{log}\rho _\phi }{d\mathrm{log}a}=\frac{3\tau ^2(5+\tau ^2)}{(3+\tau ^2)(1+2\tau ^2)}.$$
(39)
Thus, it is clearly

$$\delta _10\text{ for }\tau 0,\delta _1\frac{3}{2}\text{ for }\tau \mathrm{},$$
(40)
so that we asymptotically have two scaling regimes: $`\rho _\phi const.`$ for early times, and $`\rho _\phi a^{3/2}`$ for late times. In fact, this approximation holds for a very long time, well far from the asymptotic values. Indeed, computing the $`n`$-th derivative $`\delta _na(d\delta _{n1}/d(\tau ^2))/(da/d(\tau ^2))`$, it is possible to show that they are all asymptotically zero up to, say, $`n=10`$.

An estimate of the moment in which the regime changes can be given finding the maximum (or the minimum) of $`\delta _2`$. This is achieved at $`\tau 0.3`$, and it is remarkable that this result depends only on $`t_s`$. All this discussion shows that it is possible to generalise the concept of scaling solutions. The situation is illustrated in Figs. (2) and (3).

As a matter of fact, in the literature it is widely accepted that using an exponential potential leads to a dark energy density which scales like matter. Our results seem to be in contrast with this statement, which is a consequence of assuming $`w_\phi `$ almost perfectly constant. It is not our case, as shown in the following.

As in , we can use the function $`\mathrm{\Gamma }V^{\prime \prime }V/(V^{})^2`$. Defining

$$x\frac{\dot{\phi }^2}{2V}=\frac{1+w_\phi }{1w_\phi },\dot{x}\frac{d\mathrm{log}x}{d\mathrm{log}a},\ddot{x}\frac{d^2\mathrm{log}x}{d\mathrm{log}a^2},$$
(41)
it is possible to find

$$\mathrm{\Gamma }=1+\frac{w_mw_\phi }{2(1+w_\phi )}\frac{1+w_m2w_\phi }{2(1+w_\phi )}\frac{\dot{x}}{6+\dot{x}}\frac{2}{1+w_\phi }\frac{\ddot{x}}{(6+\dot{x})^2}.$$
(42)
If one makes the assumption that $`w_\phi `$ is nearly constant, then $`\dot{x}\ddot{x}0`$. Since it is $`\mathrm{\Gamma }=1`$ strictly in the case of our potential, Eq. (42) implies $`w_mw_\phi `$. The point is that it is not true that our exact solutions for the exponential potential lead to $`\dot{x}\ddot{x}0`$. For instance, we have

$$x=\frac{\tau ^2}{3(1+\tau ^2)},\dot{x}=\frac{3}{1+2\tau ^2},$$
(43)
so that $`\dot{x}0`$ only asymptotically. Since in our model $`\tau 1`$, $`x`$ and $`\dot{x}`$ are then far from being zero.

But let us also consider the third term in the right-hand side of Eq. (42) with $`w_m=0`$ (matter is simply dust in our model). Substituting our solution, we get

$$\frac{1+w_m2w_\phi }{2(1+w_\phi )}\frac{\dot{x}}{6+\dot{x}}=\frac{9+8\tau ^2}{12\tau ^2+16\tau ^4};$$
(44)
we see that this expression diverges for $`\tau 0`$, being always $`>0.5`$ in the useful range of $`\tau `$.

In our opinion, the main check for the solution in Eq. (27) is thus only its capability to reproduce the experimental results, which we are going to do just below.

From Eq. (35) we get

$$\tau _{0}^{}{}_{}{}^{2}=\frac{14\mathrm{\Omega }_{m}^{}{}_{0}{}^{}+\sqrt{1+8\mathrm{\Omega }_{m}^{}{}_{0}{}^{}}}{8\mathrm{\Omega }_{m}^{}{}_{0}{}^{}}.$$
(45)
Once we give an acceptable value for $`\mathrm{\Omega }_{m}^{}{}_{0}{}^{}`$, we obtain a value for $`\tau _0`$. For instance, $`\mathrm{\Omega }_{m}^{}{}_{0}{}^{}=0.3`$ gives $`\tau _0=0.82`$, and this implies $`w_{\phi }^{}{}_{0}{}^{}=0.76`$. If the value $`H_{0\text{ }}=100hkms^1Mpc^1`$ is also given, we get

$$t_s=\frac{2(1+2\tau _0^2)}{3H_0\tau _0(1+\tau _0^2)}=\frac{\mathrm{1.\hspace{0.17em}14}}{H_0};$$
(46)
assuming $`h=0.7`$, we have $`t_s=15.8\times 10^9`$ years, and $`t_0=13\times 10^9`$ years. It is also possible to obtain the relation between $`w_{\phi }^{}{}_{0}{}^{}`$ and $`\mathrm{\Omega }_{m}^{}{}_{0}{}^{}`$

$$w_{\phi }^{}{}_{0}{}^{}=\frac{1+8\mathrm{\Omega }_{m}^{}{}_{0}{}^{}3\sqrt{1+8\mathrm{\Omega }_{m}^{}{}_{0}{}^{}}}{4(1\mathrm{\Omega }_{m}^{}{}_{0}{}^{})}.$$
(47)
For $`\mathrm{\Omega }_{m}^{}{}_{0}{}^{}=0.2÷0.4`$, we get $`w_{\phi }^{}{}_{0}{}^{}=\mathrm{0.\hspace{0.17em}699}÷\mathrm{0.\hspace{0.17em}811}`$, and the value $`0.5`$ is reached only in the case of $`\mathrm{\Omega }_{m}^{}{}_{0}{}^{}=0`$. It is also possible to obtain $`w_\phi `$ as a function of the redshift

$$w_\phi =\frac{2\zeta +\sqrt{\zeta (4+\zeta )}}{\zeta +2\sqrt{\zeta (4+\zeta )}},$$
(48)
where $`\zeta (1+z)^3`$.

Figs. 4, 5, and 6 show the dependance of $`w\varphi `$ versus $`\tau `$, $`\mathrm{\Omega }_{m0}`$, and $`\zeta `$, respectively. We can see that $`w_\phi `$ varies much and has values $`0.78`$. For $`\zeta =1`$ (now) it is $`w_\phi =0.77`$, and already for $`\zeta =4`$ ($`z0.59`$) we find $`w_\phi 0.89`$. It is remarkable that the values of $`|w_\phi |`$ are greater than $`0.7`$. According to our knowledge, this feature is found only in .

Another interesting quantity is the present value of $`\phi `$. After staightforward algebra we get

$$\phi _0=\frac{1}{\sigma }\mathrm{log}\left(\frac{27H_0^2}{16B^2}(14\mathrm{\Omega }_{m}^{}{}_{0}{}^{}+\sqrt{1+8\mathrm{\Omega }_{m}^{}{}_{0}{}^{}})\right),$$
(49)
and we see that this value depends on the observed parameters and on the value of $`B^2`$, which was until now completely undetermined. Now, for $`\tau _0=0.82`$ and $`t_0=13\times 10^9`$ years, we have $`B^2=2.5\times 10^{47}\mathrm{exp}(\sigma \phi _0)GeV^4`$. Considering $`\phi _01/6M_P`$ ($`M_P`$ being the Planck mass) we see that $`\mathrm{exp}(\sigma \phi _0)1`$, and we can determine the unknown parameter for the potential

$$V(0)B^22.5\times 10^{47}GeV^4.$$
(50)
But we have also to observe that a ‘little’ change in $`\phi _0`$ entails a ‘large’ change in $`B^2`$. For instance, if $`\phi _0M_P`$, then $`\mathrm{exp}(\sigma \phi _0)0.0025`$ and $`B^2`$ changes of three orders of magnitude. Due to Eq. (26), we have that $`\phi (0)=1/\sigma \mathrm{log}(2\pi GB^2t_s^2/c^2)`$. This means therefore that a relatively wide range of initial values of $`\phi `$ ends up to a narrower set of final $`\phi _0`$’s.

Thus, everything seems to work fine, but things are more complicated. Indeed, one has to ask what is really measured in the supernovae experiment. The value $`\mathrm{\Omega }_{m}^{}{}_{0}{}^{}=0.3`$ is not a direct consequence of the data, since it depends on the model, which uses the constant $`\mathrm{\Lambda }`$-term. What we really measure is the distance modulus, so that it is this quantity that we should compare in the two situations. Here, we limit ourselves to a very qualitative discussion.

Let us recall, then, the definitions of luminosity distance (in Mpc)

$$d_L=3000(1+z)_0^z\frac{dz\stackrel{´}{}}{H(z\stackrel{´}{})},$$
(51)
and distance modulus

$$\delta mM=5\mathrm{log}_{10}d_L(z)+25.$$
(52)
We have thus to compare this last quantity in the case when $`H(z)`$ is taken from the usual model with $`\mathrm{\Lambda }`$ , that is,

$$H(z)=H_0\sqrt{(1+z)^2(1+\mathrm{\Omega }_{m}^{}{}_{0}{}^{}z)z(2+z)(1\mathrm{\Omega }_{m}^{}{}_{0}{}^{})},$$
(53)
with the one obtained eliminating $`\tau `$ from Eqs. (33) and (34).

In Fig. 7 we compare $`\delta `$ with $`\stackrel{~}{\delta }`$ (let us mark with a $``$ the values for the model with $`\mathrm{\Lambda }`$). The agreement is almost perfect, up to $`0.06\%`$. But there is a trick! $`\stackrel{~}{\delta }`$ was obtained from a value $`\stackrel{~}{\mathrm{\Omega }}_{m0}=0.37`$ . Of course, this value is still in the possible range but at its limit. If we decide to trust strongly on the value $`\stackrel{~}{\mathrm{\Omega }}_{m0}=0.3`$ and want to obtain the same good agreement, we have to change the value of $`\tau _0`$ to $`1.22`$. This gives a very different value $`\mathrm{\Omega }_{m}^{}{}_{0}{}^{}=0.16`$ in the model with $`\phi `$. This is again at the limit of possible estimates (due to other investigations on dark matter).

In conclusion, we see that this solution (with the potential in Eq. (8)) is indeed difficult to fully adapt to observed data, but for reasons which are not easy to investigate without general exact solutions. Moreover, it is not clearly incompatible (until we get better data); therefore, it seemed to us useful to present it in detail.

## III Two exponentials combined

### A Mathematical treatment

We now consider a combination of two exponentials, which will give us much better results, as expected. The procedure strictly follows the above one.

Let us consider the potential

$$V(\phi )=A^2e^{\sigma \phi }+B^2e^{\sigma \phi },$$
(54)
with $`\sigma ^2=12\pi G/c^2`$ as before, and $`A^2`$, $`B^2`$ arbitrary parameters. We use, now, the following change of variables

$$a^3=\frac{u^2v^2}{4};\phi =\frac{1}{\sigma }\mathrm{log}\frac{B(u+v)}{A(uv)},$$
(55)
which is invertible, provided that $`a0`$. This leads to

$$V(u,v)=2AB\frac{u^2+v^2}{u^2v^2}.$$
(56)
With these variables Eqs. (6) and (7) are rewritten as

$$\ddot{u}=\omega ^2u;\ddot{v}=\omega ^2v,$$
(57)
where now

$$\omega ^2=\frac{12\pi GAB}{c^2}.$$
(58)
Again, the integration is immediate, and gives the general solutions

$`u(t)`$ $`=`$ $`\alpha e^{\omega t}+\beta e^{\omega t},`$ (59)
$`v(t)`$ $`=`$ $`v_1\mathrm{sin}(\omega t+v_2),`$ (60)
with $`\alpha `$, $`\beta `$, $`v_1`$, $`v_2`$ arbitrary constants. As before, we derive $`D`$ from the constraint in Eq. (5)

$$D=\frac{c^2\omega ^2(v_1+4\alpha \beta )}{24\pi G}.$$
(61)
Being $`D>0`$, this implies $`v_1<4\alpha \beta `$. A change in the sign of $`v_1`$ has the only effect of changing the sign of $`\phi `$, and interchanging $`A^2`$ with $`B^2`$. So, we can set $`v_1>0`$ without any loss of generality (the case $`v_1=0`$ is obviously equivalent to considering a $`\mathrm{\Lambda }`$-term). As a consequence, $`\alpha `$ and $`\beta `$ must be non zero and with opposite signs.

Again, we can write some important functions in terms of $`u`$, $`v`$

$`\rho _\phi `$ $`=`$ $`{\displaystyle \frac{2((\dot{v}u\dot{u}v)^2+\omega ^2(u^4v^4))}{\sigma ^2(u^2v^2)^2}},`$ (62)
$`p_\phi `$ $`=`$ $`{\displaystyle \frac{2((\dot{v}u\dot{u}v)^2\omega ^2(u^4v^4))}{\sigma ^2(u^2v^2)^2}},`$ (63)
$`w_\phi `$ $`=`$ $`{\displaystyle \frac{(\dot{v}u\dot{u}v)^2\omega ^2(u^4v^4)}{(\dot{v}u\dot{u}v)^2+\omega ^2(u^4v^4)}}.`$ (64)
(This last expression gives $`w_\phi 1`$ when $`(\dot{v}u\dot{u}v)^2u^4v^4`$, which certainly can happen for sufficiently large times.)

For the redshift we have

$$(1+z)^3=\frac{4a_{0}^{}{}_{}{}^{3}}{u^2v^2},$$
(65)
and finally

$`H(u,v)`$ $`=`$ $`{\displaystyle \frac{2(u\dot{u}v\dot{v})}{3(u^2v^2)}},`$ (66)
$`\mathrm{\Omega }_m(u,v)`$ $`=`$ $`{\displaystyle \frac{(u^2v^2)(\dot{v}^2\dot{u}^2+\omega ^2(u^2+v^2))}{(u\dot{u}v\dot{v})^2}},`$ (67)
$`\mathrm{\Omega }_\phi (u,v)`$ $`=`$ $`{\displaystyle \frac{(\dot{v}u\dot{u}v)^2+\omega ^2(u^4v^4)}{(u\dot{u}v\dot{v})^2}},`$ (68)
with $`\mathrm{\Omega }_m+\mathrm{\Omega }_\phi =1`$, of course.

### B Physical considerations

The potential in Eq. (54) has a nonzero minimum $`V_{min}=A^2+B^2`$. This is unusual in quintessence theory, since $`V_{min}=0`$ is optimal to remove the fine-tuning problem. Indeed, $`V_{min}0`$ can be seen as a disguised $`\mathrm{\Lambda }`$-term. Anyway, in the following we find that this is actually the case when the scalar field is almost stationary near the minimum of the potential. But there is also the possibility of a slower rolling far from the minimum. The situation is then similar to that in Sec. II, but the additional term in the potential now allows to achieve a better agreement with observational data.

Let us now make a trial for the choice of the free parameters. We set again $`a(0)=0`$, and ask for nonsingular $`\phi (0)`$. The situation, and hence the interpretation, is the same as above. Thus, we pose $`\alpha =\beta =\lambda /2`$, $`v_2=0`$. It is also possible to fix an arbitrary normalization for $`a`$ and set $`v_1=1`$, obtaining at last

$$u(t)=\lambda \mathrm{sinh}(\omega t);v(t)=\mathrm{sin}(\omega t).$$
(69)
We get now

$$D=\frac{c^2\omega ^2(\lambda ^21)}{24\pi G},$$
(70)
implying $`|\lambda |>1`$, and

$$\phi (0)=\frac{1}{\sigma }\mathrm{log}\frac{B(\lambda +1)}{A(\lambda 1)}.$$
(71)
We define a dimensionless time $`\tau =\omega t`$ and get

$`a(\tau )`$ $`=`$ $`\left({\displaystyle \frac{\lambda ^2\mathrm{sinh}^2\tau \mathrm{sin}^2\tau }{4}}\right)^{1/3},`$ (72)
$`(1+z)^3`$ $`=`$ $`{\displaystyle \frac{\lambda ^2\mathrm{sinh}^2\tau _0\mathrm{sin}^2\tau _0}{\lambda ^2\mathrm{sinh}^2\tau \mathrm{sin}^2\tau }},`$ (73)
$`H(\tau )`$ $`=`$ $`{\displaystyle \frac{\omega (\mathrm{sin}(2\tau )\lambda ^2\mathrm{sinh}(2\tau ))}{3(\mathrm{sin}^2\tau \lambda ^2\mathrm{sinh}^2\tau )}},`$ (74)
$`w_\phi (\tau )`$ $`=`$ $`{\displaystyle \frac{\lambda ^2(\mathrm{cosh}\tau \mathrm{sin}\tau \mathrm{cos}\tau \mathrm{sinh}\tau )^2+(\mathrm{sin}^4\tau \lambda ^4\mathrm{sinh}^4\tau )}{\lambda ^2(\mathrm{cosh}\tau \mathrm{sin}\tau \mathrm{cos}\tau \mathrm{sinh}\tau )^2(\mathrm{sin}^4\tau \lambda ^4\mathrm{sinh}^4\tau )}},`$ (75)
$`\mathrm{\Omega }_m(\tau )`$ $`=`$ $`{\displaystyle \frac{2(\lambda ^21)(\mathrm{cos}(2\tau )+\lambda ^2\mathrm{cosh}(2\tau )1\lambda ^2)}{(\mathrm{sin}(2\tau )\lambda ^2\mathrm{sinh}(2\tau ))^2}}.`$ (76)
In comparison with the situation in Sec. II, we now have one more free parameter. This gives the possibility of a much better agreement with observational data.

We can also take into some account the final and initial values of $`\phi `$, i.e., $`\phi _0`$ and $`\phi _i=\phi (0)`$. If $`\tau _0`$ is the present dimensionless time, we get

$$\mathrm{exp}\left(\sigma (\phi _i\phi _0)\right)=\frac{\lambda ^2\mathrm{sinh}\tau _0+\lambda (\mathrm{sin}\tau _0\mathrm{sinh}\tau _0)\mathrm{sin}\tau _0}{\lambda ^2\mathrm{sinh}\tau _0\lambda (\mathrm{sin}\tau _0\mathrm{sinh}\tau _0)\mathrm{sin}\tau _0}.$$
(77)
When $`\lambda >>1`$, we have $`\phi _i\phi _0`$ . In fact, $`\phi `$ is practically constant and we have $`w_\phi 1`$, with nearly perfect emulation of a cosmological constant. On the other hand, if $`\lambda 1`$, then $`\mathrm{exp}(\sigma (\phi _i\phi _0))0`$. This can be interpreted as $`\phi _i\mathrm{}`$, or, better, as the possibility of a wide range of $`\phi _i`$’s, with nearly the same final $`\phi _0`$.

Whatever is $`\lambda `$, anyhow, it is possible to obtain a good agreement with observational data. We here give only two extreme cases, with $`\lambda =30`$ and $`\lambda =1.1`$.

In the first case (case I) we set

$$\lambda =30;\tau _0=1.2;\omega =2.8\times 10^{18}s^1=2.8\times 10^{42}GeV,$$
(78)
which gives

$$H_0=70;\mathrm{\Omega }_{m}^{}{}_{0}{}^{}=0.3;w_{\phi }^{}{}_{0}{}^{}=0.999.$$
(79)
In the second case (case II) posing

$$\lambda =1.1;\tau _0=0.44;\omega =1.07\times 10^{18}s^1=1.07\times 10^{42}GeV$$
(80)
gives

$$H_0=70;\mathrm{\Omega }_{m}^{}{}_{0}{}^{}=0.3;w_{\phi }^{}{}_{0}{}^{}=0.76.$$
(81)
Again, we find that the value of $`w_\phi `$ is less than $`0.7`$, as in Sec. II. In Figs. 8 and 9 the distance modulus $`\delta `$ for these two cases is compared with the $`\mathrm{\Lambda }`$-term case. As before, the agreement is quite good, but now we have $`\stackrel{~}{\mathrm{\Omega }}_{m0}=0.30`$, so that, with this very rough analysis, it is impossible to make a distinction.

Let us present again the plots of $`\mathrm{log}(\rho )`$ versus $`\mathrm{log}(a)`$. They are shown in Figs. 10 and 11. We see that the case I is practically indistinguishable from a $`\mathrm{\Lambda }`$-term, while the case II is more similar to the situation in Sec. II, with different scaling regimes. It is interesting to note that the regimes are now three, even if the last one seems to be important only in the remote future. We stress that also in this case, as well as in Sec. II, $`w_\phi `$ is not constant.

As a final result, we plot in Figs. 12 and 13 the equation of state for the scalar field. Now it is impossible to show an exact analytical expression, so that we only give the plots (in arbitrary units) in the two examined cases. It is interesting to note that in the case I, although so similar to the pure cosmological constant case, we nonetheless obtain a nontrivial plot for the equation of state. But, clearly, this point deserves further investigation.

## IV Conclusions

We have discussed two particular kinds of potentials which have allowed the general exact integration of Friedmann equations in presence of dust (ordinary and cold dark matter) and scalar field. This has been achieved by performing suitable transformations of variables. Such transformations have not been guessed by chance, but are the results of a well known procedure, the Nöther Symmetry approach , based on an action principle. This was not mentioned before because it was unnecessary to the main goal of the discussion. Nevertheless, it now seems appropriate to stress the power of such a procedure, which allows to solve cosmological equations, often giving also informations on the potential and/or the possible coupling between the scalar field and the curvature of spacetime, without any limitation on the validity of the solution itself. (For details, see the literature quoted above, and the references therein.)

We have seen that, with a suitable choice of integration constants for both potentials, it is possible to reproduce the main recent results from supernovae (initially interpreted in a $`\mathrm{\Lambda }`$-model), with considerable precision, especially in the case II of the second potential. (The case I does not add much to what is already known in a constant $`\mathrm{\Lambda }`$-term model.) What is interesting, in our opinion, is that such kinds of models can bring to a different evaluation of an important quantity like $`\mathrm{\Omega }_{0}^{}{}_{m}{}^{}`$. This, in a certain sense, sheds new light on the exponential potential of the first type, which is usually not considered as completely adequate for quintessence, for instance. But we have seen that, without considering a priori $`w_\phi `$ as a constant and having a general exact solution, something else can be learned. As a matter of fact, not having an almost constant $`w_\phi `$ makes it impossible to treat the tracker condition in the usual way. Also, we get appreciable values of $`w_\phi `$, i.e., surely less than $`0.7`$.

Of course, all our discussion is in part still qualitative, in that we should need to make a more punctual analysis of observational data, and verify the best fit with the various models, in order to see whether and when real differences arise. Anyway, our analysis already seems to confirm some of the considerations made in .

Another important point to note is that, to be realistic and cover the whole (or, at least, a substantially wider) range of the life of the universe, radiation (and hot dark matter) must be added into the game. This could allow to study the CMBR spectrum and the formation of structures. But it presumably destroys the possibility to integrate the system of the cosmological equations, leading to the necessity of using the results we established here only as a guide for a more complete analysis.

## Acknowledgments

We thank Prof. M. Demianski for his useful suggestions and Prof. M. Sazhin for revision of the manuscript. We also thank Prof. J. D. Barrow for his comments, and Prof. V. Sahni, who has kindly attracted our attention on some papers.

This work has been in part financially sustained by the M.U.R.S.T. grant PRIN2000 “SIN.TE.SI.”.

Figure Captions

Fig. 1 - Comparison of the state equation in Eq. (36) with an approximate quadratic equation and with an approximate straight line.

Fig. 2 - The straight line represents the dependence of $`\rho _m`$ versus $`a`$ ($`\rho _m=Da^3`$). The thick curve is for the scalar field. The two tangents show the asymptotic scaling behaviour in the period of dominating matter and scalar field, respectively. The bullets indicate present time, according to $`a=0.82`$. Units are arbitrary.

Fig. 3 - The derivative of $`\mathrm{log}(\rho )`$ with respect to $`\mathrm{log}(a)`$, shows a quick transition from one scaling regime to another one. Units are arbitrary.

Fig. 4 - The plot of $`w_\varphi `$ versus time shows that $`w_\varphi `$ is far from being constant.

Fig. 5 - $`w_{\varphi 0}`$ versus $`\mathrm{\Omega }_{m0}`$.

Fig. 6 - This plot shows the dependence of $`w_\varphi `$ from the redshift.

Fig. 7 - Comparison of the distance modulus $`\delta `$(derived from Eqs. (32), (33) and (34), and assuming $`\mathrm{\Omega }_{m0}=0.3`$) (continuous line) with $`\stackrel{~}{\delta }`$ (derived from a $`\mathrm{\Lambda }`$-term model,assuming a $`\stackrel{~}{\mathrm{\Omega }}_{m0}=0.37)`$ (dots).

Fig.8 - Comparison of the distance modulus $`\delta ]`$ (derived from Eqs. (68), (69) and (71) in case I, and assuming $`\mathrm{\Omega }_{m0}=0.3`$) (continuous line) with $`\stackrel{~}{\delta }`$ (derived from a $`\mathrm{\Lambda }`$-term model, assuming $`\stackrel{~}{\mathrm{\Omega }}_{m0}=0.3)`$) \[dots\].

Fig.9 - Comparison of the distance modulus $`\delta `$(derived from Eqs. (68), (69) and (71) in case II, and assuming $`\mathrm{\Omega }_{m0}=0.3`$) (continuous line) with a (derived from a $`\mathrm{\Lambda }`$-term model, assuming $`\stackrel{~}{\mathrm{\Omega }}_{m0}=0.3)`$) \[dots\].

Fig. 10 - Case I. The straight line with slope represents the dependance of $`\mathrm{log}\rho _m`$ versus $`\mathrm{log}a`$ ($`\rho _m=Da^3`$). The horizontal thick line is for the scalar field. It emulates a constant. The bullet indicates present time, according to $`\tau _0=1.2`$. Units are arbitrary.

Fig. 11 - Case II. The straight line with slope represents the dependance of $`\mathrm{log}\rho _m`$ versus $`\mathrm{log}a`$ ($`\rho _m=Da^3`$). The thick line is for the scalar field. There are three approximate scaling regimes. The bullet indicates present time, according to $`\tau _0==0.44`$. Units are arbitrary.

Fig.12 - Case I. State equation for the scalar field. Although it perfectly emulates a cosmological constant, the equation is nontrivial. Units are arbitrary, but coherent with Fig. 13.

Fig.13 - Case II . State equation for the scalar field. Units are arbitrary, but coherent with Fig. 12.