text:
# Chemo-spectrophotometric evolution of spiral galaxies: IV. Star formation efficiency and effective ages of spirals

## 1 Introduction

The Star Formation Rate (SFR) is the most important and the less well understood ingredient in studies of galaxy evolution. Despite more than 30 years of observational and theoretical work, the Schmidt law still remains popular among theoreticians and compatible with most available observations (e.g. Kennicutt 1998a).

It is well known that there are systematic trends in the SF history of the Hubble sequence: The ratio $`\mathrm{\Psi }/\mathrm{\Psi }`$ of the current SFR $`\mathrm{\Psi }`$ to the past average one $`\mathrm{\Psi }`$ (integrated over the galaxy’s age) increases as one goes from early to late type galaxies, albeit with a large dispersion within each morphological type. The Hubble sequence seems to be determined by the caracteristic timescale for star formation, with early type galaxies forming their stars in shorter timescales than those of late types. However, this simple description of the Hubble sequence fails to answer two important, and probably related, questions: i) what determines the caracteristic timescale for SF in galaxies? and ii) what is the role (if any at all) of a galaxy’s mass?

Gavazzi et al. (1996) found an anti-correlation between the SFR per unit total mass and galactic luminosity; part of this trend may reflect the aforementioned dependence of SFR on Hubble type, but it may also be that this trend is fundamentally related to the galactic mass. In a recent work Bell and de Jong (2000) use a large sample of spiral galaxies with resolved optical and near-infrared photometry and find that it is rather surface density that drives the star formation history of galaxies, while mass is a less important parameter.

Photometric studies alone cannot lift the age-metallicity degeneracy, namely the fact that young and metal-rich stellar populations may be redder than old and metal-poor ones. Studies of the chemical aspects of galaxy evolution (i.e. gas fractions, star formation rates, metal abundances) can help to tackle the problem from a different angle, but they are also limited by the unknown history of gaseous flows from and to the system. A study combining elements from both photometric and chemical evolution offers the best chances to understand this complex situation.

In this work we study the star formation history of spirals by using detailed models of chemical and spectro-photometric galactic evolution. The numerical code has been presented in detail elsewhere (Boissier and Prantzos 1999, hereafter Paper I) and is only briefly described in Sec. 2. It has been successfully applied to the modelling of global properties of spirals (Boissier and Prantzos 2000, hereafter Paper II), as well as to the corresponding abundance and photometric gradients (Prantzos and Boissier 2000, herefter Paper III). Based on the recent work of Boselli et al. (2000), we use here a large and homogeneous sample of data for normal spirals (presented in Sec. 3), ideally suited for the purpose of this work. Our results and the comparison to the observations are presented in Sec. 4. They suggest that galactic mass is the main driver of galactic evolution, although local surface density may also play a role. We argue that massive galaxies are, on average, older than less massive ones, based on the fact that star formation efficiencies seem to be independent of galactic mass, while gas fractions are systematically low in massive spirals. Based on the observed star formation efficiencies and gas fractions, we derive in Sec. 5 corresponding galactic ages by using simple analytical models of galactic chemical evolution, taking into account gaseous flows. We find that low mass spirals must be several Gyr younger than massive ones. Also, we find that spirals with rotational velocities in the 80-400 km/s range must not have suffered extensive mass losses during their past history. Our conclusions are summarized in Sec. 6.

## 2  Modelling the evolution of disc galaxies

### 2.1 The model and the Milky Way

In Boissier and Prantzos (1999, paper I), we presented a model for the chemical and spectrophotometric evolution of the Milky Way disc. We recall here the main ingredients of the model.

The disc is simulated as an ensemble of concentric, independently evolving rings, slowly built-up by infall of primordial composition. For each ring we solve the classical equations of chemical evolution (e.g. Pagel 1997) without the assumption of Instantaneous Recycling Approximation. We use: stellar lifetimes from Schaller et al. (1992); the yields of Woosley and Weaver (1995) for massive stars and Renzini and Voli (1981) for intermediate mass stars; SNIa producing Fe with a rate given in Matteucci and Greggio (1986); and the Initial Mass Function of Kroupa et al. (1993) between 0.1 M and 100 M. The adopted star formation rate (SFR) varies with gas surface density $`\mathrm{\Sigma }_g`$ and radius $`R`$ as:

$$\mathrm{\Psi }(t,R)=\alpha \mathrm{\Sigma }_g(t,R)^{1.5}V(R)R^1$$
(1)
where $`V(R)`$ is the rotational velocity, assumed to be 220 km/s in the largest part of the Milky Way disc. This radial dependence of the SFR is suggested by the theory of star formation induced by density waves in spiral galaxies (e.g. Wyse and Silk 1989). The efficiency $`\alpha `$ of the SFR (Eq. 1) is fixed by the requirement that the observed local gas fraction, at R=8 kpc from the Galactic centre ($`\sigma _{gas}`$0.2), is reproduced at T=13.5 Gyr (our adopted value for the age of the local disc).

The disc is built up by infall with a rate exponentially decreasing in time, and a characteristic time-scale $`\tau _{inf}(R)`$ increasing with radius (as to mimic the inside-out formation of the disc). The relation between $`\tau _{inf}(R)`$ and $`\mathrm{\Sigma }(R)`$ (total surface density) is shown on Fig. 1 and is a posteriori justified, since it is crucial in shaping the various radial profiles (of gas, SFR, abundances, etc) of the disc, which compare favourably to observations.

The spectro-photometric evolution is followed in a self-consistent way, with the metallicity dependent stellar tracks from the Geneva group (Schaller et al. 1992, Charbonnel et al. 1996) and stellar spectra from Lejeune et al. (1997). Dust absorption is included according to the prescriptions of Guiderdoni et al. (1998) and assuming a “sandwich” configuration for the stars and dust layers.

It turns out that the number of observables explained by the model is much larger than the number of free parameters. In particular the model reproduces present day “global” properties (amounts of gas, stars, SFR, and supernova rates), as well as the current disc luminosities in various wavelength bands and the corresponding radial profiles of gas, stars, SFR and metal abundances; moreover, the adopted inside-out star forming scheme leads to a scalelength of $``$4 kpc in the B-band and $``$2.6 kpc in the K-band, in agreement with observations (see paper I).

### 2.2 Scaling relations

For a simplified extension of the model to the the case of other disc galaxies we adopt the “scaling properties” derived by Mo, Mao and White (1998, hereafter MMW98) in the framework of the Cold Dark Matter (CDM) scenario for galaxy formation. For the details, the reader should refer to MMW98 and Boissier and Prantzos (2000, paper II). Discs form inside non baryonic dark matter haloes of various masses and concentration factors. Assuming constant disc to halo mass ratios (here taken to be $`m_d`$=0.05), discs are characterized by two parameters: $`V_C`$, the circular velocity measuring the mass of the disc, and $`\lambda `$, the spin parameter measuring its angular momentum. A disc is described by its scale-length $`R_d`$ and its central surface density $`\mathrm{\Sigma }_0`$, which can be related to the ones of our Galaxy (designated by $`MW`$) by:

$$\frac{R_d}{R_{d,MW}}=\frac{\lambda }{\lambda _{MW}}\frac{V_C}{V_{C,MW}}$$
(2)
and

$$\frac{\mathrm{\Sigma }_0}{\mathrm{\Sigma }_{0,MW}}=\left(\frac{\lambda }{\lambda _{MW}}\right)^2\frac{V_C}{V_{C,MW}}$$
(3)
The total mass of the disc $`M_d`$ is proportional to $`V_C^3`$, but independent of $`\lambda `$. The $`\lambda `$ distribution deduced from numerical simulations (see MMW98 and references therein) presents a maximum at $`\lambda 0.040.05`$ and extends from $``$0.01 to $``$0.20. Since the value of $`\lambda _{MW}`$ is not far from the peak value, we explored here values of $`\lambda `$ in the range 1/3 $`\lambda _{MW}`$ to 3 $`\lambda _{MW}`$. Assuming that $`\lambda `$ and $`V_C`$ are independent, we constructed models with velocities in the observed range 80 to 360 km/s. We calculated the velocity profile resulting from the disc plus an isothermal dark halo, and we used Equ. 1 to calculate self-consistently the SFR (with the same coefficient $`\alpha `$, “calibrated” on the Milky Way).

At this point it should be noticed that our scaling relations (2 and 3) are based on the assumption that angular momentum is conserved during the evolution of the disk. Numerical hydrodynamical simulations of galaxy formation (e.g. Navarro and Steinmetz 1999 and references therein) do not support this idea; indeed, it is found that as the baryonic halo gas cools down and collapses to the disk it loses most of its angular momentum. In those conditions, the final configuration of the disk cannot be related in a simple way to the initial $`\lambda `$ of the halo. However, such simulations lead , in general, to disk sizes much smaller than observed. It is possible that star formation and feedback are not properly described in those simulations. For instance, Sommer-Larsen, Gelato and Vedel (1999) found that delaying the cooling of the gas reduces the loss of angular momentum. Since the situation is not clear yet, we make here the assumption of disk evolution at constant angular momentum, a posteriori justified by the fact that the resulting disk sizes are in agreement with observations (see Paper II). For a more detailed discussion of these issues, see also Cole et al. (2000).

### 2.3 Infall timescales

The infall rate in our models is exponentially dependent on time, with a time-scale $`\tau _{inf}`$ that depends on the surface mass density and is calibrated on the Milky Way disc (see Fig. 1); as explained in paper II, we found that a dependence on the total mass of the galaxy is necessary in order to reproduce the observations of present-day discs. We consider that all the galaxies started forming their stars 13.5 Gyr ago and that the formation of their exponential discs (characterized by $`\mathrm{\Sigma }_0`$ and $`R_d`$) was completed at the present epoch. Notice that the value of T=13.5 Gyr plays no essential role in the results of this work (values in the 10-15 Gyr range would lead to quite similar conclusions); these results depend mainly on the infall timescales. The more massive galaxies are characterized by shorter formation time-scales, while less massive galaxies are formed on longer ones (Fig. 1). This assumption turned out to be a crucial ingredient of our models, allowing to reproduce an impressive amount of observed properties of spirals that depend on mass (or $`V_C`$): colours, gas fractions, abundances and integrated spectra. Those observables are thoroughly presented in paper II, while some of them are revisited in this work. In Prantzos and Boissier (2000, paper III in the series), we show that this model also reproduces fairly well the observed colour and abundance gradients in spirals.

We notice that Bell and de Jong (2000) suggested recently that the observed colour gradients of disc galaxies correlate very well with the local surface density, in the sense that inner and denser regions are older. Our assumption about the infall timescale agrees, at least qualitatively, with their findings.

We stress that infall timescales are inputs to the model, and adjusted as to reproduce observations. Star formation timescales are outputs of the model (resulting from the adopted prescriptions for infall and SFR) and are presented in Sec. 4.1.

## 3 Comparison to observations

### 3.1 The observational sample

The sample of galaxies analysed in this work, which has been extracted from the large multifrequency database of nearby galaxies of Gavazzi and Boselli, is extensively described in Boselli et al. (2000). Here we give just a brief description of the sample selection criteria: we refer the reader to Boselli et al. (2000) for the detailed references on the data and on their analysis.

Galaxies analysed in this work are taken from the Zwicky catalogue (CGCG, Zwicky et al. 1961-1968)(mpg $``$ 15.7). They are either late-type (type$`>`$S0a) members of 3 nearby (recession velocity c$`z`$ 8000 km s-1) clusters (Cancer, A1367, Coma), or located in the relatively low-density regions of the Coma-A1367 supercluster (11h30m $`<`$RA$`<`$ 13h30m; 18o$`<`$dec$`<`$32o) as defined in Gavazzi et al. (1999a). To extend the present study to lower luminosities, we include in the sample the late-type Virgo cluster galaxies brighter than mpg $``$ 14.0 listed in the Virgo Cluster Catalogue as cluster members (VCC, Binggeli et al. 1985). Furthermore VCC galaxies with 14.0 $``$ mpg $``$ 16.0 included in the “ISO” subsample described in Boselli et al. (1997a) and CGCG galaxies in the region 12h $`<`$RA$`<`$ 13h; 0o$`<`$dec$`<`$18o but outside the VCC, are considered.

To avoid systematic environmental effects we consider the subsample of late-type galaxies whose HI deficiency (defined as $`d`$ = log($`HI/HI`$), the ratio of the HI mass to the average HI mass in isolated objects of similar morphological type and linear size, see Haynes & Giovanelli 1984) is $`d`$0.3, typical of unperturbed, isolated galaxies. The final combined sample comprises 233, mainly “normal” galaxies.

We assume a distance of 17 Mpc for the members (and possible members) of Virgo cluster A, 22 Mpc for Virgo cluster B, 32 Mpc for objects in the M and W clouds (see Gavazzi et al. 1999b). Members of the Cancer, Coma and A1367 clusters are assumed at the distance of 62.6, 86.6 and 92 Mpc respectively. Isolated galaxies in the Coma supercluster are assumed at their redshift distance adopting $`H_o`$ = 75 km s-1 Mpc-1.

For the 233 optically selected galaxies, data are available in several bands as follows: 100% have HI (1420 MHz) and 99% H band (1.65 $`\mu `$m) data, while a much coarser coverage exists in the UV (2000 Å)(29 %), CO (115 GHz)(38%) and H$`\alpha `$ (6563 Å) (65%), as shown in Table 1 and 2 of Boselli et al. (2000). The distribution of the sample galaxies in the different morphological classes is given in Table 3 of Boselli et al. (2000).

As previously discussed, the present sample is optically selected and thus can be biased against low surface brightness galaxies; the inclusion of the Virgo cluster should in principle favor the presence of some low surface brightness galaxies, not easily detectable at higher distances. Being volume limited, the sample is not biased toward bright, giant spirals, but it includes also dwarfs and compact sources. Its completeness at different wavelengths makes this a unique sample suitable for statistical analysis.

### 3.2 Data analysis

H$`\alpha `$ and UV fluxes, corrected for extinction (and \[NII\] contamination) as described in Boselli et al. (2000), are used to estimate star formation rates through population synthesis models given in that work; a power-law IMF of slope -2.5 with a lower and upper mass cutoff of 0.1 M and 80 M, respectively, is adopted. Its high mass part is quite similar to the one of the IMF of Kroupa et al. (1993), adopted in our models. Given the uncertainty in the UV and H$`\alpha `$ flux determination, in the extinction correction and in the transformation of the corrected fluxes into SFRs via population synthesis models, we estimate an uncertainty of a factor of $``$3 in the determination of the SFR.

The total gas content of the target galaxies, HI$`+`$H2, has been determined from HI and CO measurements. CO (at 2.6mm) fluxes have been transformed into H2 masses assuming a standard CO to H2 conversion factor of 1.0 1020 mol cm-2 (K km s-1)-1 (Digel et al. 1996). For galaxies with no CO measurement, we assume that the molecular hydrogen content is 10% of the HI, as estimated from isolated spiral galaxies by Boselli et al. (1997b). The total gas mass has been corrected for He contribution by 30%. HI fluxes are transformed into neutral hydrogen masses with an uncertainty of $``$ 10%. The average error on CO fluxes is $``$ 20%; the error on the H2 content, however, is significantly larger (and difficult to quantify) due to the poorly known CO to H2 conversion factor (see Boselli et al. 1997b).

Galaxy colours have been determined from broad band near-IR and optical photometry. Near-IR (H band) images are available for 230 of the 233 sample galaxies, while B images or aperture photometry for 214 objects. H and B magnitudes have been corrected for extinction as in Gavazzi & Boselli (1996). No correction has been applied to galaxies of type later than Scd. The estimated error on B and H magnitudes is 15 %.

Rotational velocities have been determined from the HI line width at 21cm, and corrected for inclination as in Gavazzi (1987). To avoid large systematic errors, we estimate rotational velocities only for galaxies with inclinations $`>`$ 30 deg and with the 21cm line width accurately determined (double or single horned profile with high signal-to-noise). The uncertainty on the determination of the rotational velocity is $``$ 15 km s-1.

In summary, we have a homogeneous sample of 233 normal (non-perturbed) disc galaxies. For 96 of them, we evaluated all the quantities of interest in this work: Blue magnitude $`M_B`$, total mass $`M_T`$, total gas mass $`M_g`$, star formation rate $`\mathrm{\Psi }`$ and rotational velocity $`V_C`$. We shall see below how our models fit those properties and what kind of inferences can be made on the star formation history of those galaxies.

### 3.3 Mass-driven colours and gas fractions

In Fig. 2 we compare our model results to the colours $`BH`$ of our galaxy sample. A clear correlation is obtained between $`BH`$ and $`V_C`$, with the more massive discs being, on average, redder than their lower mass counterparts. Our models (solid curves) reproduce naturally the observed correlation, since by construction the more massive discs form their stars earlier (Fig. 1). The corresponding timescales for star formation will be discussed in Sec. 4.1. The observed dispersion in the lower panel of Fig. 2 can be accounted for by the range of $`\lambda `$ values in our models, but only partially. Discs with lower $`\lambda `$ are more compact, have higher central surface densities and evolve earlier than those with larger $`\lambda `$. Obviously, $`\lambda `$ values larger than 3 $`\lambda _{MW}`$ (the largest value used here), would lead to even smaller $`BH`$ values than the models displayed on Fig. 2, possibly accounting for the rest of the scatter.

We notice that our models do take extinction by dust into account (with the prescriptions presented in Sec. 2.1), and that extinction contributes somewhat (by $``$0.5 mag) in redenning the most massive discs in our models, which have relatively large amounts of metals. However, we insist on the fact that it is age, not extinction, which is mainly responsible for the trends of our models on Fig. 2.

In Fig. 3 we present two important pieces of data, on which the main argument of this paper is based. In the upper panel, we show the ratio $`M_g/L_H`$ (mass of gas to H-band luminosity) vs rotational velocity $`V_C`$. Since $`L_H`$ reflects the mass of the stellar population (Gavazzi et al. 1996), this ratio is a measure of the gas fraction of the system. Despite the large scatter, it is clear that a correlation does exist, the more massive galaxies being, in general, gas poor. Our model results can fairly well describe the observations, accounting for both the slope and the dispersion of the $`M_g/L_H`$ vs $`V_C`$ relation.

In the lower panel of Fig. 3 we plot the ratio $`\mathrm{\Psi }/L_H`$ (current SFR to $`H`$band luminosity) vs. $`V_C`$. As discussed in several places (e.g. Kennicutt 1998a), this ratio is a measure of the parameter $`b=\mathrm{\Psi }/\mathrm{\Psi }`$, the ratio of the current SFR to the past average one. Massive galaxies display smaller $`\mathrm{\Psi }/L_H`$, and thus lower $`b`$ values, than their low mass counterparts. This means that they formed their stars at much higher rates in the past. Again, our results compare fairly well with the data, concerning both the slope and the scatter of the correlation. For a given $`V_C`$, large $`\lambda `$ discs have larger gas fractions today and were less active in the past than their lower $`\lambda `$ counterparts. As we shall see in Sec. 4.3, and discuss in Sec. 5, when the two panels of Fig. 3 are combined, they suggest quite convincingly that low mass spirals are on average younger than massive ones.

At this point, we notice that Gavazzi and Scodeggio (1996) already advanced a similar hypothesis to explain the observed colors of galaxies as a function of their H-band luminosity (a mesure of their dynamical mass according to them). Gavazzi et al. (1996) came to similar conclusions and suggested that mass is the main parameter of galaxy evolution, on the basis of multiwavelength observations concerning a variety of disc properties (colours, gas content, star formation rate, radius, surface brightness).

## 4 Star formation histories in discs

### 4.1 Formation time-scales

The Hubble sequence of galaxies is usually interpreted in terms of different star formation timescales (e.g. Kennicutt 1998a), although such an interpretation leaves unclear the role of the galaxy’s mass (Prantzos 2000). In order to put our results in that context, we performed an exponential fit to the star formation histories of our model galaxies (excluding the first 2 Gyr, where such a fit turns out to be inadequate). The resulting timescales are shown in Fig. 4. Notice that in some models the SFR is continuously increasing with time, resulting in negative caracteristic timescales. For that reason we present $`\tau _{EXP}^1`$, which has the advantage of varying continuously when going from a SFR increasing in time to one decreasing with time. For an exponential star formation rate, $`\mathrm{\Psi }exp^{t/\tau _{EXP}}`$, $`\tau _{EXP}^1`$ is equal to $`\dot{\mathrm{\Psi }}/\mathrm{\Psi }`$ (where $`\dot{\mathrm{\Psi }}=d\mathrm{\Psi }/dt`$) and can be considered as the normalized rate of change of the SFR. Obviously, the larger is $`\tau _{EXP}^1`$, the earlier the galaxy forms its stars.

As can be seen on Fig. 4, massive galaxies in our models form their stars on shorter timescales than their lower mass counterparts. This general trend is somewhat modulated by the spin parameter $`\lambda `$: discs with smaller $`\lambda `$ (i.e. more compact) have larger $`\tau _{EXP}^1`$ than discs with larger $`\lambda `$ of similar rotational velocity. The most massive discs of our simulations ($`V_C`$=360 km/s) have decreasing SFR with $`\tau _{EXP}`$4-5 Gyr. Discs with $`V_C`$200 km/s have long timescales, of the order of $``$10 Gyr, i.e. essentially constant SFR. Finally, low mass discs ( $`V_C`$100 km/s) have SFR increasing in time, with caracteristic timescales $`\tau _{EXP}^1`$-0.2 to -0.5 Gyr-1.

### 4.2 The evolution of the star formation efficiency

What is the reason for the vastly different SF timescales obtained in our models as a function of $`V_C`$? Is the overall SF efficiency $`ϵ=\mathrm{\Psi }/M_g`$ directly affected by the mass of the galaxy? An inspection of Eq. 1 suggests that this cannot be the reason. Indeed, at the caracteristic radius of the disc $`R_d`$, the local efficiency is:

$$ϵ(R_d)=\mathrm{\Psi }(R_d)/\mathrm{\Sigma }_g(R_d)=\alpha \mathrm{\Sigma }_g^{0.5}VR_d^1\lambda ^2V^{0.5}$$
(4)
since $`R_d\lambda V`$ from Eq. 2 and $`\mathrm{\Sigma }(R_d)\mathrm{\Sigma }_0\lambda ^2V`$ from Eq. 3. In other terms, the SF efficiency at $`R_d`$ varies very little with $`V_C`$ and depends much more on $`\lambda `$ than on $`V_C`$. Since the value of any intensive quantity (like SF efficiency) at $`R_d`$ is typical of the whole disc, it is expected that the global SF efficiency also depends little on $`V_C`$. In order to show this quantitatively, we plot in Fig. 5 the evolution of the global SF efficiency $`ϵ`$ of our models as a function of time. We show the results for 3 values of $`\lambda `$ (1/3, 1 and 3 times $`\lambda _{MW}`$) and three values of $`V_C`$ (80, 220 and 360 km/s, respectively). We compare our results at T=13.5 Gyr to estimates of $`ϵ=\mathrm{\Psi }/M_g`$ in our sample of nearby spirals (presented in Sec. 4.3) . The following points should be noted concerning our models:

i) The efficiency $`ϵ`$ depends very little on $`V_C`$, especially during the last half of the history of the galaxy.

ii) $`ϵ`$ is mainly determined by $`\lambda `$: “compact” galaxies have higher efficiency because of their smaller size ($`\mathrm{\Psi }1/R`$, for a given $`V_C`$) and larger gas surface density at the caracteristic radius $`R_d`$. At T=13.5 Gyr, the variation of our $`ϵ`$ values due to $`\lambda `$ can fully account for the dispersion in $`ϵ`$ measured in nearby spirals.

iii) The SF efficiency in compact galaxies (low $`\lambda `$) presents a peak at early times and then decreases. This happens because the star formation migrates to outer regions (due to the inside-out disc formation scheme) where the local SF efficiency is lower (because of lower surface densities and of the 1/R factor). In more extended galaxies (larger $`\lambda `$) $`ϵ`$ does not present such a decrease because local properties vary little with radius.

iv) For the lowest disc velocities $`V_C`$, $`ϵ`$ may increase with time. This is due to the adopted form of infall: the gas surface density increases considerably when the gas arrives finally in the disc (which may take a very long time in the case of the largest $`\lambda `$ values and lowest $`V_C`$ values, see Fig. 4).

The main point of this section is that in our models the SF efficiency $`ϵ`$ does not depend directly on the mass of the galaxy. The range of $`ϵ`$ span by our models during galactic evolution is due to $`\lambda `$ and $`\tau _{INF}`$, not to $`V_C`$; the most important of the two parameters is $`\tau _{INF}`$. Of course, $`\tau _{INF}`$ is adjusted to $`V_C`$, so that the observations of Fig. 2 and 3 (and many others, presented in Paper II) are reproduced. This is a crucial ingredient for the success of our models, and we discuss its implications in Sec. 5.

### 4.3 Gas fraction and star formation efficiency

The analysis of Sec. 4.1 and Sec. 4.2 lead to an important conclusion: the success of our models is to be interpreted in terms of mass-dependent SF timescales; but this is not due to any explicit dependence of SF efficiency on galaxy mass. Indeed, the SF efficiency of our models is virtually independent of mass during most of galactic history, and in particular at the present time. Is this supported by observations?

In Fig. 6 we display our data of Fig. 3, this time in a more “physical” presentation, appropriate for a quantitative discussion (see Sec. 5). The gas fractions $`\sigma _{gas}`$ (upper panels) and SF efficiency $`ϵ`$ (lower panels) are presented as a function of rotational velocity $`V_C`$ (left panels), H-magnitude (middle panels) and $`BH`$ colour (right panels).

The gas fraction $`\sigma _{gas}=M_g/M_T`$ (where the total mass $`M_T=M_{}+M_g`$ is the mass of stars + gas) is obtained by converting the $`L_H`$ luminosity to star mass $`M_{}`$ through the mass to light ratio $`M_{}/L_H`$ obtained in our models. This value is in the 0.3-0.6 range in solar units ($`M_{}/L_H`$0.35 for $`V_C`$=80 km/s, $``$0.47 for $`V_C`$ 220 km/s and $``$0.56 for $`V_C`$=360 km/s). Despite the uncertainties, a clear trend is present in the upper panel: the more massive, luminous and red a galaxy is, the smaller is its gas fraction. Notice that a similar trend is obtained by McGaugh and de Blok (1997) and Bell and de Jong (2000). Our model results are in excellent agreement with the data, although we have some difficulty in reproducing blue and gas poor discs (with $`BH<`$2.5 and $`\sigma _{gas}<`$0.4).

The SF efficiency of our observed galaxies is obtained by simply dividing the SFR $`\mathrm{\Psi }`$ with the gas mass $`M_g`$. We notice that the uncertainties in deriving $`ϵ`$ are rather large: those concerning the SFR $`\mathrm{\Psi }`$ are quite large (a factor of $``$3), whereas those of the gas mass at least of $`\pm `$20 %. As a result, the observationally derived scatter of $`ϵ`$, as appears on Figs. 5 and 6, is certainly larger that the real one. It is clear, however, that the SF efficiency does not seem to be correlated with the mass or the colours of spirals. The scatter in the observed values of $`ϵ`$ (a factor of $``$10), should be compared to the range of $``$60 span by galaxy mass (since the mass of the disc $`m_dV_C^3`$) or the range of five magnitudes in H-luminosity. Our model SF efficiencies $`ϵ`$ at T=13.5 Gyr, also displayed in Fig. 6 (lower panel, solid curves), show no dependence on mass or colour. They also reproduce the observed dispersion, the more compact discs (smaller $`\lambda `$) being the more efficient in turning their gas into stars.

The upper and lower panels of Fig. 6, when combined, point to an important conclusion: since the SF efficiency is independent of galactic mass (or luminosity), the fact that low mass galaxies have larger gas fractions today may only be due to their smaller ages. This is the most straightforward interpretation, independent of any theoretical considerations (except for the implicit assumption that the SF efficiency has remained $``$constant during the galaxy’s history). This conclusion is corroborated by an independent observable, namely that low mass galaxies are, in general, bluer than more massive ones. Notice that the latter observable concerns also elliptical galaxies, but the well known problem of the age-metallicity degeneracy does not allow to conclude in that case. In the case of spirals, the situation is even worse in principle, because colours may be affected by the presence of dust (presumably more abundant in massive spirals). Because of this complication, the observed gas fraction (smaller in large spirals) is not sufficient in itself to lift the age-metallicity degeneracy. However, when combined with the fact that the SF efficiency is independent of galactic mass (as argued here), the degeneracy is lift, and the aforementionned conclusion is naturally obtained.

Before turning to a more detailed discussion of our findings, we would like to point out that similar results are obtained in other recent works. For instance, McGaugh and de Blok (1997) find a clear trend between the gas fraction and the $`B`$magnitude of their galaxy sample, that they describe by the relation $`\sigma _{gas}`$=0.12 ($`M_B`$+23). This relationship is shown in Fig. 7 (upper panel, dashed line), along with the McGaugh and de Blok (1997) data for normal spirals and our data for normal spirals. It can be seen that there is a very good agreement between the two data sets. We notice that our models fit the observations of McGaugh and de Blok (1997) fairly well (see Fig. 13 of Boissier and Prantzos 2000).

In the lower panel of Fig. 7 we show the SF efficiency, both for our data set and the one of Kennicutt (1998b). Notice that Kennicutt (1998b) gives the average SFR surface density $`\mathrm{\Psi }`$ and average gas surface density $`\mathrm{\Sigma }_{gas}`$ of normal spirals, i.e. the integrated quantities are divided by the disc surface aerea (within the optical radius). Obviously, the ratio $`\mathrm{\Psi }/\mathrm{\Sigma }_{gas}`$ gives the overall SF efficiency $`ϵ`$ (since the disc aerea cancels out). As can be seen in Fig. 7, Kennicutt’s values of the SF efficiency are slightly larger than ours (by a factor of $`2`$), and have a smaller dispersion in the adopted logarithmic scale. Taking into account the various uncertainties in estimating the SFR from the data (see Sec. 3), such a discrepancy between Kennicutt’s results and ours is not unexpected. But the important point is that Kennicutt’s values are also independent of the galaxy $`B`$luminosity and, by virtue of the Tully-Fisher relation, on the galaxy’s mass.

The anti-correlation between gas fraction and luminosity in spirals was noticed by several authors (e.g. Gavazzi et al. 1996, McGaugh and de Blok 1997, Bell and de Jong 2000, Boselli et al. 2000, etc.). Based on two independent samples, we showed here that the SF efficiency of spirals is independent of their mass. The two findings combined, point to small discs being younger, on average, than massive ones. Our model, presented in Sec. 3.2, nicely explains these features (and several other, presented in papers II and III). However, one may argue that the complex interdependence between the adopted infall and SFR prescriptions makes a straightforward interpretation difficult; he/she may also argue that other types of models could also account for the observations and give a different interpreration for the same data (e.g. by invoking outflows). For that reason, we discuss in the next section this issue on the basis of simple analytical models of galactic chemical evolution and we make a very rough evaluation of the “ages” of the galaxies in our sample. Our purpose is not to derive the exact ages of the galaxies, but rather to have an order of magnitude estimate and, in particular, to check whether there is any trend of the derived ages with galactic mass.

## 5 Galactic Ages

In the framework of simple models of galactic chemical evolution adopting the Instantaneous Recycling Appproximation (IRA), one may obtain analytical solution for various quantities. In particular, provided that the Star Formation Rate $`\mathrm{\Psi }`$ is proportional to the gas mass $`M_g`$, i.e.

$$\mathrm{\Psi }=ϵM_g$$
(5)
one may obtain a relationship between the gas fraction $`\sigma _{gas}`$ and time T (assuming that the SF efficiency $`ϵ`$ is constant in time). The form of this relationship depends on further assumptions about the evolution of the system, i.e. on the possibility of allowing for gas flows in or outside the “box” (e.g. Pagel 1997).

Assuming that the galaxies of our sample have evolved as simple, homogeneous, “boxes”, we consider three possibilities: a “closed box” (all the gas is present from the very beginning), an “infall” model (where gas mass is continuously added to the system) and an “outflow” model (with gas continuously leaving the system). In the cases of gaseous flows, further assumptions about the corresponding flow rates are required in order to obtain analytical solutions. More specifically:

a) Closed Box: In that case, we have

$$\sigma _{gas}=\mathrm{exp}[(1R)ϵT]$$
(6)
where the return fraction $`R`$ accounts for the gas returned by stars to the interstellar medium; for the IMF of Kroupa et al. (1993) adopted here we have $`R`$0.32.

b) Infall: An analytical solution is easily obtained if it is assumed that the infall rate just balances the gas depletion due to star formation. As we have shown in Paper II, with detailed numerical models reproducing a large body of observational data, this situation describes rather well the largest period in the lifetime of most spiral galaxies. In that case, we have:

$$\sigma _{gas}=[1+ϵ(1R)T]^1$$
(7)
c) Outflow: Analytical solutions may be obtained by assuming that the outflow rate is proportional to the SFR: $`f_{out}=\gamma \mathrm{\Psi }`$. In that case, we have:

$$\sigma _{gas}=\frac{R1\gamma }{(R1)\mathrm{exp}[ϵ(R1\gamma )T]\gamma }$$
(8)
We shall assume here that the outflow rate is equal to the SFR, i.e. $`\gamma `$=1 (since, for higher outflow rates, we obtain ridiculously low galactic ages).

For each of the three scenarios, we shall consider five combinations between the observed gas fraction and SF efficiency, in order to derive galactic age T through Eqs. 6, 7 and 8. For the gas fraction vs. rotational velocity, we adopt the three curves of the upper panel of Fig. 8, corresponding roughly to the mean trend and the upper and lower bounds of the observations, respectively. For each of those three curves, we derive the corresponding galactic ages by using a SF efficiency $`ϵ`$=0.3 Gyr-1 (the mean value of the observations in Fig. 6 and 7); in that way, we obtain the shaded regions in Fig. 8, with the thick curve corresponding to the mean trend, the largest ages to the lowest gas fractions and vice versa. We perform two more calculations. For the first, we adopt the high gas fractions (uppermost curve in the upper panel of Fig. 8), combined to a low SF efficiency $`ϵ`$=0.1 Gyr-1 (a reasonable lower bound for the observations of Fig. 6 and 7); in that way we obtain the dotted curves in Fig. 8. Finally, we adopt the low gas fractions (lower curve in the upper panel of Fig. 8), combined to a high SF efficiency $`ϵ`$=1 Gyr-1 (a reasonable upper bound for the observations of Fig. 6 and 7); in that way we obtain the dashed curves in Fig. 8.

An inspection of the results, plotted in the three lower panels of Fig. 8, shows that the derived “effective” galactic age is a monotonic function of rotational velocity $`V_C`$ in all cases. The absolute ages depend, of course, on the adopted gas fractions, SF efficiencies and assumptions about gaseous flows. For a given gas fraction and SF efficiency, the infall model leads to the largest ages; since the gas is constantly replenished, it takes more time to attain a given gas fraction than in a closed box (starting with $`\sigma _{gas}`$=1). Also, the outflow model produces the lowest ages; since part of the gas is constantly removed, it takes less time to reach a given gas fraction than in a closed box. The ages derived in the framework of the outflow model are, in general, too low: they are lower than 6 Gyr for all galaxies with $`V_C<`$200 km/s. It transpires that galaxies in this velocity range should not have suffered extensive mass losses, i.e. they must not have lost during their lifetime an amount of gas as important as their stellar content (since we adopted $`f_{out}=\mathrm{\Psi }`$ here).

The mean age values in the case of the infall model are 4-5 Gyr for the small discs (in the $`V_C`$100 km/s range) and $``$10 Gyr for $`V_C`$220 km/s. For discs with $`V_C>`$300 km/s, extremely large ages ($`>`$15 Gyr) are found; however, the approximation of evolution at constant gas mass is certainly not valid in that case (see Fig. 4 in Boissier and Prantzos 2000). We notice that Bell and de Jong (2000), on the basis of a different sample and with a completely different method (based on a photometric estimate of the ages) find an “effective age” of $``$10 Gyr for the more luminous galaxies ($`M_K`$-26) and $``$6 Gyr for the less luminous ones ($`M_K`$-20); the galaxies of our sample span also this luminosity range, and its interesting to see that similar ages are found with completely independent methods.

The results obtained in this section confirm the suggestion of Sec. 3.2: low mass discs are, on average, younger than massive ones. The only way to reverse this trend is by assuming that the SF efficiency is strongly correlated with disc mass. Despite the large scatter in the observational data, it is clear that such a correlation does not exist, at least at the present time; and, since the observed galaxies span a large range in masses and metallicities, it does not seem plausible that such a correlation ever existed in the past. In the case of outflow models, another possibility would be to consider that the more massive discs suffered more important mass losses (leading to low gas fractions without having to invoke large ages). However, such a hypothesis is incompatible to the fact that massive discs have deeper potential wells and are less prone to outflows than low mass ones.

## 6 Summary

In this work we investigate the properties of the star formation efficiency of spiral galaxies and study the implications for their evolution. We use a large homogeneous sample of disc galaxies, for which we measure gaseous mass $`M_g`$ (HI), rotational velocity $`V_C`$, star formation rates $`\mathrm{\Psi }`$ and luminosities in the $`B`$ and $`H`$ bands. We are then able to derive the corresponding gas fractions $`\sigma _{gas}`$ and SF efficiencies $`ϵ=\mathrm{\Psi }/M_g`$ as a function of $`V_C`$, $`H`$luminosity and $`BH`$ colour index. We find that the gas fraction is correlated to $`V_C`$, $`H`$luminosity and $`BH`$, in the sense that more massive, luminous and redder discs have smaller gas fractions. Previous work by McGaugh and de Blok (1997) reached similar conclusions. The main finding of this work is that the SF efficiency does not correlate with any of the galaxy properties; despite a rather large dispersion (within a factor of $``$10), the observed $`ϵ`$ is independent of $`V_C`$, $`L_H`$ or $`BH`$.

We interpret our data in the framework of detailed models of galactic chemical and spectrophotometric evolution, utilising metallicity dependent stellar lifetimes, yields, tracks and spectra. These models are calibrated on the Milky Way disc (paper I) and use radially dependent star formation rates, which reproduce observed gradients (paper III). They are extended to other spirals in the framework of Cold Dark Matter scenarios for galaxy formation, and are described by two parameters: rotational velocity $`V_C`$ and spin parameter $`\lambda `$. As in our previous work (paper II), we find good agreement with the observations, provided a crucial assumption is made: massive discs are formed earlier than less massive ones. With this assumption our models reproduce the observed trends of gas fractions and SF efficiencies vs. $`V_C`$, while variations due to $`\lambda `$ account for the observed dispersion in both cases.

It is important to notice that the dependence of age on galactic mass that we find is not due to any mass-dependent SF efficiency, only to disc formation timescales; in our models this is achieved by varying the infall timescales. Both observations and models suggest that the SF efficiency is independent of galaxy properties. Since the observed galaxies cover a wide range of masses, colours and metallicities, there is no reason to suppose that the SF efficiency was different in the past. The adopted SFR prescription in our models also results in a very slowly varying SF efficiency with time.

When the observed relations of gas fraction vs $`V_C`$ and SF efficiency vs $`V_C`$ are considered in combination, they convincingly suggest that low mass discs are, on average, younger than more massive ones; this conclusion is independent of any model and the only assumption is that the SF efficiency is constant in time (a quite plausible assumption, as argued above). In the framework of simple analytical models of galactic chemical evolution, we evaluate the “effective ages” for the galaxies of our sample, using the observationally derived values of $`\sigma _{gas}`$ and $`ϵ`$. We find that even models with modest outflows (with ejected masses equal to the stellar ones) lead to ridiculously low values for the galaxy ages; our conclusion is that galaxies in the range $`V_C`$80-400 km/s have not suffered extensive mass losses. Closed box models and infall models lead to more plausible values for the effective ages. In particular, infall models lead to ages of $``$4-5 Gyr for discs of $`V_C`$100 km/s and $``$8-10 Gyr for $`V_C`$200 km/s. These “chemically derived” ages are in fair agreement with those derived on the basis of our more sophisticated numerical models, which fit a much larger body of observational data for low redshift spirals. Most importantly, they are also in fair agreement with the “photometric” ages derived in a completely independent way and with a different sample by Bell and de Jong (2000).

In summary, our data, taken at face value, suggest that the bulk of stars in more massive discs are older than in less massive ones. This is supported by our detailed numerical models of galactic chemical and photometric evolution, but also by recent, independent, analysis (Boselli et al. 2000). We notice that Bell and deJong (2000) conclude that it is local surface density that mainly drives the star formation history, while mass plays a less important role. Our study suggests that mass is the main factor, while local surface density plays only a minor role (through the spin parameter $`\lambda `$: lower $`\lambda `$ values lead to higher local surface densities for a given rotational velocity $`V_C`$).

We notice that this picture is hardly compatible with the currently popular “paradigm” of hierarchical galaxy formation, which holds that large discs are formed by merging of small units at relatively late epochs. If this were the case, massive discs should have large SF efficiencies, in order to have their gas fractions reduced to lower levels than their less massive counterparts. However, such an enhanced SF efficiency is not supported by observations.


format: text

text:
# Gravitational Experiment Below 1 Millimeter and Comment on Shielded Casimir Backgrounds for Experiments in the Micron Regime

## 1 Introduction

Experimental searches for new macroscopic forces have only marginally explored the distance range under 1 mm and there is little knowledge of gravity itself in this range. The sub-millimeter region is of profound and rapidly increasing experimental interest, given a number of recent predictions of new forces from a rich variety of modern theories of fundamental interactions.

The existing experimental limits on new forces at short distances are defined by classical gravity measurements and Casimir force measurements, as shown in Fig. 2. From the figure, in which the strength $`\alpha `$ of a hypothetical new force relative to gravity is plotted versus the Yukawa range $`\lambda `$, the published experimental limits allow for forces in nature several million times stronger than gravity over ranges as great as 100 $`\mu `$m.

Also shown in Fig. 2 are theoretical predictions of new effects in this regime. Most notable is the line indicating Yukawa corrections to the inverse square law which arise from compact extra dimensions. Corrections are also predicted from massive scalars in string theories, such as Moduli and Dilatons. The other predictions shown are motivated by the Cosmological Constant Problem and the Strong CP Problem of QCD.

Since the summer of 1997, we have been constructing an experiment to explore the sub-centimeter regime. The approach uses a high-frequency technique, departing from the torsion balance method which was used in all measurements defining the current limits.

## 2 Status of the Experiment

### 2.1 Design Overview

If all dimensions of experimental test masses of a gravitational experiment are scaled by the same factor $`D`$, Newtonian attraction varies as $`D^4`$. At small separations, signals from mass–coupled forces become very weak. At the same time, background effects such as surface forces increase rapidly.

Our experimental approach is shown in Fig. 1 (see Ref. References and the last paper of Ref. References for more detail). Planar test mass geometry is chosen to concentrate as much mass as possible at the scale of interest. The test masses consist of 200 $`\mu `$m thick tungsten oscillators. The detector is driven by the source mass on resonance near 1 kHz. At this frequency it is possible to construct simple, passive vibration isolation stacks sufficient to suppress the acoustic coupling of the test masses though the apparatus. Detector oscillations are read out with a capacitive transducer and lock-in amplifier. The entire apparatus is enclosed in a vacuum chamber and operated at $`10^7`$ torr to reduce the acoustic coupling between the test masses though the residual gas in between them.

The principal backgrounds in addition to the acoustic backgrounds arise from electrostatic and magnetic forces. Electrostatic forces, as well as acoustic coupling through the residual gas, are suppressed with a stiff conducting shield, consisting of a 75 $`\mu `$m thick quartz window plated with gold, suspended between the test masses. Magnetic backgrounds have so far been avoided with the exclusive use of non-magnetic materials for the construction of the apparatus. If the need arises, in-situ imaging of magnetic contaminants on the test masses is possible.

If all of these backgrounds can be sufficiently reduced, the limiting background of this experiment will be thermal noise due to dissipation in the detector mass. The experiment is therefore being carried out in two phases. A preliminary, room–temperature version of the experiment is now in operation and nearing the practical limit of its sensitivity. The final version of the experiment will be cooled to liquid Helium temperatures for greater sensitivity.

### 2.2 Current Sensitivity

A signal from a new effect can be modeled as the force on the detector mass due to a Yukawa interaction with the source mass. This force is given by:

$$F_Y(t)=2\pi \alpha G\rho _s\rho _dA\lambda ^2\mathrm{exp}(d(t)/\lambda )[1\mathrm{exp}(t_s/\lambda )][1\mathrm{exp}(t_d/\lambda )],$$
(1)
where $`d(t)`$ is the separation between the test masses, $`\rho _s`$ and $`\rho _d`$ are the source and detector mass densities, and $`t_s`$ and $`t_d`$ are the thicknesses. This would be an exact expression if either plate had area $`A`$ and the other had infinite area, but for the real geometry there are small edge corrections, which we do not consider. For values of these parameters typical for our experiment, this force is approximately $`2\times 10^{14}`$ N, for $`\alpha =1`$ and $`\lambda `$ = 100 $`\mu `$m.

The rms thermal noise force is found from the mechanical Nyquist theorem to be:

$$F_T=\sqrt{\frac{4kT}{\tau }\left(\frac{m\omega }{Q}\right)},$$
(2)
where $`m`$ is the mass of the small rectangular section of the detector oscillator, $`\omega `$ is the resonant frequency, $`Q`$ is the detector quality factor, $`T`$ is the temperature and $`\tau `$ is the measurement integration time. For typical $`Q`$ values of our tungsten detector mass ($`2.5\times 10^4`$), a temperature of 300 K, and an integration time of 1000 s, this force (and hence the current sensitivity of our experiment) is about $`4\times 10^{14}`$ N.

In April 2000, the room temperature experiment became fully operational for the first time. No signal was observed above the detector thermal noise over the entire integration time of 1800 s. Setting the ratio of Eqs. 1 and 2 to unity and solving for $`\alpha `$, we infer the limits on the Yukawa parameters as shown in Fig. 2.

### 2.3 Recent Improvements

For the April run, thermal drifts of the detector resonance limited the integration time to 1800 s. Furthermore, a poor match between source and detector resonant frequencies limited the source mass amplitude to roughly 0.5 $`\mu `$m. Since that run we have installed temperature control, allowing for much longer integration times. We have tuned the source mass such that amplitudes greater than 50 $`\mu `$m on resonance are now possible. Also, we have fabricated a new shield from 50 $`\mu `$m thick sapphire, making smaller separations possible. With these improvements, we expect the room temperature experiment to attain the sensitivity represented by the upper dashed curve in Fig. 2.

The final cryogenic version of the experiment is under development. Assuming the backgrounds can be controlled, we expect this experiment to reach gravitational strength at 50 $`\mu `$m, as represented by the lower dashed curve in Fig. 2.

## 3 Shielded Casimir Background in the Micron Regime

Session PT6 at MG9 featured discussions of experimental searches for new forces in the sub-micron regime, where the Casimir background is expected to dominate. A series of Casimir force measurements has been carried out in this regime using an AFM, and limits on new forces down to 10 nm have been derived from these experiments. So far, the limits are many orders of magnitude stronger than gravity, motivating the design of dedicated searches for new forces. One proposal in this regime is the “iso-electronic” technique of Fischbach and Krause.

The Casimir background is not expected to be significant in our existing experiment. However, an extension of the shielding technique we employ may be useful for AFM–type searches in experiments near 1 $`\mu `$m.

Our idea is illustrated in Fig. 3. A probe at a distance $`D`$ above a flat gold sample is scanned over regions of alternating depth $`D`$ and 2$`D`$. The horizontal dimensions of all probe and sample features are taken to be much larger than $`D`$. The gold sample is backed by a low–density dielectric substrate which for the purposes of this study we take to be vacuum.

Due to the finite penetration depth of the sample metal, above a certain scale $`D`$ we expect the Casimir force between the probe and either thickness region of the sample to be essentially equal. While mass–coupled forces between probe and sample may be much more feeble, for them we do not expect a similar equalization across the thickness boundaries.

We have calculated the Casimir forces for the geometry in Fig 3. We use the method developed by A. Lambrecht and S. Reynaud, ignoring the effects of surface roughness and finite temperature. We employ a Drude model for the gold probe and sample, using a plasma frequency of $`\omega _P=1.4\times 10^{16}`$ rad/s and relaxation parameter of $`\gamma =5.3\times 10^{13}`$ rad/s in Eq. 16 of Ref. References.

Setting the plate separation $`L`$ in Eq. 9 of Ref. References equal to $`D`$, we then use Eqs. 12, 9, and 4 of that reference to calculate the difference in the Casimir force between the probe and sample regions of thickness $`D`$ and 2$`D`$: $`\mathrm{\Delta }F_C=F_C^{2D}F_C^D`$. We assume the probe to have infinite thickness. As a preliminary check, we set the sample thickness equal to infinity and note that our results for the Casimir reduction factor $`\eta _F`$ agree with those of Fig. 6 of Ref. References to within 5%.

We then compute the difference in the Yukawa force for the same geometry, using $`\alpha =1`$ and $`\lambda =D`$: $`\mathrm{\Delta }F_Y=F_Y^{2D}F_Y^D`$. The ratio $`\mathrm{\Delta }F_Y/\mathrm{\Delta }F_C`$ is shown in Fig. 4 as a function of $`D`$ for $`10^8`$ m $`<D<5\times 10^6`$ m. We note that these results change by less than 1% as we vary the limits of integration \[$`10^8`$$`10^{20}`$ rad/s\] by an order of magnitude.

The steep rise in the curve occurs when the scale $`D`$ approaches the plasma wavelength we have used for gold (about $`1.4\times 10^7`$ m), at which point the shielding increases rapidly. For the idealized geometry and conditions we have assumed, Fig. 4 implies that a gravitational–strength Yukawa force becomes distinguishable from the Casimir background as the probe is scanned across the thickness boundary at a scale $`D`$ of about 3 $`\mu `$m.

## Acknowledgments

We wish to thank Hilton Chan for his extensive work on this project. J. L. thanks A. Lambrecht and S. Reynaud for useful discussions.


text:
# Statistics of two-particle dispersion in two-dimensional turbulence

## I Introduction

Understanding the statistics of particle pairs dispersion in turbulent velocity fields is of great interest for both theoretical and practical implications. At variance with single particle dispersion which depends mainly on the large scale, energy containing eddies, pair dispersion is driven (at least at intermediate times) by velocity fluctuations at scales comparable with the pair separation. Since these small scale fluctuations have universal characteristics, independent from the details of the large scale flow, relative dispersion in fully developed turbulence is expected to show universal behavior . From an applicative point of view, a deep comprehension of relative dispersion mechanisms is of fundamental importance for a correct modelization of small scale diffusion and mixing properties.

Since the pioneering work by Richardson , many efforts have been done to confirm experimentally or numerically his law . Nevertheless, the main obstacle to a deep investigation of relative dispersion in turbulence remains the lack of sufficient statistics due to technical difficulties in laboratory experiments and to the moderate inertial range reached in direct numerical simulations.

In this Paper we present a detailed investigation of the statistics of relative dispersion from extensive direct numerical simulations of particle pairs in two-dimensional Navier-Stokes turbulence. We will see that the main ingredient of the original Richardson description, i.e. Richardson diffusion equation, is sufficient for rough description of relative dispersion in this flow. Nevertheless, our simulations show that two-particle statistics is rather sensible to finite size effects. This asks for a different analysis based on doubling time statistics which has recently introduced for the analysis of relative dispersion . Comparison of numerical results with ones based on the Richardson’s equation show that the last delivers a qualitatively good description of the doubling-time-distributions. The quantitative deviations found are attributed to the fact that the dispersion process is not purely diffusive and is influenced by ballistic (persistent) motion.

The article is organized as follows: In Section II we discuss the Richardson’s approach to the two-particle dispersion, in Section III the fixed-scale properties of dispersion process (such as doubling-time statistics) are considered. The numerical approach and the results of simulations are discussed in Section IV. Section V is devoted to conclusions. The mathematical details of calculations of doubling-time statistics for the Richardson’s case are given in Appendices A and B.

## II Statistics of relative dispersion

Relative dispersion in turbulence is often phenomenologically described in terms of a diffusion equation for the probability density function of pair separation $`p(𝐫,t)`$

$$\frac{p(𝐫,t)}{t}=\frac{}{r_j}\left(K(𝐫,t)\frac{p(𝐫,t)}{r_j}\right)$$
(1)
with a space and time dependent diffusion coefficient $`K(𝐫,t)`$ . The original Richardson proposal obtained from experimental data in the atmosphere corresponds to $`K(r,t)=K(r)=k_0\epsilon ^{1/3}r^{4/3}`$ , where $`\epsilon `$ has the dimension of energy dissipation (see below) and $`k_0`$ is an adimensional constant. In the $`d`$-dimensional isotropic case, this diffusion equation takes the form

$$\frac{p(𝐫,t)}{t}=\frac{1}{r^{d1}}\frac{}{r}r^{d1}K(r)\frac{p(𝐫,t)}{r}$$
(2)
Its solution leads to the well known non-Gaussian distribution

$$p(𝐫,t)=\frac{A}{(k_0t)^3\epsilon }exp\left(\frac{9r^{2/3}}{4k_0\epsilon ^{1/3}t}\right)$$
(3)
where $`A`$ is a normalizing factor. The growth of pair separation is described in terms of a single exponent as

$$R^{2n}(t)r^{2n}(t)=C_{2n}\epsilon ^nt^{3n}.$$
(4)
The so called Richardson constant in (4) is $`g=C_2=\frac{1280}{243}k_0^3`$.

The Richardson’s conjecture was formulated based on the scaling nature of the diffusion coefficient and in analogy with diffusion. No information about the nature of the pdf was available by that time. Any chose of the form $`K(r,t)r^{4/3\alpha }r^2(t)^{\alpha /2}`$ (i.e. $`K(r,t)r^{4/3\alpha }t^{3\alpha /2}`$) would give the same scaling law $`R^2t^3`$ but with different pdf’s (see ).

The possibility to describe the dispersion process by means of a diffusion equation is based on essentially two important physical assumptions which can be verified a posteriori. The first one is that the dispersion process is self-similar in time, which is probably true in non intermittent velocity field , the second one is that the velocity field is short correlated in time . In the limit of velocity field $`\delta `$-correlated in time the diffusion equation (1) becomes exact . As we proceed to show, the Richardson’s conjecture (2), which is exact under small values of the persistence parameter of the flow, vide infra and , still delivers a qualitatively good approximation for realistic 2d turbulent flows, whose persistence parameter of the order of 1.

Richardson scaling in turbulence is a consequence of Kolmogorov scaling for the velocity differences . Under Kolmogorov scaling, the mean square relative velocity and the correlation time in the inertial range are given by

$$\delta v(r)^2=v_0^2\left(\frac{r}{r_0}\right)^{2/3}\epsilon ^{2/3}r^{2/3}$$
(5)
and

$$\tau (r)=\tau _0\left(\frac{r}{r_0}\right)^{2/3}\epsilon ^{1/3}r^{2/3},$$
(6)
where $`r_0`$, $`\tau _0`$ and $`v_0`$ are some (large scale) characteristic length, time and velocity scale and $`\epsilon v_0^2/\tau _0`$ is the energy flux in the inertial range. The value of the dimensionless combination $`Ps=\frac{v_0\tau _0}{r_0}`$ remains however unspecified by scaling considerations. It is referred to as a persistence parameter of the flow and plays a central role in describing single particle diffusion and pair separation .

The persistence parameter gives the ratio of the velocity correlation time to the Lagrangian characteristic time. In order to see how $`Ps`$ influences Lagrangian dispersion, let us consider the following simple model, which has been used as a basis for building a stochastic model of turbulent dispersion . We take that the magnitude of the separation velocity is a function of $`r`$ only so that $`\delta v(r)=v_0(r/r_0)^{1/3}`$. The temporal changes of the flow can be accounted for by letting the particle change its velocity direction from time to time, while keeping the velocity’s magnitude constant. The probability that the relative velocity of particle separation changes its direction during time interval $`dt`$ can be assumed to be $`dp=dt/\tau (r)`$. The growth of the magnitude of the interparticle separation $`𝐫(t)`$ in $`dt`$ is $`dr\delta v(r)dt`$, thus the probability to change the direction of velocity within $`dr`$ is, using (5-6)

$$dp=p(r)dr=\frac{dr}{\delta v(r)\tau (r)}=\frac{r_0}{v_0\tau _0}\frac{dr}{r}=\frac{1}{Ps}\frac{dr}{r}.$$
(7)
The distribution of the position of turning points in the separation follows from (7) . The conditional probability density to find a next turning point at $`r_2`$ provided a previous one was at $`r_1<r_2`$ is given by

$$\mathrm{\Psi }(r_2|r_1)=\frac{1}{Psr_1}\left(\frac{r_2}{r_1}\right)^{1/Ps1}.$$
(8)
Note that the dependence of $`\mathrm{\Psi }(r_2|r_1)`$ only on the relative positions of the turning points, i.e. on $`r_2/r_1`$, is a clear consequence of scaling.

The tail of $`\mathrm{\Psi }(r_2|r_1)`$ decides about the existence of the second moment of this distribution, i.e. on the fact whether the corresponding motion is short- or long-range correlated in space. Depending on the persistence parameter $`Ps`$, the dispersion can be either diffusive ($`Ps1`$) or ballistic ($`Ps1`$) in nature. Note that the power-law tail of the distributions make the problem extremely sensitive to the finite-size effects, especially for large $`Ps`$, when the weights of ballistic events (Lévy-walks ) is considerable.

We note that the value of $`Ps`$ is not a free parameter, but is fixed for a given physical situation. The scaling nature of turbulence supposes that this parameters is a constant, depending only on general properties of the flow, e.g. on its 2d or 3d nature (in this last case also the overall geometry of the flow can be of importance). On the other hand, since the nature of dispersion process depends crucially on the value of $`Ps`$, the only way for getting quantitative information about the dispersion is through direct numerical simulations. The strong finite-size effects in relative dispersion simulations require the introduction of quantities which are less sensitive to finite resolution.

## III Exit time statistics

All the results discussed in the previous Section can be observed only in high Reynolds number flows in which the inertial range, where the scaling laws (5,6) hold, is sufficiently wide. The needs for large Reynolds is particularly severe in the case of Richardson dispersion, as a consequence of the long tails in (3). Moreover, the observation of time scaling laws as (4) requires sufficiently long times in order to forget the initial separation $`r(0)`$ .

For these reasons, the observation of Richardson scaling (i.e. (3) or (4)) is very difficult in direct numerical simulations where the Reynolds number is limited by the resolution. The same kind of limitations arise in laboratory experiments, as a consequence of the necessity to follow the Lagrangian trajectories which limits again the Reynolds number .

To partially overcome these difficulties, an alternative approach based on doubling time (or exit time) statistics has been recently proposed . Given a set of thresholds $`R_n=\rho ^nR(0)`$ within the inertial range, one computes the “doubling time” $`T_\rho (R_n)`$ defined as the time it takes for the particle pair separation to grow from threshold $`R_n`$ to the next one $`R_{n+1}`$. Averages are then performed over many dispersion experiments, i.e., particle pairs, to get the mean doubling time $`T_\rho (R)`$. The outstanding advantage of this kind of averaging at fixed scale separation, as opposite to a fixed time, is that it removes crossover effects since all sampled particle pairs belong to the same scales.

The problem of doubling time statistics is a first-passage problem for the corresponding transport process. For the Richardson case, in 2d it is given by the solution of the Richardson’s diffusion equation, Eqs. (2), with initial condition $`p(𝐫,0)=\delta (rR/\rho )/2\pi `$ and absorbing boundary at $`r=R`$ (so that $`p(R,t)=0`$). The pdf of doubling time can be obtained as the time derivative of the probability that the particle is still within the threshold

$$p_D(t)=\frac{d}{dt}_{|𝐫|<R}p(𝐫,t)𝑑𝐫.$$
(9)
Using (2) one obtains

$$p_D(t)=2\pi \epsilon ^{1/3}k_0R^{4/3}\frac{p(𝐫,t)}{r}|_{r=R}.$$
(10)
The solution using the eigenfunction decomposition is given in Appendix A and shows that the long-time asymptotic of $`p_D(t)`$ is exponential:

$$p_D(t)\mathrm{exp}(\kappa k_0\epsilon ^{1/3}R^{2/3}t)$$
(11)
where $`\kappa 2.93`$ is a number factor. This exponential nature of the tail of $`p_D(t)`$-distribution will be confirmed by direct simulations in Section IV.

Note that the combination $`\epsilon ^{1/3}R^{2/3}`$ has a dimension of time and is proportional to the average doubling time $`T_\rho (R)`$. This time can be obtained by a simple argument reported in the Appendix B. In the $`2`$-dimensional case one obtains

$$T_\rho (R)=\frac{3}{4}\frac{\rho ^{2/3}1}{\epsilon ^{1/3}\rho ^{2/3}}\frac{R^{2/3}}{k_0}$$
(12)
Prediction (12) contains the parameter $`k_0`$ which is, as shown in Section II, is equivalent to the Richardson constant $`g`$. As a consequence, the computation of average doubling time can be used for an alternative (and more robust, as we will see) estimation of the Richardson constant. It is convenient to rewrite the doubling time pdf (11) in terms of the average doubling time $`T_\rho (R)`$. Making use of (12) one obtains in 2d

$$p_D(t)\mathrm{exp}(0.252t/T(R)).$$
(13)
which is a parameterless, universal function.

## IV Direct numerical simulations

Pair dispersion statistics has been investigated by high-resolution direct numerical simulations of the inverse energy cascade in two-dimensional turbulence . There are several reasons for considering 2D turbulence. First of all, the dimensionality of the problem makes feasible direct numerical simulations at high Reynolds numbers. Moreover, the observed absence of intermittency  makes the 2D inverse energy cascade an ideal framework for the study of Richardson scaling in Kolmogorov turbulence.

The 2D Navier-Stokes equation for the vorticity $`\omega =\times 𝐯=\mathrm{\Delta }\psi `$ is:

$$_t\omega +J(\omega ,\psi )=\nu \mathrm{\Delta }\omega \alpha \omega +\varphi ,$$
(14)
where $`\psi `$ is the stream function and $`J`$ denotes the Jacobian. The friction linear term $`\alpha \omega `$ extracts energy from the system to avoid Bose-Einstein condensation at the gravest modes . The forcing is active only on a typical small scale $`l_f`$ and is $`\delta `$-correlated in time to ensure the control of the energy injection rate. The viscous term has the role of removing enstrophy at scales smaller than $`l_f`$ and, as customary, it is numerically more convenient to substitute it by a hyperviscous term (of order eight in our simulations). Numerical integration of (14) is performed by a standard pseudospectral method on a doubly periodic square domain at resolution $`N=1024`$ or $`N=2048`$. All the results presented are obtained in conditions of stationary turbulence.

In Figure 1 we plot the energy spectrum, which displays Kolmogorov scaling $`E(k)=Cϵ^{2/3}k^{5/3}`$ over about two decades with Kolmogorov constant $`C6.0`$. In the inset we plot the third order structure function $`S_3(r)=\delta v(r)^3`$ compensated with the theoretical prediction $`S_3(r)=3/2\epsilon r`$. The observation of the plateau confirms the existence of an inverse energy cascade and indicates the extension of the inertial range. Previous numerical investigation has shown that velocity differences statistics in the inverse cascade is not affected by intermittency corrections . In this case we may expect the Lagrangian statistics to be self-similar with Richardson scaling .

Lagrangian statistics is obtained by integrating the trajectories of $`64000`$ particle pairs in the turbulent velocity field, initially uniformly distributed with constant separation $`R(0)`$.

The Lagrangian data reported below are in dimensionless units in which separations are rescaled with the box size $`L=2\pi `$ and time with the large scale time $`T_0=(L^2/\epsilon )^{1/3}`$.

### A Relative dispersion analysis

In Figure 2 we plot the relative dispersion $`R^2(t)`$ for two different initial separation, $`R(0)=\delta x/2`$ and $`R(0)=\delta x`$ (where $`\delta x=2\pi /N`$ is the grid mesh). The Richardson $`t^3`$ law (4) is observed in a limited time interval, especially for the larger $`R(0)`$ run. It is remarkable that the relative separation law displays such a strong dependence on the initial conditions even in our high resolution runs.

This dependence makes the determination of the Richardson constant particularly difficult. In Figure 2 we also show the compensated plot $`R^2(t)/t^3`$ which, in the adimensional units, should directly give the constant $`g`$. It is clear that a precise determination of $`g`$ is impossible; even the Richardson scaling (4), when looked in a compensated plot, is rather poor. Figure 2 suggests that starting with an intermediate initial separation would give a wider scaling range. Of course, one would like to avoid this “fine tuning”, which is probably impossible to implement in the case of experimental data.

The probability distribution function of pair separations is plotted in Figure 3 for the $`R(0)=\delta x/2`$ run. At short time $`t=0.015`$, in the beginning of the $`t^3`$ range in Figure 2, we found that the Richardson pdf (3) fits pretty well our data, although some deviations can be detected. Of course, at time comparable with the integral time $`t=0.77`$, particle separations are of the order of the integral scale and we observe Gaussian distribution. The crossover between these two regimes is extremely broad: deviations from Richardson pdf are clearly seen already for the times well within the Richardson’s $`t^3`$ range. To observe better this transition, in Figure 4 we plot, in log-log plot, the right tail of $`\mathrm{ln}(p(r,t)/p(0,t))`$. We see that the tail slope can be fitted with an exponent $`\alpha `$ which change continuously in time, from $`2/3`$ to the Gaussian value $`2`$ (see the inset of Figure 4). Thus self similarity, if it exists, is reduced to the very short time at the beginning of dispersion. Moreover, the scaling region is strongly affected by the choice of initial separation, as shown in Figure 2.

### B Doubling time data

From the same Lagrangian trajectories discussed in the previous Section, we have computed the doubling time statistics. The average doubling time as function of the scale is plotted in Figure 5. The line represent the dimensional prediction $`T(R)R^{2/3}`$. In the inset we plot the quantity $`\frac{20}{9}\frac{(\rho ^{2/3}1)^3}{\rho ^2}\frac{R^2}{\epsilon T^3}`$ which, from (4) and (12) gives the value of the Richardson constant $`g3.8`$, which is compatible with the estimation from Figure 2. It is interesting to compare our result with previous estimations of $`g`$. The only experimental estimation of $`g`$ for 2D inverse cascade  gives a value about $`7`$ times smaller (but the Reynolds number in the experiment is even smaller than in present simulations). Other estimations of $`g`$ are based on kinematic simulations with synthetic flows. In all these cases  the reported values are even smaller. On the other hand, our numerical finding is not very far from the prediction of turbulence closure theory .

From Figure 5 we observe that at very small separations $`R10^3`$, the doubling time has a tendency to a constant value $`T(R)0.0016`$. On these scales we are below the forcing scale (see Figure1), and the velocity field can be assumed smooth. As a consequence of Lagrangian chaos we expect on these scales an exponential amplification of separations  at a rate given by the Lagrangian Lyapunov exponent $`\lambda `$. The latter can be obtained as $`\lambda =lim_{R0}\mathrm{ln}\rho /T(R)`$  and gives $`\lambda 110`$ (in adimensional units). The Lagrangian Lyapunov exponent $`\lambda `$ is a small scale quantity (i.e. depends on the Reynolds number of the simulation), and thus has to be compared with a small scale characteristic time. One can estimate the smallest characteristic time $`\tau _{min}`$ by the minimum value of $`(k^3E(k))^1`$. We obtain $`\lambda 0.23\tau _{min}^1`$.

The comparison of Figures 2 and 5 shows the advantages of fixed scale statistics with respect to fixed time statistics. In particular, the dependence of $`R(t)`$ and the Richardson constant on the initial separation are absent in fixed scale statistics $`T(R)`$.

In Figure 6 we plot the doubling time pdf $`p_D(T)`$ compensated with the mean value $`T(R)`$ at different scales in the inertial range $`0.003R0.046`$. First, we obtain a very nice collapse of the different curves, indicating that the process is really self-similar. Second, we observe the exponential tail predicted in Section III with a fitted coefficient $`0.3`$ which is indeed not far from the theoretical prediction (13) based on the Richardson’s picture. The difference between the predicted and measured values of the prefactors is not large, but perceptible: it shows that the Richardson’s equation gives a correct qualitative description of the dispersion process, but is not exact. The reasons for deviations from the diffusive picture proposed by Richardson are the long-range correlations in the particles’ motion, as seen from the analysis of the turning points of their relative trajectories.

### C Turning points statistics and the persistence parameter

A possible explanation for the deviations of our high resolution numerical data from the Richardson’s picture is the not too small value of the persistence parameter. As discussed in Section II, at large values of $`Ps`$ the contribution of ballistic events may lead to non-Richardson distributions and moreover makes the dispersion strongly sensible to finite-size effects, cutting the longer trajectories.

We have computed the persistence parameter making use of (8). We have recorded, for each pair, the set of turning points $`r_i`$ at which the pair’s relative velocity changes sign. From the set of $`r_i`$ we have then computed the pdf of the ratio $`r_{i+1}/r_i`$, accumulating for all the $`i`$ and all the pairs. The result, plotted in Figure 7, gives $`Ps0.87`$. The requirement that both $`r_1`$ and $`r_2`$ are in the inertial range, strongly limits the statistics on turning points and the numerical result is affected by rather large uncertainty. Nevertheless, it is remarkable that the power law tail in the conditional probability density $`\mathrm{\Psi }(r_2|r_1)`$ is well observed in our numerical simulations. This justifies, a posteriori, the use of models based on $`\mathrm{\Psi }(r_2|r_1)`$ for describing relative dispersion . The numerical value of the effective persistence parameter $`Ps0.87`$ is not so small, and can explain the observed deviations from Richardson pdf (which are however less pronounced in a 2d flow than in a theoretical one-dimensional model ): The transport in a 2d turbulent flow is neither purely diffusive nor ballistic.

In order to be more confident on the numerical value of $`Ps`$ obtained through the turning-points statistics, let us show that it agrees with a simple estimates based on the values of the Kolmogorov’s and the Richardson’s constants. According to the Kolmogorov’s scaling, the mean squared relative velocity of the pair is given by

$$\delta v^2(r)=C_2\epsilon ^{2/3}r^{2/3},$$
(15)
with $`C_212.9`$ . If the particles separate ballistically with the rms velocity

$$\delta v(r)=C_2^{1/2}\epsilon ^{1/3}r^{1/3}$$
(16)
the distance between them should grow as

$$R_{\mathrm{max}}^2=\left(\frac{2}{3}\right)^3C_2^{3/2}\epsilon t^3.$$
(17)
On the other hand, due to the unsteadiness of the separation velocity, the distance between the particles grows slower, namely as $`R^2=g\epsilon t^3`$, so that the factor

$$\xi ^2=R^2/R_{\mathrm{max}}^2=\left(\frac{3}{2}\right)^3\frac{g}{C_2^{3/2}}$$
(18)
serves as a measure of this unsteadiness and $`\xi ^2`$ is connected with the value of the persistence parameter. In our case, $`\xi ^20.28`$. Within the stochastic model of Ref.  this corresponds to a value of $`Ps`$ between $`1.1`$ and $`1.2`$, in reasonable agreement with the direct measurement from the turning-point statistics, and again corroborates the stochastic approach.

We also note a possibility to “tune” the $`Ps`$ value by performing simulations in which the Lagrangian trajectories are integrated according to $`\dot{𝐱}=\lambda 𝐯(𝐱,t)`$. By changing the value of parameter $`\lambda `$ one effectively changes $`v_0`$ and thus $`Ps`$. In the extreme case $`\lambda 0`$ the trajectories resemble those in a time $`\delta `$-correlated velocity field. In the opposite limit, $`\lambda >>1`$ we have dispersion on a quenched field. Of course, it is only for the standard value $`\lambda =1`$ that Lagrangian trajectories move consistently with velocity field (i.e. for $`\nu =\alpha =\varphi =0`$ (14) conserves vorticity along the Lagrangian trajectories). For other values of $`\lambda `$ such simulations suffer the typical problem of advection in synthetic field (i.e. wrong reproduction of the sweeping effect, see  for a discussion). Simulations for several values of $`\lambda `$ show that existence of the power-law tails of $`\mathrm{\Psi }(r_2|r_1)`$ is a robust effect, as supposed by the model of Refs., and that $`Ps`$ grows with $`\lambda `$. As an example, in Figure 8 we plot the probability density $`\mathrm{\Psi }(r_2|r_1)`$ obtained from a simulation with $`\lambda =0.5`$. All the Eulerian parameters are the same of Figure 7. We again observe a clear power law tail but now with $`Ps0.58`$.

## V Conclusions

In this paper we have investigated the Lagrangian relative dispersion in direct numerical simulation of two-dimensional turbulence. The inverse energy cascade of two-dimensional turbulence displays Kolmogorov scaling without intermittency and it is thus the natural framework for investigating possible deviations from the classical Richardson picture. Moreover, the large enough values of the accessible Reynolds numbers make it possible to obtain the quantitative results not inferior to those of laboratory experiments.

The analysis of the numerical data was performed by using both standard statistics at fixed time and exit time statistics at fixed scale. The latter is shown to be more robust to finite Reynolds effects. An application of exit time statistics is developed for measuring the Richardson constant. The value found is somewhat larger than one obtained experimentally and in kinematic simulations, and thus nearer to the predictions of the closure theories.

The results of numerical simulation of statistics of the interparticle distances are in good agreement with the original Richardson’s description in terms of diffusion equation. The large deviations (with respect to Richardson theory) observed in the pdf of separation at fixed times are mostly related to crossover effects due to finite Reynolds numbers and disappear when looking at exit time statistics. The deviations found here are only of the quantitative nature. Thus, the Richardson’s equation gives a good basis for qualitative description of the dispersion in turbulent flows.

The deviations found in the fixed-scale statistics and also the strong sensitivity to the finite size effects can be associated with the long-range-correlated nature of the particles’ relative motion in turbulent flow. Paying attention to the turning points of the relative trajectories allows for estimating the effective persistence parameter of the motion which is found to be of the order of unity. Thus, the motion shows a strong ballistic component and is not purely diffusive. On the other hand, the correlations are not too strong to fully destroy the Richardson’s picture.

We note that the methodology of analysis proposed here based on the fixed-scale statistics and on the analysis of the relative trajectories can be also applied to the analysis of laboratory experiments. It would be extremely interesting to see, to what extent the laboratory flows reproduce the results reported here.

###### Acknowledgements.

We gratefully acknowledge the support of the DFG through SFB428, of the Fonds der Chemischen Industrie and of INFM (PRA-TURBO). We thank M. Cencini and J. Klafter for useful discussions. Numerical simulations were partially performed at CINECA within the INFM project “Fully developed two-dimensional turbulence”.

## A The PDF of Doubling Times

Let us discuss the probability $`p_D\left(t\right)`$ of the time when the a pair of particles initially at distance $`R/\rho `$ separates up to the distance $`R`$, and obtain its asymptotic decay of this probability for $`t`$ large.

Changing to a variable $`\xi =(k_0\epsilon ^{1/3})^{1/2}r^{1/3}`$ reduces (2) to a radial part of a spherically symmetric diffusion equation in $`3d`$ dimensions with constant diffusion coefficient. In 2d one has

$$\frac{p}{t}=\frac{1}{9\xi ^5}\frac{}{\xi }\xi ^5\frac{}{\xi }p$$
(A1)
with the initial condition $`p(\xi ,0)=\delta (\xi _{\mathrm{min}}\xi )`$ with $`\xi _{\mathrm{min}}=(k_0\epsilon ^{1/3})^{1/2}(R/\rho )^{1/3}`$ and with the boundary condition $`p(\xi _{\mathrm{max}},t)=0`$, with $`\xi _{\mathrm{max}}=(k_0\epsilon ^{1/3})^{1/2}R^{1/3}`$.

The solution of a boundary-value problem for (A1) can be obtained by means of eigenfunction decomposition. Assuming the variable separation we get the solution in the form $`p(\xi ,t)=_ie^{\lambda _i^2t}\psi _i(\xi )`$, where $`\psi _i(\xi )`$ is an eigenfunction of the equation

$$\frac{1}{9\xi ^5}\frac{}{\xi }\xi ^5\frac{}{\xi }\psi _i=\lambda _i^2\psi _i$$
(A2)
satisfying the boundary condition $`\psi (\xi _{\mathrm{max}})=0`$. The corresponding solution which is nonsingular in zero is $`\psi _i=\xi ^2J_2(3\lambda _i\xi )`$ ($`J_2`$ is the Bessel function ). The fact that $`\psi `$ vanishes at $`\xi _{\mathrm{max}}`$ gives $`3\lambda _i\xi _{max}=j_{2,i}`$, where $`j_{2,i}`$ is the $`i`$-the real zero of $`J_2(x)`$. For example, the smallest eigenvalue is $`\lambda _1^2=j_{2,1}^2/9\xi _{\mathrm{max}}^22.93k_0\epsilon ^{1/3}R^{2/3}`$. Thus, the long-time asymptotic of the doubling-time distribution is $`\mathrm{exp}(2.93k_0\epsilon ^{1/3}R^{2/3}t)`$.

## B Average doubling time

The mean doubling time can be obtained from a stationary solution of the Richardson diffusion equation. Imagine that one particle per unit time is introduced at $`r=R/\rho `$ and there are respectively a reflecting and absorbing boundaries at $`r=0`$ and $`r=R`$. The stationary solution of (2) in 2d with the appropriate boundary conditions and continuity are $`R/\rho `$ is

$$p(r)=\{\begin{array}{ccc}C\left[\rho ^{4/3}1\right]\hfill & \text{for}\hfill & 0<r<R/\rho \hfill \\ C\left[\left(\frac{r}{R}\right)^{4/3}1\right]\hfill & \text{for}\hfill & R/\rho <r<R\hfill \end{array}$$
(B1)
The number of particle in $`r<R`$ is

$$N=_{|𝐫|<R}p(r)𝑑𝐫=2\pi _0^Rrp(r)𝑑r$$
(B2)
where $`S_d`$ is a surface of a unit sphere in $`d`$ dimensions. By using (B1) one obtains

$$N=2\pi C(1\rho ^{2/3})R^d.$$
(B3)
The current at $`r=R`$, i.e. the number of particle exiting from the boundary $`R`$ per unit time, is given, as in (10), as

$$J=2\pi \epsilon ^{1/3}k_0R^{7/3}\frac{p(r)}{r}|_{r=R}=\frac{8\pi }{3}C\epsilon ^{1/3}k_0R^{4/3}$$
(B4)
The mean doubling time is the average time spent by a particle at $`r<R`$. It is given by the ratio $`N/J`$ and thus

$$T_\rho (R)=\frac{3}{4}\frac{\rho ^{2/3}1}{\epsilon ^{1/3}k_0\rho ^{2/3}}R^{2/3}$$
(B5)
which is Eq.(12).