text:
# Helioseismic limit on heavy element abundance

## 1 Introduction

The precisely measured frequencies of solar oscillations have been used to probe the solar interior. The primary inversions of these observed frequencies yield the sound speed and density profiles inside the Sun. In order to infer the temperature and chemical composition profiles, we also need to know the input physics such as opacities, equation of state and nuclear energy generation rates (Gough & Kosovichev dog88 (1988); Kosovichev kos96 (1996); Shibahashi & Takata st96 (1996); Takata & Shibahashi tak98 (1998); Antia & Chitre ac98 (1998)). In all these works the heavy element abundance profile is assumed to be known; attempts to determine heavy element abundance profile from helioseismic data have not been particularly successful (Antia & Chitre ac99 (1999); Takata & Shibahashi tak01 (2001)) as the resulting inverse problem becomes extremely ill-conditioned. Fukugita & Hata (fuk98 (1998)) have obtained limits on heavy element abundance in the solar core using observed solar neutrino fluxes. It would be interesting to enquire if such limits can be independently obtained from helioseismic data.

In general, the computed luminosity in a seismically computed solar model is not expected to match the observed solar luminosity. By applying the observed luminosity constraint it is possible to constrain the input physics, particularly, the cross-section of proton-proton (pp) nuclear reaction. Antia & Chitre (ac98 (1998)) estimated this cross-section to be $`S_{11}=(4.15\pm 0.25)\times 10^{25}`$ MeV barns. Similar values have been obtained by comparing the computed solar models with helioseismic data (Degl’Innocenti, Fiorentini & Ricci inn98 (1998); Schlattl, Bonanno & Paterno bon99 (1999)). The main source of error in these estimates is the uncertainty in the $`Z`$ profile and, therefore, Antia & Chitre (ac99 (1999)) attempted to find the pp reaction rate as a function of $`Z`$ in the solar core. In all these works the plasma screening of nuclear reaction cross-sections was calculated using intermediate screening formulation of Graboske et al. (gra73 (1973)). The treatment of screening in stellar nuclear reaction rates is not yet adequately understood (Dzitko et al. dzi95 (1995); Gruzinov & Bahcall gru98 (1998)). Wilets et al. (wil00 (2000)) have done a sophisticated treatment of plasma screening and compared their results with earlier prescriptions. Their results indicate that for the solar core the intermediate screening treatment due to Mitler (mit77 (1977)) is better than that due to Graboske et al. (gra73 (1973)). It would thus be interesting to study the effect of different treatment of plasma screening on the helioseismically estimated pp reaction cross-section.

Antia & Chitre (ac99 (1999)) included the effect of heavy element abundance $`Z`$ only on the opacity of the solar material. If we make the reasonable assumption that the abundances of C, N, O also increase with $`Z`$, then the CNO cycle will become more effective in contributing to the nuclear energy generation in the solar core. At normally accepted values of $`Z`$ it is estimated that less than 2% of energy generated in the central region is produced by the CNO cycle (Bahcall et al. bp01 (2001)). But if $`Z`$ value is increased, this proportion will clearly increase and consequently, the pp reaction rate needs to be reduced to maintain the observed solar luminosity. In this work we demonstrate that this effect can be exploited to set an upper limit on $`Z`$ in the solar core.

## 2 The technique

With the use of accurately measured p-mode frequencies from the first year of operation of MDI instrument (Rhodes et al. mdi97 (1997)) we infer the sound speed and density profiles adopting a Regularised Least Squares technique (Antia a96 (1996)). The determination of thermal and chemical composition profiles, however, necessitates the use of equations of thermal equilibrium (Antia & Chitre ac98 (1998)) with the supplementary input of the equation of state, opacities and nuclear energy generation rates, provided we have a knowledge of the heavy element abundance profile inside the Sun. There is no guarantee for the resultant seismic model to yield the observed solar luminosity, $`L_{}=3.846\times 10^{33}`$ ergs s-1, unless we adjust the nuclear reaction rates slightly. This turns out to provide a valuable tool for constraining the nuclear reaction rates, in particular, the cross-section of the pp nuclear reaction which has not been measured in the laboratory because of the slow weak interaction rate. The pp reaction cross-section has only been calculated theoretically and it would be instructive to test the validity of calculated results using the helioseismic constraints.

There is an uncertainty of about 2% in evaluating the luminosity of seismic models because of possible errors in primary inversions (including an estimate of systematic errors), solar radius, equation of state, opacity, nuclear reaction rates for other reactions. Possible error due to uncertainties in treatment of plasma screening of nuclear reactions is not included in this error estimate. This can be estimated by using different treatments of plasma screening. The uncertainty arising from errors in $`Z`$ profiles is, however, much larger (cf., Antia & Chitre ac98 (1998)). We, therefore, estimate for each $`Z`$ profile, the range of cross-section of pp nuclear reaction, which reproduces the luminosity to within 2% of the observed value. The computed luminosity in the seismic models will, of course, depend on the metallicity, $`Z_c`$ in the solar core; we therefore, attempt to delineate the region in the $`S_{11}`$-$`Z_c`$ plane which gives the correct solar luminosity. In order to obtain the thermal structure we adopt the OPAL opacities (Iglesias & Rogers igl96 (1996)), the OPAL equation of state (Rogers, Swenson & Iglesias rog96 (1996)) and the nuclear reaction rates from Adelberger et al. (fusion (1998)). The plasma screening effects are calculated using either weak screening (Salpeter sal54 (1954)) or intermediate screening (Mitler mit77 (1977)). The effect of variation of $`Z_c`$ on nuclear energy generation rate is incorporated in the computations by assuming the abundances of all heavy elements to increase in the same ratio for obtaining the abundances of C, N, O.

## 3 Results

With the help of the inverted profiles for sound speed and density, along with a homogeneous $`Z`$ profile, covering a wide range of $`Z`$ values, we first calculate a seismic model by employing the equations of thermal equilibrium. For each central value of $`Z`$ we estimate the range of cross-section of pp nuclear reaction, which reproduces the luminosity to within 2% of the observed value. The numerical results are shown in Fig. 1, which delineates the region in $`Z_c`$$`S_{11}`$ plane that is consistent with helioseismic and luminosity constraints. This figure shows the results obtained using intermediate screening (Mitler mit77 (1977)). For comparison the central value obtained using weak screening is also shown in this figure. In order to compare these results with the earlier work of Antia & Chitre (ac99 (1999)), in these calculations the effect of $`Z`$ on nuclear energy generation rate is not included. Thus the difference is mainly due to treatment of plasma screening. It is clear that with these treatments of screening the current estimate of pp reaction cross-section (Adelberger et al. fusion (1998)) is consistent with helioseismic estimate within the expected error bars, though the best value is about 1.5% higher than the theoretical estimates. Further, there is not much difference between results obtained using the two different treatments of plasma screening. In general, the estimated value of $`S_{11}`$ using weak screening is smaller than that obtained using intermediate screening by about 1%.

With a view to study the effect of $`Z`$ on energy generation through CNO cycle, we repeat the calculations by varying the abundances of all heavy elements in the same proportion as $`Z`$ in the nuclear reaction network. These results are shown by long-dashed line in Fig. 1. We do not expect much difference at low values of $`Z_c`$, since CNO cycle produces a tiny fraction of energy in the Sun. But with augmented $`Z_c`$, this fraction becomes larger due to increased $`Z`$ as well as due to the resulting increase in the temperature, since the CNO reactions are more sensitive to temperature. As a consequence, in order to maintain the solar luminosity constraint we need to reduce $`S_{11}`$, with the result that around $`Z_c=0.035`$ the estimated value of $`S_{11}`$ decreases very sharply. This is because a significant fraction of the luminosity is then accounted for by the CNO reactions with seismically inferred temperature and composition profiles. For $`Z_c=0.035`$, CNO cycle accounts for about 8% of solar luminosity and this fraction, in fact, increases rapidly with $`Z_c`$. Thus, we can consider this as the upper limit on $`Z_c`$, which is somewhat larger than that inferred by Fukugita and Hata (fuk98 (1998)). A lower limit on $`Z_c`$ can also be obtained from this figure if we accept the range of theoretically computed values of $`S_{11}`$ (Bahcall & Pinsonneault bp95 (1995); Turck-Chiéze & Lopes tc93 (1993)) as representing the acceptable range of $`S_{11}`$. Assuming an upper limit of $`4.3\times 10^{25}`$ MeV barns for $`S_{11}`$ we get a lower limit of $`Z_c=0.006`$, which is again comparable to the lower limit obtained by Fukugita and Hata (fuk98 (1998)). If we adopt a lower limit of $`3.8\times 10^{25}`$ MeV barns for $`S_{11}`$, the resulting upper limit on $`Z_c=0.032`$. The exact limits on $`Z_c`$ will of course, depend on the assumed limits on $`S_{11}`$, but because of the sharp fall in the long-dashed curve in Fig. 1, the upper limit is not particularly sensitive to the assumed lower limit on $`S_{11}`$ and any value $`Z_c>0.035`$ would certainly be difficult to reconcile with seismic models.

In the foregoing calculations we have adopted a homogeneous $`Z`$ profile, which is perhaps not too realistic. Since most of the energy generation occurs in the solar core ($`r0.25R_{}`$) our results are unaffected by the $`Z`$ profile in the outer radiative region. To demonstrate this we also use the following $`Z`$ profiles :

$$Z=\{\begin{array}{cc}Z_c\hfill & \text{if }r<0.25R_{}\text{,}\hfill \\ Z_c+\frac{(Z_sZ_c)}{0.463}(r/R_{}0.25)\hfill & \text{otherwise,}\hfill \end{array}$$
(1)
and

$$Z=Z_c+(Z_sZ_c)r/(0.713R_{}),$$
(2)
where $`Z_s=0.018`$ is the known value of $`Z`$ at the solar surface. The seismically estimated value of $`S_{11}`$ for these $`Z`$ profiles is compared with those for the homogeneous profile in Fig. 2. As to be expected, for $`Z`$ profile given by Eq. 1 the results are almost identical to those obtained with the homogeneous profile, while for the linear $`Z`$ profile given by Eq. 2 the average value of $`Z`$ in the solar core would be somewhat lower and as a result the estimated $`S_{11}`$ values are somewhat larger, but the difference is comparable to the estimated errors. If instead of central value of $`Z`$ we had used the value of $`Z`$ at around $`0.1R_{}`$ for the linear profile, the two curves would have been almost identical. Thus it is clear that exact form of the $`Z`$ profile is not very important and the upper limit on $`Z_c`$ is only weakly dependent on the $`Z`$ profile.

Apart from uncertainties in $`Z`$ profile the nuclear reaction rate for the CNO cycle reactions are also rather uncertain and these may affect the inferred upper limit. To estimate this uncertainty we repeat the calculations using the profile given by Eq. 1 for $`Z`$, but with the reaction rate of the $`{}_{}{}^{14}\mathrm{N}+p`$ reaction reduced by 40%, which is the estimated uncertainty in this reaction rate (Adelberger et al. fusion (1998)). These results are also shown by the dotted line in Fig. 2. From the figure it is clear the uncertainties in CNO reaction rates will not particularly affect the upper limit on $`Z_c`$ in a significant way.

Table 1 lists the neutrino fluxes in a seismic model obtained using the Z-profile of a model including tachocline mixing of Brun et al. (bru02 (2002)), with the use of intermediate screening, for calculating the nuclear energy generation rates. In this case we need to increase $`S_{11}`$ by 1.6% over the currently accepted value. As noted by Brun et al. (bru02 (2002)), an increase in $`S_{11}`$ leads to a better agreement between solar models and the seismically inferred sound speed and density profiles. Thus a small increase in $`S_{11}`$ may still be required, but it is within the errors in helioseismic estimates and those in theoretically computed values. These neutrino fluxes are generally larger than those in seismic model of Antia & Chitre (ac98 (1998)), presumably because of reduction in $`S_{11}`$. The 8B neutrino flux in the seismic model is consistent with the current estimate of $`(5.09\pm 0.65)\times 10^6`$ cm-2 s-1 using the neutral current channel of Sudbury Neutrino Observatory (SNO) (Ahmad et al. sno02 (2002)).

## 4 Discussion and Conclusions

With the help of inverted sound speed and density profiles, it is possible to infer the $`T,X`$ profiles in the solar interior, provided the $`Z`$ profile and the input physics are known. The resulting seismic models have the correct solar luminosity, provided the heavy element abundance $`Z_c`$ in the solar core and the cross-section for pp nuclear reaction are within the shaded region shown in Fig. 1. It appears that the currently accepted values of $`Z_c`$ or $`S_{11}`$ need to be increased marginally to make them consistent with helioseismic constraints. The required increase is within the error estimates. The higher estimates for $`S_{11}`$ obtained earlier were due to differences in treatment of plasma screening. With the use of weak (Salpeter sal54 (1954)) or intermediate screening due to Mitler (mit77 (1977)) the theoretically estimated value of $`S_{11}`$ is in reasonable agreement with seismically estimated value. With a $`Z`$ profile in a standard solar model N0 of Brun et al. (bru02 (2002)), the seismically estimated value of $`S_{11}`$ is $`4.07\times 10^{25}`$ MeV barns with the use of intermediate screening and $`4.02\times 10^{25}`$ MeV barns for weak screening used while calculating the nuclear energy generation rate.

If the value of heavy element abundance in the solar core, $`Z_c`$, is increased beyond 0.035 the CNO cycle generates a good fraction of solar luminosity and it is not possible to get any consistent seismic model unless $`S_{11}`$ is decreased substantially below the accepted value. This puts a clear upper limit on the heavy element abundance in the solar core, which is comparable to that independently obtained by Fukugita and Hata (fuk98 (1998)). This upper limit is not very sensitive to the $`Z`$ profile or to the uncertainties in the CNO reaction rates. Even if the additional heavy elements in solar core do not include CNO, the effect of opacity alone will also put an upper limit on $`Z_c`$, but in that case the limit will depend on the assumed lower limit on $`S_{11}`$. For currently accepted theoretical limits, the upper limit in this case turns out to be around 0.035.

Note that the seismic model satisfies the normal stellar structure equations, though the inferred $`X`$ profile may not match that given by an evolutionary solar model. Further, since the seismic model is confined to the radiative interior, it will not be possible to match it to an acceptable convection zone model, unless the heavy element abundance at the top of the radiative region is close to the known surface value. Thus a relatively high value of $`Z`$ in the solar core can be realised only if the $`Z`$ profile has a significant gradient in the radiative region. From experiments with varying $`Z`$ or $`S_{11}`$ in evolutionary solar models also it is known that an increase in $`Z`$ can be compensated for by a reduction in $`S_{11}`$ to match the seismically inferred sound speed and density profiles (Brun et al. bru02 (2002)). Thus it is quite possible that a comparable upper limit on $`Z_c`$ may be obtained from these models if we restrict the range of $`S_{11}`$, provided the initial $`Z`$ profile at the zero age main sequence stage is not homogeneous.

###### Acknowledgements.

SMC is grateful to DAE-BRNS for support under the Senior Research Scientist Scheme and to Ian Roxburgh for supporting his visit to Queen Mary, University of London under the Leverhulme Trust Visiting Professorship scheme. We thank the Referee, S. Degl’Innocenti for valuable suggestions, which have led to an improvement in the presentation of our results.


format: text

text:
# Generalized (m,k)-Zipf law for fractional Brownian motion-like time series with or without effect of an additional linear trend

## 1 Introduction

Many phenomena which contain discernible events which can be counted can be ranked according to their frequency, and a so called Zipf plot can be drawn.  Very often a $`quasi`$ linear relationship is found on a log-log plot. The slope corresponds to an exponent $`s`$ describing the frequency $`P`$ of the $`cumulative`$ occurrence of the events according to their rank $`R`$ through, e.g. $`P(>R)R^s`$.

Such a ($`Zipf`$) $`power`$ $`law`$ is found in many cases, with $`s1`$: see the distribution of $`income`$ of individuals or companies in countries (Pareto distribution) , in economy with the $`size`$ of companies , in $`earthquakes`$ (Gutenberg-Richter law) , in $`city`$ distribution , etc. … The Zipf law universal feature is thought to originate from stochastic processes , in particular when they can be modeled as random walks in a log scale , - though it is still often said in a lay language that the law is a description (or a result ?) of uniformity and diversity.

A simple extension of the Zipf analysis is to consider $`m`$-letter words, i.e. the words strictly made of $`m`$ characters without considering the white spaces. The available number of different letters or characters $`k`$ in the alphabet should also be specified. A power law for the word frequency $`f(R)`$ is expected to be observed

$$fR^\zeta .$$
(1)
The Zipf exponent can be estimated through the derivative of the best linear fit on a log-log graph. There is no theory at this time predicting whether the exponent is a function of ($`m,k`$). Elsewhere we had already shown that the Zipf exponent could be different, call it $`\zeta ^{}`$ if the frequency $`f`$ of occurrence is normalized with respect to the theoretical frequency $`f^{}`$ which should be that expected for pure, or unbiased, (stochastic) Brownian processes, thus

$$f/f^{}R^\zeta ^{}.$$
(2)
E.g. suppose a binary alphabet, i.e. made up of two characters, $`u,d`$ (or 0 and 1 in electronics) with occurrence $`p_u`$ and $`p_d`$; the theoretical frequency $`f^{}`$ for the number $`n`$ of characters, say of the type $`d`$ in a word of length $`m`$ is $`f^{}=p_u^{(mn)}p_d^{(n)}`$.

Previous Zipf-like analyses have often neglected the possibility that the ($`u,d`$) distribution could be biased, over a finite size time interval, thus a finite tendency might exist, - and in fact could have a quite varied structure, as when in finance a trend is considered to be mimicked by some moving average . It seems clear that if a tendency is positive the number of positive volatility values is larger than the number of negative ones (and conversely) . Therefore the Zipf law exponent might be affected because of some bias in the ranking. Whether or not the $`\zeta `$ exponent depends on the bias is briefly examined here.

The final goal of this series of investigations is to apply the idea in the study of financial temporal series, or generally translating a time series into a sentence, a particular letter corresponding to a particular variation of a signal. This is in line with previous studies in econophysics e.g. in order to search whether some investment strategy can be derived in particular from a Zipf law observation. This would fall into the same type of studies as those implying the detrended fluctuation analysis (DFA).

Here we analyze time series based on a one dimensional fractional Brownian motion  characterized by the so called Hurst exponent $`H`$.

Why such a series? Because the FBM is not a Markovian process since its value at a given time depends on all past points, whence we consider that it can be a useful model for modeling financial data time series. In fact, Peters  has shown that a FBM is a good model for describing returns (but it does not work for $`options`$) in financial series.

## 2 Data

In order to develop a FBM we have used Rambaldia and Pinazza algorithm . According to the latter a FBM can be obtained from

$$B_H(j)=\underset{i=1}{\overset{j}{}}\omega _{ji+1}\xi _i,$$
(3)
where $`\xi `$ represents the ”walker” position during a time interval $`\mathrm{\Delta }t`$; $`\xi `$ is a random variable to be extracted from a Gaussian distribution with zero mean and a given variance. Thereby the signal is a stochastic one, with a diffusion growing with time as $`2H`$. The weight function $`\omega _{ji+1}`$ is given by

$$\omega _{ji+1}=\frac{\gamma }{H+\frac{1}{2}}[(ji+1)^{H+\frac{1}{2}}(ji)^{H+\frac{1}{2}}],$$
(4)
where $`\gamma `$ is such that $`<B_H^2(1)>=1`$. If $`H=1/2`$ one has the usual Brownian motion. The signal is said to be persistent for $`H>1/2`$, and antipersistent otherwise. There has been some conjecture  , sometimes thought to be proven like a theorem that

$$\zeta ^{}=\zeta =|2H1|.$$
(5)
We have created six different FBM signals with $`H`$ values respectively equal to 0.17, 0.41, 0.47, 0.60, 0.67, 0.82. The series have a 16 384 ( = 214) length. They are normalized such that

$$B_H(t)=\frac{B_H(t)1.1B_{H,min}(t)}{1.1B_{H,Max}(t)}.$$
(6)
The coefficient 1.1 in the denominator and numerator is used in order to avoid zero and unity values. The series are shown in Fig.1. Their characteristic is summarized in Table 1. We have recalculated the Hurst exponent  by the box counting method. The error bar is given in Table 1, as $`\mathrm{\Delta }H`$. The error bars are those resulting from a root mean square analysis. The linear trend has been measured and is reported to be of the order of $`10^5`$, obviously due to the finite size of the system.

In order to apply the Zipf method and extract the $`\zeta ^{}`$ exponent relating the word frequency and its rank, we have translated the FBM signals into a text based on two $`letters`$, $`u`$ and $`d`$, occurring with a frequency $`p_u`$ and $`p_d`$ respectively. The bias defined as $`ϵ=p_u0.5`$ has been measured. The bias and also the signal tendency have been observed as a function of time but are not shown here for lack of space. They decay quickly, can be positive or negative and are of the order or smaller than $`1\%`$ after 4000 time steps.

The partial distribution functions (pdf) of the logarithm of the signal volatility , i.e.

$$Z(t)=\mathrm{ln}(y(t+\mathrm{\Delta }t))\mathrm{ln}(y(t))$$
(7)
have been fitted to Gauss and stretched exponential distributions. Instead of ranking the $`Z(t)`$ values in constant size box histograms, one can as suggested by Adamic  use binned histograms with bin size increasing exponentially, thereby obtaining other best fit parameters. The stretched exponential distribution seems to well describe the signal fluctuations. A Kolomogorov-Smirnov test has also been made.  Notice that even though the KS distance increases when the exponential size box is used, this scheme reduces the uncertainty values and is found to lead to better fits.

Another test of the stochasticity (or not) of the data is based on the surrogate data method  in which one randomizes either the sign of the fluctuations or shuffles their amplitude and Fourier transform the resulting signal in order to observe whether a white noise signal is so obtained. Finally we have observed whether the error bars (or confidence intervals) of the raw signal and the surrogate data signal (not shown here) overlap. The characteristic of the spectral functions $`S(f)f^\beta `$ so obtained are available from the authors if necessary. The results allow us to conclude that the above FBM signals are satisfactory for further treatment in presence of a to be pre-imposed bias.

Notice that similar histograms of such ”words” were already published for ($`m=3`$, $`k=2`$) and ($`m=5`$, $`k=3`$) respectively, but the authors were more interested in deviations from randomness than in the Zipf exponent.

## 3 Zipf Analysis of Fractional Brownian Motion Raw Data

Consider the ($`m`$,2) Zipf method, thus for an only two character alphabet, and words of arbitrary length $`m`$. Therefore there are $`2^m`$ different possible words. We searched whether these words exist in the series of Fig.1, counted them and ranked them in decreasing order, as shown in Fig. 2, for $`2m8`$. The $`\zeta ^{}`$ (and $`\zeta `$, also but not shown) exponents seem to increase with $`m`$ (Fig. 3). The result may be attributed to the finite size of the series, if one realizes that for such series the number of long words necessarily occurs more rarely than the number of short ones.

As a test of finite size effects, we have successively removed one by one the less frequently occurring words, and recalculated the $`\zeta ^{}`$ (and $`\zeta `$) exponents, in some sense taking a rank $``$ zero limit, or first derivative. The exponents resulting from the average of the latter values are shown in Fig. 4. It is found that for large $`H`$, i.e. a persistent signal, the exponents are rather constant, but still markedly vary for the antipersistent signals.

As a further step, we can also consider whether there is a Zipf law $`evolution`$, i.e. search for the $`\zeta ^{}`$ exponent evolutions, as a function of time $`t`$ or as a function of the number of points in the series. We have calculated values of $`\zeta ^{}`$ and $`\zeta `$ for the (first) box containing 100 points, then for a box containing 200, 300, … etc points up to 16 000. The evolution (Fig. 5) is rather drastic for the first, say, 6000 points but is moderately stable thereafter.

Next recall the so called local (or better $`instantaneous`$) DFA method.  In order to probe the existence of correlated and/or decorrelated sequences, a so-called observation box of finite size can be constructed, and is placed at the beginning of the data. A DFA is performed on the data contained in that box. The box is then moved along the historical time axis by a few points toward the right along the sequence. Iterating this procedure for the sequence, a ”local measurement” of the ”degree of correlations” is obtained, i.e. a local measure of the Hurst exponent in the DFA case. The results indicate that the $`H`$ exponent value varies with time. This is similar in finance to what is observed along $`DNA`$ sequences in biology where the $`H`$ exponent drops below 1/2 in so-called non-coding regions. Doing the same here, we obtain a $`local`$ Zipf law and $`local`$ Zipf exponent. To show a full list of figures or data as a function of $`m`$ $``$  would generate a quite aversive stimulus in the reader, - therefore only the case $`m`$ =5 is illustrated here. This value is so chosen within the financial idea background having motivated this study, i.e. $`m=5`$ is the (true) length of a week !

Results for the three FBM signals, with $`H`$ close to 0.5, like in financial time series signals, are illustrated in Figs. 6-8, considering windows (boxes) of size 250, 500 and 1000 respectively moved along the signal. These values correspond to 1, 2 and 4 year type investment window in finance. Notice that the $`local`$ exponents are usually larger than the corresponding average one, in some sense corroborating the previous finite size effect analysis results. In turns it seems that the method is also of interest in order to observe short range correlation in fractional Brownian motions, and non ergodic properties of finite size series.

## 4 Zipf Analysis of Fractional Brownian Motion + Linear Trend Data

Next consider a FBM on which a linear trend with amplitude $`A_L`$ is $`additionally`$ superposed (and is equal to zero at the origin), i.e. we add

$$y_L(i)=A_Li$$
(8)
to Eq. (4). We have taken values of $`A_L`$ ranging from $`2^{16}`$ to $`2^7`$ depending on the FBM considered. The case $`k`$=2 is only considered at this stage, again as a first step and having in mind scales used in other works in econophysics and electronics. The Zipf plots are shown in Fig.9. For small $`H`$ a sharp drop is still seen at large $`R`$ values. The amplitude (and trend) effect is hardly noticed at small $`H`$, even with a large trend, while the linearity (on a log-log plot) is interestingly observed for larger $`H`$.

The effect of the trend when $`m`$ is larger than 2, thus up to $`m`$=8, as in the preceding section is illustrated in Fig. 10. As should be expected the exponent $`\zeta ^{}`$ is rather stable (= $`\zeta _0^{}`$ ) for a small trend (or slope) $`A_L`$ and small $`m`$, but markedly increases when $`m`$ increases. The $`\zeta ^{}`$ stability is observed apparently below some sort of crossover amplitude $`A_{L,cr}(m)`$, … depending on $`m`$. The variation above such a crossover amplitude follows a power law ($`\zeta ^{}A_L^{\theta _L}`$ ) which has been determined. The values of the exponent and of the power law variation on both sides of the $`A_{L,cr}`$ are given in Table 2, with their error bar $`\mathrm{\Delta }\theta _L`$.

We should point out that there seems to be some difference whether the FBM signal is persistent or not. For antipersistent signals, $`\zeta ^{}`$ is quasi $`A_L`$ independent. The Brownian motion is clearly an $`in`$ $`between`$ case.

## 5 Conclusions

We have considered fractional Brownian motions on which a linear trend is superimposed. We have translated the signals each into a text based of two letters $`u`$ and $`d`$, (or 0 and 1 bits) according to the fluctuations in the corresponding random walk. We have studied the variation of the Zipf exponent(s) giving the relationship between the frequency of occurrence of words of length $`m8`$ made of such ”letters” for a binary alphabet. We have searched how the $`\zeta ^{}`$ exponent of the Zipf law is influenced by the trend and its amplitude. We can distinguish finite size effects, and results depending whether the starting FBM is persistent or not, depending on the Hurst exponent. It seems that $`\zeta ^{}`$ varies as a power law in terms of $`m`$. This seems to be due to the fact that due to the trend a marked bias occurs in the relative fluctuations, thus in the word occurrences. It seems proven, even though it might have been expected so, that a persistent signal short range correlations as analyzed through the Zipf method is more influenced by a trend than an antipersistent signal.

In the spirit of the local DFA, we have introduced considerations based on the effect of finite sizes of texts, and on the notion of a local (or ”time” dependent) Zipf law, in order to show their influence on the exponent values.

Finally, coming back to the Czirok et al. conjecture, it seems of interest to display the relationship for both studied classes of cases, i.e. (i) the raw FBM, and (ii) the FBM + linear trend. The results are shown in Figs. 11-12. Even within considerations taking into account finite size effects and a poor man statistical analysis, it appears that such a law only holds near $`H=0.5`$.

Other considerations are in order showing that many other cases can still be considered: first one could wonder more about signal stationarity effects. Next, either a non linear (thus like a power law or a moving average) trend or a periodic background could be superposed on the raw signal. Also multiplicative rather than additive trends could be used. This should be put in line with the remarkable studies of Hu et al.  on DFA but is a work an order of magnitude higher here because there are mk parameters to consider. Applications of the above to financial data will be presented elsewhere .

Acknowledgments

We don’t thank anyone for stimulating discussions and comments, but the referee for imposing precision and conciseness.

Figure Captions

Figure 1 – Six fractional Brownian motions studied in the text with characteristics summarized in Table 1 characterized by an $`H`$ exponent

Figure 2 – ($`m`$,2) Zipf plots of the six FBM for $`2m8`$

Figure 3 – Evolution of the $`\zeta ^{}`$ exponent with $`m`$ for different $`H`$ values

Figure 4 – Evolution of the $`\zeta ^{}`$ exponent calculated when successively removing one by one the less frequently occurring words as a function of $`m`$ for different $`H`$ values

Figure 5 – Zipf law $`\zeta ^{}`$ exponent $`evolution`$ as a function of time $`t`$ starting from the (first) box containing 100 points, then for a box containing 200, 300, … etc. points up to 16 000

Figure 6 – Time dependence of the $`local`$ Zipf exponent in a box of size $`T`$ displaced along the FBM with $`H`$= 0.41; three size boxes are illustrated : 205, 500 and 1000

Figure 7 – Time dependence of the $`local`$ Zipf exponent in a box of size $`T`$ displaced along the FBM with $`H`$= 0.47; three size boxes are illustrated : 205, 500 and 1000

Figure 8 – Time dependence of the $`local`$ Zipf exponent in a box of size $`T`$ displaced along the FBM with $`H`$= 0.60; three size boxes are illustrated : 205, 500 and 1000

Figure 9 – ($`m`$,2) Zipf plots of the six FBM+linear trend for $`2m8`$ with slope $`A_L`$ of the linear trend given in insert

Figure 10 – Variation of the $`\zeta ^{}`$ exponent as a function of the slope trend when $`m`$ is larger than 2 and less than 8

Figure 11 – ”Verification” of the relationship $`\zeta =|2H1|`$ for the six FBM

Figure 12 – ”Verification” of the relationship $`\zeta =|2H1|`$ for the six FBM + linear trend for different trends


text:
# Multifractal analysis of DNA walks and trails

## I Introduction

The search for patterns in nature and their interpretation in terms of general principles is one of the main purposes of science. The explosive accumulation of DNA sequence data in the last two decades has provided a rich source of raw material that hides valuable hints about the evolutionary mechanisms of the genome organization. Unveiling the patterns in those sequences has become an exciting challenge to the present generation of statistical physicists and information scientists. In that vein, a truly remarkable result was the finding that intron-containing DNA coding regions exhibit long-range power-law correlations extending across more than $`10^4`$ nucleotides, whereas intronless coding regions display only short-range correlations Peng92 ; Peng94 .

Characterizing the long-range correlations in DNA sequences is a highly nontrivial task because of the mosaic structure of DNA consisting of patches with different nucleotide composition Karlin93 . In fact, in order to use the standard techniques (e.g., power spectrum and R/S fluctuation analysis) to study the correlations of DNA sequences, it is necessary first to eliminate the local biases in nucleotide composition (trends), thus avoiding spurious effects due to the mosaic character of the sequence. In the context of time series or records there are two main techniques to eliminate these trends, namely, the detrended fluctuation analysis (DFA) Peng94  and the wavelet transform (WT) Chui92 . The former is a physically appealing ad hoc technique of easy implementation and the latter is a mathematically well-established transform used to study the regularity of arbitrary functions via the systematic elimination of local polynomial behavior. Both techniques can be readily applied to the analysis of ordered linear sequences, such as DNA, by considering the reading direction as the time axis (see e.g. Peng94 ; Galvan96  for the DFA and Arneodo95 ; Tsonis96  for the WT analyses of DNA sequences).

A common strategy to represent graphically a given DNA sequence consists of transforming it into a $`d`$-dimensional random walk, so-called DNA walk, by associating a space direction to each nucleotide (A,C,T,G) or class of nucleotides (e.g. purine and pyrimidine) and the time direction to the reading direction of the sequence Peng92 . Here, A, C, T and G are the bases adenine, cytosine, thymine, and guanine, respectively. In particular, a very popular representation is the one-dimensional walk in which a purine (A or G) at position $`i`$ is associated to one step down ($`y_i=1`$) while a pyrimidine (T or C) is associated to one step up ($`y_i=1`$). We note that the purine-pyrimidine classification is related to the hydrophobic-hydrophilic characteristics of the encoded amino acids. Of course, there are many alternative DNA walks and the most complete one that considers the four bases equivalently in base space is the three-dimensional DNA walk Luo98 , in the sense that its projections on appropriate axes or planes recover all possible one and two-dimensional walks.

For the sake of concreteness, in this section we will focus on the fluctuations of the cumulative variable of the purine-pyrimidine random walk only, defined as

$$s(n)=\underset{i=1}{\overset{n}{}}y_in=1,\mathrm{},N,$$
(1)
where $`N`$ is the length of the sequence. This variable is displayed in Fig. 1 for the intron-rich human beta-globin region (GenBank name HUMHBB) for which $`N=73326`$ as well as for a shuffled sequence with same length and nucleotide composition. The statistical quantity that characterizes the fluctuations of $`s(n)`$ is the root mean square fluctuation $`F_2(l)`$ defined as

$$F_2\left(l\right)=\left\{\overline{\left[\mathrm{\Delta }s\left(l\right)\right]^2}\left[\overline{\mathrm{\Delta }s\left(l\right)}\right]^2\right\}^{1/2}$$
(2)
where $`\mathrm{\Delta }s(l)=s(n+l)s(n)`$ and the bars indicate an average over all positions $`n`$ in the record. We expect $`F(l)`$ to increase with increasing $`l`$, the length of the window considered. Explicitly, in the case of stationary series the fluctuations are described by a power law

$$F_2\left(l\right)l^H$$
(3)
where the scaling exponent (Hurst exponent) $`H`$ equals $`1/2`$ in the case of random and short-ranged correlated sequences. Any value different from $`1/2`$ is evidence of the existence of long-range correlations in the sequence Mandelbrot82 ; Voss89 ; Feder88 . In addition, for self-affine records (e.g. DNA walks) there is a simple relation between $`H`$ and the local fractal dimension of the record, $`D_r=2H`$. We note that the power law behavior implies that there is no characteristic length scale, i.e., very large fluctuations are likely due to the same mechanisms as smaller ones.

Complex records are unlikely to be fully characterized by a single scaling exponent as evidenced by the case where the scaling behavior is different in distinct parts of the series so that $`H`$ is actually dependent on the part of the record being analyzed. A complete framework to describe the situation where a multitude of scaling exponents is required is provided by the multifractal formalism Feder88 ; Halsey86  and, in this sense, the multifractal spectrum may be viewed as the ultimate tool to characterize a stationary time series. However the direct calculation of $`F_2(l)`$ using Eq. (2) or the application of the standard multifractal formalism yield wrong results for nonstationary time series that are affected by local trends. In this contribution we apply a generalization of the DFA recently proposed by Kantelhardt et al. for the characterization of nonstationary time series to study the multifractal properties of DNA walks Kantelhardt02 . In agreement with a previous multifractality analysis based on the wavelet transform modulus maxima (WTMM) method Arneodo95 , we find that the DNA walks exhibit a simple monofractal scaling behavior. The multifractality analysis of one-dimensional DNA-walks is the subject of Sec. II.

A seemingly alternative graphical representation of DNA sequences that has received less attention and so less criticism than the aforementioned DNA walks is the so-called DNA pseudo random-walks Gates86 ; Berthelsen92 ; Berthelsen94 ; Glazier95 ; Oiwa02 . They are simply the projection of the DNA walks into their $`d`$ spatial components, i.e., the trails of the DNA walks. Of course, drawing these trails makes sense for $`d2`$ only, and Fig. 2 illustrates such a trail for a $`d=2`$ representation of the HUMHBB sequence where right, left, up and down steps correspond to the presence of nucleotides T, A, G and C, respectively. This representation preserves some basic symmetries of the DNA such as complementarity, reflection, substitution and compatibility Berthelsen92 .

Despite the clear visual evidence of local trends in the trails exhibited in Fig. 2, it has been a common practice to apply standard fractal and multifractal methods (e.g., box counting Mandelbrot82  and sandbox Vicsek92 ) to characterize their shapes without much concern about the mosaic structure of the underlying DNA sequence. Actually the efforts focused on the understanding of the disturbing effects of the finite length of the sequences on the estimate of the fractal dimension $`D_t`$ as well as on the multifractal characterization of the trail. For instance, although it is well-known that in two-dimensions an infinite length true random walk (i.e., a sequence of bases generated at random) is space filling and so $`D_t=2`$, the box-counting and sandbox methods yield estimates significantly lower than 2 even for walks as large as $`10^5`$ steps Berthelsen92 . Moreover, finite-length random walks exhibit an effective multifractal spectrum due mainly to crossover effects between the usually distinct bulk and surface properties of the trail Berthelsen94 . These caveats, however, have not prevented claims that relatively short (typically $`2\times 10^4`$ bp) mitochondrial DNA genomes have a definite multifractal structure Oiwa02 . In Sec. III of this contribution we invoke a classical result of fractal theory to link the Hurst exponent of the DNA walk to the fractal dimension of the corresponding trail in space (see, e.g., Mandelbrot82 ; Voss89 ), then arguing in favor of the monofractality of DNA trails too.

## II DNA walks

Consider a record, such as that given in Eq. (1), where the variables $`y_i`$ are not necessarily Ising variables but must form a compact support, i.e., $`y_i=0`$ for a very small fraction of the elements only. The generalized multifractal DFA (MF-DFA) involves the following steps (we refer the reader to ref. Kantelhardt02  for a thorough presentation of the method). First, divide the entire sequence into $`N_l=\mathrm{int}(N/l)`$ nonoverlapping segments of length $`l`$. Second, for each segment $`\nu =1,\mathrm{},N_l`$ of length $`l`$ calculate the local trend by a least-square fit of the record in the segment, and denote by $`s_\nu (n)`$ the fitting polynomial in segment $`\nu `$. Then evaluate the variance

$$F^2(l,\nu )=\frac{1}{l}\underset{n=1}{\overset{l}{}}\left\{s\left[\left(\nu 1\right)l+n\right]s_\nu (n)\right\}^2$$
(4)
for each segment. An important parameter is the order $`m=0,1,2,\mathrm{}`$ of the polynomial $`s_\nu (n)`$ used in the fitting procedure. The choices of $`m`$ correspond to different orders of DFA, denoted by DFA$`m`$, which differ in their capability of eliminating trends in the series. For instance, DFA$`0`$ can eliminate only constant trends in the series, DFA$`1`$ eliminates constant as well as linear trends, and so on, so that a comparison of the results for different orders of DFA yields information on the type of trend in the series. Third, determine the order $`q`$ fluctuation function averaging over all segments

$$F_q\left(l\right)=\left\{\frac{1}{N_l}\underset{\nu =1}{\overset{N_l}{}}\left[F^2(l,\nu )\right]^{q/2}\right\}^{1/q},$$
(5)
where in general the index variable $`q`$ takes on any real value, except zero. As the final step, analyze log-log plots of $`F_q(l)`$ versus $`l`$ for each value of $`q`$ in order to determine the scaling behavior of the fluctuation functions. For large values of $`l`$, we expect these functions to increase with increasing $`l`$ as a power law

$$F_q\left(l\right)l^{h(q)}.$$
(6)
Clearly, for $`q=2`$ one has $`h(2)=H`$ by construction, and in the case that the scaling behavior of $`F^2(l,\nu )`$ is identical for all segments $`\nu `$, i.e., the record is a monofractal, the exponent $`h(q)`$ is independent of $`q`$, as expected. The definition (5) is well-suited to detect discrepancies in the scaling behavior of the large and small fluctuations. In particular, for positive $`q`$, $`h(q)`$ yields the scaling behavior of the segments with large fluctuations while for negative $`q`$, $`h(q)`$ yields the scaling of the segments with small variances. In addition, there is a simple relation between $`h(q)`$ and the scaling exponents $`\tau (q)`$ defined by the standard partition function-based multifractal formalism Kantelhardt02 ,

$$\tau (q)=qh(q)1.$$
(7)
The importance of this relation is that it allows comparison of the results obtained in the MF-DFA scheme with those of the standard multifractal analysis in the case of stationary series, and with those of the WTMM in the case of nonstationary series.

In Fig. 3 we show the multifractal spectrum $`\tau (q)`$ obtained by the application of different orders of the MF-DFA method to the intron-rich HUMHBB sequence. The results for orders larger than $`5`$ are indistinguishable in the scale of the figure, so the data for $`m=5()`$ yields the trendless spectrum. Analysis of this figure gives valuable information on the nature of the correlations of the DNA walk. In particular, the linearity of the trendless spectrum indicates that the record $`s(n)`$ depicted in Fig. 1 is monofractal. Actually, this conclusion is not affected by the presence of trends in the series, since the $`q`$ dependence of $`\tau (q)`$ is linear regardless of the value of $`m`$. However, the value of the scaling exponent $`H=h(q)q`$, given by the slope of the straight lines that fit the data, is sensitive to local trends. For example, we find $`H=0.72\pm 0.02`$ for $`m=0`$, $`H=0.68\pm 0.02`$ for $`m=1`$, and $`H=0.60\pm 0.02`$ for $`m=5`$. As expected, $`H=0.50\pm 0.02`$ for the shuffled sequence. All these numerical values are in perfect agreement with the values obtained using the WTMM method Arneodo95 . Since the trendless scaling exponent is definitively different from $`1/2`$ we can conclude for the existence of long-range correlations in this particular intron-containing sequence. The long-range correlation induced by the local trends is reflected in the larger value of the scaling exponent calculated with the MF-DFA$`0`$, which overestimates the value of $`H`$ in about $`20\%`$. It is interesting to note that the trends influencing the large fluctuations ($`q>0`$) are practically unaffected by the application of DFA$`1`$.

The results of the application of the MF-DFA method to the characterization of a sequence composed predominantly of coding regions is summarized in Fig. 4. The sequence considered is a portion of the E. Coli K12 genome (GenBank name ECO110K) for which $`N=111401`$. We find $`H=0.60\pm 0.02`$ for $`m=0`$, and $`H=0.51\pm 0.02`$ for $`m1`$ as well as for the shuffled sequence. As before, the monofractal structure of the record is confirmed by the linearity of the $`\tau (q)`$ spectrum. Moreover, failure to eliminate the local trends leads to a wrong conclusion about the existence of long-range correlations in this DNA walk, as evidenced by the overestimate of $`H`$ that results from the application of MF-DFA$`0`$. The local trends in this sequence seem to have a particularly simple linear nature since application of MF-DFA$`1`$ is sufficient to eliminate them altogether.

To conclude this section we note that we have considered several alternative one-dimensional representations, such as the strong-weak bond classification where C or G are associated to one step up ($`y_i=1`$) and A or T to one step down ($`y_i=1`$). All the alternative schemes analysed have produced within the numerical precision exactly the same exponents of the purine-pyrimidine classification. Our results could also be obtained within the WTMM framework and, in fact, the wavelet analysis of the HUMHBB sequence can be found in ref. Arneodo95 . The aim of this section was to give additional support to the claim of Kantelhardt et al. that, despite the operational and conceptual simplicity of the MF-DFA method, it yields results that are identical to those derived by the sophisticated WTMM method Kantelhardt02 . For instance, this easiness of implementation has allowed us to study higher-dimensional DNA walks, an awkward task to be carried out using the WTMM method, for several representations proposed in the literature (see, e.g. Luo98 ; Berthelsen92 ; Glazier95 ; Oiwa02 ). As we will show in the sequel, the results were essentially the same as those reported here for the one-dimensional DNA walks.

## III DNA trails

As pointed out before, a DNA trail is the projection of a two-dimensional DNA walk into the space plane. Fig. 5 shows the time record of the walk that gives rise to the trail illustrated in Fig. 2 for the HUMHBB sequence. We begin by noting that the generalization of the MF-DFA method to two or higher-dimensional records is straightforward as one needs only to change the basic Eq. (4), which is rewritten as

$$\widehat{F}^2(l,\nu )=\frac{1}{l}\underset{n=1}{\overset{l}{}}\left\{𝐬\left[\left(\nu 1\right)l+n\right]𝐬_\nu (n)\right\}^2$$
(8)
where $`𝐬(n)=x(n)𝐢+y(n)𝐣`$ and similarly for the fitting polynomial vector $`𝐬_\nu `$. Over an interval of $`l`$ steps the position vector in the trail will vary by typically $`\delta 𝐫F_2(l)=l^H`$ while the ‘mass’ $`M`$ of the trail, i.e., the number of points generated by these steps is $`Ml`$, provided the contribution of overlapping points is negligible. Hence $`M\delta 𝐫^{1/H}`$, from where one concludes that $`D_t=1/H`$. This relation is valid only if $`1/H<d`$ where $`d`$ is the space dimension. When $`1/H>d`$ overlap cannot be neglected and so $`D_t=d`$ (see Mandelbrot82 ; Voss89  for mathematical and pictorial details of this argument). Thus knowledge of the Hurst exponent $`H`$ for the $`d`$-dimensional DNA walks determines uniquely the fractal dimension of the DNA trail. Calculation of $`H`$ or, more generally, of $`h(q)`$ follows the same procedure sketched before for the one-dimensional case, except that the variance given in Eq. (4) is now replaced by $`\widehat{F}^2`$.

Application of the MF-DFA method to two-dimensional DNA walks representing the HUMHBB (see Fig. 5) and the ECO110K sequences yields, within the numerical precision, the same results as for the one-dimensional case. Actually, so far as the value of the Hurst exponent $`H`$ is concerned this agreement is expected since using Eqs. (8) and (5) for $`q=2`$ we have simply $`F_2l^{H_x}+l^{H_y}`$, where $`H_x`$ and $`H_y`$ are the scaling exponents of the one-dimensional DNA walks in the plane $`(n,x)`$ and $`(n,y)`$. Hence for large $`l`$ one has $`H=\text{max}(H_x,H_y)`$. In Table 1 we present the estimate of the fractal dimension $`D_t`$ of the HUMHBB and ECO110K trails, as well as of a random trail characterized by the same length and base frequency as the ECO110K sequence. The errors in the estimates of $`D_t`$ are statistical errors. For the box-counting method the systematic errors are probably much larger than those shown in Table 1, due mainly to finite-length effects Berthelsen94 .

Interestingly, although the $`D_t`$ estimates of the MF-DFA$`m`$ method are based on the analysis of finite-length records, this has practically no effect on the results, as evidenced by the correct calculation of $`D_t`$ for the random sequence, a test the box-counting method clearly fails. Moreover, these data corroborate the intuition, stemming from the visual inspection of Fig. 2, that the effect of the local trends is to decrease the fractal dimension of the trail. The excellent agreement between the box-counting estimate, which is affected by both local trends and finite-length effects, and the MF-DFA$`0`$ estimate, which is affected by local trends only, indicates that the surface effect is negligible for the HUMHBB sequence. In fact, this is expected since the box-counting method is very reliable to calculate the fractal dimensions of structures that do indeed exhibit fractality, predictably failing in the case of space-filling structures. The results for the ECO110K sequence illustrate the bad performance of box-counting when the two sources of systematic errors are present. The aforementioned agreement between the box-counting and the MF-DFA$`0`$ methods demonstrates inequivocally that the spurious effects of local trends must be filtered out for a correct account of the fractal or multifractal properties of DNA trails, a procedure that has been ignored in all previous analyses presented in the literature Luo98 ; Gates86 ; Berthelsen92 ; Berthelsen94 ; Glazier95 ; Oiwa02 .

Since the MF-DFA analysis has indicated that even the higher-dimensional DNA walks are monofractals, one concludes that the variance $`\widehat{F}^2(l,\nu )`$ defined in Eq. (8) is the same for all segments $`\nu `$ of the record. As a result, the variation of the position vector $`\delta 𝐫`$ in the trail will independ of the segment $`\nu `$ and so the trail must have a monofractal structure too. This conclusion seems at odds with those of multifractal analyses based on the sandbox and box-counting methods Berthelsen94 ; Glazier95 ; Oiwa02 . However, we note that the original claim of ref. Berthelsen94  was that finite-length DNA trails, as well as randomly generated trails, are characterized by effective multifractal spectra. Clearly, these spectra are artifacts of the sandbox and box-counting methods which cannot deal adequately with the surface effects. The unjustified dropping of the adjective ‘effective’ Glazier95 ; Oiwa02 , led then to the apparent disagreement between the two approaches. These effective multifractal spectra might indeed be an useful tool for comparing random or biological sequences of similar lengths Glazier95 , especially if one manages to free them from the local trends effects.

## IV Conclusion

The main purpose of this contribution was to point out that the statistical properties of two apparently distinct and pictorially appealing representations of DNA sequences, namely, DNA walks (Figs. 1 and 5) and DNA trails (Fig. 2), are in fact closely related. In particular, we argued that the patchiness of the DNA sequences, which has greatly hindered the characterization of the long-range correlations of DNA walks Peng94 ; Karlin93 , also plagues the fractal and multifractal analyses of the DNA trails, a fact that has been largely overlooked in previous investigations Berthelsen94 ; Glazier95 ; Oiwa02 . Using the recently proposed multifractal detrended fluctuation analysis (MF-DFA) method Kantelhardt02  to filter the local trends of DNA walks in $`d=1`$ and $`2`$ dimensions, we were able to show that these walks are monofractals. This conclusion holds true for the non-detrended walks as well, since the multifractal spectrum $`\tau (q)`$ is linear regardless of the order $`m0`$ of the fitting polynomials used to eliminate the local trends (Figs. 3 and 4). More importantly, the independence of the variance $`F^2(l,\nu )`$ \[see Eqs. (4) and (8)\] on the segment $`\nu `$ of the DNA walk carries over to the DNA trail which is thus monofractal too. In particular, the fractal dimension of the trail is $`D_t=1/H`$ where $`H`$ is the Hurst exponent of the record Mandelbrot82 ; Voss89 . In the case $`1/H>d`$, as for the one-dimensional walks, one has $`D_t=d`$. Although we have presented data for the HUMHBB and the ECO110K sequences only, which probably are, respectively, the most popular examples of intron-rich and intron-poor sequences, we have verified that our main results hold true for many other sequences in the GenBank.

Once one has discarded the simple local variations in nucleotide composition along DNA sequences as the cause of the long-range correlations observed in intron-rich sequences (though these variations do cause the spurious long-range correlations in intron-poor sequences), it is natural to ask then what are the sources for these correlations, which, as shown above, are also responsible for the nontrivial (fractal) geometry of the DNA trails. The answer is that neither the internal structure of patches nor their order in the sequence are relevant: it is the power-law distribution of patch lengths that determines the true long-range correlations Galvan96 . This as well as several other findings concerning the statistical properties of DNA sequences have prompted the proposal of minimal evolutionary models, based on biologically motivated mechanisms, to account for those features (see, e.g., Li91 ; Oliveira99 ; Almirantis01 ). In fact, rather than providing measures to characterize or distinguish classes of sequences, we think the thrust of the research on large-scale statistical properties of genomes is to provide quantitative standards for modeling the inherently stochastic process of molecular evolution Kimura71 , bearing thus on fundamental issues such as the origin of life and the evolution of complexity Bonner88 .

###### Acknowledgements.

This research was supported by Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP), Proj. No. 99/09644-9. The work of J.F.F. was supported in part by Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq).


text:
# 1 Introduction

## 1 Introduction

Gamma–ray bursts(GRBs), which are still mysterious, have very complex temporal structure. Their temporal profiles are enormously varied — no two bursts have ever been found to have exactly the same temporal and spectral development. The temporal activity is suggestive of a stochastic process (Nemiroff et al. 1993). The diversity of the bursts seems to be due to random realization of the same process that is self-similar over the whole range of timescale. Attempts to quantify these structures have not been successful (e.g., Fishman 1999).

Most of the observed profiles of GRBs are composed of pulses, each comprising a fast rise and an exponential decay (a FRED; e. g., Desai 1981, Fishman et al. 1994). Many methods for pulse analysis have been developed, e. g., the parametric analysis in model fitting (Nemiroff et al. 1993, Norris et al. 1996), the auto-correlation method (Fenimore et al. 1995), the nonparametric method (Li & Fenimore 1996), the peak alignment and normalized flux averaging method (Mitrofanov et al. 1996, Mitrofanov et al. 1998, Ramirez–Ruiz & Fenimore 1999, Ramirez–Ruiz & Fenimore 2000), and the pulse decomposition analysis method (Lee et al. 2000), etc. These statistical studies have revealed many observed temporal signatures of pulses. The pulses are hypothesized to have the same shape at all energies, differing only by scale factors in time and amplitude (“pulse scale conjecture”). And, the pulses are hypothesized to start at the same time, independent of energy (“pulse start conjecture”). The two conjectures were confirmed by Nemiroff (2000). In general, higher energy channels show shorter temporal scale factors (e.g., Norris et al. 1996., Nemiroff 2000). It is found that the temporal scale factors between a pulse measured at different energies are related to that energy by a power law, possibly indicating a simple relativistic mechanism is at work (Fenimore et al. 1995, Norris et al. 1996., Nemiroff 2000).

The statistical pulse shape has been well studied by the peak alignment and normalized flux averaging method. The peak aligned averaging pulse is spiky. A succinct pulse model, which well describes many pulse shapes, was proposed by Norris et al. (1996):

$$I\left(t\right)=I_0e^{\left(\left|tt_{\mathrm{max}}\right|/t_{r,d}\right)^\upsilon }$$
(1)
where $`t_{max}`$ is the time of the pulse’s maximum intensity ($`I_0`$); $`t_r`$ and $`t_d`$ are the rise and decay time constants, respectively; and $`\upsilon `$ is a measure of the pulse sharpness, which was referred to “peakness” by Norris et al. (1996).

However, both the duration and total count of pulses vary significantly. Statistical properties of the pulses revealed by the peak alignment and flux–normalized averaging method are limited. In this paper, this method is developed to study the shape of pulses in a more detailed manner.

To reach a result of high quality, we concern in this paper only the highest pulse of bursts, where one finds the highest level of signal–to–noise. We make the noise decomposition for the time profile of bursts by performing the wavelet analysis (which is described in section 2), then modify the pulse–finding algorithm proposed by Li & Fenimore (Li & Fenimore 1996) to identify the highest pulse in a burst profile (see section 3). In sections 4 and 5, we employ and develop the pulse aligned method to study the flux–normalized aligned averaging pulse shape and the count–and–duration–normalized aligned averaging pulse shape, respectively. Conclusions and discussion are presented in section 6.

## 2 Data Analysis

The data used for analyzing is the 64 ms temporal resolution and four–channel spectral resolution GRB data observed by BATSE from 1991 April 21 to 1999 January 26. There are 1738 bursts included. It is a concatenation of three standard BATSE data types, DISCLA, PREB, and DISCSC. All these data types are derived from the on-board data stream of BATSE’s eight Large Area Detectors (LADs). There are four energy channels observed, with the following approximate channel boundaries: 25-55 keV, 55-110 keV, 110-320 keV, and $`>`$320 keV. The DISCLA data are a continuous stream with 1.024 second resolution. They are independent of burst occurrence and taken as the background. The PREB data cover the interval 2.048 second just prior to a burst trigger.

We make the noise decomposition for the time profiles by the wavelet package analysis technique. The technique is suitable to treat those signals which cannot be analyzed by the traditional Fourier method. It was successful in de–noising the original signal and identifying the structure within a burst (e.g., Hurley et al. 1997, Quilligan et al. 1999, Lee et al. 2000). We use DB3 wavelet to make the first–class decomposition with the MATLAB software. The profile is decomposed into the signal component and the noise component. Figure 1 illustrates an example of the decomposition.

The method of the background treatment used here is similar to that in Li & Fenimore (1996). Since the DISCLA data are a continuous stream prior to and independent of the burst occurrence, they are always taken as the background of bursts. The data of the background is obtained by a linear fitting to the DISCLA data.

## 3 Pulse–Finding Algorithm and Sample Selection

Many burst time profiles appear to be composed of a series of overlapping pulses, mingling with noises. It is not easy to determine their actual light curves and to distinguish a pulse from the time profile. The result of pulse analysis strongly relies on the algorithm of pulse–finding and sample selection. Several pulse–finding algorithms have been proposed (e.g., , Norris et al. 1996,Mitrofanov et al. 1998). Li & Fenimore (1996) suggested an efficacious algorithm to identify “true peaks” in a profile. A “true peak” is not necessarily to be regarded as a pulse. If a profile is composed of only one “true peak”, the “true peak” can be regarded as a pulse. However, most of the profiles are composed of many overlapping “true peaks”. To identify a pulse in such situation is not easy. Norris et al. (1996) introduced a definition of “inseparable pulse”. We adopt this concept and regard an “inseparable pulse” as a true pulse. Our pulse–finding algorithm is described in the following.

(1) The peak–finding criterion proposed by Li & Fenimore is $`C_pC_{1,2}N_{var}\sqrt{C_p}`$ , where $`C_i`$ ($`i=1,2`$) is the photon count at time bin $`t_i`$, $`C_p`$ (at $`t_p`$) is the maximum count of a candidate peak, and $`N_{var}`$ is an adjustable parameter, typically $`3N_{var}5`$. This criterion strongly relies on absolute photon count of the candidate peak. We follow Norris et al. (1996) to employ the concept of “inseparable pulse”, and then the pulse–finding criterion becomes $`1C_{1,2}/C_p0.5`$. It means that a candidate peak is a true peak only when $`C_1`$ (at $`t_1`$) and $`C_2`$ (at $`t_2`$) are lower than the half of the $`C_p`$ on both sides of $`t_p`$. With this method, one might find more than one true peaks within a burst. For the reason mentioned in section 1, we select only the highest one.

(2) In order to maintain a high level of signal–to–noise, we adopt the intensity criterion as $`C_{max}>10\sigma `$, where $`\sigma `$ is the standard deviation of the background.

(3) Only those pulses with at least 10 bins of time are selected. Those with less bins do not provide enough structure information and thus are ignored.

We apply the above pulse–finding algorithm to select highest pulses in the profiles of bursts. There are 760, 885, 885, and 334 bursts, for which the highest pulses can be identified, in channels 1 to 4, respectively. The number of bursts for which all of the highest pulses in four channels can be identified is 275. The flux in this sample ranges widely, from 0.513 photons.cm-1.s-1 to 183.370 photons.cm-1.s-1. We select the sample to study the statistical properties of the pulse morphology.

## 4 The Flux–Normalized aligned Averaging Pulse

It was found that the peak–aligned averaging pulse well illustrates how the average pulse evolves with energy. To do that, one averages the time profile of individual events by the normalized peak–alignment technique, where each time profile is normalized by the peak number of counts Cmax, aligned at the peak time bin tmax, and then averaged for all bins along the timescale(e. g., Norris et al. 1996, Mitrofanov et al. 1996, Ramirez–Ruiz & Fenimore 2000). Though the pulse–finding method and the sample adopted in this paper are somewhat different from the previous ones, the peak–aligned averaging pulses we obtain are quite similar to that in Norris et al. (1996) and Ramirez-Ruiz & Fenimore (2000) (see Fig. 2). Figure 2 shows the same result that the higher the energy, the narrower the pulse. One can also find this from the flux–normalized–and–beginning–aligned averaging pulses shown in Fig. 3.

We find that the shapes of pulses in Fig. 2 are quite different from that in Fig. 3. Though both figures come from the sum of the same normalized pulses, but the ways of the alignment are different. The difference between Figs. 2 and 3 must come from the diversity of the duration and the asymmetry of the shape of the normalized pulses. Obviously, the peak–aligned method would lead to a spiky shape. The statistical result must conceal most of the diversity of the duration and the asymmetry of the pulses. This explains why the averaging pulses in Fig. 2 are very different from observation. Also, the flux–normalized–beginning–aligned method does not take into account the diversity of the duration of pulses. Thus, neither Fig. 2 or 3 truly embodies the temporal structure of pulses.

## 5 The Count-and-Duration-Normalized aligned Averaging Pulse

We observe that, the above peak–aligned method allows those pulses with longer durations to possess larger total counts and then to contribute more significantly to the averaging pulse. To investigate the averaged shape of pulses, one needs to get rid of the effects from both the total count and the duration of the selected pulses. This leads to the count-and-duration-normalized aligned averaging method employed bellow.

First, both the total count and the duration for each selected pulse are normalized. Then the averaging pulses for the four channels are obtained with the same way used in last section. The results are shown in Figs. 4 and 5. They are different respectively from that in Figs. 2 and 3.

Figure 4 shows that the count–and–duration–normalized–peak–aligned averaging pulses in the four channels are almost the same. One can not tell any difference between the count–and–duration–normalized–beginning–aligned averaging pulses of the four channels. It shows that the shapes of these pulses are independent of energy bands. The pulses of the four channels have the same temporal structure.

The pulse model in Norris et al. (1996) is employed to fit them, and it shows that these pulses can be well described by the model. The fitting parameters are: $`t_r`$=0.28, $`t_d`$=0.38, $`\upsilon =2.29`$ for the count–and–duration–normalized–peak–aligned averaging pulse (in channel 3); $`t_r`$=0.12, $`t_d`$=0.16, $`\upsilon =1.09`$ for the count–and–duration–normalized–beginning–aligned averaging pulse (in channel 3). The ratios of $`t_r`$ to $`t_d`$ for both pulses are about 0.75.

## 6 Conclusions and Discussion

In this paper we study the flux–normalized–peak(beginning)–aligned averaging pulses and the count-and-duration-normalized–peak(beginning)–aligned averaging pulses with the highest pulse in the profiles of GRBs. We apply the wavelet package analysis technique and a developed pulse–finding algorithm to select the highest pulse. The sample we get includes 275 bursts which fluxes range from 0.513 erg.cm-1.s-1 to 183.370 photons.cm-1.s-1.

The wavelet package analysis technique is suitable to treat those signals which cannot be analyzed by the traditional Fourier method. It is successful in de–noising the original signal and identifying the structure within a burst. The pulse–finding algorithm used in Li & Fenimore (1996) is developed in this paper so that the selection of pulses depends on the relative value of counts rather than the absolute value.

The number of bursts concerned in this paper is the largest one of that used for pulse analysis so far, and the sample adopted here covers the biggest flux range. Many samples of very bright bursts have been employed to study statistical properties of pulses (e.g., Norris et al. 1996). Though the numbers of bursts concerned are much smaller, the authors were able to get more pulses by selecting not only the highest pulse but also other pulses in a burst (e.g., Norris et al. 1996). Their results refer only to bright bursts, and the pulses so selected might include possible evolutionary effects of pulses.

Figures 2 and 3 support the well-known conclusion that the higher the energy, the narrower the pulse.

Different from the previous works, we make in this paper not only the normalization of the total count but also the normalization of the duration. In this way, all the pulses (strong or weak) contribute equally to the averaging pulses. And, in this way, the averaging pulses obtained stand only for the statistical shape of pulses. The effects from both the duration and the total count are removed. The results so obtained are quite different from that got by the previous method. One can find this by comparing Figs. 2 and 4.

For the pulses shown in Figs. 4 and 5, we prefer those in the latter. Since the pulses in Fig. 4 come from aligning the normalized pulses at the moment of maximum count of the pulses, those asymmetric normalized pulses would contribute differently to different sides of the averaging pulses. Thus the distribution of the peak count position in the shape of selected pulses must be at work. It would lead to a spiky shape of the averaging pulses. Fig. 4 must conceal most of the diversities of the duration and the asymmetry of pulses. Differently, the pulses shown in Fig. 5 stand only for the average shape of the original pulses. These pulses are well described by the model in Norris et al. (1996). The fitting parameters for the pulse in channel 3 in Fig. 5 are: $`t_r`$=0.12, $`t_d`$=0.16, $`\upsilon =1.09`$. The ratios of $`t_r`$ to $`t_d`$ for the pulse is 0.75.

We find that the count-and-duration-normalized–peak(beginning)–aligned averaging pulses are the same for different channels. Our results strongly support the previous conclusion that the temporal profiles in different channels are self–similar. The averaging pulse shape is independent of energy bands. Due to these results, we believe that the mechanisms of bursts in different gamma-ray bands must be the same.

The mechanism generating the bursts is still unknown. Many models for interpreting the origin and emission of the event have been proposed (e. g., Rees & M$`\stackrel{´}{e}`$sz$`\stackrel{´}{a}`$ros 1992, Vietri et al. 1998, Fuller et al. 1998, Dai & Lu 1998, Daigne et al. 1998, , etc.), mostly in the context of two major scenarios involving relativistic shells. An approach frequently used in these models is to identify each pulse in the light curve with a single event. Depending on the model chosen, this event could be the collision between inhomogeneities in a relativistic wind in the internal models or the “activation” of a region on a single external shell. Our study shows that, the shock, either an internal or an external one, producing the pulse, might produce photons over the four energy channels in the same way.


text:
# Long-range interactions in Sznajd consensus model

e-mail: ab127@uni-koeln.de

Abstract: The traditional Sznajd model, as well as its Ochrombel simplification, for opinion spreading are modified to have a convincing strength proportional to a negative power of the spatial distance. We find the usual phase transition in the full Sznajd model, but not in the Ochrombel simplification. We also mix the two rules, which favours a phase transition.

Keywords: Sociophysics, phase transition, distance dependence, quenched disorder

PACS: 05.50 +q, 89.65 -s

Ising models have been studied since nearly one century, and simulated on computers since more than four decades. A new version is the Sznajd model  where again each lattice site carries a spin $`\pm 1`$. If two randomly selected neighbouring spins have the same value, they force their neighbours to accept this value; otherwise nothing is changed and a new pair is selected. In the Ochrombel simplification, instead of a pair, already a single site “convinces ” its neighbours . This model can be interpreted as the spreading of opinions until a consensus is reached. Instead of two values $`\pm 1`$ we also can work with $`q`$ values: 1,2, …, $`q`$.

The Sznajd model on the square lattice shows a phase transition: If initially one of the two opinions has a slight majority in a random distribution, then at the end all spins have that value and the dynamics stops. The Ochrombel modification lost this transition . We now check for this transition in the case of long-range interactions, decaying with a power law of the distance, and with a mixture of Sznajd and Ochrombel rule. The program is similar to the published one .

We made 100 or 1000 simulations for $`L\times L`$ square lattices with $`L=7,13,26`$ and 53, sometimes also 73, usually allowing $`q=5`$ values. A spin convinces, alone (Ochrombel) or together with an equally-minded neighbour (Sznajd), a neighbour at Euclidean distance $`R`$ with probability $`1/R^{2x}`$. Initially the spins are distributed randomly among the $`q`$ opinions except that with a bias probability $`p`$ the just initialized spins are set to +1. A quenched fraction $`r`$ of the sites follows the Sznajd pair rule, the remaining fraction $`1r`$ the Ochrombel single-site rule. A success is a sample where at the end all spins had the bias value +1.

Figs.1 and 2 show for the Sznajd case $`r=1`$ the phase transition: For large $`L`$ a small bias $`p`$ suffices to make nearly all samples successes. It does not matter much whether the interactions decay slowly ($`x=1/2`$) or fast ($`x=2`$) with distance. For the $`r=0`$ Ochrombel case, however, analogous simulations (not shown) give no phase transition, and this situation persists even if we take a very small $`x=0.1`$ (Fig.3) and reduce $`q`$ from 5 to 2 (Fig.4).

Thus we mixed the two rules in Figs. 5 ($`r=0.5`$) and 6 ($`r=0.1`$) which show a phase transition, for $`x=1/2`$, in both cases. With a faster decay, $`x=2`$ instead of $`x=1/2`$, the phase transition for $`r=0.1`$ becomes more pronounced, Fig.7.

In summary, contrary to our expectation from thermal phase transitions, the introduction of long-range interactions instead of nearest-neighbour interactions did not create a phase transition.

We thank Deutsche Forschungsgemeinschaft for support, and D. Stauffer, to whom this note is dedicated because of his senility, for help.


text:
# Iron Emission Lines from Extended X-ray Jets in SS 433: Reheating of Atomic Nuclei

> Powerful relativistic jets are among the most ubiquitous and energetic observational consequences of accretion around supermassive black holes in active galactic nuclei and neutron stars and stellar-mass black holes in x-ray binary (XRB) systems. But despite more than three decades of study, the structure and composition of these jets remain unknown. Here we present spatially resolved x-ray spectroscopy of arc second-scale x-ray jets from the XRB SS 433 analized with the Chandra advanced charge-coupled device imaging spectrometer. These observations reveal evidence for a hot continuum and Doppler-shifted iron emission lines from spatially resolved regions. Apparently, in situ reheating of the baryonic component of the jets takes place in a flow that moves with relativistic bulk velocity even more than 100 days after launch from the binary core.

The jets of XRB SS 433 are well established, quasi-persistent and the best example of regular precession of the jet axis (1–3). Extensive observing campaigns have placed constraints on the emission physics, geometry and scales associated with the baryonic component of the jets of SS 433. The generally accepted model for the thermal component of the jet is an adiabatic cooling model (4, 5). This model assumes that the emitting matter is moving along ballistic trajectories and expanding adiabatically in a conical outflow, and that the temperature and matter density are only a function of the distance from the core. The matter is ejected (the physical processes of ejection are still uncertain) with a temperature of $`>10^8`$ K (mean thermal energy per particle kT$`>10`$ keV) (4–6), comparable to that measured for the inner regions of accretion disks. Close to the binary core, at distances less than $`10^{11}`$ cm, this jet plasma emits high-temperature continuum and Doppler-shifted lines in the soft x-ray band (6). The jet subsequently emits Doppler-shifted optical emission lines that are constrained to originate at distances of less than $`10^{15}`$ cm because of a lack of measurable offsets between blue- and red-shifted beams. Neither line-emitting region has been spatially resolved from the core of SS 433. Furthermore, the temperature of the optical line-emitting region is limited to less than $`10^4`$ K by the absence of higher excitation lines such as He II and C III/N III (7). In the adiabatic model the temperature decreases with the radial distance $`R`$ from the core as temperature T$``$R-4/3; thus, at distances greater than $`10^{16}`$ cm, the matter is expected to be at a temperature of $`<100`$ K, too cool to thermally emit x-ray and optical radiation. Besides this thermal radiation, there is nonthermal synchrotron radiation, which traces the distribution of relativistic electrons and the magnetic field in the jet, observed in the radio band at distances between $`10^{15}`$ and $`10^{17}`$ cm (2). Further out, the jet is not detected again until degree-scale (i.e., at distances of $`10^{19}`$ to $`10^{20}`$ cm from the binary core) termination shocks are observed, deforming the surrounding W50 nebula (Fig. 1, top), which is thought to be a supernova remnant associated with the formation of the neutron star or black hole in SS 433 (8, 9).

Marshall et al. (6) have recently analyzed the x-ray line emission from SS 433 with the Chandra high energy transmission grating spectrometer (HETGS). They refined the adiabatic cooling model, limiting the thermal x-ray emission region to within $`10^{11}`$ cm of the jet base. They also discovered faint arc sec-scale x-ray emission on either side of the core, but they were not able to measure any emission lines and concluded that the jet gas would be too cool at such a distance to emit thermal x-rays.

We have analysed a 9.6-ks observation from 27 June 2000 with the Chandra advanced charge-coupled device (CCD) imaging spectrometer (ACIS)-S. The reduced image (Fig. 1) resolves arc sec-scale x-ray jets approximately symmetrically placed on either side of a weaker core and aligned with both the sub-arc sec radio jets and degree-scale structures in W50. Analysis of the image indicates that photon pile-up (10) in the core artificially reduced the measured count rate there, although according to the orbital ephemeris of Dolan et al. (11), our observation was close to x-ray eclipse, phase $`0.9`$ (12), which may also have reduced the core flux. At a distance of 5 kpc, these x-ray jets correspond to physical scales of $`10^{17}`$ cm and soft x-ray (2 to 10 keV) luminosities of 3$`\times 10^{33}`$ to 4$`\times 10^{33}`$ erg s-1 each (about 3% of the average source total luminosity).

We were able to spatially isolate and take x-ray spectra of the two jets. The spectra of the two lobes between $`0.8`$ and 10 keV (Fig. 2) can be fit by a bremsstrahlung continuum with a temperature of $`5`$ keV and Gaussian emission lines at $`7.3`$ keV for the approaching (east) jet and at $`6.4`$ keV for the receding (west) jet (Table 1). In both spectra there is an indication of additional emission, possibly also lines, at higher energies, but this is near to the limit of the Chandra sensitivity, and given our relatively low count rates, we did not attempt to fit this component. A hydrogen column density of $`10^{22}`$ cm-2 was found for both lobes, consistent with the column density found by Marshall et al. (6) for the core of SS 433. A power law with a photon index of $`2.1\pm 0.2`$ fits the continuum equally well, and on the basis of our current data, we cannot distinguish between a power law and bremsstrahlung for the continuum emission.

The significantly different line centroid energies measured for the east and west jets suggest that they may originate in the moving jet material. To test this hypothesis, we developed a model that sums the emission from ten equally spaced phase bins (13) around the precession cycle, calculating for each bin the Doppler shift $`\delta `$ (14) and Doppler boosting of the line peak intensity $`\delta ^2`$ (assuming that the jets are intrinsically symmetric about the precession axis and east to west). A simulation based on this model for a 10-ks observation, assuming a rest energy for the line of 7.06 keV, corresponding to Fe XXV k$`\beta `$ and using the same normalizations for the Gaussians as found in our observation, corresponds well to our data (Fig. 3 and Table 1).

The line energies and widths we obtain in our simulation are a good match to those measured in our observation (Fig. 2 and 3). Thus, we infer that the jet is emitting around the precession cone at a large distance from the core. Our model predicts twin-peaked lines, due to the distribution of Doppler shifts as a function of time, with a stronger blue peak (due to Doppler boosting); however, this twin-peak is not seen in the 10-ks simulation, which have a relatively low signal-to-noise ratio. However, if our model is correct, then additional observations should reveal this structure. Although we cannot use our model to conclusively establish that the Doppler-shifted line is Fe XXV k$`\beta `$, we do need a transition of hot, highly ionized Fe in order to achieve a high enough energy for the line in the east jet. We cannot exclude an origin for the $`6.4`$ keV line in the west jet from neutral Fe; however, we cannot envisage a plausible geometry whereby we are observing neutral Fe at 6.4 keV in the west jet and the same line extremely blue-shifted in the east jet. A further prediction of this model is that the east jet should be stronger than the west jet as a result of the overall larger Doppler boosting. The ratio of 2 to 10 keV fluxes (line and continuum) from the east to west jets is $`1.3`$, approximately what would be expected for the model considered here. Although we cannot explicitly show that the source of the continuum emission arises in the moving jet gas, the measured bremsstrahlung temperature is consistent with an origin in the same plasma as the high-excitation Fe lines. We therefore conclude that we observed lines, and probably also continuum, from a hot ($`>10^7`$K) plasma which is moving with a similarly high velocity to that measured for the “inner” jets (6), $`0.26c`$ (where $`c3\times 10^{10}`$ cm s-1 is the speed of light), on scales of $`10^{17}`$ cm from the binary core.

This is contrary to the adiabatic cooling model (4–6) which at that distance predicts that the baryonic component will have cooled to the ambient interstellar medium temperature; therefore, the jet has to re-heat itself. The most plausible energy source is the kinetic energy ($`10^{38}`$ ergs/s) associated with the bulk motion of the jet (6). Because we observed optical emission lines coming from matter with temperatures of $`<10^4`$ K at a distance of $`10^{15}`$ cm, this re-heating must occur between $`10^{15}`$ and $`10^{17}`$ cm from the core.

Alternatives to this scenario can be considered. A power law fits the continuum well, and by analogy with the extended x-ray jets seen in active galactic nuclei (AGN) (15), we cannot exclude an inverse-Compton or synchrotron origin for the continuum component. However, because both of these emission mechanisms require a population of hot electrons, they also imply reheating at those large distances from the central source. Regardless of which form of reheating is associated with the continuum spectrum, the emission lines require a baryonic component. Reflection of the core emission does not seem possible because the line flux we observed at $`10^{17}`$ cm is only a factor of $`10`$ less than that observed to originate at $`10^{12}`$ cm (6), and yet our images indicate it does not subtend a significantly larger solid angle than the jets on smaller scales. Furthermore, the line to continuum flux ratio we observed is an order of magnitude greater than in the core \[$`F_{\mathrm{line}}/F_{\mathrm{continuum}}`$ is $`1`$ for our east jet, but is $`0.1`$ for the 7.3-keV line in Chandra HETGS observations (6)\]. The lines from these extended jets are still an order of magnitude weaker than the moving core x-ray lines (5, 6), and would not be expected to show measurable changes in their Doppler shifts (unless, e.g., a major outburst from SS 433 resulted in enhanced emission at some particular precession phase). Nevertheless, they should be present as weak features in spectra that do not spatially resolve the core and jets.

Reheating (in situ acceleration) of a leptonic (electron and possibly also positron) component of jet plasma, presumably by conversion of the jet kinetic energy, is commonly invoked to explain synchrotron or inverse Compton energy emission on large scales in the jets of AGN (15-17). The reheating of the baryonic plasma that we observed in SS 433 requires a similar tapping of the bulk kinetic energy of the flow, by processes that are able to act on atomic nuclei. In the internal shock model for blazars (16), a slightly varying jet speed results in shocks downstream in the jet flow. In the context of this model, random discrepancies of about 15% in the velocity of the jets from SS 433, inferred from long-term optical monitoring (3), could produce shocks a few hundred days downstream (for a precession period of 162 days and a mean jet velocity of $`0.26c`$), comparable to what we observed (18) .

To conclude, these observations spatially resolve a line-emitting region in a relativistic jet. They reveal that the hot x-ray emitting region is still moving relativistically on physical scales that are orders of magnitude larger than previously inferred, with no evidence for substantial deceleration, supporting models for extended x-ray jets in AGN that require bulk relativistic motion on large physical scales (15). Furthermore, this observation demonstrates that particle reacceleration in a relativistic jet, previously inferred only to act on the leptonic component, can act also on atomic nuclei. Comparison of this highly super-Eddington, mass-loaded jet with those of other x-ray binary systems and AGN can provide unique insights into the matter and energy content of relativistic jets and hence the coupling of accretion and outflow in conditions of extreme gravity, pressure, and energy density.

> References and Notes
>
> 1. G. O. Abell and B. Margon, Nature 279, 701 (1979).
> 2. R. M. Hjellming and K. J. Johnston, Nature 290, 100 (1981).
> 3. S. S. Eikenberry, et al., Astrophys. J. 561, 1027 (2001).
> 4. W. Brinkmann, N. Kawai, M. Matsuoka, H. H. Fink, Astron. Astrophys. 241, 112 (1991).
> 5. T. Kotani, N. Kawai, M. Matsuoka, W. Brinkmann, Publ. Astron. Soc. Pac. 48, 619 (1996).
> 6. H. L. Marshall, C. R. Canizares, N. S. Schulz, Astrophys. J. 564, 941 (2002).
> 7. J. Shaham, Vistas Astron. 25, 217 (1981).
> 8. W. Brinkmann, B. Aschenbach, N. Kawai, Astron. Astrophys. 312, 306 (1996).
> 9. G. M. Dubner, M. Holdaway, W. M. Goss, I. F. Mirabel, Astron. J. 116, 1842 (1998).
> 10. Pile-up is an effect by which two or more photons that hit the same pixel on the CCD in one integration are counted as a single photon. If the summed energy is greater than the ACIS threhold ($`15`$ keV), the detection is rejected, and no counts are recorded.
> 11. J. F. Dolan, et al., Astrophys. J. 327, 648 (1997).
> 12. The orbital phase $`0\varphi 1`$ indicates the fraction of one entire orbit completed. Eclipse of the x-ray source (i.e., superior conjunction of the neutron star or black hole) corresponds to $`\varphi =0`$. The orbital period is 13.08 days.
> 13. The precession phase $`0\psi 1`$ indicates the fraction of one entire jet precession cycle. $`\psi =0`$ is, by convention (2), one of the two instances when the beams of SS 433 are perpendicular to the line of sight. The precessional period is 162.4 days.
> 14. The Doppler factor $`\delta =[\gamma (1\beta \mathrm{cos}\theta )]^1`$ where the velocity $`v=\beta c`$, $`c`$ is the speed of light, $`\gamma =(1\beta ^2)^{1/2}`$ is the Lorentz factor, and $`\theta `$ is the angle to the line of sight.
> 15. D. E. Harris and H. Krawczynski, Astrophys. J. 565, 244 (2002).
> 16. M. Spada, G. Ghisellini, D. Lazzati, A. Celotti, Mon. Not. R. Astron. Soc. 325, 1559 (2001).
> 17. G. V. Bicknell and D. B. Melrose, Astrophys. J. 262, 511 (1982).
> 18. The relationship $`\mathrm{\Delta }t=P(1\delta )/2\delta `$, where $`P`$ is the precession period and $`\delta `$ is the fractional deviation from the mean velocity $`v_m`$ of the ejected blobs, gives the time interval $`\mathrm{\Delta }t`$ over which a blob with velocity $`(1+\delta )v_m`$ reaches the blob launched with velocity $`(1\delta )v_m`$ at the same precessional phase one cycle earlier. Using $`P=162`$ days and $`\delta =0.15`$ we obtain $`t450`$ days.
> 19. Fig. 1 (top) is reproduced by permission of the American Astronomical Society (Fig. 1 in (9)).
> 20. We thank Dan Harris and Michiel van der Klis for comments on a draft of this manuscript, and Wolfgang Brinkmann and Mike Watson for earlier useful discussions. MM and RF are grateful to Max-Planck-Institut für Astrophysik for their hospitality during the initial analysis of these data.