text:
# Long-term persistence and multifractality of river runoff records: Detrended fluctuation studies

## 1 Introduction

The analysis of river flows has a long history. Already more than half a century ago Hurst found by means of his $`R/S`$ analysis that annual runoff records from various rivers (including the Nile river) exhibit ”long-range statistical dependencies” (Hurst, 1951), indicating that the fluctuations in water storage and runoff processes are self-similar over a wide range of time scales, with no single characteristic scale. Hurst’s finding is now recognized as the first example for self-affine fractal behaviour in empirical time series, see e.g. Feder (1988). In the 1960s, the ”Hurst phenomenon” was investigated on a broader empirical basis for many other natural phenomena (Hurst et al., 1965; Mandelbrot and Wallis, 1969).

The scaling of the fluctuations with time is reflected by the scaling of the power spectrum $`E(f)`$ with frequency $`f`$, $`E(f)f^\beta `$. For stationary time series the exponent $`\beta `$ is related to the decay of the corresponding autocorrelation function $`C(s)`$ of the runoffs (see Eq. (1)). For $`\beta `$ between 0 and 1, $`C(s)`$ decays by a power law, $`C(s)s^\gamma `$, with $`\gamma =1\beta `$ being restricted to the interval between 0 and 1. In this case, the mean correlation time diverges, and the system is regarded as long-term correlated. For $`\beta =0`$, the runoff data are uncorrelated on large time scales (”white noise”). The exponents $`\beta `$ and $`\gamma `$ can also be determined from a fluctuation analysis, where the departures from the mean daily runoffs are considered as increments of a random walk process. If the runoffs are uncorrelated, the fluctuation function $`F_2(s)`$, which is equivalent to the root-mean-square displacement of the random walk, increases as the square root of the time scale $`s`$, $`F_2(s)\sqrt{s}`$. For long-term correlated data, the random walk becomes anomalous, and $`F_2(s)s^H`$. The fluctuation exponent $`H`$ is related to the exponents $`\beta `$ and $`\gamma `$ via $`\beta =1\gamma =2H1`$. For monofractal data, $`H`$ is identical to the classical Hurst exponent. Recently, many studies using these kinds of methods have dealt with scaling properties of hydrological records and the underlying statistics, see e. g. Lovejoy and Schertzer (1991); Turcotte and Greene (1993); Gupta et al. (1994); Tessier et al. (1996); Davis et al. (1996); Rodriguez-Iturbe and Rinaldo (1997); Pandey et al. (1998); Matsoukas et al. (2000); Montanari et al. (2000); Peters et al. (2002); Livina et al. (2003a, b).

However, the conventional methods discussed above may fail when trends are present in the system. Trends are systematic deviations from the average runoff that are caused by external processes, e.g. the construction of a water regulation device, the seasonal cycle, or a changing climate (e.g. global warming). Monotonous trends may lead to an overestimation of the Hurst exponent and thus to an underestimation of $`\gamma `$. It is even possible that uncorrelated data, under the influence of a trend, look like long-term correlated ones when using the above analysis methods. In addition, long-term correlated data cannot simply be detrended by the common technique of moving averages, since this method destroys the correlations on long time scales (above the window size used). Furthermore, it is difficult to distinguish trends from long-term correlations, because stationary long-term correlated time series exhibit persistent behaviour and a tendency to stay close to the momentary value. This causes positive or negative deviations from the average value for long periods of time that might look like a trend.

In the last years, several methods such as wavelet techniques (WT) and detrended fluctuation analysis (DFA), have been developed that are able to determine long-term correlations in the presence of trends. For details and applications of the methods to a large number of meteorological, climatological and biological records we refer to Peng et al. (1994); Taqqu et al. (1995); Bunde et al. (2000); Kantelhardt et al. (2001); Arneodo et al. (2002); Bunde et al. (2002). The methods, described in Section 2, consider fluctuations in the cumulated runoffs (often called the ”profile” or ”landscape” of the record). They differ in the way the fluctuations are determined and in the type of polynomial trend that is eliminated in each time window of size $`s`$.

In this paper, we apply these detrending methods to study the scaling of the fluctuations $`F_2(s)`$ of river flows with time $`s`$. We focus on 23 runoff records from international river stations spread around the globe and compare the results with those of 18 river stations from southern Germany. We find that above some crossover time (typically several weeks) $`F_2(s)`$ scales as $`s^H`$ with $`H`$ varying from river to river between 0.55 and 0.95 in a nonuniversal manner independent of the size of the basin. The lowest exponent $`H=0.55`$ was obtained for rivers on permafrost ground. Our finding is not consistent with the hypothesis that the scaling is universal with an exponent close to 0.75 (Hurst et al., 1965; Feder, 1988) with the same power law being applicable for all time scales from minutes till centuries.

The above detrending approaches, however, are not sufficient to fully characterize the complex dynamics of river flows, since they exclusively focus on the variance which can be regarded as the second moment $`F_2(s)`$ of the full distribution of the fluctuations. Note that the Hurst method actually focuses on the first moment $`F_1(s)`$. To further characterize a hydrological record, we extend the study to include all moments $`F_q(s)`$. A detailed description of the method, which is a multifractal generalization of the detrended fluctuation analysis (Kantelhardt et al., 2002) and equivalent to the Wavelet Transform Modulus Maxima (WTMM) method (Arneodo et al., 2002), is given in Section 3. Our approach differs from the multifractal approach introduced into hydrology by Lovejoy and Schertzer (see e.g. Schertzer and Lovejoy (1987); Lovejoy and Schertzer (1991); Lavallee et al. (1993); Pandey et al. (1998)) that was based on the concept of structure functions (Frisch and Parisi, 1985) and on the assumption of the existence of a universal cascade model. Here we perform the multifractal analysis by studying how the different moments of the fluctuations $`F_q(s)`$ scale with time $`s`$, see also Rodriguez-Iturbe and Rinaldo (1997). We find that at large time scales, $`F_q(s)`$ scales as $`s^{h(q)}`$, and a simple functional form with two parameters ($`a`$ and $`b`$), $`h(q)=(1/q)[\mathrm{ln}(a^q+b^q)]/[q\mathrm{ln}(2)]`$ describes the scaling exponent $`h(q)`$ of all moments. On small time scales, however, a stronger multifractality is observed that may be partly related to the seasonal trend. The mean position of the crossover between the two regimes is of the order of weeks and increases with $`q`$.

## 2 Correlation Analysis

Consider a record of daily water runoff values $`W_i`$ measured at a certain hydrological station. The index $`i`$ counts the days in the record, $`i=1`$,2,…,$`N`$. To eliminate the periodic seasonal trends, we concentrate on the departures $`\varphi _i=W_i\overline{W_i}`$ from the mean daily runoff $`\overline{W_i}`$. $`\overline{W_i}`$ is calculated for each calendar date $`i`$ (e.g. April 1st) by averaging over all years in the runoff series. In addition, we checked that our actual results remained unchanged when also seasonal trends in the variance have been eliminated by analysing $`\varphi _i^{}=(W_i\overline{W_i})/(\overline{W_i^2}\overline{W_i}^2)^{1/2}`$ instead of $`\varphi _i`$.

The runoff autocorrelation function $`C(s)`$ describes, how the persistence decays in time. If the $`\varphi _i`$ are uncorrelated, $`C(s)`$ is zero for all $`s`$. If correlations exist only up to a certain number of days $`s_\times `$, the correlation function will vanish above $`s_\times `$. For long-term correlations, $`C(s)`$ decays by a power-law

$$C(s)=\varphi _i\varphi _{i+s}s^\gamma ,0<\gamma <1,$$
(1)
where the average $`\mathrm{}`$ is over all pairs with the same time lag $`s`$. For large values of $`s`$, a direct calculation of $`C(s)`$ is hindered by the level of noise present in the finite hydrological records, and by nonstationarities in the data. There are several alternative methods for calculating the correlation function in the presence of long-term correlations, which we describe in the following sections.

### 2.1 Power Spectrum Analysis

If the time series is stationary, we can apply standard spectral analysis techniques and calculate the power spectrum $`E(f)`$ of the time series $`W_i`$ as a function of the frequency $`f`$. For long-term correlated data, we have $`E(f)f^\beta ,`$ where $`\beta `$ is related to the correlation exponent $`\gamma `$ by $`\beta =1\gamma `$. This relation can be derived from the Wiener-Khinchin theorem. If, instead of $`W_i`$ the integrated runoff time series $`z_n=_{i=1}^n\varphi _i`$ is Fourier transformed, the resulting power spectrum scales as $`\stackrel{~}{E}(f)f^{2\beta }`$.

### 2.2 Standard Fluctuation Analysis (FA)

In the standard fluctuation analysis, we consider the “runoff profile”

$$z_n=\underset{i=1}{\overset{n}{}}\varphi _i,n=1,2,\mathrm{},N,$$
(2)
and study how the fluctuations of the profile, in a given time window of size $`s`$, increase with $`s`$. We can consider the profile $`z_n`$ as the position of a random walker on a linear chain after $`n`$ steps. The random walker starts at the origin and performs, in the $`i`$th step, a jump of length $`\varphi _i`$ to the right, if $`\varphi _i`$ is positive, and to the left, if $`\varphi _i`$ is negative.

To find how the square-fluctuations of the profile scale with $`s`$, we first divide each record of $`N`$ elements into $`N_s=\mathrm{int}(N/s)`$ nonoverlapping segments of size $`s`$ starting from the beginning and $`N_s`$ nonoverlapping segments of size $`s`$ starting from the end of the considered runoff series. Then we determine the fluctuations in each segment $`\nu `$.

In the standard fluctuation analysis, we obtain the fluctuations just from the values of the profile at both endpoints of each segment $`\nu `$, $`F^2(\nu ,s)=[z_{\nu s}z_{(\nu 1)s}]^2`$, and average $`F^2(\nu ,s)`$ over the $`2N_s`$ subsequences to obtain the mean fluctuation $`F_2(s)`$,

$$F_2(s)\left\{\frac{1}{2N_s}\underset{\nu =1}{\overset{2N_s}{}}F^2(\nu ,s)\right\}^{1/2}.$$
(3)
By definition, $`F_2(s)`$ can be viewed as the root-mean-square displacement of the random walker on the chain, after $`s`$ steps. For uncorrelated $`\varphi _i`$ values, we obtain Fick’s diffusion law $`F_2(s)s^{1/2}`$. For the relevant case of long-term correlations, where $`C(s)`$ follows the power-law behaviour of Eq. (1), $`F_2(s)`$ increases by a power law (see, e.g., Bunde et al. (2002)),

$$F_2(s)s^H,$$
(4)
where the fluctuation exponent $`H`$ is related to the correlation exponent $`\gamma `$ and the power-spectrum exponent $`\beta `$ by

$$H=1\gamma /2=(1+\beta )/2.$$
(5)
For power-law correlations decaying faster than $`1/s`$, we have $`H=1/2`$ for large $`s`$ values, like for uncorrelated data.

We like to note that the standard fluctuation analysis is somewhat similar to the rescaled range analysis introduced by Hurst (for a review see, e. g., Feder (1988)), except that it focusses on the second moment $`F_2(s)`$ while Hurst considered the first moment $`F_1(s)`$. For monofractal data, $`H`$ is identical to the Hurst exponent.

### 2.3 The Detrended Fluctuation Analysis (DFA)

There are different orders of DFA that are distinguished by the way the trends in the data are eliminated. In lowest order (DFA1) we determine, for each segment $`\nu `$, the best linear fit of the profile, and identify the fluctuations by the variance $`F^2(\nu ,s)`$ of the profile from this straight line. This way, we eliminate the influence of possible linear trends on scales larger than the segment. Note that linear trends in the profile correspond to patch-like trends in the original record. DFA1 has been proposed originally by Peng et al. (1994) when analyzing correlations in DNA. It can be generalized straightforwardly to eliminate higher order trends (Bunde et al., 2000; Kantelhardt et al., 2001).

In second order DFA (DFA2) one calculates the variances $`F^2(\nu ,s)`$ of the profile from best quadratic fits of the profile, this way eliminating the influence of possible linear and parabolic trends on scales larger than the segment considered. In general, in $`n`$th-order DFA, we calculate the variances of the profile from the best $`n`$th-order polynomial fit, this way eliminating the influence of possible $`(n1)`$th-order trends on scales larger than the segment size.

Explicitly, we calculate the best polynomial fit $`y_\nu (i)`$ of the profile in each of the $`2N_s`$ segments $`\nu `$ and determine the variance

$$F^2(\nu ,s)\frac{1}{s}\underset{i=1}{\overset{s}{}}\left[z_{(\nu 1)s+i}y_\nu (i)\right]^2.$$
(6)
Then we employ Eq. (3) to determine the mean fluctuation $`F_2(s)`$.

Since FA and the various stages of the DFA have different detrending capabilities, a comparison of the fluctuation functions obtained by FA and DFA$`n`$ can yield insight into both long-term correlations and types of trends. This cannot be achieved by the conventional methods, like the spectral analysis.

### 2.4 Wavelet Transform (WT)

The wavelet methods we employ here are based on the determination of the mean values $`\overline{z}_\nu (s)`$ of the profile in each segment $`\nu `$ (of length $`s`$), and the calculation of the fluctuations between neighbouring segments. The different order techniques we have used in analyzing runoff fluctuations differ in the way the fluctuations between the average profiles are treated and possible nonstationarities are eliminated. The first-, second- and third-order wavelet method are described below.

(i) In the first-order wavelet method (WT1), one simply determines the fluctuations from the first derivative $`F^2(\nu ,s)=[\overline{z}_\nu (s)\overline{z}_{\nu +1}(s)]^2`$. WT1 corresponds to FA where constant trends in the profile of a hydrological station are eliminated, while linear trends are not eliminated.

(ii) In the second-order wavelet method (WT2), one determines the fluctuations from the second derivative $`F^2(\nu ,s)=[\overline{z}_\nu (s)2\overline{z}_{\nu +1}(s)+\overline{z}_{\nu +2}(s)]^2`$. So, if the profile consists of a trend term linear in $`s`$ and a fluctuating term, the trend term is eliminated. Regarding trend-elimination, WT2 corresponds to DFA1.

(iii) In the third-order wavelet method (WT3), one determines the fluctuations from the third derivative $`F^2(\nu ,s)=[\overline{z}_\nu (s)3\overline{z}_{\nu +1}(s)+3\overline{z}_{\nu +2}(s)\overline{z}_{\nu +3}(s)]^2`$. By definition, WT3 eliminates linear and parabolic trend terms in the profile. In general, in WT$`n`$ we determine the fluctuations from the $`n`$th derivative, this way eliminating trends described by $`(n1)`$st-order polynomials in the data.

Methods (i-iii) are called wavelet methods, since they can be interpreted as transforming the profile by discrete wavelets representing first-, second- and third-order cumulative derivatives of the profile. The first-order wavelets are known in the literature as Haar wavelets. One can also use different shapes of the wavelets (e.g. Gaussian wavelets with width $`s`$), which have been used by Arneodo et al. (2002) to study, for example, long-range correlations in DNA. Since the various stages of the wavelet methods WT1, WT2, WT3, etc. have different detrending capabilities, a comparison of their fluctuation functions can yield insight into both long-term correlations and types of trends.

At the end of this section, before describing the results of the FA, DFA, and WT analysis, we note that for very large $`s`$ values, $`s>N/4`$ for DFA and $`s>N/10`$ for FA and WT, the fluctuation function becomes inaccurate due to statistical errors. The difference in the statistics is due to the fact that the number of independent segments of length $`s`$ is larger in DFA than in WT, and the fluctuations in FA are larger than in DFA. Hence, in the analysis we will concentrate on $`s`$ values lower than $`s_{\mathrm{max}}=N/4`$ for DFA and $`s_{\mathrm{max}}=N/10`$ for FA and WT. When determining the scaling exponents $`H`$ using Eq. (4), we manually chose an appropriate (shorter) fitting range of typically two orders of magnitude.

### 2.5 Results

In our study we analyzed 41 runoff records, 18 of them are from the southern part of Germany, and the rest is from North and South America, Africa, Australia, Asia and Europe (see Table 1). We begin the analysis with the runoff record for the river Weser in the northern part of Germany, which has the longest record (171 years) in this study. Figure 1(a) shows the fluctuation functions $`F_2(s)`$ obtained from FA and DFA1-DFA5. In the log-log plot, the curves are approximately straight lines for $`s`$ above 30 days, with a slope $`H0.75`$. This result for the Weser suggests that there exists long-term persistence expressed by the power-law decay of the correlation function, with an exponent $`\gamma 0.5`$ \[see Eq. (4)\].

To show that the slope $`H0.75`$ is due to long-term correlations and not due to a broad probability distribution (Joseph- versus Noah-Phenomenon, see Mandelbrot and Wallis (1968)), we have eliminated the correlations by randomly shuffling the $`\varphi _i`$. This shuffling has no effect on the probability distribution function of $`\varphi _i`$. Figure 1(b) shows $`F_2(s)`$ for the shuffled data. We obtain $`H=1/2`$, showing that the exponent $`H0.75`$ is due to long-term correlations.

To show that the slope $`H0.75`$ is not an artefact of the seasonal dependence of the variance and skew, we also considered records where $`\varphi _i`$ was divided by the variance of each calendary day and applied further detrending techniques that take into account the skew (Livina et al., 2003b). In both cases, we found no change in the scaling behaviour for large times (see also Sect. 3.4). This can be understood easily, since kinds of seasonal trends cannot effect the fluctuation behaviour on time scales well above one year. It is likely, however, that the seasonal dependencies of the variance and possibly also of the skew contribute to the behaviour at small times, where the slope $`H`$ is much larger than 0.75 in most cases (see also Sect. 3.4, where a seasonal trend in the variance is indeed used for modelling the crossover).

Figure 2 shows the fluctuation functions $`F_2(s)`$ of 4 more rivers, from Africa, South America, Australia, and Europe. The panels on the left-hand side show the FA and DFA1-4 curves, while the panels on the right-hand side show the results from the analogous wavelet analysis WT1-WT5. Most curves show a crossover at small time scales; a similar crossover has been reported by Tessier et al. (1996) for small French rivers without artificial dams or reservoirs. Above the crossover time, the fluctuation functions (from DFA1-4 and WT2-5) show power-law behaviour, with exponents $`H0.95`$ for the Zaire, $`H0.73`$ for the Orinoco, $`H0.60`$ for the Mary river, and $`H0.55`$ for the Gaula river.

Accordingly, there is no universal scaling behaviour since the long-term exponents vary strongly from river to river and reflect the fact that there exist different mechanisms for floods where each may induce different scaling. This is in contrast to climate data, where universal long-term persistence of temperature records at land stations was observed (Koscielny-Bunde et al., 1998; Talkner and Weber, 2000; Weber and Talkner, 2001; Eichner et al., 2003).

The Mary river in Australia is rather dry in the summer and the Gaula river in Norway is frozen in the winter. For the Mary river, the long-term exponent $`H0.60`$ is well below the average value. For the Gaula river, the long-term correlations are not pronounced ($`H=0.55`$) and even hard to distinguish from the uncorrelated case $`H=0.5`$. We obtained similar results for the other two ”frozen” rivers (Tana from Norway and Dvina from Russia) that we analysed. For interpreting this distinguished behaviour of the frozen rivers we like to note that on permafrost ground the lateral inflow (and hence the indirect contribution of the water storage in the catchment basin) contributes to the runoffs in a different way than on normal ground, see also Gupta and Dawdy . Our results (based on 3 rivers only) suggest that the contribution of snow melting leads to less correlated runoffs than the contribution of rainfall, but more comprehensive studies will be needed to confirm this interesting result.

Figure 3(a) and Table 1 summarize our results for $`H`$. One can see clearly that the exponents $`H`$ do not depend systematically on the basin area $`A`$. This is in line with the conclusions of Gupta et al. (1994) for the flood peaks, where a systematic dependence on $`A`$ could also not be found.

There is also no pronounced regional dependence: the rivers within a localized area (such as South Germany) tend to have nearly the same range of exponents as the international rivers. The three ”frozen” rivers in our study have the lowest values of $`H`$. As can be seen in the figure, the exponents spread from 0.55 to 0.95. Since the correlation exponent $`\gamma `$ is related to $`H`$ by $`\gamma =22H`$, the exponent $`\gamma `$ spreads from almost 0 to almost 1, covering the whole range from very weak to very strong correlations.

## 3 Multifractal Analysis

### 3.1 Method

For a further characterization of hydrological records it is meaningful to extend Eq. (3) by considering the more general fluctuation functions (Barabasi and Vicsek (1991), see also Davis et al. (1996)),

$$F_q(s)=\left\{\frac{1}{2N_s}\underset{\nu =1}{\overset{2N_s}{}}|z_{\nu s}z_{(\nu 1)s}|^q\right\}^{1/q}.$$
(7)
where the variable $`q`$ can take any real value except zero. For $`q=2`$, the standard fluctuation analysis is retrieved. The question is, how the fluctuation functions depend on $`q`$ and how this dependence is related to multifractal features of the record.

In general, the multifractal approach is introduced by the partition function

$$Z_q(s)\underset{\nu =1}{\overset{N_s}{}}|z_{\nu s}z_{(\nu 1)s}|^qs^{\tau (q)},$$
(8)
where $`\tau (q)`$ is the Renyi scaling exponent. A record is called ’monofractal’, when $`\tau (q)`$ depends linearly on $`q`$; otherwise it is called multifractal.

It is easy to verify that $`Z_q(s)`$ is related to $`F_q(s)`$ by

$$F_q(s)=\left\{\frac{1}{N_s}Z_q(s)\right\}^{1/q}.$$
(9)
Accordingly, Eq. (8) implies

$$F_q(s)s^{h(q)},$$
(10)
where

$$h(q)=[\tau (q)+1]/q.$$
(11)
Thus, $`h(q)`$ defined in Eq. (10) is directly related to the classical multifractal scaling exponents $`\tau (q)`$.

In general, the exponent $`h(q)`$ may depend on $`q`$. Since for stationary records, $`h(1)`$ is identical to the well-known Hurst exponent (see e. g. Feder (1988)), we will call the function $`h(q)`$ the generalized Hurst exponent. For monofractal self-affine time series, $`h(q)`$ is independent of $`q`$, since the scaling behaviour of the variances $`F^2(\nu ,s)`$ is identical for all segments $`\nu `$, and the averaging procedure in Eq. (7) will give just this identical scaling behaviour for all values of $`q`$. If small and large fluctuations scale differently, there will be a significant dependence of $`h(q)`$ on $`q`$: If we consider positive values of $`q`$, the segments $`\nu `$ with large variance $`F^2(\nu ,s)`$ (i. e. large deviations from the corresponding fit) will dominate the average $`F_q(s)`$. Thus, for positive values of $`q`$, $`h(q)`$ describes the scaling behaviour of the segments with large fluctuations. Usually the large fluctuations are characterized by a smaller scaling exponent $`h(q)`$ for multifractal series. On the contrary, for negative values of $`q`$, the segments $`\nu `$ with small variance $`F^2(\nu ,s)`$ will dominate the average $`F_q(s)`$. Hence, for negative values of $`q`$, $`h(q)`$ describes the scaling behaviour of the segments with small fluctuations, which are usually characterized by a larger scaling exponent.

In the hydrological literature (Rodriguez-Iturbe and Rinaldo, 1997; Lavallee et al., 1993) one often considers the generalized mass variogram $`C_q(\lambda )`$ (see also Davis et al. (1996)),

$$C_q(\lambda )|z_{i+\lambda }z_i|^q\lambda ^{K(q)}.$$
(12)
Comparing Eqs. (7), (10), and (12) one can verify easily that $`K(q)`$ and $`h(q)`$ are related by

$$h(q)=K(q)/q.$$
(13)
Another way to characterize a multifractal series is the singularity spectrum $`f(\alpha )`$, that is related to $`\tau (q)`$ via a Legendre transform (e.g. Feder (1988); Rodriguez-Iturbe and Rinaldo (1997)),

$$\alpha =\frac{d\tau (q)}{dq}\mathrm{and}f(\alpha )=q\alpha \tau (q).$$
(14)
Here, $`\alpha `$ is the singularity strength or Hölder exponent, while $`f(\alpha )`$ denotes the dimension of the subset of the series that is characterized by $`\alpha `$. Using Eq. (11), we can directly relate $`\alpha `$ and $`f(\alpha )`$ to $`h(q)`$,

$$\alpha =h(q)+q\frac{dh(q)}{dq}\mathrm{and}f(\alpha )=q[\alpha h(q)]+1.$$
(15)
The strength of the multifractality of a time series can be characterized by the difference between the maximum and minimum values of $`\alpha `$, $`\alpha _{\mathrm{max}}\alpha _{\mathrm{min}}`$. When $`q\frac{dh(q)}{dq}`$ approaches zero for $`q`$ approaching $`\pm \mathrm{}`$, then $`\mathrm{\Delta }\alpha =\alpha _{\mathrm{max}}\alpha _{\mathrm{min}}`$ is simply given by $`\mathrm{\Delta }\alpha =h(\mathrm{})h(\mathrm{})`$.

The multifractal analysis described above is a straightforward generalization of the fluctuation analysis and therefore has the same problems: (i) monotonous trends in the record may lead to spurious results for the fluctuation exponent $`h(q)`$ which in turn leads to spurious results for the correlation exponent $`\gamma `$, and (ii) nonstationary behaviour characterized by exponents $`h(q)>1`$ cannot be detected by the simple method since the method cannot distinguish between exponents $`>1`$, and always will yield $`F_2(s)s`$ in this case (see above). To overcome these drawbacks the multifractal detrended fluctuation analysis (MF-DFA) has been introduced recently (Kantelhardt et al. (2002), see also Koscielny-Bunde et al. (1998); Weber and Talkner (2001)). According to Kantelhardt et al. (2002, 2003), the method is as accurate as the wavelet methods. Thus, we have used MF-DFA for the multifractal analysis here. In the MF-DFA, one starts with the DFA-fluctuations $`F_\nu ^2(s)`$ as obtained in Eq. (6). Then, we define in close analogy to Eqs. (3) and (7) the generalized fluctuation function,

$$F_q(s)\left\{\frac{1}{2N_s}\underset{\nu =1}{\overset{2N_s}{}}\left[F^2(\nu ,s)\right]^{q/2}\right\}^{1/q}.$$
(16)
Again, we can distinguish MF-DFA1, MF-DFA2, etc., according to the order of the polynomial fits involved.

### 3.2 Multifractal Scaling Plots

We have performed a large scale multifractal analysis on all 41 rivers. We found that MF-DFA2-5 yield similar results for the fluctuation function $`F_q(s)`$. We have also cross checked the results using the Wavelet Transform Modulus Maxima (WTMM) method (Arneodo et al., 2002), and always find agreement within the error bars (Kantelhardt et al., 2003). Therefore, we present here only the results of MF-DFA4. Figure 4(a,b) shows two representative examples for the fluctuation functions $`F_q(s)`$, for (a) the Weser river and (b) the Danube river. The standard fluctuation function $`F_2(s)`$ is plotted in full symbols. The crossover in $`F_2(s)`$ that was discussed in Sect. 2.5 can be also seen in the other moments. The position of the crossover increases monotonously with decreasing $`q`$ and the crossover becomes more pronounced. We are only interested in the asymptotic behaviour of $`F_q(s)`$ at large times $`s`$. One can see clearly that above the crossover, the $`F_q(s)`$ functions are straight lines in the double logarithmic plot, and the slopes increase slightly when going from high positive moments towards high negative moments (from the bottom to the top). For the Weser, for example, the slope changes from 0.65 for $`q=10`$ to 0.9 for $`q=10`$ (see also Fig. 6(b)). The monotonous increase of the slopes, $`h(q)`$, is the signature of multifractality.

When the data are shuffled (see Figs. 4(c,d)), all functions $`F_q(s)`$ increase asymptomatically as $`F_q(s)s^{1/2}`$. This indicates that the multifractality vanishes under shuffling. Accordingly the observed multifractality originates in the long-term correlations of the record and is not caused by singularities in the distribution of the daily runoffs (see also Mandelbrot and Wallis (1968)). A reshuffling-resistent multifractality would indicate a ’statistical’ type of non-linearity (Sivapalan et al., 2002). We obtain similar patterns for all rivers. Figure 5 shows four more examples; Figs. 5(a,b) are for two rivers (Amper and Wertach) from southern Germany, while Figs. 5(c,d) are for Niger and Susquehama (Koulikoro, Mali and Harrisburg, USA).

From the asymptotic slopes of the curves in Figs. 4(a,b) and 5(a-d), we obtain the generalized Hurst-exponents $`h(q)`$, which are plotted in Fig. 6 (circles). One can see that in the whole $`q`$-range the exponents can be fitted well by the formula

$$h(q)=\frac{1}{q}\frac{\mathrm{ln}[a^q+b^q]}{q\mathrm{ln}2},$$
(17)
or

$$K(q)=1+\tau (q)=1\frac{\mathrm{ln}[a^q+b^q]}{\mathrm{ln}2}.$$
(18)
The formula can be derived from a modification of the multiplicative cascade model that we describe in Sect. 3.3. Here we use the formula only to show that the infinite number of exponents $`h(q)`$ can be described by only two independent parameters, $`a`$ and $`b`$. These two parameters can then be regarded as multifractal finger print for a considered river. This is particularly important when checking models for river flows. Again, we like to emphasize, that these parameters have been obtained from the asymptotic part of the generalized fluctuation function, and are therefore not affected by seasonal dependencies.

We have fitted the $`h(q)`$ spectra in the range $`10q10`$ for all 41 runoff series by Eq. (17). Representative examples are shown in Fig. 6. The continuous lines in Fig. 6 are obtained by best fits of $`h(q)`$ (obtained from Figs. 4 and 5 as described above) by Eq. (17). The respective parameters $`a`$ and $`b`$ are listed inside the panels of each figure. Our results for the 41 rivers are shown in Table 1. It is remarkable that in each single case, the $`q`$ dependence of $`h(q)`$ for positive and negative values of $`q`$ can be characterized well by the two parameters, and all fits remain within the error bars of the $`h(q)`$ values.

From $`h(q)`$ we obtain $`\tau (q)`$ (Eq. (11)) and the singularity spectrum $`f(\alpha )`$ (Eq. (15)). Figure 7 shows two typical examples for the Danube and the Niger. The width of $`f(\alpha )`$ taken at $`f=0`$ characterizes the strength of the multifractality. Since both widths are very different, the strength of the multifractality of river runoffs appears to be not universal.

In order to characterize and to compare the strength of the multifractality for several time series we use as a parameter the width of the singularity spectrum $`f(\alpha )`$ \[see Eqs. (14) and (15)\] at $`f=0`$, which corresponds to the difference of the maximum and the minimum value of $`\alpha `$. In the multiplicative cascade model, this parameter is given by

$$\mathrm{\Delta }\alpha =\frac{\mathrm{ln}a\mathrm{ln}b}{\mathrm{ln}2}.$$
(19)
The distribution of the $`\mathrm{\Delta }\alpha `$ values we obtained from Eq. (19) is shown in Fig. 3(b), where we plot $`\mathrm{\Delta }\alpha `$ versus the basin area. The figure shows that there are rivers with quite strong multifractal fluctuations, i.e. large $`\mathrm{\Delta }\alpha `$, and one with almost vanishing multifractality, i.e. $`\mathrm{\Delta }\alpha 0`$. Two observation can be made from the figure: (1) There is no pronounced difference between the width of the distribution of the multifractality strength for the runoffs within the local area of southern Germany (open symbols in Fig. 3(b)) and for the international runoff records from all rivers around the globe (full symbols). In fact, without rivers Zaire and Grand River the widths would be the same. (2) There is a tendency towards smaller multifractality strengths at larger basin areas. This means that the river flows become less nonlinear with increasing basin area. We consider this as possible indication of river regulations that are more pronounced for large river basins.

Our results for $`K(q)=1+\tau (q)`$ (see Eq. (18) and Table 1) may be compared with the functional form

$$K(q)=(H^{}+1)q\frac{C_1}{\alpha ^{}1}(q^\alpha ^{}q)q0$$
(20)
with the three parameters $`H^{}`$, $`C_1`$, and $`\alpha ^{}`$, that have been used by Lovejoy, Schertzer, and coworkers (Schertzer and Lovejoy, 1987; Lovejoy and Schertzer, 1991; Lavallee et al., 1993; Tessier et al., 1996; Pandey et al., 1998) successfully to describe the multifractal behaviour of rainfall and runoff records. The definition of $`K(q)`$ we used in this paper is taken from Rodriguez-Iturbe and Rinaldo (1997) and differs slightly from their definition. We like to note that Eq. (18) for $`K(q)`$ is not only valid for positive $`q`$ values, but also for negative $`q`$ values. This feature allows us to determine numerically the full singularity spectrum $`f(\alpha )`$. In the analysis we focused on long time scales, excluding the crossover regime, and used detrending methods. We consider it as particularly interesting that only two parameters $`a`$ and $`b`$ or, equivalently, $`H`$ and $`\mathrm{\Delta }\alpha `$, are sufficient to describe $`\tau (q)`$ and $`K(q)`$ for positive as well as negative $`q`$ values. This strongly supports the idea of ’universal’ multifractal behaviour of river runoffs as suggested (in different context) by Lovejoy and Schertzer.

It is interesting to note that the generalized fluctuation functions we studied do not show any kind of multifractal phase transition at some critical value $`q_D`$ in the $`q`$-regime ($`10q10`$) we analysed. Instead, our analysis shows a crossover at a specific time scale $`s_\times `$ (typically weeks) that weakly increases with decreasing moment $`q`$. In this paper, we concentrated on the large-time regime ($`ss_\times `$), where we obtained coherent multifractal behaviour and did not see any indication for a multifractal phase transition. But this does not exclude the possibility that at small scales a breakdown of mulifractality at a critical $`q`$-value may occur, as has been emphasized by Tessier et al. (1996) and Pandey et al. (1998).

### 3.3 Extended Multiplicative Cascade Model

In the following, we like to motivate the 2-parameter formula Eq. (17) and show how it can be obtained from the well known multifractal cascade model (Feder, 1988; Barabasi and Vicsek, 1991; Kantelhardt et al., 2002). In the model, a record $`\varphi _k`$ of length $`N=2^{n_{\mathrm{max}}}`$ is constructed recursively as follows: In generation $`n=0`$, the record elements are constant, i.e. $`\varphi _k=1`$ for all $`k=1,\mathrm{},N`$. In the first step of the cascade (generation $`n=1`$), the first half of the series is multiplied by a factor $`a`$ and the second half of the series is multiplied by a factor $`b`$. This yields $`\varphi _k=a`$ for $`k=1,\mathrm{},N/2`$ and $`\varphi _k=b`$ for $`k=N/2+1,\mathrm{},N`$. The parameters $`a`$ and $`b`$ are between zero and one, $`0<a<b<1`$. Note that we do not restrict the model to $`b=1a`$ as is often done in the literature (Feder, 1988). In the second step (generation $`n=2`$), we apply the process of step 1 to the two subseries, yielding $`\varphi _k=a^2`$ for $`k=1,\mathrm{},N/4`$, $`\varphi _k=ab`$ for $`k=N/4+1,\mathrm{},N/2`$, $`\varphi _k=ba=ab`$ for $`k=N/2+1,\mathrm{},3N/4`$, and $`\varphi _k=b^2`$ for $`k=3N/4+1,\mathrm{},N`$. In general, in step $`n+1`$, each subseries of step $`n`$ is divided into two subseries of equal length, and the first half of the $`\varphi _k`$ is multiplied by $`a`$ while the second half is multiplied by $`b`$. For example, in generation $`n=3`$ the values in the eight subseries are $`a^3,a^2b,a^2b,ab^2,a^2b,ab^2,ab^2,b^3`$. After $`n_{\mathrm{max}}`$ steps, the final generation has been reached, where all subseries have length 1 and no more splitting is possible. We note that the final record can be written as $`\varphi _k=a^{n_{\mathrm{max}}n(k1)}b^{n(k1)}`$, where $`n(k)`$ is the number of digits 1 in the binary representation of the index $`k`$, e. g. $`n(13)=3`$, since 13 corresponds to binary 1101.

For this multiplicative cascade model, the formula for $`\tau (q)`$ has been derived earlier (Feder, 1988; Barabasi and Vicsek, 1991; Kantelhardt et al., 2002). The result is $`\tau (q)=[\mathrm{ln}(a^q+b^q)+q\mathrm{ln}(a+b)]/\mathrm{ln}2`$ or

$$h(q)=\frac{1}{q}\frac{\mathrm{ln}(a^q+b^q)}{q\mathrm{ln}2}+\frac{\mathrm{ln}(a+b)}{\mathrm{ln}2}.$$
(21)
It is easy to see that $`h(1)=1`$ for all values of $`a`$ and $`b`$. Thus, in this form the model is limited to cases where $`h(1)`$, which is the exponent Hurst defined originally in the $`R/S`$ method, is equal to one. In order to generalize this multifractal cascade process such that any value of $`h(1)`$ is possible, we have subtracted the offset $`\mathrm{\Delta }h=\mathrm{ln}(a+b)/\mathrm{ln}(2)`$ from $`h(q)`$. The constant offset $`\mathrm{\Delta }h`$ corresponds to additional long-term correlations incorporated in the multiplicative cascade model. For generating records without this offset, we rescale the power spectrum. First, we fast-Fourier transform (FFT) the simple multiplicative cascade data into the frequency domain. Then, we multiply all Fourier coefficients by $`f^{\mathrm{\Delta }h}`$, where $`f`$ is the frequency. This way, the slope $`\beta `$ of the power spectra $`E(f)f^\beta `$ (the squares of the Fourier coefficients) is decreased from $`\beta =2h(2)1=[2\mathrm{ln}(a+b)\mathrm{ln}(a^2+b^2)]/\mathrm{ln}2`$ into $`\beta ^{}=2[h(2)\mathrm{\Delta }h]1=\mathrm{ln}(a^2+b^2)/\mathrm{ln}2`$, which is consistent with Eq. (17). Finally, backward FFT is employed to transform the signal back into the time domain. A similar Fourier filtering technique has been used by Tessier et al. (1996) when generating surrogate runoff data.

### 3.4 Comparison with Model Data

In order to see how well the extended multiplicative cascade model fits to the real data (for a given river), we generate the model data as follows: (i) we determine $`a`$ and $`b`$ for the given river (by best fit of Eq. (17)), (ii) we generate the simple multiplicative cascade model with the obtained $`a`$ and $`b`$ values, and (iii) we implement the proper long-term correlations as described above.

Figure 8(a) shows the DFA analysis of the model data with parameters $`a`$ and $`b`$ determined for the river Weser. By comparing with Fig. 4(a) we see that the extended model gives the correct scaling of the fluctuation functions $`F_q(s)`$ on time scales above the crossover. By comparing Fig. 8(b) with Fig. 4(c) we see that the shuffled model series becomes uncorrelated without multifractality similar to the shuffled data. Below the crossover, however, the model does not yield the observed $`F_q(s)`$ in the original data. In the following we show that in order to obtain the proper behaviour below the crossover, either seasonal trends that cannot be completely eliminated from the data or a different type of multifractality below the crossover, represented by different values of $`a`$ and $`b`$, have to be introduced.

To show the effect of seasonal trends, we have multiplied the elements $`\varphi _i`$ of the extended cascade model by $`0.1+\mathrm{sin}^2(\pi i/365)`$, this way generating a seasonal trend of period 365 in the variance. Figure 8(c) shows the DFA4 result for the generalized fluctuation functions, which now better resembles the real data than Fig. 8(a). Finally, in Fig. 8(d) we show the effect of a different multifractality below the crossover, where different parameters $`a`$ and $`b`$ characterize this regime. The results also show better agreement with the real data. When comparing Figs. 8(a,c,d) with Figs. 4 and 5, it seems that the Danube, Amper and Wertach fit better to Fig. 8(d), i.e. suggesting different multifractality for short and large time scales, while the Weser, Susquehanna and Niger fit better to Fig. 8(c) where seasonal trends in the variance (and possibly in the skew) are responsible for the behaviour below the crossover.

## 4 Conclusion

In this study, we analyzed the scaling behaviour of daily runoff time series of 18 representative rivers in southern Germany and 23 international rivers using both Detrended Fluctuation Analysis and wavelet techniques. In all cases we found that the fluctuations exhibit self-affine scaling behaviour and long-term persistence on time scales ranging from weeks to decades. The fluctuation exponent $`H`$ varies from river to river in a wide range between 0.55 and 0.95, showing non-universal scaling behaviour.

We also studied the multifractal properties of the runoff time series using a multifractal generalization of the DFA method that was crosschecked with the WTMM technique. We found that the multifractal spectra of all 41 records can be described by a ’universal’ function $`\tau (q)=\mathrm{ln}(a^q+b^q)/\mathrm{ln}2`$, which can be obtained from a generalization of the multiplicative cascade model and has solely two parameters $`a`$ and $`b`$ or, equivalently, the fluctuation exponent $`H=\frac{1}{2}\mathrm{ln}(a^2+b^2)/\mathrm{ln}4`$ and the width $`\mathrm{\Delta }\alpha =\mathrm{ln}\frac{a}{b}/\mathrm{ln}2`$ of the singularity spectrum. Since our function for $`\tau (q)`$ applies also for negative $`q`$ values, we could derive the singularity spectra $`f(\alpha )`$ from the fits. We have calculated and listed the values of $`H`$, $`a`$, $`b`$, and $`\mathrm{\Delta }\alpha `$ for all records considered. There are no significant differences between their distributions for rivers in South Germany and for international rivers. We also found that there is no significant dependence of these parameters on the size of the basin area, but there is a slight decrease of the multifractal width $`\mathrm{\Delta }\alpha `$ with increasing basin area. We suggest that the values of $`H`$ and $`\mathrm{\Delta }\alpha `$ can be regarded as ’fingerprints’ for each station or river, which can serve as an efficient non-trivial test bed for the state-of-the-art precipitation-runoff models.

Apart from the practical use of Eq. (17) with the parameters $`a`$ and $`b`$ that was derived by extending the multiplicative cascade model and that can serve as finger prints for the river flows, we presently are lacking a physical model for this behaviour. It will be interesting to see, if physically based models, e.g. the random tree-structure model presented in Gupta et al. (1996), can be related to the multiplicative cascade model presented here. If so, this would give a physical explanation for how the multiplicative cascade model is able to simulate river flows.

We have also investigated the origin of the multifractal scaling behaviour by comparison with the corresponding shuffled data. We found that the multifractality is removed by shuffling that destroys the time correlations in the series while the distribution of the runoff values is not altered. After shuffling, we obtain $`h(q)1/2`$ for all values of $`q`$, indicating monofractal behaviour. Hence, our results suggest that the multifractality is not due to the existence of a broad, asymmetric (singular) probability density distribution (Anderson and Meerschaert, 2002), but due to a specific dynamical arrangement of the values in the time series, i.e. a self-similar ’clustering’ of time patterns of values on different time scales. We believe that our results will be useful also to improve the understanding of extreme values (singularities) in the presence of multifractal long-term correlations and trends.

Finally, for an optimal water management in a given basin, it is essential to know, whether an observed long-term fluctuation in discharge data is due to systematic variations (trends) or the results of long-term correlation. Our approach is also a step forward in this directions.

Acknowledgments: We would like to thank Daniel Schertzer and Diego Rybski for valuable discussions. This work was supported by the BMBF, the DAAD, and the DFG. We also would like to thank the Water Management Authorities of Bavaria and Baden-Württemberg (Germany), and the Global Runoff Center (GRDC) in Koblenz (Germany) for providing the observational data.


format: text

text:
# The Librating Companions in HD 37124, HD 12661, HD 82943, 47 Uma and GJ 876: Alignment or Antialignment?

## 1 Introduction

The number of known extrasolar planetary systems are quickly growing (Butler et al. 2003,Paper I) in recent years. The analysis and orbital fitting of the curves of radial velocity can reveal the fact that the main-sequence stars may host one or more planets. Recently, Fischer et al.(2003)(Paper II) reported the properties of ten multiple-planet systems (see Table 8 in Paper II) and pointed out that most of the systems with the ratio of orbital periods less than 5:1 are characterized by mean motion resonance (MMR), such as GJ 876(2:1) and HD 82943(2:1). The other important feature for these systems is the apsidal lock of the orbiting companions, which indicates that the relative apsidal longitudes of two orbits librate about a constant, such that two planets have common time-averaged rate of apsidal precession. A number of researchers studied this dynamical mechanism in the planetary systems. The work of Ji et al.(2003a) confirmed that the inner two companions of 55 Cnc (Marcy et al. 2002) are in a 3:1 resonance and still experience the apsidal phase-locking. In fact, one may be familiar to the presence of the apsidal lock in our solar system, such as the Galilean satellites of Io-Europa where the difference of their apsidal longitudes $`\mathrm{\Delta }\varpi `$ are anti-aligned (Lee & Peale 2002). In the case of the pulsar PSR 1257+12, Wolszczan & Frail (1992) pointed out that there exists almost a $`\mathrm{\Delta }\varpi `$-libration about $`180^0`$ for the two planets in a 3:2 MMR. Then a question arises: are all discovered multiple-planet systems in apsidal lock? If so, how frequently does this mechanism take place for these systems?

According to Laplace-Lagrange (L-L) linear secular theory (Murray & Dermott 1999, MD99), if $`\mathrm{\Delta }\varpi `$ is a libration point, then $`360^{}\mathrm{\Delta }\varpi `$ could also be the other equilibrium due to the symmetry of planetary configurations (also see Beauge et al. 2003). So the second question is that if two planets undergo alignment or antialignment, could the mirror geometry occur? In this Letter, our goal is to examine the apsidal motion for HD 37124, HD 12661, HD 82943, 47 Uma, GJ 876.

## 2 Results

In the present work, we performed a preliminary study of various systems by using N-body codes(Ji, Li & Liu 2002). We take the masses of the parent star and companions(with sini = 1) from Table 1 for each system. The time step is adopted to be 1% to 2.5% of the orbital period of the innermost planet. Here we simply account for coplanar configurations for these systems. We remain the values of the semi-major axes unchanged for all orbits. The eccentricities and arguments of periapse are generated in the orbital parameters web based on the best-fit orbital solutions given the nominal observation errors (see Table 1)111Here we study the possible motion near the best-fit orbital solutions. This may reveal some important dynamical features for a certain investigated system.. And the remaining arguments of nodal longitudes and mean anomalies are randomly made between $`0^{}`$ and $`360^{}`$. Finally,five two-planet systems were individually produced 100 pairs of coplanar orbits for integration from 1 Myr to 10 Myr.

Throughout the paper, $`\lambda `$ and $`\varpi `$ denote the mean and periastron longitudes of a planet respectively. The subscripts $`b`$ and $`c`$ separately represent the inner and outer companions222For the sake of convenience, we exchange the order of the subscripts for two planets of HD 82943 such that the inner planet can be always described by the symbol $`b`$.. Consequently, the relative apsidal longitude is $`\mathrm{\Delta }\varpi `$ = $`\varpi _b\varpi _c`$.

### 2.1 HD 37124

Paper I revealed that two companions orbit HD 37124, but indicated that with only 30 observations the eccentricity of the outer planet was poorly determined,with any eccentricity between 0.3 and 0.8 fitting the observations within measurement errors (Butler 2003, private communication). In our simulations, we let $`e_c`$ initially be in the interval \[0.20, 0.60\] in an attempt to find likely dynamical constraints for $`e_c`$. However, all orbits can survive for $`10^6`$ yr, indicating the system is extremely stable for any initial $`e_c`$ no greater than 0.60. Moreover, we noticed that $`e_c`$ ranges from 0 to 0.62 in secular evolution for all cases, coinciding with the outcomes by Paper I and Kiseleva-Eggleton et al.(2002). They pointed out that the system self-destructed for the critical initial value of $`e_c`$ above 0.65. In addition, there are two typical secular variations of both eccentricities for stable orbits:(a)Case of large oscillations:$`e_b`$ is up to 0.85 and $`e_c`$ is less than 0.62.(b)Case of small modulations, both $`e_b`$ and $`e_c`$ are less than 0.40 (see Figure 1). The dynamical origin of high eccentricity for Case (a) may result from Kozai mechanism (Kozai 1962;Wu & Murray 2003) that implies the coupled relationship of the maximum eccentricity and minimum inclination.

Furthermore, it is noteworthy that 43%333The percentage of the aligned and antialigned cases may be affected by the uncertainty of $`e_c`$, but imply that two planets of this system likely tend to be aligned. and 1% of the stable orbits are separately in the periapse alignment and antialignment. To the best of our knowledge, such results are a new finding for this system. In Figure 1, the upper and low panels separately exhibit the variations of eccentricities and relative apsidal longitudes. The results of L-L theory (dash lines) confirm the numerical results (thick lines) in the amplitudes and periods of libration for both eccentricities and $`\mathrm{\Delta }\varpi `$.

### 2.2 HD 12661

The orbital fits revealed that the HD 12661 system is close to 11:2 MMR (Paper II;Gozdziewski 2003a) or 6:1 MMR (Gozdziewski & Maciejewski 2003b). With the best-fit data from paper II, our simulations exhibit all the orbits are stable for 1 Myr. We observe that 46% and 49% of the orbits are in $`\mathrm{\Delta }\varpi `$-libration about $`0^{}`$ and $`180^{}`$ respectively. Moreover, Gozdziewski(2003a) indicated that the possibility of $`\mathrm{\Delta }\varpi `$-circulation is much less likely than $`\mathrm{\Delta }\varpi `$-libration about $`0^{}`$ or $`180^{}`$ in the numerical investigations. Here the two types of the apsidal librations almost share an equal opportunity to be acted as one of the dynamical mechanisms to remain the stability of HD 12661, because the two planets with moderate eccentricities are always far from each other with isolated orbits and protected from close encounters. Still, the fact of $`\mathrm{\Delta }\varpi `$-libration indicates that one planet travels through the periapse or apoapse when the other moves to the opposite side at the same time.

The critical arguments for 11:2 MMR $`\varphi _k`$ can be written:

$$\varphi _k=2\lambda _b11\lambda _c+k\varpi _b+(9k)\varpi _c,k=0\mathrm{}9.$$
(1)
Here we define $`\theta _1=\varphi _3=2\lambda _b11\lambda _c+3\varpi _b+6\varpi _c`$, $`\theta _2=\varphi _2=2\lambda _b11\lambda _c+2\varpi _b+7\varpi _c`$ and $`\theta _3=\theta _1\theta _2=\mathrm{\Delta }\varpi `$. Most systems in a 11:2 MMR are only in behavior of $`\theta _1`$-libration. A particular case is exhibited in Figure 2, which not only ten resonant arguments all librate but $`\theta _3`$ librates about $`180^{}`$ for 1 Myr. Using octupole-level secular theory, Lee & Peale (2003) found that if the ratio of the maximum orbital angular momenta is approximately equal to a critical value,then the $`\mathrm{\Delta }\varpi `$-libration about $`0^{}`$ or $`180^{}`$ could happen with large amplitude variations of the eccentricities. Additionally, we further found that the status of $`\mathrm{\Delta }\varpi `$-libration are relevant to the ratio of the initial eccentricities and the initial relative apsidal longitudes of two planets (see §3).

### 2.3 47 Uma

The two companions of 47 Uma is a close analog to the duo of Jupiter and Saturn in our solar system in the ratios of the planet masses and orbital periods. Laughlin et al.(2002) found the system is in $`\mathrm{\Delta }\varpi `$-libration about $`0^{}`$ with low eccentricities of the companions. Further, Gozdziewski (2002) used MEGNO technique to investigate in wide ranges of the orbital parameters for the system. Using the updated parameters by Paper II, we still found this system can experience not only a $`\mathrm{\Delta }\varpi `$-libration about $`0^{}`$, but about $`180^{}`$. We note that all the experiments can last for 10 Myr. The fractions in alignment, antialignment and circulation are respectively 20%, 10% and 70%. As well known, Jupiter and Saturn are near 5:2 MMR, but the $`\mathrm{\Delta }\varpi `$ circulates at present-day. As for 47 Uma, we did find one of the critical arguments $`3\lambda _b8\lambda _c+2\varpi _b+3\varpi _c`$ temporarily librates for hundreds of years but $`\mathrm{\Delta }\varpi `$ simultaneously circulates in the secular evolution. However, we argue that 47 Uma may be captured into 8:3 MMR in the past during the process of migration due to the planet-disk interaction.

### 2.4 HD 82943 and GJ 876

The celebrated nature for HD 82943 and GJ 876 is that the ratio of the orbital periods of two companions is respectively close to 2:1. The critical arguments for 2:1 resonance are $`\theta _1=\lambda _b2\lambda _c+\varpi _b`$, $`\theta _2=\lambda _b2\lambda _c+\varpi _c`$, and the relative apsidal longitude is $`\theta _3=\theta _1\theta _2`$. In the simulations, we found three types of stable orbits for HD 82943:(I)Only $`\theta _1`$ librates about $`0^{}`$ (II)Case of alignment, $`\theta _1\theta _2\theta _30^{}`$ (III)Case of antialignment, $`\theta _1180^{}`$, $`\theta _20^{}`$, $`\theta _3180^{}`$. In addition, we noticed that 2%, 3% and 2% of the cases belong to Classes I, II and III respectively. The geometry of Class I was also discovered by Gozdziewski & Maciejewski (2001). But Classes II and III are quite a new discovery that exhibits the HD 82943 system can undertake apsidal lock (see Figure 3). Recently, Lee M. H.(2003, private communication) pointed out that the best-fit solution444The velocities for HD 82943 have not yet been published, and the best-fit orbital parameters taken from http://obswww.unige.ch. makes the orbits closer to be aligned than anti-aligned in the study of the 2:1 resonance. We used a semi-analytical (Ji et al. 2003b) model to investigate the apsidal motion for this system and found that although the two planets have high eccentricities up to 0.54 and 0.41, the system remain stable for 10 Myr because of both dynamical mechanisms of the 2:1 MMR and apsidal lock.

Several best-fit orbital solutions for GJ 876 appeared in literature (Marcy et al. 2001; Laughlin & Chambers 2001;Rivera & Lissauer 2001). Here we presented numerical results based on Laughlin-Chambers fit. Two families of orbits are found to be stable:Classes I (14%) and II (7%). In summary, for above two systems, we stress that the 2:1 MMR always takes place for stable cases and the apsidal alignment(Laughlin & Chambers 2001;Lee & Peale 2002) seems to have more likely chances.

On the other hand, as the two orbits of GJ 876 (or HD 82943) are not well separated, repeated interplanetary close approaches may lead to disruption of the systems:79% for GJ 876 rapidly become unstable, and 93% for HD 82943.

### 2.5 Other systems

Butler et al.(1999) reported the first triple-planet system of $`\upsilon `$ And, consisting of a 51 Peg-like planet with a 4.6-day orbit and two librating outer planets (Rivera & Lissauer 2000;Chiang, Tabachnik & Tremaine 2001;Chiang & Murray 2002). With the updated data(Paper II), we again examined several tests for integration of $`10^6`$ yr and discovered that two outer components can still remain in the apsidal lock about $`0^{}`$. As for 55 Cnc, we revealed that the inner two companions perform the asymmetric $`\mathrm{\Delta }\varpi `$-libration about $`250^{}`$ or $`110^{}`$ in the secular evolution (Ji et al. 2003a).

Nagasawa et al.(2003) showed that the two planets of HD 168443 and HD 74156 are both in circulation. Moreover,we performed additional experiments for HD 38529, but simply found the planets circulate. The systems of HD 168443, HD 74156 and HD 38529 are similar (Table 2), because each of them has a more massive outer planet and a close inner companion. The reason for circulation will be explained in §3. Nevertheless, if one of three systems can host an extra middle planet with new observations, there may exist librating planets similar to those of 55 Cnc and $`\upsilon `$ And.

## 3 Discussion on apsidal motion of $`\mathrm{\Delta }\varpi `$

However, one should make clear the conditions for libration or circulation of $`\mathrm{\Delta }\varpi `$ for two planets. Several theoretical works (Malhotra 2002;Lee & Peale 2003;Nagasawa et al. 2003;Zhou & Sun 2003)were concentrated on the problem. Beauge et al.(2003) presented analytical and numerical outcomes on the existence and location of stable equilibrium solutions about the systems in mean motion resonance. Here we present our results based on L-L theory. Following MD99, we have:

$$\frac{d\varpi _b}{dt}=\frac{e_c}{e_b}A_{12}\mathrm{cos}(\varpi _b\varpi _c)A_{11}$$
(2)
$$\frac{d\varpi _c}{dt}=\frac{e_b}{e_c}A_{21}\mathrm{cos}(\varpi _b\varpi _c)A_{22}$$
(3)
$$\frac{d\mathrm{\Delta }\varpi }{dt}=A_{22}A_{11}+(\frac{e_b}{e_c}A_{21}\frac{e_c}{e_b}A_{12})\mathrm{cos}\mathrm{\Delta }\varpi $$
(4)
where $`A_{11}`$, $`A_{12}`$, $`A_{21}`$ and $`A_{22}`$ are initially determined constants (see MD99) and $`e_b`$, $`e_c`$ are the eccentricities of each planet. Here we rewrite the coefficients in equation (4), such that $`B_1=A_{22}A_{11}`$ and $`B_2=(e_b/e_c)A_{21}(e_c/e_b)A_{12}`$. For hierarchical planetary systems (HD 168443, HD 74156 and HD 38529, the ratio of the semi-major axes $`\alpha =a_b/a_c1`$), we have $`|B_1||B_2|`$(see Table 2), thus the part of $`|B_1|`$ is dominant and $`\mathrm{\Delta }\varpi `$ circulates(see also Nagasawa et al.2003). According to (4), the critical values of $`e_b/e_c`$ for $`\mathrm{\Delta }\varpi `$-libration:

$$(\frac{e_b}{e_c})_\pm =\frac{\pm (\gamma 1)\sqrt{(\gamma 1)^2+4\gamma \beta ^2}}{2\beta }$$
(5)
where the plus sign in $`\pm `$ corresponds to the case for $`\mathrm{\Delta }\varpi `$-libration about $`0^{}`$, and the minus sign to that of $`\mathrm{\Delta }\varpi `$-libration about $`180^{}`$, with $`\beta =b_{3/2}^{(2)}(\alpha )/b_{3/2}^{(1)}(\alpha )`$ ($`b_{3/2}^{(j)}(\alpha )`$ are Laplace coefficients, $`j=1,2`$) and $`\gamma =(1/\sqrt{\alpha })(\mu _c/\mu _b)`$$`\sqrt{(1+\mu _c)/(1+\mu _b)}`$ ($`\mu _{b,c}`$ are respectively the ratio of the masses of Companion b (c) and host star). Hence, there are two requirements for apsidal libration of a pair of planets in the absence of mean motion commensurabilities. Firstly, the ratio of initial eccentricities of two planets should not be far from the critical values;secondly, the starting relative apsidal longitude should be satisfied to approach $`0^{}`$ or $`180^{}`$. As a paradigm of HD 37124 (see Figure 1), the initial ratio of the eccentricities $`(e_b/e_c)_00.30`$ is close to $`(e_b/e_c)_+=0.34`$, and $`\mathrm{\Delta }\varpi _07^{}0^{}`$, thus the two planets are aligned. We can apply the criteria to other systems.

## 4 Summary

In this work, we investigated the apsidal motion for the multi-planet systems:we find that the two planets of HD 37124, HD 12661, 47 Uma and HD 82943 separately undergo the apsidal alignment or antialignment. But the companions of GJ 876 and $`\upsilon `$ And are only in apsidal lock about $`0^{}`$, which means that no mirror configurations occur for them. Moreover, we again used the L-L theory to discern whether two planets of a system librate or not. The status of $`\mathrm{\Delta }\varpi `$-libration depends on the ratio of their initial eccentricities,semi-major axes and planetary masses and original relative apsidal longitudes. And it seems to be a selection effect for the companions in the multi-planet systems that are involved in the apsidal libration. However, this should be further testified by additional observations of the planetary systems.

We thank G. Laughlin, R.P. Butler, D.A. Fischer, M.H. Lee for informative discussions and an anonymous referee for good comments to improve the manuscript.This work is financially supported by National Natural Science Foundation of China(Grants 10203005, 10173006, 10233020) and Foundation of Minor Planets of Purple Mountain Observatory.


text:
# Dynamical ordering induced by preferential transitions in Planar Arrays of Superheated Superconducting granules

## ACKNOWLEDGMENTS

We acknowledge T. Girard for helpful discussions and a careful criticism of the manuscript.. We acknowledge financial support from Dirección General de Investigación Científica y Técnica (Spain) (Project BFM2000-0624-C03-02, BFM2002-02629) and Comissionat per a Universitats i Recerca (Spain) Projects (2001SGR97-0021). We also acknowledge computing support from Fundació Catalana per a la Recerca-Centre de Supercomputació de Catalunya (Spain).


text:
# Structure of a large social network

## Abstract

We study a social network consisting of over $`10^4`$ individuals, with a degree distribution exhibiting two power scaling regimes separated by a critical degree $`k_{\mathrm{crit}}`$, and a power law relation between degree and local clustering. We introduce a growing random model based on a local interaction mechanism that reproduces all of the observed scaling features and their exponents. Our results lend strong support to the idea that several very different networks are simultenously present in the human social network, and these need to be taken into account for successful modeling.

The ubiquity of networks has long been appreciated: complex systems in the social and physical sciences can often be modelled on a graph of nodes connected by edges. Recently it has also been realized that many networks arising in nature and society, such as neural networks watts\_strog , food webs food , cellular networks jeong2000 , networks of sexual relationships sex , collaborations between film actors watts\_strog ; ba  and scientists new\_pnas ; barab2001 , power grids watts\_strog ; classes , Internet routers faloutsos  and links between pages of the World Wide Web ajb ; kumar ; broder  all share certain universal characteristics very poorly modelled by random graphs er : they are highly clustered “small worlds” milgram ; watts\_strog ; watts  with small average path length between nodes, and they have many highly connected nodes with a degree distribution often following a power law ajb ; ba . The network of humans with links given by acquaintance ties is one of the most intriguing of such networks milgram ; pool\_kochen ; was\_faust ; watts , but its study has been hindered by the absence of large reliable data sets.

— The new data set  The WIW project was started by a small group of young professionals in Budapest, Hungary in April 2002 on the web site www.wiw.hu with the aim to record social acquaintance. The network is invitation–only and new members join by an initial link connecting to the person who invited them. New links are recorded between members after mutual agreement. This scheme results in members preferring to use their real names and effectively prevents proliferation of multiple pseudonyms. Because of the relatively short age of the network, links formed between people newly acquainted through the web site have a minimal structural effect; thus the majority of the links represent genuine social acquaintance and the WIW develops as a growing subgraph of the underlying social acquaintance network. Indeed, the growth process of the WIW network is essentially equivalent to the “snowball sampling” method well known to sociologists was\_faust , and to the crawling methods used to investigate the World Wide Web and other computer networks. We study the WIW network using two anonymous snapshots taken in October 2002 (with 12388 nodes and 74495 links) and January 2003 (with 17496 nodes and 127190 links).

The degree distribution of the WIW network is plotted on Figure 1. The graph shows two power law regimes

$$P(k)\{\begin{array}{cc}k^{1.0}& \text{ if }k<k_{\mathrm{crit}}\hfill \\ k^{2.0}& \text{ if }k>k_{\mathrm{crit}}\hfill \end{array}$$

The two regimes are separated by a critical degree $`k_{\mathrm{crit}}25`$. The exponent $`\gamma _22`$ of the large-$`k`$ power law falls in a range that has often been observed before in a variety of contexts sex ; ajb ; kumar ; broder ; ba ; faloutsos ; new\_pnas ; barab2001 ; classes ; jeong2000 . The value $`\gamma _11`$ of the small-$`k`$ power law exponent is much less common, observed before only in some scientific collaboration networks new\_pnas  and food webs food . The possibility of a double power law was discussed in barab2001 ; new\_coll , but the WIW network is the first data set which conclusively demonstrates the existence of double power law behaviour. The two snapshots give essentially identical distributions. Since the network grew by about 50% during this period, the described distribution can be regarded as essentially stationary in time.

The two scaling regimes in the degree distribution of the WIW graph are indicative of two distinct growth processes: the invitation of new members, and the recording of acquaintance between already registered members. The degree distribution of the invitation tree graph is shown on Figure 2, where a power law is observed for large degrees with an exponent $`\gamma 3`$. Since this distribution is qualitatively different from the total degree distribution, it is reasonable to conclude that there are at least two different types of social linking at play here: the network of friends defined by ties strong enough to warrant an invitation is different from the network of acquaintances that drives the mutual recognition, once both parties are registered.

The density of edges in the neighbourhood of a node is measured by the local clustering coefficient. For a node $`v`$ of degree $`k`$, the local clustering coefficient $`C(v)`$ is the number of acquaintance triangles of which $`v`$ is a vertex, divided by $`k(k1)/2`$, the number of all possible triangles. Figure 3 plots $`C(k)`$, the average of $`C(v)`$ over nodes of degree $`k`$, against the degree, showing the existence of a power law

$$C(k)k^{0.33}.$$

A relationship $`C(k)k^\alpha `$ was observed before in vazvesp ; rav\_barab , but with significantly larger exponents. Such power laws hint at the presence of hierarchical architecture rav\_barab : when small groups organize into increasingly larger groups in a hierarchical manner, the local clustering decreases on different scales according to such a power law.

The average clustering coefficient $`C0.2`$ is obtained as the average of $`C(v)`$ over all nodes. The average diameter between two members of the WIW network is about $`4.5`$. These two measures indicate the “small–world” nature of the WIW network in the sense of watts\_strog .

—Time development  As Figure 4 shows, the number of nodes of the WIW network grew approximately linearly with time. This appears to be related to the fact that the WIW network develops as a subgraph of the underlying social network, and thus the availability of new members is constrained by high clustering of the existing social links: a substantial proportion of the aquaintances of a newly invited member will have been invited already.

On the other hand, the number of edges also grew linearly with time, and thus the edge/node ratio only grew moderately during the existence of the network. This observation is in contradiction with any purely local time-independent edge creation mechanism. If every member of the network actively participates in edge creation indepently of its age in the network, the edge/node ratio would also increase linearly with time. This linear growth of the edge/node ratio was not observed in the network, and hence we conclude that the edge creatioon activity of members necessarily decreased with time.

The fact that the edge/node ratio changes little over time is consistent with the observed stationary nature of the degree distribution. To see this, consider a growing network with $`V(t)`$ nodes at time $`t`$ and a time-independent degree distribution $`P(k)`$ with $`_kP(k)=1`$ and finite first moment $`_kkP(k)`$. At time $`t`$, there are

$$n(k,t)=V(t)P(k)$$

nodes of degree $`k`$, and hence the number of edges is

$$E(t)=\frac{1}{2}\underset{k}{}kn(k,t)=\frac{V(t)}{2}\underset{k}{}kP(k).$$

Consequently the edge/node ratio $`E(t)/V(t)`$ is essentially constant, and it only changes because the maximal degree increases. This argument applies to the WIW network with stationary distribution

$$P(k)=\{\begin{array}{cc}\frac{c_1}{k^{\gamma _1}}& \text{ if }k<k_{\mathrm{crit}}\hfill \\ \frac{c_2}{k^{\gamma _2}}& \text{ if }k>k_{\mathrm{crit}}\hfill \end{array}$$

with $`\gamma _11`$, $`\gamma _22`$. This distribution is on the boundary of distributions with finite first moment: the first moment exists for $`\gamma _2<2`$ but not otherwise.

—Network modelling  A random graph process based on linear preferential attachment for the creation of new edges was proposed in ba  to account for the observed power laws in natural networks. Such a process leads indeed to a graph with a power law degree distribution ba ; bollob . However, this model is by definition macroscopic, requiring information about the entire network in every step. This assumption is realistic for the World Wide Web or some collaboration networks, where all nodes are “visible” from all others. For human social networks however, it is reasonable to assume some degree of locality in the interactions. Also, the original scale-free models are not applicable to networks with high, degree dependent clustering coefficients. These problems motivated the introduction of new models which use local triangle creation mechanisms watts ; hk ; vazquez ; deb ; jin , which increase clustering in the network. These models have degree-dependent local clustering, and can also lead to power law degree distributions, though no existing model of this kind shows a double power law.

We now present a new model to account for the observed properties of the WIW network. As mentioned above, the WIW can be viewed as a growing subgraph of the underlying social acquiantance graph. This suggests a model obtained by a two-step process, first modelling the underlying graph, and then implementing a growth process. The lack of available data on the underlying graph however prevents us from following this programme directly. We build instead a growing graph in one single process, choosing the local triangle mechanism as our basic edge creation method. This models the social introduction of two members of the WIW network by a common friend some time in the past, such edges being gradually recorded in the WIW network itself. The invitation of new members is modelled by sub-linear preferential attachment ba , motivated by experimental results on scientific collaboration networks barab2001 ; new\_timedev , where the data permits the analysis of initial edge formation. We also impose the constant edge/node ratio to be consistent with the observed stationary distribution. Note that this constant has to be tuned from the shape of observed distributions, and cannot be inferred from the WIW data directly. The reason for this is that the WIW has a disproportionate number of nodes of degree one (Figure 1), representing people who once responded to the invitation but never returned, which distorts the edge/node ratio without invalidating our other conclusions.

The precise description of our process is as follows.

* We begin with a small regular graph.
* New nodes arrive at a rate of one per unit time and are attached to an earlier node chosen with a probability distribution giving weight $`k^q`$ to a node of degree $`k`$, where $`q0`$ is a parameter.
* Internal edges are created as follows: we select two random neighbours of a randomly chosen node $`v`$, and if they are unconnected, we create an edge between them. Otherwise, we select two new neighbours of the same node $`v`$ and try again.
* A constant number of internal edges is created per unit time, so that the edge/node ratio equals a constant $`m`$ after each time step.

The degree distribution of graphs generated by our process is shown on Figure 5. We found a very robust large $`k`$ power law of exponent $`\gamma _22`$, essentially independently of the invitation mechanism. We measured the joint probability distribution of the degrees $`k`$,$`k^{}`$ of nodes connected by new internal edges, and found that for large values, it was proportional to $`kk^{}`$. This directly leads to a power law exponent of $`2`$ via standard mean-field arguments barab2001 . The small $`k`$ behaviour was found to be sensitive to the invitation mechanism; Figure 5 shows that a second power law only appears at a critical $`q`$. The critial value of $`q`$ depends on the edge/node ratio.

To test the hypothesis that the low degree power law is indeed related to the invitation mechanism, we plot on Figure 2 the degree distribution of the invitation tree of the model network for various values of the parameter $`q`$. At the critical value, we obtain a scale-free distribution with exponent $`3`$. Decreasing $`q`$ leads to a much sharper drop in the curve, with an exponential tail for $`q=0`$, whereas increasing $`q`$ above the critical value leads to a gelation-type behaviour: new nodes connect only to very large degree nodes.

Figure 1 shows that, for appropriate choice of parameters, the degree distribution of our model reproduces that of the real network extremely well. Figure 3 plots the dependence of the local clustering cofficient $`C(k)`$ as a function of the degree $`k`$ in our model network. While no simple power law can be observed, there is a clear trend of decreasing clustering with increasing degree, indicative of the presence of hierarchy in the model network rav\_barab .

—Conclusion  We have presented and analyzed a large new data set of a human acquintance network with a stable degree distribution which exhibits a new feature: two power law regimes with different exponents. The observed approximately constant edge/node ratio is a consequence of the stability of the degree distribution, and implies that the average activity of members is time dependent, whereas the growth of the number of nodes is constrained by the underlying social network. We also introduced a model which reproduces the observed degree distribution extremely well, and concluded that the small-$`k`$ power law is a related to the scale-free nature of the invitation tree, whereas the large-$`k`$ power law is a result of the triangle mechanism of social introductions. Our results show that human social networks are likely to be composed of several networks with different characteristics, and directly observable processes will exhibit a mixture of features resulting from distinct underlying mechanisms.

###### Acknowledgements.

We thank Zsolt Várady and Dániel Varga for access to the data of the WIW network, and Risi Kondor for helpful discussions.


text:
# Untitled Document

Can we observe galaxies that recede faster than light ? — A more clear-cut answer

T. Kiang

Dunsink Observatory, Dublin Institute for Advanced Studies,
Dublin 15, Ireland

> Abstract A more clear-cut answer to the title question is: Yes, if the expanding universe started with a big bang; No, if it started infinitely slowly.
>
> Key words: cosmology—cosmological models

1. Formulations of the Question The present Letter is a follow-up of a previous, much longer paper (Kiang 1997) on the same question stated in the title. The verbs in the question were left deliberately vague as to tense, because the previous paper dealt with both of the two formulations, (i) “Can we now observe galaxies that were receding faster than light at the time of emission ?” and (ii) “Can we eventually observe galaxies that are now receding faster than light ?”

In the present Letter I shall concentrate on the first formulation, which seems to be of more practical interest, for it can be re-stated in a sharper form: “Among the objects we now observe can we identify those that were receding faster than light at the time of emission ?”

It appears that, when so formulated, the question is capabale of quite a clear-cut answer. That is the main thrust of this Letter. But before I present it, I shall (1) summarize the relevant results so far obtained and (2) consider the question for the currently favored ($`\mathrm{\Omega }_M=0.3,\mathrm{\Omega }_\mathrm{\Lambda }=0.7,\mathrm{\Omega }_K=0`$) model (the “30/70 model” hereinafter). It was the result obtained with the 30/70 model that prompted me to seek, and eventually to find, the more general result.

2. Previous Results

The main relevant results of my (1997) paper are as follows.

1. This question of observability has nothing to do with horizon.

2. It depends solely on the form of the scale factor, $`R(t)`$, of the cosmological model. More precisely, it simply depends on whether or not two well-defined curves in the $`rt`$ plane ($`r`$ the coordinate distance, $`t`$ the cosmic time) intersect at some point in the past. The two curves are, 1) the curve that separates the regions of subluminal and superluminal recession velocities,

$$r_{v=c}(t)=1/\dot{R}(t),$$
(1)
and 2) the equation for our past light cone:

$$r_{\mathrm{PLC}}(t)=_t^{t_0}\frac{\mathrm{d}t}{R(t)},$$
(2)
$`t_0`$ being the present epoch. Here and throughout this Letter, we always use $`c=1`$ units.

3. It then followed that for the steady-state model, the answer is “no”; on the other hand, for all the three varieties ($`k=0,\pm 1`$) of the big bang ($`\mathrm{\Lambda }=0`$) model, the answer is “yes”; in particular, for the $`k=0`$ model (the ”standard” big bang), we have the result that all quasars with redshifts greater than 1.25 (and we now know thousands of such objects) had superluminal recession velocities at the time of emission.

3. The 30/70 Model The equations appear simplest when we use the Hubble units (Kiang 1987), i.e., when all times are in units of the current Hubble time $`H_0^1`$ and all distances, the current Hubble radius $`cH_0^1`$. Then, $`c`$ and $`H_0`$ will not appear in the equations (being equal to 1), and all the physical quantities are non-dimensional. In these units, the Friedmann equation for the 30/70 model is, (cf. e.g., Davis et al., 2003, Equation (11); Bondi 1960, p.80),

$$\dot{R}^2=0.3R^1+0.7R^2,$$
(3)
and we have the explicit solution (Bondi 1960, p.82, putting $`C=0.3,\lambda =3\times 0.7=2.1`$ in the Bondi formula for his case 2 (i)),

$$R^3(t^{})=0.2143[\mathrm{cosh}(2.510t^{})1],$$
(4)
where $`t^{}`$ is time reckoned from the big bang. Since for $`t^{}=0.964`$ we have $`R=1`$, the age of the universe in this model is $`t_0=0.964H_0^1`$. In the figures that follow, I shall be using time reckoned from the present epoch, $`t`$, rather than $`t^{}`$; we have $`t^{}=t+0.964`$.

From these expressions it is easy to calculate the two curves (1) and (2). These are shown in Fig. 1, along with the scale factor $`R`$. The scale factor for the standard big bang model is included for interest;—it is labelled $`R_1`$:

The $`v=c`$ locus (1) and our past light cone (2) intersect at point A $`(t=0.681,r=1.063)`$. At $`t=0.681`$, we have $`R=0.383`$ (Point B). This last value corresponds to $`z=1.61`$. Hence we may assert that according to the 30/70 model, all objects with redshifts greater than 1.61 have the property that, at the time of emission, they were receding from us faster than light.

This result is qualitatively the same as for the standard big bang model; the only difference is that, in the latter, the dividing redshift is even lower, at 1.25.

4. Towards a More Clear-Cut Answer. It thus appears that the answer to the stated question may not depend on whether or not there is a non-zero cosmological constant $`\mathrm{\Lambda }`$, that the answer is always “yes”, as long as the universal expansion started with a big bang, that is, as long as there was a point in the finite past, $`t_{\mathrm{BB}}`$, when $`R=0`$ and $`\dot{R}=+\mathrm{}`$. This conjecture is easily proved to be true as follows. Consider the behavior of the two curves (1) and (2) as we move backward in time from $`t=0`$ to $`t=t_{\mathrm{BB}}`$. Consider curve (2) first. Since $`R(t)`$ is always positive, $`r_{\mathrm{PLC}}`$ must start from 0 and increase monotonically to some finite value at $`t_{\mathrm{BB}}`$ (The integral in (2) converges at $`t_{\mathrm{BB}}`$, because generally, $`R(t^{})^{2/3}`$ as $`t^{}0`$). On the other hand, the $`r_{v=c}`$ curve (1) always starts at unity (when the Hubble units are used) at $`t=0`$, and ends at 0 at $`t=t_{\mathrm{BB}}`$, whether the decrease is monotonic, as in $`\mathrm{\Lambda }=0`$ cases, or whether it first increases then decreases, as may happen in a $`\mathrm{\Lambda }0`$ case. In all cases, the two curves must intersect at some time between the present epoch and $`t_{\mathrm{BB}}`$, and we have proved

> Lemma 1 If the universal expansion started with a big bang, then of the objects we now observe, those with redshifts above a certain calculable value have the property that, at emission, they were receding from us at speeds greater than $`c`$.

Now, what happens if the universe did not start with a bang, but very, very gently from time minus infinity ? The steady-state model is one such model and I found in my 1997 paper that here, the $`r_{v=c}`$ curve lies consistently one Hubble radius above the $`r_{\mathrm{PLC}}`$ curve; in fact, for all $`t<0`$,

$$r_{v=c}(t)=1+\mathrm{e}^{|t|},r_{\mathrm{PLC}}(t)=\mathrm{e}^{|t|}.$$
(5)
Hence, as already stated, the answer is “no” here.

Now, the steady-state model is a non-relativistic model. Of the relativistic models (more precisely, Friedmann models), only those with a non-negative mass parameter ($`\mathrm{\Omega }_M0`$) and a non-negative curvature parameter ($`\mathrm{\Omega }_K0`$) and the corresponding critical dark-energy parameter ($`\mathrm{\Omega }_\mathrm{\Lambda }=\mathrm{\Omega }_{\mathrm{\Lambda }_c}=(4\mathrm{\Omega }_K^3)/(27\mathrm{\Omega }_M^2)`$) start expansion infinitely slowly from a value of $`R`$ equal to $`(3\mathrm{\Omega }_M)/(2\mathrm{\Omega }_K)`$ (Bondi p.82). Felten and Isaacman (1986) have comprehensively reviewed all Friedmann models with non-negative mass density $`\mathrm{\Omega }_M`$ for all values of the cosmological constant $`\mathrm{\Lambda }`$.

Fig. 2 shows the $`R(t)`$ for the model ($`\mathrm{\Omega }_M=0.1,\mathrm{\Omega }_K=0.45,\mathrm{\Omega }_\mathrm{\Lambda }=1.35`$), as a typical representative of $`\mathrm{\Lambda }=\mathrm{\Lambda }_c`$ models. Also shown is the scale factor, $`R^{}(t)`$, for the limiting case ($`\mathrm{\Omega }_M=\mathrm{\Omega }_K=0,\mathrm{\Omega }_{\mathrm{\Lambda }_c}=1`$).

The limiting (0,0,1) model is just the de-Sitter model. Its scale factor is the same as that of the steady-state model: $`R^{}(t)=\mathrm{exp}(t)`$ for all $`t`$ between $`\mathrm{}`$ and $`+\mathrm{}`$. Hence, always using the asterisk to mark the limiting case, we have, from (5),

$$r_{v=c}^{}(t)>r_{\mathrm{PLC}}^{}(t)$$
(6)
for all $`t<0`$. Now, from Fig. 2 we see that, for all $`t<0`$,

$$R(t)>R^{}(t),\dot{R}(t)<\dot{R}^{}(t),$$
(7)
hence, from the definitions (1) and (2), we have, for all $`t<0`$,

$$r_{v=c}(t)>r_{v=c}^{}(t),r_{\mathrm{PLC}}(t)<r_{\mathrm{PLC}}^{}(t).$$
(8)
Combining the inequalities (8) and (6), we have, for all $`t<0`$,

$$r_{v=c}(t)>r_{\mathrm{PLC}}(t).$$
(9)
This result is valid for all ($`\mathrm{\Lambda }_c,\mathrm{\Omega }_M>0`$) models. Combining this with (6), we have that for all ($`\mathrm{\Lambda }_c,\mathrm{\Omega }_M0`$) models, the two curves (1) and (2) never intersected at any time in the past. Hence we have

> Lemma 2 If the universe started expansion infinitely slowly from either zero or a finite size, then all the objects we observe now had subluminal recession velocities at the time of emission

Combining the two Lemmas stated above, I now conclude with the following statement: “to the question whether among the objects we now see, there are some that were receding faster than light at the time of emission, the answer is ‘Yes’, if the universe started the expansion with a bang; ‘No’, if the expansion started infinitely slowly.”

Acknowledgement I thank Tamara M. Davis for stimulating correspondence.

References

Bondi H., 1960, Cosmology. Cambridge University Press. pp. 80, 82

Davis, Tamara M. et al., 2003, Amer. J. Phys., 71(4), 358-364

Felten J.E., Isaacman R., 1986, Rev. Mod. Phys., 58/3, 689-698

Kiang T., 1997, Chin. Astron. Astrophys., 21/1, 1–18. Chinese version in Acta Astrophys. Sinica (Tianti Wuli Xuebao) 1997, 17/3, 225-238

Kiang T., 1987, Quaterly J, R. astron. Soc., 28, 456-471