text:
# Change of Variable for Multi-dimensional Integral

(4 March 2003)

## Abstract

The change of variable theorem is proved under the sole hypothesis of differentiability of the transformation. Specifically, it is shown under this hypothesis that the transformed integral equals the given one over every measurable subset on which the transformation is injective; that countably many of these subsets cover the domain of invertibility; and that its complement – the domain of non-invertibility – is measurable and so may be broken off and handled separately.

Isidore Fleischer

A.M.S. Subj. Classification 26B12, 26B15, 28A75.

The change of variable formula in multi-dimensional calculus, $`_{F(E)}\phi 𝑑y=_E\phi (F)|J_F|𝑑z`$, is usually established on the assumption that $`F`$ is a one-to-one selfmap of $`n`$-space continuously differentiable with non-vanishing Jacobian $`J_F`$ in an open domain containing $`E`$, $`F(E)`$ is a closed bounded domain and $`\phi `$ is real-valued integrable. The integrals make sense however – and thus the equality can be formulated – under more general conditions, the natural most general being just measurability of the domains and (Lebesgue) integrability of the integrands. Given the standard properties of the Lebesgue integral, it suffices to obtain the equality for the integrand $`\phi =1`$:

$$|F(E)|=_E\mathrm{\Delta }F^{}𝑑z,$$

a formula which will be shown valid on measurable subsets $`E`$ on which $`F`$ is injective. The absolute value of the Jacobian of $`F`$ has been rewritten $`\mathrm{\Delta }F^{}`$ (for reasons which will appear below).

If one wishes to cumulate these subsets (which will be shown to cover $`\mathrm{\Delta }F^{}>0`$ countably) one will need to replace the value of the measure on the left with the integral of the Banach indicatrix $`N(y)`$, which counts the number of elements of $`E`$ sent by $`F`$ on $`y`$. Also, the zero-set of $`\mathrm{\Delta }F^{}`$ will be shown to have measure zero, which permits establishing the equality over it seperately.

This form has been attained gradually – see the successive editions of Rudin – and appears finally in his most recent edition as well as in Smith with, however, some superfluous restrictions and – what is perhaps more irritating – an appeal in both cases to Brouwer theory. A purely analytic proof can be extracted by specialization from a proof of the “area formula” in Evans *et al*. but as given there contains some gaps and unnecessary detours. Thus the proof below, accessible to anyone with a modest command of Lebesgue theory, should offer some appeal.

This paper was submitted to the Monthly in September 1998, thus was contemporaneous with Lax’s which (although far less complete) appeared in the Monthly just when this one was rejected; a couple of years later, Lax published an extension (II in the References) which still does not tell the full story.

The sequel is organized as follows: The next five paragraphs motivate and introduce the “scale factor”, which is the ratio by which a linear selfmap augments (or reduces) volume, and verifies its continuity for the operator norm. The following three paragraphs recall the definition of the derivative for vector-valued functions and show its Borel measurability as a function from the domain of differentibility to the operator normed linear selfmaps. The remaining three paragraphs then justify the integrability of the integrand, the break-off of its zero set, and the desired equality.

For real-valued functions of a single real variable, the derivative equals simultaneously two a priori distinct characteristics of a function $`F`$ in the neighborhood of an argument $`x`$: it is the coefficient $`F^{}(x)`$ of the closest linear approximation to $`F`$ near $`x`$ and it is the slope of the tangent to the graph at $`x.`$ These entities need to be distinguished in higher dimension; let’s start with “slope.”

For a linear function on the line, the absolute value of its slope is the ratio of the length of the image of an interval by its length: it is the same for all

intervals and so may as well be defined as the length of the image of the unit semi-open interval $`[0,1).`$

For a linear selfmap $`L`$ of $`n`$-space, define analogously the “slope” or “scale factor” $`\mathrm{\Delta }L`$ as the ratio of the volume of the image of a box (i.e. a product of intervals) by its volume. This will be shown in a moment to be independent of the box and so may as well be defined as the volume of the image of the unit semi-open box (s.o.b.), i.e. the set of $`x`$’s whose co-ordinates satisfy $`0x_i<1.`$ (In the plane the image of the unit square is the parallelogram spanned by the images of the co-ordinate vectors, whose cross product gives its area; to recognize $`\mathrm{\Delta }L`$ in general as the absolute value of the determinant of the image vectors, one should arrange to obtain these as non–negative multiples of an orthonormal basis. However, this identification is not used in what follows.)

Divide the unit s.o.b. into disjoint equal size s.o.b.’s of side $`\frac{1}{k}:`$ their volume is $`\frac{1}{k^n}`$ and there are $`k^n`$ of them, sent into each other by translation. If $`L`$ is invertible it sends them into disjoint parallelepipeds, also sent into each other by translation, hence of equal volume, whose union has volume $`\mathrm{\Delta }L`$ – thus each image has volume $`\left(\frac{1}{k^n}\right)\mathrm{\Delta }L.`$ (If $`L`$ is not invertible it maps into a lower dimensional space and all images have $`n`$-volume zero.) Any open $`U`$ is a disjoint union of countably many s.o.b.’s, so its Lebesgue $`n`$-volume is also increased $`\mathrm{\Delta }L`$ times by $`L`$ and thus finally $`|L(E)|=(\mathrm{\Delta }L)|E|`$ for every set $`E`$.

$`\mathrm{\Delta }`$ is a continuous function on the $`L`$, topologized by uniform convergence on bounded sets – i.e. on any bounded with non-void interior. The image (by every $`L`$) of the closed unit box is compact, hence its epsilon neighborhoods are a base around the image: so a uniformly convergent sequence has image measures with lim sup $``$ the image measure of the limit $`L`$. To obtain the dual inequality, apply the same reasoning to the complement of the open unit box in some larger closed box.

A vector-valued function of a real variable can be differentiated in the classical way by taking the vector limit of the difference quotient, but this is already impossible even for a scalar-valued function of a vector variable. What one can do is to fix a vector $`u`$ and take the limit of $`\frac{F(z+hu)F(z)}{h}`$ as $`h0.`$ This limit would then exist for, and be positive homogeneous in, every positive multiple $`pu`$ in place of $`u`$. It thus suffices to postulate these limits for unit vectors $`u`$; they are called the directional derivatives (of $`F`$ at $`z`$) in the directions $`u`$ (the unit vectors are construed to point in the direction from the origin to their endpoint on the unit sphere; for a function of a single real variable this is the one-sided derivative and for co-ordinate $`u`$, the one-sided partials).

Assuming the directional derivatives exist in all directions $`u`$ at $`z,`$ one could define the “derivative” $`F^{}(z)`$ of $`F`$ at $`z`$ as the function which takes each $`u`$ to the directional derivative $`F^{}(z)u.`$ One will however require the approach of $`\frac{F(z+hu)F(z)}{h}`$ to $`F^{}(z)u`$ to be uniform in the unit vectors $`u`$; and if $`F^{}(z)`$ is bounded this will result in $`\frac{F(x)F(z)}{|xz|}`$ bounded for $`x`$ sufficiently close to $`z.`$ This last may be shown to entail that $`F^{}(z)pu:=pF^{}(z)u,p0`$, acts linearly on vectors $`pu`$ for a.e. $`z`$: accordingly, we now strengthen the definition by making $`F^{}(z)`$ a linear operator whose value at every $`v`$ is the limit as $`h0`$, uniform for bounded vector arguments $`v`$, of $`\frac{F(z+hv)F(z)}{h}`$.

The derivative is thus a function from the domain, assumed measurable (follows from $`F`$ continuous), of differentiability of $`F`$ to the space of linear selfmaps of $`n`$-space topologized by uniform convergence on the unit ball111This topology is defined by the norm $`L`$, the radius of the smallest ball centered at $`0`$ which contains the image of the centered ball of radius one; by linearity this $`L`$ is also the Lipschitz constant. – a base at $`0`$ is obtained by taking for every centered open ball the subset of those which send the unit ball into it. This function will next be shown to be (Borel) measurable: i.e., the inverse image of every open is measurable. Let $`r_mr,`$ the radius of centered open ball $`B`$ in $`R^n,`$ which is an increasing union of closed balls $`B_m`$ of radius $`r_m`$ with center zero. Then $`F^{}(z)uB`$ entails (by compactness of the $`u`$’s) that all but finitely many (in $`k`$) of the continuous $`k\left[F\left(z+\frac{u}{k}\right)F(z)\right]`$ belong to some $`B_m`$ (recall) uniformly in $`u;`$ the converse already holds if the belonging is uniform in some countable dense subset of the unit ball. Hence $`\{z:F^{}(z)uB\}`$ is the countable union (over $`m`$) of the sets where a countable intersection (over $`u`$) of lim sup (over $`k`$) of the modulus of continuous functions of $`z`$ and $`u`$ is $`r_m.`$

Since $`\mathrm{\Delta }L`$ has been shown continuous for linear operators $`L`$, it follows that the composite, $`\mathrm{\Delta }F^{}(z)`$, is a non-negative measurable function on a measurable subset $`E`$ of $`n`$-space, hence has a well-defined Lebesgue integral $`_E\mathrm{\Delta }F^{}(z)𝑑z.`$

(It also follows that the set where $`\mathrm{\Delta }F^{}(z)=0`$ is measurable; since the equality for this set can be established by a simple calculation, we shall restrict to $`\mathrm{\Delta }F^{}(z)>0`$ in the sequel.)

If $`F^{}:=F^{}(z)`$ is invertible at $`z`$ then for $`|v|<`$ some $`\delta ,|F(z+v)F(z)F^{}v|\left[\epsilon /F^1\right]|v|\epsilon |F^{}v|`$; thus $`(1\epsilon )|F^{}xF^{}z||FxFz|(1+\epsilon )|F^{}xF^{}z|`$ for $`|xz|<\delta `$. Observe that for fixed $`\delta `$ the set of $`z`$’s satisfying (the extremes of) the first inequality is measurable: both terms are continuous in $`v`$ and measurable in $`z`$ so imposing it at a countable dense subset of the ball of radius $`\delta `$ reduces the implication to a countable set of inequalities between measurable functions. Since $`E`$ is a countable union of these sets for a sequence of $`\delta `$’s, it suffices to establish the theorem on each of them. Then one can replace the fixed argument $`z`$ of $`F^{}`$ with any $`y`$ in this subset and still obtain (for this $`\delta `$) the extremes of the first, hence also the second, inequality with $`y`$ in place of $`z`$. Restricting further to measurable subsets of diameter $`<`$ $`\delta `$ on which the values of $`|F^{}|`$ lie between $`(1\epsilon )`$ and $`(1+\epsilon )`$ of its value at $`z`$ for every vector in the ball of radius $`\delta `$, pass from the inequality with $`y`$ in place of $`z`$ back to $`F^{}`$ at z, to attain $`(1\epsilon )^2|F^{}xF^{}y||FxFy|(1+\epsilon )^2|F^{}xF^{}y|`$ even for $`x,yz`$. (Since there is a countable cover by such sets, it suffices to restrict attention to each of them.) By the left inequality and invertibility of $`F^{}`$ at $`z`$, $`F`$ is one-to-one on this subset. From $`|FxFy||GxGy|`$ for any one-to-one $`F`$ and all $`x,y`$ in some measurable set follows: $`FG^1`$ is non-expansive – Lip(1) – hence does not increase measure; so $`|FE||FG^1GE||GE|`$ for each of its subsets $`E`$: thus $`(1\epsilon )^{2n}\mathrm{\Delta }F^{}|E|=(1\epsilon )^{2n}|F^{}(E)||F(E)|(1+\epsilon )^{2n}|F^{}(E)|=(1+\epsilon )^{2n}\mathrm{\Delta }F^{}|E|`$. Since $`F`$ is one-to-one, it preserves disjointness of subsets and since differentiable, measurability of subsets; hence $`|F()|`$ is additive on measurable decompositions. By passing to the limit in the definition of the integral one is led at last to

$$(1\epsilon )^{2n}_E\mathrm{\Delta }F^{}𝑑z|F(E)|(1+\epsilon )^{2n}_E\mathrm{\Delta }F^{}𝑑z$$

with arbitrary $`\epsilon `$.

References

T. M. Apostol, *Mathematical Analysis*, Addison-Wesley (1957).

P. Billingsley, Probability and Measure, John Wiley & Sons (1979).

L. C. Evans & R. F. Gariepy, *Measure Theory and Fine Properties of Functions*, CRC Press (1992).

P. D. Lax, *Change of variables in multiple integrals*, Amer. Math. Monthly, 106 (1999), 497–501; II, Ibid 108 (2001) 115–119.

I. P. Natanson, *Theory of Functions of a Real Variable*, F. Ungar (1955).

W. F. Pfeffer, The Riemann Approach to Integration, Cambridge U. Press (1993).

W. Rudin, Real and Complex Analysis, 3rd edition, McGraw-Hill (1987).

K. T. Smith, *Primer of Modern Analysis*, Springer (1983).

T. Traynor, *Change of Variable for Hausdorff Measure (from the beginning)*, Rendiconti dell’Istit di Mat., Trieste 26 (1994), suppl.,328–347/

Centre de recherches mathématiques

Université de Montréal

C.P. 6128, succursale centre-ville

Montréal, Qc H3C 3J7

email:fleischi@crm.umontreal.ca


format: text

text:
# Long-Term Timing of 374 Pulsars

## 1. Introduction

More than 500 pulsars are being regularly observed using the 76-m Lovell radio telescope at Jodrell Bank Observatory at frequencies between 408 and 1630 MHz. This archive contains more than 5600 years of pulsar rotational history which we supplement, for 18 pulsars, with early observations from the Jet Propulsion Laboratory (Downs & Reichley 1983) to provide individual data spans of up to 34 years. In a series of three papers we plan to carry out a full analysis of the data for more than 350 pulsars that have been observed for longer than six years. The first paper provides accurate timing solutions including proper motion measurements. In the second paper we will use these proper motion values to improve our understanding of pulsar velocities. In the third paper we will discuss the remnant structures in the timing residuals after fitting a timing solution for rotational frequency and its first derivative. Some qualitative results, which will be discussed in more detail in the third paper, are highlighted here.

## 2. Structure in the timing residuals

Many different structures are observed in the timing residuals, including cubic features, cusps, flat residuals, pseudo-sinusoids and more irregular features. Cusps (corresponding to a local maxima) indicate the pulsar undergoing a glitch. Pseudo-sinusoidal variations in the timing residuals may be caused by the pulsar precessing or from unmodelled binary companions.

The amplitudes of the observed cubic structures are many orders of magnitude greater than those predicted from magnetic dipole radiation; measured ‘braking indices’ range between $`\pm 3\times 10^8`$.

The structures observed in the majority of the timing residuals are not time resolved. For example, about 30% are dominated by a large scale cubic feature and a further 40% currently have no clearly defined features. However, we have a sample of 29 pulsars for which it is possible to study individual features in detail. We find that the distribution of the residuals around the mean value has no general asymmetry (with 15 having negative skew and 14 positive). We do find that 70% have, on average, sharper (larger second derivative) local maxima than local minima. No explanation is given for this effect at the present time.

The time corresponding to the first local minima in the auto-correlation function of the timing residuals provides an estimate of the time–scale of the structures present. However, we see no correlation with this parameter and the observed or derived pulsar parameters. We do find strong correlations between the amplitude of the timing noise and the first and second rotational frequency derivatives (Figure 1).

Our new sample of timing residuals that span up to $``$30 years will allow the first large scale analysis of timing noise to be carried out. Already it is clear that many of the structures observed in the timing residuals can not have been caused by a simple ‘noise’ process.

## 3. Acknowledgements

Many people have been involved in timing pulsars at Jodrell Bank. We thank, in particular, Christine Jordan for her work in developing much of the software used in the analysis of the observations.

## References

Downs, G. S. & Reichley, P. E., ApJS, 53, 169. 1983


text:
# The potential of the ground state of NaRb

## I Introduction

The heteronuclear alkali dimers attract interest of both experimental and theoretical researchers involved in collision dynamics, photoassociative spectroscopy, laser cooling and trapping of alkali atoms Shaffer:99 ; Ferrari:02 ; Hadzibabic:02 ; Weiss:03 . Special interest is put on the study of the ground states and especially near the atomic asymptote, since the precise knowledge of the long-range interactions between two different types of alkali atoms is necessary for understanding and realization of such cold collision processes as sympathetic cooling, formation of two species BEC (TBEC) and ultracold heteronuclear molecules.

Although the NaRb molecule is one of the promising candidates for TBEC Weiss:03 , the experimental information on the ground singlet state is still limited.

Accurate term values of the X$`{}_{}{}^{1}\mathrm{\Sigma }_{}^{+}`$ state were determined by the Doppler-free polarization spectroscopy Wang2:91 . They cover a range of low vibrational quantum numbers $`v^{\prime \prime }[0,6]`$ with limited sets of rotational quantum numbers (e.g. $`J^{\prime \prime }[1,80]`$ for $`v^{\prime \prime }=0`$ and $`J^{\prime \prime }[18,36]`$ for $`v^{\prime \prime }=6`$). Later Kasahara et al. Kasahara:96  extended the range of $`v^{\prime \prime }`$ up to $`v^{\prime \prime }=30`$, but only for $`J^{\prime \prime }=10,12`$. On the other hand the recent experiment by the Riga group Docenko:02  collected about 300 transition frequencies to energy levels with a much wider spectrum of vibrational and rotational quantum numbers ($`v^{\prime \prime }[15,76]`$, $`J^{\prime \prime }[12,64]`$) but with relatively moderate accuracy ( 0.1 cm-1 ). A hybrid potential for the NaRb ground state based on extrapolation between regions determined from the early experimental data Wang2:91 ; Kasahara:96  and the theoretical dispersion coefficients from Ref. Marinescu:99  was constructed by Zemke et al. Zemke:NaRb  giving an improved estimate on the dissociation energy of the ground state - $`D_\mathrm{e}=5030.75`$ cm-1.

On the theoretical side two papers have recently appeared: by Korek et al. Korek:00  and by Zaitsevskii et al. Zaitsevskii:01 . In both of them theoretical potential energy curves for the ground and several excited states are given. Results from Ref. Korek:00  on selected low electronic states are shown also in Figure 1, which will help in following the electronic assignment of the new observations, reported below.

The goal of the present investigation is to study extensively the ground electronic state X$`{}_{}{}^{1}\mathrm{\Sigma }_{}^{+}`$ of the NaRb molecule. We chose the Fourier transform spectroscopy which is able to provide abundant and precise spectroscopic data, needed for construction of an accurate potential energy curve. In addition, this study is considered as the necessary background for further investigations of this molecule at long internuclear distances required for modeling cold collisions at 1mK and below.

## II Experiment

The NaRb molecules were produced in a single section heat-pipe oven similar to that described in Ref. Allard:02  by heating 5 g of Na (purity 99.95 %) and 10 g of Rb (purity 99.75 %, natural isotopic composition) from Alfa Aesar. The oven was operated at temperatures between 560 K and 600 K and typically with 2 mbar of Ar as buffer gas. At these conditions apart from the atomic vapors all three types of molecules are formed - Na2, Rb2 and NaRb. This mixture was illuminated by three different laser sources: a single mode, frequency doubled Nd:YAG laser, Ar+ ion laser and Rhodamine 6G dye laser. Their wavelengths fall in a spectral window with weak absorption of Rb2. Setting the working temperatures as low as possible we reduced the density of the Na2 molecules reaching a relatively intense NaRb spectra free of Na2 emission. The heat pipe oven was in operation more than 60 hours without refilling and at the end of the experiment was still in good working condition.

The fluorescence from the oven was collected in the direction opposite to the propagation of the laser beam and then focused into the input aperture of Bruker IFS 120 HR Fourier transform spectrometer. The signal was detected with a Hamamatsu R928 photomultiplier tube or a silicon diode. The scanning path of the spectrometer was set to reach a typical resolution of 0.0115 cm-1 $`÷`$ 0.03 cm-1 and the number of scans for averaging varied between 10 and 20. In order to avoid the illumination of the detector by the He-Ne laser, used for calibration of the spectrometer, a notch filter with 8 nm FWHM was introduced. For better signal-to-noise ratio some spectra were recorded by limiting the desired spectral window with colored glass or interference filters.

The single mode, frequency doubled Nd:YAG laser with a typical output power of 70 mW at 532.1 nm excited the C$`{}_{}{}^{1}\mathrm{\Sigma }_{}^{+}`$ X$`{}_{}{}^{1}\mathrm{\Sigma }_{}^{+}`$ and D$`{}_{}{}^{1}\mathrm{\Pi }`$ X$`{}_{}{}^{1}\mathrm{\Sigma }_{}^{+}`$ systems of NaRb, compare Fig. 1. The frequency was varied between 18787.25 cm-1 and 18788.44 cm-1 and spectra were recorded at frequencies which excite strong fluorescence.

The Ar+ ion laser was operated both in single (typical power 100$`÷`$500 mW) and multi mode (typical power 0.5$`÷`$3 W) regimes. The 514.5 nm, 501.7 nm, 496.5 nm, 488.0 nm and 476.5 nm Ar+ ion laser lines induced fluorescence mainly due to the D$`{}_{}{}^{1}\mathrm{\Pi }`$ X$`{}_{}{}^{1}\mathrm{\Sigma }_{}^{+}`$ system of NaRb. The Na2 B$`{}_{}{}^{1}\mathrm{\Pi }_{\mathrm{u}}^{}`$ X$`{}_{}{}^{1}\mathrm{\Sigma }_{\mathrm{g}}^{+}`$ band was also observed, especially exciting with the bluer Ar+ laser lines, but its intensity was reduced by decreasing the working temperature down to 560 K. Along with the D $``$ X system, the 514.5 nm line excites also transitions in the C$`{}_{}{}^{1}\mathrm{\Sigma }_{}^{+}`$ X$`{}_{}{}^{1}\mathrm{\Sigma }_{}^{+}`$ system of NaRb Docenko:02 . The special form of potential energy curve of the C$`{}_{}{}^{1}\mathrm{\Sigma }_{}^{+}`$ state Korek:00  results in favorable Franck-Condon factors for transitions from some C state levels both to low and high lying energy levels of the ground state. Such C state levels are excited by the 514.5 nm line and the green florescence of the D $``$ X system is accompanied by a red fluorescence due to C $``$ X transitions, which reach $`v^{\prime \prime }=76`$ of the ground state (Fig. 2).

In order to bring the frequency of the single mode Ar+ laser exactly in resonance with desired transitions we did the following. As a first step we recorded the LIF for a given frequency simultaneously with the Fourier spectrometer and a GCA/McPherson Instruments scanning monochromator (1m focal length). The high resolution Fourier spectrum helped us to assign unambiguously the record of the monochromator. Then we set the monochromator on a strong and pronounced line of the progression of interest and we tuned the Ar+ laser frequency until the maximum signal on the monochromator was reached. The frequency stability of the laser with a temperature stabilized intracavity etalon was sufficient in order to perform 10 to 20 scans with the Fourier spectrometer.

The Rhodamine 6G laser excites the B$`{}_{}{}^{1}\mathrm{\Pi }`$ X$`{}_{}{}^{1}\mathrm{\Sigma }_{}^{+}`$ transition and also weak transitions in the C$`{}_{}{}^{1}\mathrm{\Sigma }_{}^{+}`$ X$`{}_{}{}^{1}\mathrm{\Sigma }_{}^{+}`$ system and was used for two reasons. The first was to enrich the information on the ground state levels with intermediate vibrational quantum numbers because of the gap between the levels from the D$`{}_{}{}^{1}\mathrm{\Pi }`$ X$`{}_{}{}^{1}\mathrm{\Sigma }_{}^{+}`$ system (mainly with low $`v^{\prime \prime }`$) and those from the C$`{}_{}{}^{1}\mathrm{\Sigma }_{}^{+}`$ X$`{}_{}{}^{1}\mathrm{\Sigma }_{}^{+}`$ system (mainly with high $`v^{\prime \prime }`$). Secondly, we wanted to excite levels of the B state with significant triplet admixture due to perturbations from the neighboring b$`{}_{}{}^{3}\mathrm{\Pi }`$ and c$`{}_{}{}^{3}\mathrm{\Sigma }_{}^{+}`$ states Wang2:91 ; Wang:92 . Indeed, scanning the frequency of the dye laser from 16729 cm-1 to 16965 cm-1 we encountered a large number of excitations where a second fluorescence band around 12000 cm-1 appeared along with the B $``$ X system (see Fig. 3). The analysis of this band (a report is in preparation) confirmed that we have observed transitions to the ground triplet state a$`{}_{}{}^{3}\mathrm{\Sigma }_{}^{+}`$.

All transitions excited in Na85Rb and Na87Rb by the single mode Nd:YAG and Rhodamine 6G lasers and by the Ar+ laser in single- and multi mode, assigned and used in this analysis are listed in Table 1 of the supplementary material (EPAPS No. ….).

## III Analysis

The assignment of the recorded spectra was simplified by the published potential energy curve in Ref. Docenko:02 . After we identified the strongest progressions, which were already observed in previous studies Takahashi:81 ; Nikolayeva:00 ; Zaitsevskii:01 ; Docenko:02 , a new potential was fitted. This preliminary potential was further improved filling additionally collected experimental data; this sequential procedure allows continuous checking of the assignment, especially for high $`J^{\prime \prime }`$ and in cases of large gaps in $`v^{\prime \prime }`$. The total data set (Fig. 4) consists of more than 6150 transitions in Na85Rb and 2650 in Na87Rb given in Table 2 of the supplementary material (EPAPS No…..).

For fitting the experimental data we used the approach adopted in Ref. Allard:02 . In order to exclude any fitting parameters concerning the excited states, we fitted the ground state potential directly to the observed differences between ground state levels. In our case of 8800 transitions we selected about 43300 differences for representing the different cases of vibrational intervals (for details see also Allard:02 ). The uncertainties of the differences were estimated taking into account the instrumental accuracy and resolution, and the signal-to-noise ratio for both transitions forming the difference. The uncertainty of strong lines was set to 0.002 $`÷`$ 0.003 cm-1 increasing to 0.007 cm-1 for weaker ones with SNR $`2`$. Differences which were measured several times were averaged with weighting factors determined from their uncertainties.

As already discussed in Allard:02 , for our experimental set up a fluorescence progression excited by a single mode laser will be shifted in frequency if the laser is not tuned to the center of the Doppler profile. This problem is partially overcome by fitting differences instead of transition frequencies. Nevertheless there can be still a residual shift of the differences due to the frequency dependence of the Doppler shift. For example if the laser frequency around 20000 cm-1 is detuned from the line center by one Doppler width (about 0.03 cm-1 for NaRb at 600 K) the difference between two transitions $`\nu _1`$ and $`\nu _2`$ separated by $`\nu _1\nu _2`$= 4000 cm-1 will be shifted from the “true” one by 0.006 cm-1. This shift is comparable with the typically expected uncertainty ($`0.004`$ cm-1) but it cannot be determined in our experiment. Therefore, initially we did not correct the uncertainties of the large differences (a possible solution could be to increase the uncertainty of the difference $`\nu _1\nu _2`$ proportionally to the ratio $`(\nu _1\nu _2)/\nu _1`$). Large differences however were checked carefully during the fitting procedure and finally no need for such correction was found.

## IV Construction of PEC

The ground singlet states of both isotopomers of NaRb are described in the adiabatic approximation with a single potential energy curve. The potential was constructed as a set of points $`\{R_i,U(R_i)\}`$ (see Ref. ipaasen ) connected with cubic spline function. The pointwise representation of the potential, however, is not convenient for long internuclear distances. On the other hand, in order to ensure proper boundary conditions for solving the Schrödinger equation for high vibrational energy levels, distances up to $`R20`$ Å are needed. Therefore for $`RR_{\mathrm{out}}`$ we adopted the usual dispersion form:

$$U(R)=D_eC_6R^6C_8R^8C_{10}R^{10}\text{ ,}$$
(1)
with coefficients $`C_6,C_8`$ taken from Ref. Derevianko:01a ; Porsev:03 . The connecting point $`R_{\mathrm{out}}`$ and the parameters $`D_e`$ and $`C_{10}`$ were varied in order to ensure a smooth connection with the pointwise potential as follows. Initially, the PEC was constructed in a pointwise form up to 16 Å. After each fitting iteration $`D_e`$ and $`C_{10}`$ were adjusted to fit the shape of the pointwise potential between 13.0 Å and 16.0 Å to the analytic formula (1). The crossing point of the pointwise and the long range curves was taken as $`R_{\mathrm{out}}`$. Since the experimental data are almost insensitive to details of the potential shape in this region, we found this procedure sufficient and did not fit the long range parameters directly to the experimental data. Once the final form of the PEC was achieved, we found that it is possible to extend the validity of the long range expression down to approximately 11.8 Å. At the same time the quality of the derived potential did not change if the analytic formula (1) was fitted to the shape of the pointwise potential between 11.8 Å and 13.24 Å. A discussion on the advantage of using this matching procedure can be found in Ref. MC:03 .

The final set of potential parameters consists of 51 points. In order to calculate the value of the potential for $`R<R_{\mathrm{out}}`$ a natural cubic spline Numer  through *all* 51 points and for $`RR_{\mathrm{out}}`$ the long range expression (1) should be used with parameters listed in Table 1. This final PEC gives a standard deviation of $`\sigma =0.0031`$ cm-1 and a normalized standard deviation of $`\overline{\sigma }=0.70`$ showing the internal consistency of the data set and the quality of the derived potential.

In order to find the region where the PEC is unambiguously characterized by the experimental data we analyzed the uncertainty of the fitting parameters as described in detail in Ref. Li2F:00 ; Allard:02 . By minimizing the merit function $`\chi ^2`$ using the Singular Value Decomposition (SVD) method Wilkinson:SVD ; Numer  it is possible to determine the parameters on which $`\chi ^2`$ and consequently the quality of the fit depend only weakly. The fitting parameters in our case are the values of the potential itself in a preset grid of points. In Fig. 5 the uncertainties of the fitted potential for intermediate internuclear distances are presented for three different values of the singularity parameter $`\zeta `$ (see Li2F:00 ). Although the outer classical turning point of the last observed energy level ($`v^{\prime \prime }=76`$, $`J^{\prime \prime }=27`$) is 12.4 Å we see that the shape of the potential energy curve is unambiguously fixed (within the experimental uncertainty) approximately only between 3.0 Å and 9.8 Å since the uncertainties there almost do not depend on $`\zeta `$. Therefore, although the fitted potential energy curve describes all the experimental data up to $`v^{\prime \prime }=76`$, there is no rigorous way to estimate the uncertainty of its shape beyond 9.8 Å (see also Allard:02 ). The result of Fig. 5 might indicate that the correlation between the left (R $`<`$ 3Å) and right (R $`>`$ 9.8 Å) branch of the potential becomes significant. In this respect the reported potential in Table 1 for these regions is only one possibility which describes our observations within the error limits. From Table 1 we derive some parameters which will be useful when applying the potential for other calculations and spectroscopic studies. Their values are not rounded in order to keep consistency with the data from Table 1.

* equilibrium distance of the potential,
$`R_e=3.643415`$ Å,
* position of the level with $`v^{\prime \prime }=0`$ and $`J^{\prime \prime }=0`$,
$`E_{00}=53.3117`$ cm-1 with respect to potential minimum,
* the dissociation energy of the ground state, calculated with respect to $`E_{00}`$,
$`D_0^\mathrm{X}`$=4977.5363 cm-1.

In Fig. 6 the potential energy curves published in Ref. Zemke:NaRb ; Docenko:02  are compared with the pointwise potential (PW) of this study. The reason for the difference between the PW potential and the Modified Lennard-Jones (MLJ) potential Docenko:02  around $`R_e`$ is probably that for the construction of the later potential the rotational quantum numbers for levels with $`v^{\prime \prime }<30`$ were limited only to the narrow interval from Kasahara:96  ($`J^{\prime \prime }=10,12`$). With such data set the rotational constant $`B_v`$ can be determined only with large uncertainty which corresponds to shifts of the whole potential along the internuclear axes which is clearly indicated in Fig. 6. As already discussed in Section I the hybrid potential from Ref. Zemke:NaRb  was constructed using the whole range of available rotational quantum numbers from the literature (Wang2:91 ; Kasahara:96 ) and therefore the agreement with the present study around $`R_e`$ is much better. The deviations between the PW and the hybrid potential reaching 20 cm-1 for intermediate internuclear distances are caused by the extrapolation between the theoretical long range part of the potential and the experimental short range part. In the same region the experimentally determined shape of the MLJ potential is in a much better agreement with the present study.

Since the primary data from Wang2:91  were not available to us, we were able to check our potential only against the data published in Kasahara:96 . Forming differences between transition frequencies exactly as for our LIF data we found an agreement with the differences predicted by the present PEC from Table 1 with a standard deviation of 0.0026 cm-1, well within the expected experimental error of 100 MHz. Only two transitions ((12,11)-(17,10) and (12,11)-(19,10)) were excluded from this comparison since differences with them showed systematically large deviations reaching 0.01 cm-1.

In Table 2 a set of Dunham coefficients for 23Na85Rb is listed which describes the whole set of the measured differences for both isotopomers with a deviation of 0.0030 cm-1 and a normalized standard deviation of 0.84. For calculating the eigenvalues for 23Na87Rb the usual scaling rules were used. The distribution of the indexes of the Dunham coefficients indicates that for description of such a large set of vibrational and rotational quantum numbers the parameters of the Dunahm expansion lose their original meaning of “spectroscopic constants”. The extremely small values of some coefficients ($`10^{40}`$ cm-1 !) indicate that they should be used with caution due to roundoff errors. At last, the extrapolation properties of such set of Dunham coefficients are doubtful. Therefore, we prefer to give the full description of the experimental data by a potential energy curve, which is obviously a much better physical model, and to present Dunham coefficients only for convenience since there are still applications where their use is easier and faster.

## V Results and Discussion

In this paper a spectroscopic study of the NaRb ground state is presented. A variety of laser sources was used to excite several band systems in both isotopomers of the NaRb molecule and the induced fluorescence was analyzed with a high resolution Fourier transform spectrometer. More than 8800 transitions were identified resulting in a set of more than 4000 ground state energy levels with wide range of vibrational and rotational quantum numbers. An accurate potential energy curve was fitted to the experimental data reproducing differences between ground state energy levels with $`\sigma =0.0031`$ cm-1 and $`\overline{\sigma }=0.70`$.

In order to check the consistency of the derived potential energy curve we performed additional fits in several different ways. Instead of fitting differences, we applied the commonly used approach of fitting directly transition frequencies, getting the term energies of the excited state levels as free parameters. In this case two representations of the potential were tried: a pointwise and an analytic form described in detail in Ref. Samuelis:00 . Both representations gave the same quality of the fit compared to the pointwise potential determined from differences and the potentials were practically identical in the studied range. Thus no additional parameters sets are given here in order to avoid the confusion of the reader what potential form should be preferred.

The classical turning point of the last observed energy level ($`v^{\prime \prime }=76`$, $`J^{\prime \prime }=27`$) is around 12.4 Å. Although this point lies beyond the Le Roy radius for the NaRb 32S+52S asymptote ($`R_{\mathrm{LeRoy}}=11.2`$ Å), the analysis of the PEC uncertainties indicates that with the present body of experimental data it is unsafe to determine $`C_6`$ and the other dispersion coefficients. For this reason the experimental PEC was smoothly connected to a long range potential formed by dispersion terms taken from the literature Derevianko:01a ; Porsev:03 . For better description of the near asymptotic region of the PEC transition frequencies to weakly bound energy levels are needed (see for example MC:03 ). Some improvement might be expected already from the combined analysis of the X$`{}_{}{}^{1}\mathrm{\Sigma }_{}^{+}`$ and a$`{}_{}{}^{3}\mathrm{\Sigma }_{}^{+}`$ states at long internuclear distances which is in progress.

Weiss et al. Weiss:03  applied the MLJ potential from Docenko:02  and extended it for $`R>11`$ Å by the long range parameters from Marinescu:99 ; Derevianko:01a  and the exchange interaction taken from Smirnov:65  to determine the singlet and triplet scattering lengths for 32S Na + 52S Rb. Using the newly derived potential for the inner part we constructed the full potential with the same long range parameters as in Weiss:03  and derived a singlet scattering length of +305a0 ($`a_00.52918`$ Å, the Bohr radius) for 23Na85Rb and +103a0 for 23Na87Rb. These results do not agree with the ones from Ref. Weiss:03 : +167a0 (+50a0/-30a0) for 23Na85Rb and 55a0 (+3a0/-3a0) for 23Na87Rb. We used exactly the same approach for the long range as in Weiss:03  to see clearly the influence of the inner part of the potential. By the new values we will not claim that we have determined more realistic scattering lengths. Rather we want to point out that the existing data do not allow to specify the scattering length in an useful uncertainty interval. Note also that slight changes in the exchange energy still allow to alter the sign of the scattering length in 23Na85Rb. Thus energy levels closer to the asymptote ($`<4.5`$ cm-1) must be studied.

## VI Acknowledgments

This work is supported by DFG through SFB 407. The authors appreciate the assistance of M. Klug and Ch. Samuelis during the experiments. A.P. gratefully acknowledges the research stipend from the Alexander von Humboldt Foundation. O. D., M. T. and R. F. acknowledge the support by the University of Latvia, by the Latvian Science Council grant 01.0264, the Latvian Ministry of Education and Science (grant TOP 02-45) as well as by NATO SfP 978029 - Optical Field Mapping grant.


text:
# Stochastic resonance and noise delayed extinction in a model of two competing species

## 1 Introduction

Recently the influence of noise on population dynamics has been the subject of intense theoretical investigations . The problem of the stability of the biological complex systems in the presence of noise for example is one of the more discussed in research activity of these last years . The noise through its interaction with the nonlinearity of the systems can give rise to new counterintuitive phenomena like stochastic resonance , noise enhanced stability , noise delayed extinction, etc… The principal aim of this work is the investigation of a generalized Lotka-Volterra model with a random interaction parameter between two competing species and in the presence of a multiplicative noise . We find that the noise has a constructive role. In fact it is responsible for the generation of quasi-periodic temporal oscillations and the enhancement of the response of the system to a driving force producing stochastic resonance. Moreover we find noise delayed extinction, i. e. a nonmonotonic behavior of the average extinction time of one of the two species as a function of the noise intensity.

## 2 The model

Time evolution of two competing species is obtained within the formalism of the Lotka-Volterra equations  in the presence of a multiplicative noise

$`{\displaystyle \frac{dx}{dt}}=\mu _1x(\alpha _1x\beta _1(t)y)+x\xi _x(t)`$ (1)
$`{\displaystyle \frac{dy}{dt}}=\mu _2y(\alpha _2y\beta _2(t)x)+y\xi _y(t),`$ (2)
where $`\xi _x(t)`$ and $`\xi _y(t)`$ are statistically independent Gaussian white noises with zero mean and correlation function $`\xi _i(t)\xi _j(t^{})=\sigma \delta (tt^{})\delta _{ij}`$ ($`i,j=x,y`$).

The time series for the two populations are obtained setting $`\alpha _1=\alpha _2=\alpha `$, $`\beta _1(t)=\beta _2(t)=\beta (t)`$. It is known that for $`\beta <1`$ a coexistence regime takes place, that is both species survives, while for $`\beta >1`$ an exclusion regime is established, that is one of the two species vanishes after a certain time. Coexistence and exclusion of one of the two species correspond to stable states of the Lotka-Volterra’s deterministic model . Real ecosystems are immersed in a noisy nonstationary environment, so also the interaction parameter is affected by the noise and some other deterministic periodical driving such as the temperature. The change in the competition rate between exclusion and coexistence e. g. occurs randomly because of the coupling between the limiting resources and the noisy environment. A random variation of limiting resources produces a random competition between the species. The noise therefore together with the periodic force determines the crossing from a dynamical regime ($`\beta <1`$, coexistence) to the other one ($`\beta >1`$, exclusion). To describe this continuous and noisy behaviour of the interaction parameter $`\beta (t)`$ we consider a stochastic differential equation with a bistable potential and a periodical driving force

$$\frac{d\beta (t)}{dt}=\frac{dU(\beta )}{d\beta }+\gamma cos(\omega _0t)+\xi _\beta (t),$$
(3)
where $`U(\beta )`$ is a bistable potential (see Fig.1)

$$U(\beta )=h(\beta (1+\rho ))^4/\eta ^42h(\beta (1+\rho ))^2/\eta ^2,$$
(4)
and $`h`$ is the height of the potential barrier. The periodic term takes into account for the environment temperature variation. Here $`\gamma =10^1`$ and $`\omega _0/(2\pi )=10^3`$. In Eq.(3) $`\xi _\beta (t)`$ is a Gaussian white noise with the usual statistical properties: $`\xi _\beta (t)=0`$ and $`\xi _\beta (t)\xi _\beta (t^{})=\sigma \delta (tt^{})`$. Due to the shape of $`U(\beta )`$ it is reasonable to expect a coexistence regime for $`\beta (0)<1`$, when deterministic case ($`\xi _\beta (t)=0`$) is considered.

### 2.1 Stochastic resonance

First we investigate the effect of the noise on the time behavior of the species. Since the dynamics of the species strongly depends on the value of the interaction parameter, we initially analyze the time evolution of $`\beta (t)`$ for different levels of the additive noise $`\sigma _\beta `$.

Specifically for $`\sigma _\beta =0`$ we obtain a periodical behavior of $`\beta (t)`$ in the coexistence region (see Fig.2a). For low noise intensity ($`\sigma _\beta h`$) we get the same periodical behavior of deterministic case, slightly perturbed by the noise (Fig.2b). For higher noise intensity ($`\sigma _\beta h`$) the periodical behavior of the interaction parameter jumps between the two values $`\beta =0.94<1`$ and $`\beta =1.04>1`$, which characterize the coexistence and the exclusion regimes. We see in Fig.(2c) the typical picture of stochastic resonance. A further increase of the noise intensity produces a loss of coherence and the dynamical behavior is strongly controlled by the noise (Fig.2d). There is therefore a cooperation between periodical driving of the temperature, due to some geological cause, and the environmental noise. The noise can synchronize with the periodical driving giving rise to the stochastic resonance effect, which affect considerably the dynamics of the ecosystem. The exact value of $`\sigma _\beta =1.7810^3`$ used in Fig.2c is obtained from the statistical synchronization typical of the SR phenomenon

$$\tau _k=T_0/2$$
(5)
where $`\tau _k`$ is the Kramers time given by

$$\tau _k=\frac{2\pi }{\sqrt{|U^{\prime \prime }(0.99)|U^{\prime \prime }(0.94)}}\mathrm{exp}[2h/\sigma _\beta ].$$
(6)
and $`T_0`$ is the period of the driving force. In Eq.(6) $`U^{\prime \prime }(0.99)`$ and $`U^{\prime \prime }(0.94)`$ are the second derivative respectively calculated in the unstable and stable states of the potential. To analyze the dynamics of the two species we fix the additive noise intensity at the value $`\sigma _\beta =1.7810^3`$, corresponding to a competition regime between the two species periodically switched from coexistence to exclusion. The temporal series of the two species are obtained for different values of the multiplicative noise intensity $`\sigma =\sigma _x=\sigma _y`$. The initial values of the two species are x(0) = y(0) = 1. After a short transient in the coexistence regime, which is not visible in Fig.3 because of the scale used for the x axes, both species reach the value: $`x_{st}=y_{st}=\alpha /(1+\beta )1/2`$. From then on the densities of both species oscillate quasi regularly around this value with amplitudes dependent on the multiplicative noise intensity. In particular for $`\sigma =0`$ (see Fig.3a) and for very low levels of multiplicative noise ($`\sigma =10^{12}`$) the two species coexist and the populations show correlated oscillations. For higher level of the multiplicative noise ($`\sigma =10^4`$) the amplitude of the oscillations increases, anti-correlated oscillations appear and random periodical inversions of populations occur (see Fig.3c). A further increase of the noise ($`\sigma =10^2`$) produces a degradation of the signal and a loss of coherence of the temporal series for the species (see Fig.3d). In the absence of multiplicative noise the Lotka-Volterra equations are symmetric. The choice of the same values for the initial conditions of the two species densities maintains this symmetry, so the species oscillate in phase around the stationary value $`\alpha /(1+\beta )`$, even if the additive noise that controls the SR effect of the interaction parameter was fixed at a value that produces many transitions to the exclusion regime (see Fig.3a). In the presence of very weak multiplicative noise ($`\sigma =10^{12}`$) the species oscillate in phase almost regularly with some little anticorrelated spikes. By increasing the amplitude of the multiplicative noise we obtain an enhancement of the amplitude of both species, as we can see from Figs. 3b, 3c and 3d. These figures show this transient behaviour from a regime of quasi correlated oscillations (Fig.3b) to a regime of anticorrelated oscillations (Fig.3d). The presence of the multiplicative noise breaks the symmetric dynamical behaviour of the ecosystem. For $`\beta >1`$, i. e. in the exclusion regime, this symmetry breaking determines the ”divergent” behaviour of the trajectories of the two species: one tends to survive, while the other one tends to extinguish. From the inspection of Fig.3 we note that in the presence of a bistable potential modulated by a weak periodic force the response of the system may be enhanced by the presence of the noise. The periodicity of the noise-induced oscillations in the time behavior of the population densities shown in Fig.3e is the same of the driving periodic term of Eq.(3). This is the signature of stochastic resonance (SR). In order to underline the presence of SR we analyze the squared difference of population densities $`(xy)^2`$. In Fig.4 it is shown the SNR of this quantity as a function of the multiplicative noise intensity $`\sigma `$, for $`\sigma _\beta =1.7810^3`$. We calculate the SNR of $`(xy)^2`$ because of the random periodical inversion of population densities. We note that dynamics of $`(xy)`$ is mainly affected by the multiplicative noise, as we can see from Eqs.(2). The SNR has been obtained by performing $`190`$ realizations of Eqs. (2). A maximum at $`\sigma =10^4`$ is present. From the above analysis it is clear the role of the two noise sources: the additive noise determines the conditions for the different dynamical regimes of the two species, the multiplicative noise produces a coherent response of the system by a mechanism of symmetry breaking of the dynamical evolution of the ecosystem.

### 2.2 Noise delayed extinction

Now we investigate the average extinction time of one species as a function of the noise intensity $`\sigma _\beta `$. To do this we fix the multiplicative noise at a low level intensity in such a way that the system is far enough from the SR regime and the dynamics is not strongly perturbed by the noise.

So we choose $`\sigma =10^9`$ because we are not interested in the coherent behavior of the ecological system but we are focused on the effect of the additive noise on the average extinction time of the species. We do simulations of Eqs. (1) and (2) by performing 200 realizations. The results are shown in Fig.(5). We use the bistable potential $`U(\beta )`$ of Fig.1 with the initial condition $`\beta (0)=0.94`$. In this condition the ecosystem is in the coexistence regime, that is the deterministic extinction time of both species is infinite. By introducing noise, exclusion takes place and a finite mean extinction time (MET) appears. By varying the intensity of the additive noise in Eq.3 we obtain, of course, a variation of the average extinction time. The delayed extinction is obtained for noise intensities ranging from the intermediate regime (2 in Fig.5a) to the coexistence regime obtained with higher values of $`\sigma `$ (3 in Fig.5a). In real ecosystems we have always a given noise intensity, which corresponds to a finite mean extinction time. Due to some environmental cause the noise intensity can considerably change, as it is observed in experimental data of populations in a very long time interval . Therefore the dynamical behaviour shown in Fig.5 should explain such physical situations, where the variation of the environmental noise produces a delayed extinction of some population. By increasing the noise intensity we obtain noise delayed extinction and the average extinction time grows reaching a saturation value, which corresponds to a situation where the potential barrier is absent. We find nonmonotonic behaviour of the MET as a function of the noise intensity $`\sigma _\beta `$, with a minimum value $`\tau _{min}=40.47`$ at $`\sigma _\beta =2.7510^3`$, which is of the same order of magnitude of the barrier height $`h`$ (see Fig.5a). The Kramers time corresponding to this noise intensity is $`\tau _k=41.6`$, that is approximately equal to $`\tau _{min}`$. This result is due to the noise driven dynamics. In fact for low value of the noise intensities the average time to overcome the potential barrier is very high, i.e. long Kramers times. The ecosystem remains in the coexistence regime for a long time and the extinction time is very large. For noise intensity of the same order of magnitude of the barrier height the system goes towards the exclusion regime of one of two species and the average extinction time is approximately equal to the Kramers time. We get the minimum value of MET.

For higher values of noise intensity the Kramers time becomes very small and the representative point of the $`\beta `$ parameter moves between the two minima in a very short time. In this condition the system ”sees” the average value of the interaction parameter ($`\beta =0.99`$), which gives a coexistence regime. In Figs.5b,c,d we show the time evolution of the ecosystem corresponding to the points $`1`$, $`2`$ and $`3`$ of Fig.5a. We have coexistence regime in points $`1`$ and $`3`$ and an exclusion regime in point $`2`$. We report finally in the following Fig.6 the behavior of the average extinction time as a function of the Kramers time. This figure confirms our physical description of the extinction dynamics of the system.

## 3 Conclusions

We report a study on the role of the noise in the dynamics of two competing species. We consider two noise sources: a multiplicative noise and an additive noise, which produces a random interaction parameter between the species. The noise induces a coherent time behavior of two species giving rise to temporal oscillations and enhancement of the response of the system through stochastic resonance phenomenon. Specifically the additive noise controls the switching between the coexistence and the exclusion dynamical regimes, the multiplicative noise is responsible for coherent oscillations of the two species. The SR in the dynamics of interaction parameter $`\beta `$ induces SR phenomenon in two competing species. These time behaviors are absent in the deterministic dynamics. We note that our model is useful to describe physical situations in which the amplitude of the periodical driving force, due to the temperature variations, is weak and therefore unable to produce considerable variations of the dynamical regime of the ecosystem. The synergetic cooperation between the nonlinearity of the system and the random and periodical environmental driving forces produces therefore a coherent time behaviour of the ecosystem investigated. The noise is also responsible for a delayed extinction which gives rise to a nonmonotonic behavior of the average extinction time as a function of the additive noise intensity. We note finally that these noise induced effects should be useful to explain the time evolution of ecological species, whose dynamics is strongly affected by the noisy environment .

## 4 Acknowledgments

The authors are grateful to Dr. A. La Barbera for useful discussions concerning numerical simulations. This work was supported by INFM, MIUR and INTAS Grant 01-450.


text:
# Lagrangian Velocity Statistics in Turbulent Flows: Effects of Dissipation

## Abstract

We use the multifractal formalism to describe the effects of dissipation on Lagrangian velocity statistics in turbulent flows. We analyze high Reynolds number experiments and direct numerical simulation (DNS) data. We show that this approach reproduces the shape evolution of velocity increment probability density functions (PDF) from Gaussian to stretched exponentials as the time lag decreases from integral to dissipative time scales. A quantitative understanding of the departure from scaling exhibited by the magnitude cumulants, early in the inertial range, is obtained with a free parameter function $`D(h)`$ which plays the role of the singularity spectrum in the asymptotic limit of infinite Reynolds number. We observe that numerical and experimental data are accurately described by a unique quadratic $`D(h)`$ spectrum which is found to extend from $`h_{min}0.18`$ to $`h_{max}1`$, as the signature of the highly intermittent nature of Lagrangian velocity fluctuations.

Statistical properties of homogeneous three dimensional turbulence have been studied for a long time in the Eulerian framework FrischBook . Recently a growing interest in studying intermittency from a dynamical point of view has been motivated by high precision Lagrangian experiments. Essentially two experimental groups have performed particle tracking in highly turbulent flows. The group at Cornell JFMBoden  reports measurements of Lagrangian acceleration in a turbulent water flow between two counter-rotating disks for Taylor-based Reynolds numbers $`200<R_\lambda <900`$. The experiment carried out at ENS-Lyon MordantPint , in a similar von Kármán flow, is based on acoustic tracking. It provides Lagrangian velocity records covering the inertial range of turbulent motion, up to several integral time scales. In addition to these complementary experiments, DNS of the Navier-Stokes equations MordantPint ; PopeYeung  have produced comparative numerical results in the range $`75<R_\lambda <380`$. Both experimental and numerical studies have revealed the existence of a very strong intermittency in the Lagrangian dynamics: the PDF of the velocity is Gaussian, while the PDF of the acceleration exhibits extremely large tails. In between, the PDFs of velocity increments changes continuoulsy between these two functional forms as the time lag is decreased from integral times (the Lagrangian auto-correlation time) to dissipative ones (below the Kolmogorov time scale). Several stochastic models have been proposed that reproduce the behavior of the Lagrangian acceleration SawfordPOF  and of the Lagrangian velocity increments Rey03 . They rely on different physical assumptions and all include several ad-hoc hypothesis and parameters in order to fit the experimental observations — a rather comprehensive review can be found in AriMaz03 . The aim of the present work is not to add to these models but to provide a comprehensive description of the Lagrangian intermittency, using a formalism that describes both the inertial and dissipative range of time scales. It is motivated by the desire to analyze globally the Lagrangian scales, which are likely to be dynamically connected (for example it is shown in MordantPOF  that the Lagrangian acceleration is influenced by the large scale dynamics, a feature that is also pointed out by models based on rapid distortion ideas LavDubNaz01 ). In addition we require that the description of intermittency be independent on model assumptions. The multifractal description FrischParisi , already widely used in Eulerian studies of turbulence, is a natural choice. To encompass inertial and dissipative features, we recast the multifractal picture of the intermediate dissipative range, originally proposed for Eulerian velocity fluctuations PalVul87 ; Men96 , into the context of Lagrangian velocity. The resulting analysis provides a synthetic and comprehensive description of the various sets of experimental and numerical data.

In the present description, a first-order Lagrangian velocity increment over a time scale $`\tau `$ is written as :

$$\delta _\tau v(t)=v(t+\tau )v(t)=\beta \left(\tau /T\right)\delta _Tv,$$
(1)
where all the time scale dependence is contained in the independent random function $`\beta (\tau /T)`$ ($`>0`$ since Lagrangian velocity increment PDFs are symmetric), $`T`$ being the integral (Lagrangian) time scale above which velocity increments become uncorrelated. The PDF of integral time scale increments $`\delta _Tv`$ is thus assumed to be Gaussian ($`𝒢`$) — a result of a central limit argument, also in agreement with Eulerian observations. Once the distribution of $`𝒫(\beta )`$ is known, the PDF of increments at any time scale $`\tau `$ is computed as

$$𝒫(\delta _\tau v)=\frac{d\beta }{\beta }𝒢\left(\frac{\delta _\tau v}{\beta }\right)𝒫(\beta ).$$
(2)
In the standard mutifractal formalism FrischParisi , $`\beta `$ is assumed to have a power law scale dependence in the inertial range, $`\beta (\tau /T)^h`$, with a spectrum $`𝒟(h)`$ (meaning that the PDF of observing an exponent $`h`$ at scale $`\tau `$ is proportional to $`(\tau /T)^{1𝒟(h)}`$). In the inviscid limit, $`h`$ and $`𝒟(h)`$ acquire the mathematical status of H$`\ddot{\text{o}}`$lder exponent and singularity spectrum, respectively. Note that this description based on first order increments is restricted to exponents $`h<1`$ MuzBacArn93 , a limitation which will be addressed later. To describe the entire range of scales covered in experimental measurements and computer simulations, one must take into account the effects of viscosity (finite $`R_\lambda `$). In the dissipative range, velocity fluctuations are smoothed by viscous damping (or by measurement filtering) and the velocity increments become proportional to the time scale $`\delta _\tau v(t)=\tau a(t)`$, where $`a(t)`$ is the Lagrangian acceleration. We shall consider that the cross-over between inertial and dissipative statistics occurs when the local Reynolds number is of order unity,

$$Re(\tau /T)=\frac{\tau }{T}\beta ^2\left(\frac{\tau }{T}\right)Re=1,$$
(3)
where $`Re`$ is the integral scale Reynolds number. This defines a local Kolmogorov dissipative time $`\tau _\eta (h)=TRe^{\frac{1}{2h+1}}`$ where the local velocity increments change from inertial scale invariance to dissipative scaling (this implies $`h1/2`$, for time scales to be shorter than $`T`$). Changing the integration variable from $`\beta `$ to $`h`$ in Eq.(2), the PDF of velocity increments at any scale $`\tau /T`$ can be written as the sum of two contributions

$`𝒫(\delta _\tau v)={\displaystyle _{1/2}^{h^{}(\frac{\tau }{T},Re)}}𝑑h{\displaystyle \frac{𝒫_i(h,\frac{\tau }{T},𝒟(h))}{\beta _i(\frac{\tau }{T},h)}}𝒢\left({\displaystyle \frac{\delta _\tau v}{\beta _i(\frac{\tau }{T},h)}}\right)`$
$`+{\displaystyle _{h^{}(\frac{\tau }{T},Re)}^+\mathrm{}}𝑑h{\displaystyle \frac{𝒫_d(h,Re,𝒟(h))}{\beta _d(\frac{\tau }{T},h,Re)}}𝒢\left({\displaystyle \frac{\delta _\tau v}{\beta _d(\frac{\tau }{T},h,Re)}}\right)`$ (4)
where the functions $`\beta _{i,d}`$ and $`𝒫_{i,d}`$ have the proper inertial ($`\beta _i(\tau /T)^h,𝒫_i(\tau /T)^{1𝒟(h)}`$) and dissipative ($`\beta _d\tau /T,𝒫_d(\tau _\eta (h)/T)^{1𝒟(h)}`$) scalings. The change occurs at the critical value $`h^{}`$ for which the local Reynolds number is unity:

$$h^{}(\frac{\tau }{T},Re)=\frac{1}{2}\left(1+\frac{\mathrm{ln}Re}{\mathrm{ln}\frac{\tau }{T}}\right)\text{ .}$$
(5)
For $`h<h^{}(\tau /T,R_e)`$, the increments $`\delta _\tau v`$ are in the inertial range, while they lie in the dissipative range for $`h>h^{}`$. Finally, we impose that the function $`\beta (\tau /T)`$ be continuous and differentiable at the transition, following a strategy used in the Eulerian domain Men96 , and inspired from an elegant interpolation formula proposed by Batchelor Bat51 . In this framework, a single function $`\beta (h;\tau /T,Re)`$ covers the entire range of scale

$$\beta (h;\frac{\tau }{T},Re)=\frac{\left(\frac{\tau }{T}\right)^h}{\left(1+\left(\frac{\tau }{\tau _\eta (h)}\right)^\delta \right)^{\frac{1h}{\delta }}}\text{ ,}$$
(6)
with probability density function

$$𝒫(h;\frac{\tau }{T},Re,𝒟(h))\frac{\left(\frac{\tau }{T}\right)^{1𝒟(h)}}{\left(1+\left(\frac{\tau }{\tau _\eta (h)}\right)^\delta \right)^{\frac{𝒟(h)1}{\delta }}}.$$
(7)
The PDFs of Lagrangian velocity increments are then computed using Eqs.(4-7), and compared to the experimental / numerical data. The entire range of time scales is covered when the four parameters $`T,Re,\delta `$ and $`𝒟(h)`$ are prescribed. $`T`$ and $`Re`$, both global parameters imposed by the flow production, are explicitly incorporated. $`𝒟(h)`$ is a parameter function which must be derived from the experimental or numerical data. In this work, and as a posteriori justified, we assume a quadratic form $`𝒟(h)=1(hc_1)^2/2c_2`$ (a lower order approximation customary in intermittency studies). The constants $`c_1`$ and $`c_2`$ are constrained to satisfy $`c_1=1/2+c_2`$ in order for Kolmogorov scaling to be observed in the inertial range ($`\delta _\tau v^2\tau `$). $`\delta `$ is a free parameter that accounts for the smoothing of the transition from inertial to dissipative scales. In his modeling of Eulerian velocity structure function, Batchelor Bat51  sets $`\delta =2`$, a value used in other Eulerian studies Men96  and supported by theoretical arguments SirYakSmi94 . In practice $`T`$ and $`Re`$ are computed using the experimental data (up to some empirical multiplicative constants), and the free parameters $`c_2`$ and $`\delta `$ are estimated using a least square minimization scheme.

In Fig.1(a) we show the results for the experimental data from the ENS-Lyon group ($`R_\lambda =740`$) MordantPint . The velocity increment PDFs are well reproduced at all scales for the set of parameter values $`\delta =1.08`$, $`c_1=0.575`$ and $`c_2=0.075`$. This agreement is emphasized in Fig.1(c), where the fourth order moment $`(\overline{\delta _\tau v})^4𝒫(\overline{\delta _\tau v})`$ is shown to peak at the same values as the experimental curves and this for all considered time scales. In Figs.1(a,c) we have also included the Cornell acceleration data at $`R_\lambda =690`$ JFMBoden . What is remarkable is that with quite consistent parameter values, namely $`\delta =1.3`$ and still a parabolic $`𝒟(h)`$ curve with $`c_1=0.579`$ and $`c_2=0.079`$, one again reproduces the experimental acceleration PDF with a great accuracy in the tails; in particular, the flatness $`F=a^4/a^2^2=56.1`$ is very close to the experimental value $`F=55\pm 8`$ JFMBoden . For our DNS data at $`R_\lambda =140`$ (Figs.1(b,d)), the agreement is even more pronounced than for the experimental data. In particular, it is hard to distinguish the fit from the data points in Fig.1(d). Note that the parameter values obtained to get the best fit of the numerical data are not significantly different from those used in Figs.1(a,c): $`𝒟(h)`$ corresponds to $`c_1=0.586`$ and $`c_2=0.086`$, while $`\delta =1.98`$ is now quite close to the Batchelor value $`\delta =2`$ Bat51 . The fact that $`\delta `$ is smaller in experimental measurements is due to the filtering induced by the finite size of the tracer particle and signal processing algorithms.

A first finding of our analysis is that almost identical functions $`𝒟(h)`$ are obtained for the three sets of data (Cornell, Lyon and DNS), although they cover a wide range of scales and of turbulent Reynolds numbers: the symbols in Fig.2 are undistiguishable, certainly within error bars. This suggests that such a parabolic $`𝒟(h)`$ curve should also be relevant in the limit of infinite Reynolds number. A parabolic singularity spectrum is the hallmark of log-normal statistics as originally proposed by Kolmogorov and Obukhov Kol62  for Eulerian velocity statistics. Our findings for Lagrangian velocity statistics is a parabolic $`𝒟(h)`$ spectrum centered at $`c_1=0.58\pm 0.01`$, significantly larger than the K41 value $`1/2`$, and of width $`c_2=0.08\pm 0.01`$ (commonly called the intermittency exponent). This is significantly larger than the corresponding value $`c_2^E=0.025\pm 0.003`$ derived for Eulerian velocity data DelMuzArn01 . It corroborates the fact that Lagrangian velocity statistics are more intermittent than Eulerian velocity statistics MordantPint . Let us remark that the ratio of the Lagrangian to Eulerian intermittency exponents $`c_2/c_2^E=0.08/0.025=3.2\pm 0.2`$ is very close to the value $`\left(3/2\right)^3=3.375`$ which can be predicted using a Kolmogorov-Richardson argument MordantPint ; Bor93 .

To emphasize further the quality of the proposed description of Lagrangian intermittency, we now consider the evolution across the scales of the moments of the velocity increments, the so-called structure functions FrischBook , $`S(p,\tau )=<|\delta _\tau v|^p>`$. As advocated in Ref.DelMuzArn01  for Eulerian velocity data analysis, the magnitude cumulant analysis provides a more reliable alternative to the structure function method. Indeed, it is straightforward to derive the relationship between the moments of $`𝒫(|\delta _\tau v|)`$ and the cumulants $`C_n(\tau )`$ of $`𝒫(\mathrm{ln}|\delta _\tau v|)`$:

$$<|\delta _\tau v|^p>=\mathrm{exp}\left(\underset{n}{}C_n(\tau )\frac{p^n}{n!}\right).$$
(8)
In the inertial range, where multifractal power law scaling is expected to be observed for the structure functions, the cumulants should behave like $`C_n(\tau )c_n\mathrm{ln}\tau `$. In Fig.3, we report the results of the computation of the cumulant $`C_n(\tau )`$ for both the ENS-Lyon experimental data (Figs 3(a-c)) and our DNS numerical data (Figs 3(d-f)). We have explicitly subtracted the cumulant of the $`\delta _Tv`$ Gaussian PDF and the curves are computed using the same parameter set as in Fig.1. For both the experimental and numerical cumulants, it is quite convincing that our multifractal decription provides a comprehensive understanding of the observed departure from scaling when going from large $`\tau /T1`$ to small $`\tau /T10^2`$ time scales. In fact, an unambiguous inertial range is observed for the experimental second order cumulant only (Fig.3(b)). For the first order cumulant, both the experimental (Fig.3(a)) and numerical (Fig.3(d)) curves display some curvature over the entire range of scale, a feature that is well reproduced by our extended multifractal description. For $`C_2(\tau /T)`$, the model describes quite well the cross-over behavior observed in the experimental (Fig.3(b)) and numerical (Fig.3(e)) data down to time scales $`\tau /T`$ of the order of the smoothing filtering scale (including the finite size of the tracer particles). The predicted plateau $`C_2(\tau /T)\overline{C}_2(R_\lambda )`$ in the limit $`\tau /T0`$, is reached in the numerical but not in the experimental data. For $`C_3(\tau /T)`$, we note that the DNS data in the inertial range are quite well reproduced taking $`c_3=0`$, at odds with the experiment.

We now return to the assumption that the scaling exponent $`h`$ should remain smaller than one. It may be seen in Fig.2 that the largest measured $`h`$ values reach 1. In this case, it is known that first order increments are not adapted to detect singularities with exponents $`h1`$ MuzBacArn93 . We have thus repeated our analysis using second order velocity increments $`\delta _\tau ^{(2)}v(t)=v(t+\tau )2v(t+\tau /2)+v(t)`$. Eqs.(4,6,7) must be slightly modified because (i) the dissipative range should now scale as $`(\tau /T)^2`$ and (ii) $`h^{}`$ may be larger than one and the integrals must be computed in the order of increasing $`h`$’s. When using the same parameters as before, we observe in Fig.3 that the experimental and numerical cumulants are robustly and even better reproduced. Let us point out that as seen in Figs. 3(b,c,e,f), because the $`C_n(\tau /T)`$ for $`n2`$ are predicted to reach a plateau $`\overline{C}_n(R_\lambda ,N)(1+2N)^n\mathrm{ln}R_\lambda `$ in the limit of vanishing $`\tau `$’s, then the experimental and numerical cumulants for the second-order increments are much more affected by the smoothing filtering process in the intermediate dissipative range than previously observed for the first-order increments.

To conclude, we return to our observation that a unique $`𝒟(h)`$ spectrum yields an accurate description of the Lagrangian velocity statistics at all scales. Such a spectrum $`𝒟^E(h)`$ has been extensively studied in the Eulerian domain FrischBook . Two widely used forms (corresponding to log-normal and log-Poisson statistics) are shown in Fig.2. They can be mapped into the Lagrangian domain

$$𝒟(h)=h+(1+h)𝒟^E\left(h/(1+h)\right),$$
(9)
using a Kolmogorv Refined Similarity argument in the spirit of the work done by Borgas Bor93  . The resulting curves are shown in Fig.2; we note that the agreement with the measured Lagrangian $`𝒟(h)`$ functions is excellent on the left-hand side of the curves, i.e. for values $`h<c_1`$ corresponding to intense velocity increments. On the right-hand side ($`h>c_1`$) there is a noticeable difference. Whether this difference is significant deserves more investigation. It may be of importance since the above relationship clearly shows that the Eulerien and Lagrangian singularity spectra cannot be both log-normal.

Numerical simulations were performed at CINES (France) using an IBM SP computer. We wish to acknowledge Bernard Castaing for his critical comments. We are grateful to the Cornell team for making their data available.


text:
# The spectra of short Gamma–Ray Bursts

## 1 Introduction

The possible existence of different classes of GRBs was considered since their discovery, and the strongest evidence for different populations is their bimodal duration distribution, with $``$1/3 of ‘short’ events with a mean duration of $``$0.3 s, and the majority of ‘long’ events with mean duration of $``$ 20 s (Kouveliotou et al. Kouveliotou1993 (1993), Norris et al. Norris2000 (2000)). Further support to such a bimodal behavior emerges from the analysis of their spectral and temporal properties: short bursts seem to be harder than long ones (Kouveliotou et al. Kouveliotou1993 (1993), Hurley et al. Hurley1992 (1992)) and their distributions of pulse width, separation and number of pulses per bursts also indicate that the two classes might be physically distinct (Norris et al. Norris2000 (2000), Nakar & Piran Nakar2002 (2002)). The distinction between short and long bursts has been also considered as indication of the existence of two distinct progenitors. If the duration of the GRB emission (as predicted by the internal shock model - see e.g. Piran Piran1999 (1999) for a review) is linked to the duration of the inner engine activity, short bursts might be produced by the merger of compact objects (Ruffert & Janka Ruffert1999 (1999)) while the core collapse of massive stars would give raise to long duration GRBs. While the properties of long events (e.g. redshifts, broad band spectral emission and evolution, environment etc., see Hurley et al. Hurley 2003 (2003) for a recent review) have been unveiled with increasing details, the understanding of short bursts is still limited. Recently, Schmidt (Schmidt2001 (2001)) suggested that short bursts have a similar luminosity to long events. So far the characterization of the spectral properties of short bursts detected by BATSE has been based on the comparison of the ratio of the fluxes emitted in different (broad) energy bands (Cline et al. 1999; Yi–ping Qui 2001). The spectrum of long GRBs, typically represented by smoothly connected power laws, is different for different bursts (Band et al. Band1993 (1993)) and it may also considerably evolve in time within the same burst (Ford et al. Ford1995 (1995), Crider et al. Crider1997 (1997)). This complex behaviour compels to consider the complete spectrum of any GRB with high time and spectral resolution in order to describe and compare the emission properties of long and short bursts. Clearly the main difficulty when fitting short burst spectra is their low signal to noise ratio due to their small duration. Paciesas et al. (Paciesas2001 (2001)) compared the spectral parameters of short and long bursts obtained from spectral fits: they found that in short bursts the low energy spectral index and the peak energy are harder than in long events. However, the time resolution (2 s) of the spectra used to describe the class of short bursts was much larger than their typical duration (0.3 s) and also the spectral resolution of their data was low compared to that of the data presented in this work. Also the distance scale of short GRBs is still a matter of debate. Their spatial distribution seems to be consistent with low redshift sources (e.g. Magliocchetti et al. Magliocchetti2003 (2003), Che et al. Che1997 (1997)), but nothing is directly known due to the lack of any redshift measurement. On the other hand, possible correlations among the spectral properties of long bursts have been recently claimed (Amati et al. Amati2002 (2002), see also Yonetoku et al. Yonetoku2003 (2003)) and confirmed by the Hete–II long GRBs and X–Ray Flashes (Lamb 2003a ). These relations might be key in explaining the still obscure energy conversion mechanism operating in long GRBs and it is thus interesting to investigate whether similar correlations also hold for short bursts.

## 2 Sample selection and spectral analysis

From the BATSE on–line catalog111http://cossc.gsfc.nasa.gov/cossc/batse/ we have selected bright short ($`T_{90}2`$ s) events as those with a peak flux (computed on the 64 ms timescale and integrated over the energy range 50–300 keV) exceeding 10 phot cm-2 s-1. The hardness ratios of these 36 selected short GRBs, computed over the energy ranges 1=25–50 keV, 2=50–100 keV, 3=100–300 keV, are $`HR_{21}=1.65`$ and $`HR_{32}=5.9`$. These values are larger than the corresponding HRs of the population of long bursts ($`HR_{21}=1.56`$ and $`HR_{32}=3.8`$), in agreement with the hardness–duration relation (Kouveliotou et al. Kouveliotou1993 (1993)). The average fluences (for energies $`25`$ keV) of the selected short GRBs is $`6.2\times 10^6`$ erg s-1 (i.e., only a factor $``$ 2.5 lower than for long bursts). We analyzed their spectra using the Large Area Detector (LAD) data and applied the standard spectral analysis technique (e.g. Preece et al. Preece2000 (2000), hereafter P00). Each spectrum was fitted over the energy range $`30`$ keV – 1.8 MeV. Due to their short duration, the minimum integration time (S/N limited) of the LAD spectra is typically $`0.20.4`$ s, so that we could analyze in most cases one single spectrum per short GRB. In 7 cases out of 36, the spectrum had a low S/N over the entire energy range which resulted in poorly constrained best fit parameters. These cases were not included in the final sample. In one case the analysis was not possible because of missing data. The remaining 28 sources (with 100-300 keV fluence $`2.4\times 10^7`$ erg/cm2), although belonging to a complete peak flux limited sample, do not form themselves a complete sample.

The spectral properties of short bursts have been compared with the results of the spectral analysis of a sample of bright long BATSE bursts (Ghirlanda et al. 2002 – hereafter G02) whose spectra (time averaged and time resolved) were fitted with different spectral models. The sample of G02 was selected on the basis of the burst peak flux, similarly to the criterion applied for the sample of short events presented in this work.

### 2.1 Spectral results

The spectrum of the selected short bursts is in most cases well fitted by a single power law with an exponential cutoff at high energies. The spectral parameters are the low energy power law photon index $`\alpha `$ and the peak energy $`E_{\mathrm{peak}}`$ (in keV) of the $`EF_E`$ spectrum. This model fit resulted in a lower reduced $`\chi ^2`$ and in better constrained spectral parameters compared with the Band function, single power law and broken power law models. In most cases the statistics in the high energy channels of the spectra is too poor to constrain the high energy power law component of the Band model (Band et al. Band1993 (1993)). In 5 cases out of the 28 analyzed the reduced $`\chi ^2`$ of the powerlaw cutoff model is high and excludes the fit at 99% confidence level. In these cases, however, the Band model resulted in even higher reduced $`\chi ^2`$ and we adopted the parameters of the powerlaw cutoff model for homogeneity with the rest of the sample. Tab. 1 222This is available in electronic form at ……..  lists our selected bursts, together with the spectral results and their errors at 99% confidence level. As also found by Paciesas et al. (Paciesas2001 (2001)), there is no evidence for any correlation between $`\alpha `$, $`E_{\mathrm{peak}}`$ and the burst duration represented by the $`T_{90}`$ parameter given in the BATSE catalog of Paciesas et al. (Paciesas1999 (1999)). The short duration of most bursts ($`T_{90}[0.104,1.8]`$ s) does not allow a time resolved spectral analysis, but in 5 cases we could extract at least two (in one case even three) time resolved spectra within the $`T_{90}`$ interval. Nonetheless, the low statistics of these time resolved spectra results in large uncertainties on the spectral parameters, with only a weak indication of a hard–to–soft spectral evolution similar to what found in long GRBs (e.g. Ford et al. Ford1995 (1995)). This supports the trend found from the analysis of the time resolved hardness ratio (e.g. Cline et al. Cline1999 (1999)).

In 6 short GRBs the low energy spectrum is harder than the optically thin synchrotron limit $`\alpha 2/3`$ (Katz 1994). In one case $`\alpha >0`$. A similar fraction of long bursts violating these limits has been found in the population of long GRBs (Crider et al. Crider1997 (1997), P00, Ghirlanda et al. Ghirlanda2003 (2003)).

## 3 Short vs Long GRBs

We compared the spectral properties of our short GRBs with those of the long–bright events of G02. For homogenity we considered the spectral parameters of long bright bursts obtained from the fits of the same model, i.e. a cutoff–powerlaw, adopted for the short ones.

Firstly we compare the time integrated spectral parameters. The distributions of $`\alpha `$ and $`E_{\mathrm{peak}}`$ for long and short bursts are reported in Fig. 1 (top panels). Long bursts have an average peak energy ($`E_{\mathrm{peak}}=520\pm 90`$ keV) slightly larger than that of short events ($`E_{\mathrm{peak}}=355\pm 30`$ keV). The distribution of the latter ones spreads between a few keV and a few MeV. A Kolmogorov–Smirnov (KS) test gives a probability of $``$0.8% that the two samples are drawn from the same parent population. A more significant diversity results from the comparison of the distributions of $`\alpha `$: the average values are $`\alpha =1.05\pm 0.14`$ and $`\alpha =0.58\pm 0.10`$ for long and short GRBs, respectively (KS probability of 0.04%).

Then we conclude that the hardness–duration relation, discussed for the BATSE bursts (Kouveliotou et al. Kouveliotou1993 (1993), Cline et al. 1999, Yi-ping Qui 2001), is caused more by short bursts having a harder low energy spectral slope rather than a higher peak energy. Indeed, if $`E_{\mathrm{peak}}`$ was the only parameter characterizing the spectral hardness, short events would be softer than long ones. This conclusion is further supported by the search of correlations between the hardness ratios and the spectral parameters of our short bursts. The only significant trend is of harder $`\alpha `$ for larger HR32 (Spearman’s correlation coefficient $`r0.6`$ with null hypothesis probability of $`10^3`$). A larger sample of bright BATSE long bursts was studied by Preece et al. (Preece2000 (2000)). Although their sample has not been fitted homogeneously with the same model their spectral results also confirm our finding that short bursts are harder than long GRBs because of a harder index $`\alpha `$ rather than a higher $`E_{\mathrm{peak}}`$.

The average spectral parameters give only an indication of the spectrum, which is likely to evolve in time, in all its spectral parameters (e.g., $`\alpha `$, $`E_{\mathrm{peak}}`$ etc., Crider et al. Crider1997 (1997)). Therefore, in order to test the tantalizing hypothesis that the spectrum of short bursts is similar to the initial emission phases of long events, we considered the time resolved spectral parameters reported for the G02 sample relative to the first seconds since the burst onset.

The spectrum of short bursts, as described by $`\alpha `$ and $`E_{\mathrm{peak}}`$, is more similar to the first 2 s than the integrated spectrum of long events, as shown in Fig. 1 (bottom panels): the low energy spectral index distributions are similar (with a KS probability of 83%) while the peak energy distributions still indicate that short events are softer than long ones (with a KS probability of 10%).

Another appealing possibility is that short bursts might be similar to the peak spectra of long bursts (i.e. they represent the ‘tip of the iceberg’ in the long GRB light curve). To this aim we extracted from the catalog of Preece et al. (2000) only those spectra accumulated around the peak with an integration time (centered around the peak time) at least as long as the average duration of short bursts (0.3 s). From the comparison of the spectral parameters distributions we conclude that short bursts are still different – especially for the $`E_{\mathrm{peak}}`$ value – from the peak of long GRBs as also indicated by the small KS probability (0.4%). The low energy spectral index instead, similarly to what found from the comparison with the first 2 seconds of long events, presents a distribution similar to that for the peak spectra of long events (with a KS probability of 23%). Moreover, we stress that if the fits were performed with the Band model (regardless of the indetermination of the high energy spectral slope) we should find a systematically lower peak energy than what found with the powerlaw cutoff model (e.g. Preece et al. Preece2000 (2000), Ghirlanda et al. Ghirlanda2002 (2002)). This would strengthen our conclusions.

## 4 The energy and luminosity of short bursts

To further test the relationship between short and long GRBs, we considered the recently found spectral correlations for long GRBs between $`E_{\mathrm{peak}}`$ and the equivalent isotropic energy $`E_{\mathrm{iso}}`$ in the source reference frame (Amati et al. Amati2002 (2002)). A similar result was found from a sample of BATSE bursts (Lloyd et al. Lloyd2000 (2000), Bloom et al. Bloom2001 (2001)) and also by the Hete–II long bursts, with the inclusion of 2 X–Ray Flashes, extending this relation to low energies (Lamb et al. 2003a ). Moreover Yonetoku et al. (Yonetoku2003 (2003)) found a similar correlation between the isotropic luminosity $`L_{\mathrm{iso}}`$ and $`E_{\mathrm{peak}}`$ of the BATSE and BeppoSAX samples. If short GRBs are similar to long events they might satisfy the above correlations. However, the redshift of short bursts is unknown. Thus we can only verify whether the observed spectral properties of short bursts, scaled in the source rest frame (for any $`z`$), are consistent with the above correlations. In other words, from the spectral fits of short bursts we can derive the peak energy ($`E_{\mathrm{peak}}^{\mathrm{obs}}`$) and the fluence ($`F^{\mathrm{obs}}`$) integrated over the same band used by Amati et al. (Amati2002 (2002)) and Bloom et al. (2001), i.e. 1 keV – 10 MeV. These can be converted in the source reference frame for any redshift: $`E_{\mathrm{iso}}(z)=F^{\mathrm{obs}}(4\pi D_L^2)/(1+z)`$, $`E_{\mathrm{peak}}(z)=E_{\mathrm{peak}}^{\mathrm{obs}}(1+z)`$ and $`L_{\mathrm{iso}}(z)=F^{\mathrm{obs}}\times (4\pi D_\mathrm{L}^2)/T_{90}`$ \[for cosmological parameters $`(H_0,\mathrm{\Omega }_\mathrm{\Lambda },\mathrm{\Omega }_\mathrm{m})=(65,0.7,0.3)`$\]. These were compared with the relation $`E_{\mathrm{iso}}E_{\mathrm{peak}}^{1.93}`$ found by Amati et al. (Amati2002 (2002)) for a redshift range 0.001-10. The same was done for the luminosity $`L_{\mathrm{iso}}(z)`$. We find that for any given $`z`$ (up to 10) short bursts (except for 2 cases) populate a region below the $`E_{\mathrm{iso}}`$$`E_{\mathrm{peak}}`$ correlation of long GRBs. The luminosity of short events 333 Note that $`L_{\mathrm{iso}}`$ for short bursts is obtained from the spectrum integrated over their typical duration ($``$0.5 sec). Similarly Yonetoku et al. (2003) derive this luminosity for long bursts either from the spectrum integrated around the GRB peak on a comparable timescale (for BATSE GRBs) or by dividing for the burst duration (for SAX GRBs).  is instead consistent with the correlation proposed for long bursts (Yonetoku et al. Yonetoku2003 (2003)). Following the spectral results suggesting that short GRB might be similar to the first $`1`$ s of long ones, we scaled the $`L_{\mathrm{iso}}`$$`E_{\mathrm{peak}}`$ correlation (computed for the peak spectra of long events) to the first second of their emission by considering the typical ratios between fluxes and $`E_{\mathrm{peak}}`$ estimated at the peak and those integrated over the first second. The curves of $`L_{\mathrm{iso}}(z)`$ of short bursts are still consistent with this relation. Although neither obvious nor unique to interpret, these results are at least consistent with the hypothesis that the luminosity of short GRBs is comparable to that of long bursts within the first seconds, with instead a lower isotropic equivalent energy. In this case the different duration might be responsible for a lower total energy emitted in short bursts. Under such hypothesis it is also possible to tentatively infer the redshift distribution of (our) short bursts by assuming that they satisfy the $`L_{\mathrm{iso}}`$$`E_{\mathrm{peak}}`$ correlation. The resulting $`z`$ distribution (Fig. 2) is compared to that of a sample of long bursts for which the redshift has been inferred from the claimed correlation between the variability properties of their prompt emission light curve and their luminosity (Fenimore & Ramirez–Ruiz 2001). Short bursts have a slightly smaller average redshift compared to the sample of long events.

We stress that this result is heavily dependent on the assumptions that short GRBs i) have the same luminosity of long ones, and ii) obey the same $`L_{\mathrm{iso}}`$$`E_{\mathrm{peak}}`$ correlation followed by long GRBs. Both these assumptions will be tested when a reasonable number of redshift of short GBRs will be known.

## 5 Conclusions

We selected a sample of short bright bursts from the BATSE catalog and performed a standard spectral analysis in order to characterize their spectrum. No correlation was found between the spectral parameters and the global properties (duration and flux) of these bursts. The low energy spectral index $`\alpha `$ is distributed between –2 and small positive values and the peak energy $`E_{\mathrm{peak}}`$ is between 20 keV and 2 MeV. Similarly to long bursts (e.g. Preece et al. 2000), some short bursts have a low energy spectrum harder than the optically thin synchrotron limit.

The comparison of the spectral properties of short and long GRBs (G02 sample) revealed that: (i) the higher hardness of short GRBs with respect to long events (typically described in terms of fluence hardness ratio) is the effect of a harder low energy spectrum rather than of a marginally different peak energy; (ii) the spectral properties of short bursts are similar to those of the first 2 s of long GRBs.

Short bursts are then harder than the time–average spectra of long GRBs, but their properties are compatible with a similar mechanism operating at the beginning of all bursts, independently of their duration.

Short bursts cannot obey the energy–peak frequency correlation found for long bursts (Amati et al. Amati2002 (2002), Lamb et al. 2003a ), but they could obey the similar correlation found by Yonetoku et al. (Yonetoku2003 (2003)) between the luminosity and the peak frequency. This may suggest that short bursts have the same luminosity, but lower total energy, than long bursts. If this is the case the redshift distribution inferred for the short bursts of our sample is similar to that of long events, with a slightly smaller average redshift.

###### Acknowledgements.

This research has made use of data obtained through the High Energy Astrophysics Science Archive Research Center Online Service, provided by the NASA/Goddard Space Flight Center. We thank F. Tavecchio for discussion and the referee for useful comments.