text:
# Gompertz law in simple computer model of ageing of biological population

## 1 Introduction

Most computer simulations on biological ageing use the Penna model . This model, motivated by the accumulation mutation hypothesis, provides a population of bit-strings that exhibits many features known for real senescence . Recently Stauffer  in a general review of the Penna model has suggested a simpler alternative with less parameters. This new attempt is based on the postulate of a minimum reproduction age and a maximal genetic lifespan. Only these two numbers are transmitted from generation to generation, with certain mutations, by asexual reproduction. This simple model shows the basic features required for senescence. The age distribution shows an increase of mortality with age. The model reproduces successfully the catastrophic senescence of salmon  and some specific human aspects like the social structures of technological societies  or a transition between extinction and survival if a minimal population for survival is required .

In the following we examine the other modification to the basic model. The privilege conditions are introduced to earlier and then to late reproduction. By doing this, in case of late reproduction, we achieve the exponential growth of mortality in agreement with the Gompertz law. The increase of fertility with age exists in nature. Some fish, lobsters and trees share this feature  and therefore this life strategy has been considered by the Penna model also .

## 2 The model

The population consists of $`N`$ individuals $`i`$ $`(i1,\mathrm{},N)`$ initially. Each individual is characterized by three integers: its age $`a(i)`$, its minimum reproduction age $`a_m(i)`$ and its maximal genetic lifespan $`a_d(i)`$ with $`0a_m(i)<a_d(i)32`$. The maximal lifetime is restricted to 32 time units (called years), the minimum reproduction age may be chosen between zero and $`a_d(i)1`$. Within these constraints the values of $`a_m(i)`$ and $`a_d(i)`$ are randomly mutated for an offspring by $`\pm 1`$, away from the maternal values, and for each child separately. After an $`i`$th individual has reached its minimum reproduction age, it gives birth to one offspring with probability $`b`$, chosen as

$$b=\frac{1}{a_d(i)a_m(i)}=\frac{1}{\mathrm{\Delta }_i}$$
(1)
Thus during its whole mature life any individual gives one offspring at the cumulated probability equal to 1.

Independently of the genetic death, which happens automatically and unavoidably if $`a(i)=a_d(i)`$, at each time interval an individual can also die “accidentally”, with the Verhulst probability $`N/N_{max}`$. $`N_{max}`$ is called the carrying capacity to account for the fact that any given environment can only support populations up to some maximal size $`N_{max}`$. Otherwise the individuals die because of food and space limitations.

In this paper differently than in the standard model we consider populations in which the cumulated birth rate is not uniformly distributed among years of reproduction: $`a_m,a_m+1,\mathrm{},a_d1`$. Instead we propose to consider:

1. privilege to younger third part of a reproduction period:

$$b^{young}=\{\begin{array}{cc}2b\hfill & \text{ if }a_m(i)a(i)a_m(i)+\frac{1}{3}\mathrm{\Delta }_i\hfill \\ \frac{1}{2}b\hfill & \text{ if }a_m(i)+\frac{1}{3}\mathrm{\Delta }_i<a(i)<a_d(i)\hfill \end{array}$$
(2)
2. privilege to older third part of a reproduction period:

$$b^{old}=\{\begin{array}{cc}\frac{1}{2}b\hfill & \text{ if }a_m(i)a(i)a_m(i)+\frac{2}{3}\mathrm{\Delta }_i\hfill \\ 2b\hfill & \text{ if }a_m(i)+\frac{2}{3}\mathrm{\Delta }_i<a(i)<a_d(i)\hfill \end{array}$$
(3)

In the following figures we present characteristics of populations developed when three ways of reproduction are considered: standard, (A) and (B) models. To avoid possible divergence we follow the Stauffer recipe  and shift the birth rate accordingly:

$$b=\frac{1+ϵ}{\mathrm{\Delta }_i+ϵ}$$
(4)
with $`ϵ=0.08`$. In figure 1 we show the population volumes in percentage of the $`N_{max}`$ capacity. The population where the reproduction of older individuals is privileged, i.e., model (B), is the smallest one and varies in time by $`33\%`$.

Figures 2, 3 compare the distribution of $`a_m`$ the minimal reproduction age and $`a_d`$ the maximal genetic lifespan, respectively, in the three models considered. A significant change in the length of life of individuals is observed when the reproduction at old age is favorable, see figure 4 for the distribution of age $`n(a)`$ in a stationary population.

Finally figure 5 shows the mortality in the populations simulated. The mortality is calculated according to the formula $`q(a)=1n(a+1)/n(a)`$ and plotted on log-scale to extract the Gompertz law.

## 3 Summary

In Stauffer’s new proposition of a model for biological ageing, by considering the privilege conditions for time of reproduction we can easily manipulate the mortality in a population. When reproduction at old age is preferred then the model provides the mortality in agreement with the Gompertz law.

In the simple ageing model discussed here the maximal genetic lifespan for each individual is limited to 32 years. This limitation was introduced to the model to make results better comparable with many Penna model results. But one could ask if discarding this assumption would help in getting a nice Gompertz law within the standard model. This useful suggestion pointed out to us by the referee has motivated us for investigations on the simple ageing model. Specially we were interested in where the model (B) would drive a population. Preliminary studies show that this model simulated without the limitation for the maximal genetic lifespan leads to the stationary state such that the distribution of $`a_d`$ is concentrated around age of 50. The mortality in this population is presented in Figure 6. One can notice that the Gompertz law of the exponential increase of mortality holds for adult individuals at age lower than about 35.

Aknowledgement
D. M. thanks KBN for financial support: Project PB0273/PO3/99/16.

List of figures

Figure 1. Usage of environment capacity in time in standard, (A) and (B) models.

Figure 2. Distribution of minimal reproduction age in standard, (A) and (B) models.

Figure 3. Distribution of lifespan in standard, (A) and (B) models.

Figure 4. Distribution of age in standard, (A) and (B) models.

Figure 5. Mortality in populations of standard, (A) and (B) models.

Figure 6. Mortality in populations of model (B) without the limitation to the maximal genetic lifespan.


format: text

text:
# Is the old nova RW UMi (1956) an intermediate polar?

## 1 Introduction

Classical novae (CNe) belong to the class of Cataclysmic Variables (CVs) and are formed by a white dwarf (WD) primary accreting matter from a Roche lobe filling companion, usually a lower main sequence star. The WD is surrounded by an accretion disk, unless its magnetic field is strong enough to partially (Intermediate Polar systems) or totally (Polar systems) control the accretion geometry. Stable mass transfer is ensured by angular momentum losses from the binary system due to the magnetic braking mechanism and gravitational radiation. Thus, CVs evolve towards shorter and shorter orbital periods and lower and lower masses of the secondary.

The orbital periods of CVs show a bimodal distribution characterized by a minimum period of $`P80`$ min and a gap in the range 2 h $`P_{\mathrm{orb}}`$ 3 h (see Warner 1995 for a comprehensive review). The 3 h upper limit of the gap corresponds to secondary stars which are just becoming fully convective ($`M_20.3M_{}`$) and in which the dynamo mechanism is partially suppressed leading to a drastic decrease of the efficiency of the magnetic braking and, consequently, of the accretion rate . The secondary then shrinks inside its Roche lobe until gravitational radiation will bring them again in contact, which will happen when the binary has reached $`P_{\mathrm{orb}}`$ 2 h.

The accretion luminosities of old novae are at any given orbital period systematically brighter than those of dwarf novae (DNe), suggesting that the former ones have higher mass transfer rates and/or hotter white dwarfs and thermally stable accretion discs. There seems to be no systematic difference between the pre- and the post- nova luminosities, although a few old novae appear 2–4 magnitudes brighter than in their pre-nova state. Amongst these, GQ Mus, CP Pup and V1974 Cyg have orbital periods below the period gap, while V1500 Cyg is just above it . However, while strongly magnetized systems like GQ Mus and V1500 Cyg, as ‘polars’, are expected to alternate between high and low states of their accretion rates, the explanation for the luminosity standstills of disc accreting systems like CP Pup and V1974 Cyg suggested by Retter & Naylor  is that, after the outburst, they experience a constant superhump phase.

RW UMi is a high galactic latitude, probably slow, nova that brightened in 1956 reaching $`m6`$. It was first considered as a supernova by Kukarkin  due to its unusually large outburst amplitude and because it was apparently located out of the Galactic halo. The pre-nova magnitude was estimated $`21`$  but its post-outburst magnitude still is about 2.2 mag brighter. Esenoglu et al.  estimated the nebular parallax of the nova as $`5000\pm 800`$ pc. This yields $`M_V^{\mathrm{max}}7.7`$ at light maximum, and $`M_V^{\mathrm{min}}7.3`$ or $`5.1`$ before and after the outburst, respectively. Thus the nova is located in the Galactic halo.

The first spectroscopic observation of RW UMi was obtained by Zwicky  seven years after the maximum. The spectrum shows a bluish continuum with emissions of hydrogen, HeI and HeII. Subsequent spectra obtained by Kaluzny & Chlebowski , Szkody & Howell  and Ringwald et al.  confirmed the steep blue continuum and the relatively high excitation state of the object.

The photometric behaviour of RW UMi at light minimum has been object of debate. No light variations were found by Kaluzny & Chlebowski ; Szkody et al.  suggested a periodicity of the order of 3 h, while Howell et al.  found a period of 1.88 h. The possibility that the orbital period of RW UMi falls below the 2–3 h gap of CVs  triggered more systematic observations until Retter & Lipkin  were able to unambiguosly determine a photometric period of $`0.05912\pm 0.00015`$ d with a 0.025 mag amplitude. In the following we will refer to this period as the main period $`P_{\mathrm{main}}`$, as its precise nature is not really clear, although it is probably closely related to the orbital motion.

An orbital period of the order of 85 min definitely places RW UMi well below the 2 h lower limit of the period gap. This region of the Period-Luminosity diagram of CVs is populated by SU UMa systems. Retter & Lipkin  suggested that RW UMi is a non-magnetic CV which is experiencing a long post-outburst superhump phase like two more short-orbital-period old novae: CP Pup and V1974 Cyg.

In this paper we present the results of a photometric study of the old nova aimed to identify possible sub-orbital periodicities in its light curve. We have obtained CCD images in 1999 and in 2000 using the Wyoming Infrared Observatory (WIRO) and the Ekar (Asiago) telescope, respectively. In section 2 we present the observation and the data reduction, and present the photometric periods discovered. In section 3 we discuss the general photometric behaviour of the old nova and suggest a physical explanation for the sub-orbital periods observed. The possible intermediate polar state of the post-nova is also discussed.

## 2 The observations and data reduction

The $`2.5^{\prime \prime }\times 2.5^{\prime \prime }`$ field centered on the nova is shown in Fig. 1. No object is seen at the position of the nova because the POSS image was taken before the 1956 outburst. However, the luminosity of the post-nova is higher than the pre-nova, and its image, together with the magnitudes and colours of two field stars, labelled b and c in Fig. 1, were provided by Kaluzny & Chlebowski . These field stars allowed an absolute calibration of the Ekar $`V`$ magnitudes. Stars a and b in Fig. 1 were used to obtain the mean relative fluxes of the nova from the WIRO observations. The journal of observations is given in Table 1.

### 2.1  WIRO observations

WIRO observations were secured in June 1999 using the Michigan State University visual camera attached at the prime focus of the 2.3 m telescope (WIRO) on Jelm Mountain (Laramie, WY). The space resolution was 2.3 arcsec/pixel. We observed with a time resolution of 27 s for 2.72 h on June 18, and for 2.60 h on June 19. In order to obtain a higher S/N ratio we did not use filters. We selected a small region for the CCD readout in order to improve the time resolution of our observations. The selected field therefore contained only RW UMi and field stars a and b (see Fig. 1). We then measured the mean flux ratios between the nova and the two nearby stars following standard photometric procedure. The light curves of the nova, expressed as flux ratio variations, for the two observing nights at WIRO are shown in Figs. 2 and 3. We note that, in spite of the fact that each run was about 1.9 times the main period, its modulation is not clearly seen in our data. The reason is that the light curve of the nova is intrinsically rather irregular which also explains why only the systematic and extended observations performed by Retter & Lipkin  were able to provide a precise period.

### 2.2 Ekar observations

Additional observations in $`V`$ were obtained in May 2000 at Asiago, Italy, using the 1.82 m Mt. Ekar telescope with the CCD of the AFOSC system as an imager. The time resolution was 430 s. Data reduction was performed with IRAF software packages. Field stars b and c of Fig. 1 allowed the calibration of the $`V`$ magnitudes . The light curves of the old nova from the two Ekar observing nights are shown in Figs. 4 and 5. The length of each observation is about 2.3 orbital periods, and in this case the orbital modulation is more evident than in the WIRO data.

### 2.3 The photometric periods.

As anticipated, the pathological irregularity of the light curve of this old nova makes the determination of the main period  difficult, when only two or three cycles are observed. For example, the Fourier anlysis of the WIRO data, probably due to the rather different appearance of the two data sets, would suggest a period of 0.098 d. The Ekar data, instead, show a best period at 0.06614 d, while a peak at 0.05856 d, the closest one to $`P_{\mathrm{main}}`$, is only the fifth highest alias in the power spectrum.

In order to look for the presence of higher frequency secondary modulations, we have corrected our data sets subtracting the observed orbital modulation and long term trends. We performed harmonic analysis using the clean program which eliminates aliases . Power spectra are presented in Figs. 68. Fig. 6 shows that all the significant peaks in WIRO data are at frequencies below 200 d-1. Fig. 7 shows the results for individual nights. Night 1 shows three main periods, $`P_154`$ min, $`P_233`$ min, and $`P_417`$ min; night 2 still shows $`P_1`$ and $`P_4`$, but also $`P_324`$ min.

Fourier analysis of the two Ekar runs is shown in Fig. 8. A year after WIRO observations, $`P_1`$ is still the more stable period, while $`P_2`$ and a modulation $`P_521`$ min are observed only in the first night.

Each of the detected periodicities is represented by a folded light curve in Figs. 9 and 10. The periods adopted are those of Table 2 and do not coincide exactly with those of Figs. 7 and 8. They were found by performing Fourier analysis on our data sets after the subtraction of all the other periods except the one under examination. In this way, each modulation is no longer distorted by the simultaneous presence of the other periods. The small error bars suggest that the modulations are rather well defined.

## 3 Discussion

### 3.1  The long term light curve

Table 3 collects some historical magnitude determinations of RW UMi, and reports, whenever possible, the range of the light modulation observed and the colours. While this list suggests a slow negative trend of the luminosity of the old nova, the reported (mean) magnitudes all lie within the observed short-term variations of the individual runs. Thus, only future observations will be able to confirm such a possible decline.

### 3.2 Physical parameters of the binary

Assuming that the period derived by Retter & Lipkin  is very close to the actual orbital period (see next section), and taking the recent correlation between the orbital period and the mass of the secondary star presented by Howell , we may estimate the mass of the secondary in RW UMi as $`M_20.115M_{}`$.

In order to constrain the mass of the primary we make use of the velocity of decline from the nova outburst. Warner  reports $`t_3`$ = 114 d, which corresponds to a decline velocity $`v_{\mathrm{dec}}`$ = 0.02 mag/d. Using the calibration given by della Valle  we obtain a maximum absolute magnitude $`M_{\mathrm{max}}7`$ mag. According to Livio  this is the luminosity produced by a white dwarf primary with $`M_10.7M_{}`$. With this, we obtain a mass ratio $`q=M_2/M_10.16`$.

The above calculated masses and main geometrical parameters of the binary system are collected in Table 4.

### 3.3 Superhumps and the orbital period

CVs can show several photometric periods, the most important ones are related to the orbital period. The most obvious case is that of CVs at sufficiently high inclination showing eclipses and/or luminosity bumps due to the hot spot. Other important photometric modulations, which do not seem to be strongly related with the inclination of the binary system, are called ‘superhumps’ and may differ by a few percent from the orbital period. Actually, two opposite cases are observed:

a) Positive superhumps. They mostly characterize SU UMa DN systems during super outbursts and are thought to be produced at relatively large mass transfer rates, when the outer edge of the accretion disc expands reaching the 3:1 resonance radius. Positive superhumps then represent the beat period between the apsidal prograde precession of the eccentric outer edge of the disc and the orbital motion of the secondary.

Taking equations (3.34) and (3.42) from Warner  we find a relation between the positve superhump period $`P_{\mathrm{SH}_+}`$ and the orbital one $`P_{\mathrm{orb}}`$:

$$\frac{P_{\mathrm{SH}_+}}{P_{\mathrm{orb}}}=\frac{1+q}{1+0.74q}.$$
(1)
With $`q`$ = 0.16, we obtain a positive superhump period 4% larger than the orbital one.

This situation corresponds to the scenario suggested by Retter & Lipkin . They also suggested that, similarly to the other two short-orbital-period old novae CP Pup and V1974 Cyg, the presence of positive superhumps might be a consequence of the post-outburst higher accretion luminosity.

b) Negative superhumps. A few CVs are known to show cyclical brightness modulations with periods shorter than the orbital one. This phenomenon is generally interpreted as the result of a precessing disc. The precession of the outer anulus of a tilted disc in this case is retrograde. From Warner’s  equation (2.62c) we obtain for the negative superhump period $`P_{\mathrm{SH}_{}}`$:

$$\frac{P_{\mathrm{SH}_{}}}{P_{\mathrm{orb}}}=\left(1+0.35\frac{q}{(1+q)^2}\right)^1.$$
(2)
For RW UMi this results in $`P_{\mathrm{SH}_{}}`$ being 4% shorter than $`P_{\mathrm{orb}}`$.

We conclude that the period observed by Retter & Lipkin  should represent the orbital period within 4% uncertainty. We furthermore note that our derived mass ratio $`q`$ = 0.16 well satisfies the condition $`q`$ 0.22 for the formation of positive superhumps. This might favour the interpretation by Retter & Lipkin  that $`P_{\mathrm{main}}`$ represents a positive superhump period.

### 3.4  The suborbital photometric periods

Besides the superhumps, which periods are close to the orbital one, CVs may show other periodic or quasi-period oscillations (QPOs) in their light curves. These phenomena might be similar to those observed in LMXBs. However, probably due to the longer timescales in CVs, they have not been studied as thoroughly. The secondary, suborbital, periodicities found in RW UMi are listed in Table 2.

Lai  has shown that the inner region of the accretion disc around a rotating magnetized collapsed object (neutron star, white dwarf) is subject to magnetic torques that induce warping and precession of the disc. The result is that a number of magnetically driven resonances between the disc and the rotating central dipole will arise. In some intermediate polar systems (IPs) the spin period of the WD is the same as that of the Keplerian orbit at the magnetosphere radius, so the two rotational motions are equal. In this case, the magnetic field will easily produce a warped inner disc which rotates with the spin period of the WD. However, in order to produce observable interferences with the orbital motion, the warped structure must extend to a large fraction of the disc, so that its variable geometry can interact with the accretion stream from the secondary. This possibility will depend on the viscosity inside the disc and the intensity of the magnetic field of the primary. In general, outside some critical radius, viscous forces will overcome the magnetic torque and warping will be destroyed .

However, since RW UMi is a short orbital period system, which can host only a very small accretion disc, even a modest magnetic field might produce a warping that extends to a considerable portion of the disc. Still, most IPs are not synchronous rotators, so that the warping will be a much more complicated process, in which the outer parts of the disc might switch between prograde and retrograde precession. Thus, RW UMi might display both positive and negative superhumps, perhaps in a rather complex and unstable fashion.

A key problem in the analysis of the suborbital periods is to detect the spin period of the WD. From our data we note that the period $`P_1`$ (Table 2) is the most stable one over at least one year. In terms of Keplerian velocities, it would correspond to an orbit with a radius $`r_1>R_{L1}`$, which is obviously impossible. We however note that $`P_1`$ fits quite well with the beat period between $`P_2`$ and $`P_{\mathrm{main}}`$. Rather curiously, this match is better than the one obtained using the two possible alternative values for the orbital period $`P_{\mathrm{orb}}=P_{\mathrm{main}}\pm 4\%`$ (see Section 3.3), although our precision does not really allow to conclusively decide between those possibilities.

$`P_2`$ was observed only during one night at WIRO and during one at Ekar. Although it is not the shortest of the detected periods, for the above reasons we tentatively assume it as the spin period of the magnetic WD. In the most general case of an inclined rotator, the inner disc region truncated by the magnetosphere of the primary should be warped and precess with the WD spin period. However, the most easily observable effect is the direct interaction of the mass stream from the secondary with the rotating magnetosphere of the WD. The changing inclination of the magnetic axis with respect to the stream produces a periodic modulation of the accretion rate onto the polar caps. In this way $`P_1`$ results as the beat period between the orbital and the rotational one ($`P_2`$).

Following this scenario, we should also look for other optical orbital sidebands like $`\omega _{\mathrm{rot}}+\mathrm{\Omega }_{\mathrm{orb}}`$ and $`\omega _{\mathrm{rot}}\pm 2\mathrm{\Omega }_{\mathrm{orb}}`$, where $`\omega _{\mathrm{rot}}=P_2^1`$ and $`\mathrm{\Omega }_{\mathrm{orb}}=P_{\mathrm{main}}^1`$. We find that

* $`\omega _{\mathrm{rot}}+\mathrm{\Omega }_{\mathrm{orb}}=60P_3^1`$;
* $`\omega _{\mathrm{rot}}+2\mathrm{\Omega }_{\mathrm{orb}}=77`$, between $`P_4^1`$ and $`P_5^1`$;
* $`\omega _{\mathrm{rot}}2\mathrm{\Omega }_{\mathrm{orb}}=9`$.

The latter, low frequency value corresponds to the time length of our observations, and is therefore not detectable. However, we note that a peak of minor significance at roughly 8.6 d-1 is present in the power spectrum published by Retter & Lipkin .

Taking a different approach, we might try to interprete the observed secondary periods in terms of inhomogeneities in the disc. We therefore compute their Keplerian radii (Table 5) and look for correlations with those radii where matter streams are in resonance with the tidal action by the secondary. Table 6 lists the main critical radii of the accretion disc for dynamical stability of orbiting matter streams. Radii labeled as $`r_{ij}`$ refer to resonances between the radial excursion of a test particle and the position of the secondary star .

Periods $`P_3`$ and $`P_5`$ would indicate Keplerian radii $`r_3`$ and $`r_5`$ which are close to resonance radii $`r_{32}`$, where superhumps are thought to be produced, and $`r_{43}`$, respectively. Finally, we note that $`P_4`$ is about $`0.5\times P_2`$ and could be its second harmonic.

If we assume that, in its steady quiescent state, the nova is a synchronous rotator and take $`P_2`$ as the spin period of the WD, then $`r_2`$ should coincide with the Alfven radius. However, this would suggest Keplerian orbits beyond the last stable orbit of the disc $`r_{\mathrm{max}}`$, and no disc could form at all. As a result, if we assume that some disc does presently exist, then RW UMi, in its post-outburst state, cannot be a syncronous rotator.

In order to have a standard disc, the stream from the secondary should be able to penetrate the magnetosphere of the primary down to the circularization radius $`r_{\mathrm{circ}}`$. Our data do not give us any clear indication about what the inner radius of the disc can be. However, we might reasonably assume that some material is presently orbiting around the WD at least between radii $`r_3`$ and $`r_{32}`$ which accounts for the probable presence of superhumps.

### 3.5 Pre-nova and post-nova state of RW UMi

In the previous subsection we have tentatively assumed $`P_2`$ as the spin period of the WD. If its Keplerian radius, $`r_2`$, is also the Alfven radius for the stream, $`r_{\mu ,\mathrm{stream}}`$, then, as $`r_{\mu ,\mathrm{stream}}>r_{\mathrm{circ}}`$, no disc can form.

However, it appears plausible that a sort of an equilibrium status, with no disc formation, might have been achieved by the fainter pre-nova, while the stream was carrying only $`\dot{M}10^9M_{}\mathrm{yr}^1`$. Using equations (6.5) and (7.6) from Warner  we may write

$$r_{\mu ,\mathrm{stream}}=3.66\times 10^{10}\mu _{34}^{4/7}M_1^{1/7}\dot{M}_{16}^{2/7}\mathrm{cm}.$$
(3)
Thus, within our hypothesis we derive $`\mu _{34}0.84`$, or, consequently, $`B=1.5\times 10^7`$ G.

The post-nova is however brighter than the pre-nova. Again we use equation (3) to derive $`r_{\mu ,stream}`$ for the high state of the nova, characterized by $`\dot{M}10^8M_{}\mathrm{yr}^1`$, which would account for yielding $`M_{\mathrm{min}}5.1`$ mag. We find that the stream can now penetrate the magnetosphere down to a radius

$$r_{\mu ,\mathrm{stream}}1\times 10^{10}\mathrm{cm}$$

This radius is probably not close enough to the circularization radius to allow for the formation of a standard disc, but material can finally penetrate the magnetosphere ensuring stable accretion and also the formation of a probably chaotic disc. Note that the resonance radii of Table 6 are now populated by streaming gas. All intermediate polars have $`r_{\mathrm{inner}}<r_{\mathrm{corotation}}`$, so they are asynchronous rotators. The gas between the inner disc radius and the corotation radius might be rather turbulent and instability phenomena at the resonance radii could be enhanced.

In conclusion, if we assume $`P_2`$ as the spin period of the WD, we find that RW UMi in its steady quiescent state should be a discless intermediate polar, while in its post-outburst state it appears as an intermediate polar with a somewhat unstable truncated accretion disc. The magnetic moment of the primary $`\mu _{34}0.84`$ is rather large for short orbital period intermediate polars. RW UMi should then tend with time to go back to a pure discless state.

## 4 Summary

1. We have presented new photometric time-series of the old nova RW Umi.
2. The data revealed the presence of several suborbital periods. The most stable one of these is $`P_154`$ min, which is detected in all data sets. We find that $`P_1`$ agrees well with a beat period between our period $`P_243`$ min and the period found by Retter & Lipkin  $`P_{\mathrm{main}}`$ = 85 min, which is assumed to be very close to the orbital period.
3. The most probable explanation for the presence of $`P_2`$ is that it represents the spin period of the white dwarf. The other suborbital periods can then be found as being optical orbital sidebands with respect to $`P_2`$ and/or correspond to Keplerian resonance radii.
4. Assuming a mixture of ‘disc’ and ‘magnetic’ periods we find that RW UMi, in its post-nova state, must be an asynchronous intermediate polar, as the Keplerian radius of the assumed spin period $`P_2`$, i.e. the Alfven radius, is larger than the circularization radius.
5. Comparing the mass-transfer rates from the pre- and the post-nova, we conclude that RW UMi in its quiescent, i.e. pre-nova, state is a discless intermediate polar.

## ACKNOWLEDGMENTS

One of us (F.T.) wishes to thank the Department of Astronomy of the University of Padova for the kind hospitality. This work has been supported by NSF REU Grant AST 9732039, and by the Italian MURST.


text:
# Ge(001)-(("2 1"¦"0 3"))-Pb ↔ (("2 1"¦"0 6"))-Pb: Low-temperature two-dimensional phase transition

## 1 Introduction

Surfaces and interfaces determine many chemical and physical properties of materials and play a key role in nanoscale electronic devices. A well-known example of the influence of an interface on electronic properties is the Schottky barrier at metal-semiconductor interfaces. In principle the height of the barrier is determined by the difference between the work function of the metal and the electron affinity of the semiconductor, but experimental investigations showed that the atomic structure directly at the interface may influence the barrier height considerably . Model systems chosen to study this phenomenon were, e.g., lead adsorbed on different silicon and germanium surfaces, since for these systems abrupt interfaces with planar layers free of alloying were expected. Extensive scanning tunneling microscopy (STM) investigations of the different phases of lead on Ge(001) surfaces for different lead coverages and preparation temperatures partially confirmed and partially disproved this simple idea . For lead coverages above $``$1.5 monolayer (ML) quite flat and uniform lead overlayers forming a commensurate $`\left(\genfrac{}{}{0pt}{}{\text{2\hspace{0.33em}1}}{\text{0\hspace{0.33em}3}}\right)`$ and at higher coverage an incommensurate c(8$`\times `$4)i reconstruction are observed. That the latter overlayer exhibits metallic properties is evident from the low STM tip bias of 20 mV used to image this surface. Below 1.5 ML lead coverage and at preparation temperatures below 300C commensurate reconstructions with (2$`\times `$2) and c(8$`\times `$4) unit-cells are formed. Both reconstructions are described by lead dimers covalently bonded to the germanium substrate. For the Ge(001)-(2$`\times `$2)-Pb reconstruction with a nominal lead coverage of 0.5 ML the germanium substrate is still dimerized whereas in the c(8$`\times `$4) phase with a nominal lead coverage of 0.75 ML half the dimer bonds are broken  – as confirmed recently by first-principle total energy calculations . The lead dimers in these two reconstructions are tilted to avoid the energetically unfavourable configuration of two approximately sp3 hybridized lead atoms each with a half-filled “dangling bond”. This buckling is static as opposed to the clean Si(001)- and Ge(001)-(2$`\times `$1) surfaces where the dimers flip dynamically at room temperature . The lead-induced reconstructions of the Ge(001) surface described so far more or less exhibit the expected flat overlayer necessary for an ideal Schottky barrier model system – although only the phases with a sufficiently high lead coverage fulfill the condition of being metallic. At low lead coverages below 1 ML and preparation temperatures above $``$300C the surface exhibits a completely different morphology. This phase is not well ordered and has been described as the “rough phase”  with height differences of more than five atomic layers. The formation of this phase with an unusually high roughness for lead overlayers may be explained by the significant stress induced by the incorporation of lead atoms with a covalent radius of 1.47 Å on germanium lattice sites (covalent radius 1.22 Å), i.e., an unexpected alloying may occur at low lead coverages and high preparation temperatures.

The lead-induced surface reconstructions of the Ge(001) surface described so far possess a wide diversity of properties. But the Ge(001)-$`\left(\genfrac{}{}{0pt}{}{\text{2\hspace{0.33em}1}}{\text{0\hspace{0.33em}3}}\right)`$-Pb reconstruction, formed at an intermediate lead coverage of around 5/3 ML seems to be even more interesting. From the point of view of coverage it resides between the low-coverage covalently-bonded and the high-coverage metallic phases and may share properties common to both.

Based on thorough STM investigations and some diffraction results a structural model was proposed for the Ge(001)-$`\left(\genfrac{}{}{0pt}{}{\text{2\hspace{0.33em}1}}{\text{0\hspace{0.33em}3}}\right)`$-Pb reconstruction with 5/3 ML lead in a distorted Pb(111) layer on top of a completely dimerized Ge(001) substrate . The complete dimerization of the germanium substrate is surprising since at lower lead coverage in the Ge(001)-c(8$`\times `$4)-Pb reconstruction the dimerization of the (even lower coverage) Ge(001)-(2$`\times `$2)-Pb surface and the clean Ge(001)-(2$`\times `$1) is already partially removed. To verify the proposed structural model for the room-temperature phase and to elucidate the properties of the novel Ge(001)-$`\left(\genfrac{}{}{0pt}{}{\text{2\hspace{0.33em}1}}{\text{0\hspace{0.33em}6}}\right)`$-Pb low-temperature phase  we investigated both systems in a surface x-ray diffraction (SXRD) experiment using synchrotron radiation.

## 2 Experimental

### 2.1 Preparation and measurements

The sample was prepared in an ultra-high vacuum (UHV) system equipped with reflection high-energy electron diffraction (RHEED) and low-energy electron diffraction (LEED) facilities. A clean (2$`\times `$1) reconstructed Ge(001) sample was prepared by the standard procedure of cycles of sputtering at 350-450C for 30-60 min using 500 eV Ar+ ions and annealing at 600-650C for 15 min. Lead was deposited while observing the RHEED-patterns: Starting from the (2$`\times `$1) reconstruction of the clean surface first the c(8$`\times `$4)-Pb and finally the Ge(001)-$`\left(\genfrac{}{}{0pt}{}{\text{2\hspace{0.33em}1}}{\text{0\hspace{0.33em}3}}\right)`$-Pb phase was observed. The change from c(8$`\times `$4) to $`\left(\genfrac{}{}{0pt}{}{\text{2\hspace{0.33em}1}}{\text{0\hspace{0.33em}3}}\right)`$ periodicity is accompanied by the disappearance of the characteristic c(8$`\times `$4)-streaks and the emergence of half-order reflections in the $`(1\overline{1})`$ azimuth. The long-range order of the surface reconstruction was improved by annealing the sample at 150C.

The sample was then transferred to a portable UHV chamber equipped with a helium closed-cycle two-stage sample cooling system and brought to the BW2 wiggler beamline at the Hamburg Synchrotron Radiation Laboratory (HASYLAB) for the x-ray diffraction measurements. The x-ray photon energy was set to 10 keV ($`\lambda =1.24`$ Å) and the glancing angle of incidence to 0.65 (i.e. above the critical angle to reduce the uncertainties in the measured intensities arising from mechanical displacements). A data-set consisting of 182 in-plane reflections, 236 reflections along nine fractional-order rods and 81 reflections along three integer-order rods (so-called crystal truncation rods, CTRs) was recorded for the Ge(001)-$`\left(\genfrac{}{}{0pt}{}{\text{2\hspace{0.33em}1}}{\text{0\hspace{0.33em}3}}\right)`$-Pb structure determination. After completing the room-temperature measurements the sample was cooled until the temperature of the sample holder reached 20 K. For the low-temperature $`\left(\genfrac{}{}{0pt}{}{\text{2\hspace{0.33em}1}}{\text{0\hspace{0.33em}6}}\right)`$-phase 313 in-plane reflections, 371 reflections along 18 fractional-order rods and 106 reflections along four CTRs were recorded.

### 2.2 Data-analysis

The integrated intensities were corrected for the Lorentz factor, polarization factor, active sample area and the rod interception appropriate for the z-axis geometry . The in-plane data were averaged assuming that all four domains which are possible within the p2-symmetry are equally populated. In this way the room-temperature in-plane data set was reduced to 96 and the low-temperature in-plane data set to 223 symmetry-inequivalent reflections. The recording of extensive data sets with 413 symmetry-inequivalent data points for the room-temperature and 700 data points for the low-temperature phase was necessary to determine all the relevant parameters in these complex reconstructions with large unit cells.

The conventional surface coordinate system with $`𝐚=1/2[110]_{\text{cubic}}`$, $`𝐛=1/2[\overline{1}10]_{\text{cubic}}`$ and $`𝐜=[001]_{\text{cubic}}`$ is used. The cubic coordinates are in units of the germanium lattice parameter of 5.66 Å at 300 K and therefore $`|𝐚|=|𝐛|=4.00`$ Å and $`|𝐜|=5.66`$ Å.

Only seven layers were included in the data-analysis to describe the substrate relaxations, i.e., less than the ten layers used for the analysis of the clean surface  to keep the number of free parameters reasonably low. In addition the elastic energy calculated using the Keating algorithm  was used as a constraint on the substrate parameters in the data analysis. The Keating term restricts the substrate parameters if the optimization runs into a shallow $`\chi ^2`$ minimum. Nevertheless the data analysis was only possible using a computer program with the ability to handle efficiently the 152 parameters necessary to describe the atomic positions in the low-temperature phase together with further parameters describing thermal motion and disorder and the Keating constraints for 105 of the parameters.

The overall good agreement between the result of the structural refinement and the measured data is apparent from figure 1 for the room-temperature results and figure 2 for the low-temperature data. Using in total 82 parameters a reduced $`\chi ^2(|F|)`$ value of 2.2 was obtained for the room-temperature phase. For the low-temperature phase in total 155 parameters were determined and the reduced $`\chi ^2(|F|)`$ was 3.3. From a technical point of view the number of parameters is exceptionally large for a surface structure determination.

## 3 Results and discussion

In the following two subsections we focus on the most important results and suggest a possible interpretation.

### 3.1 Room-temperature phase

The analysis of the room-temperature data basically confirmed the structural model with a distorted Pb(111) layer on top of the dimerized germanium substrate proposed previously by Falkenberg et al (see figure 3 and table 1). We could confirm that the germanium dimers of the clean surface broken when forming the Ge(001)-c(8$`\times `$4)-Pb reconstruction are reestablished at a higher lead coverage of 5/3 ML. In addition to this remarkable result some interesting details are revealed.

The atomic displacement parameter (ADP) used in the analysis of diffraction data describes both thermal motion and static disorder . The ADP listed in table 1 for the lead atoms 1-9 (see figure 3) is anisotropic – larger in the surface plane than perpendicular to it. This suggests the picture of a metallic overlayer with static disorder and a significant mobility of the individual atoms around their equilibrium position. An extra ADP was necessary to describe the behavior of lead atom 10. This atom is on a high-symmetry site and has a larger distance to the germanium substrate than the other lead atoms. It is therefore more weakly bonded. This is reflected in its higher ADP.

There is a striking anisotropy in the ADP of the first layer germanium atoms. The amplitude of $`0.4\pm 0.1`$ Å in the $`z`$-direction (perpendicular to the surface plane) matches well to the values of 0.50 Å  and 0.42 Å  for the clean surface. This may indicate that, similar to the clean surface, the germanium dimers are dynamically flipping beneath the lead overlayer. This analogy to the clean surface seems to be obvious, although static disorder and dynamic movement cannot be distinguished by SXRD. The bond length of $`2.45\pm 0.08`$ Å for two germanium dimers (12-13 and 14-15) is in excellent agreement with the value of $`2.44\pm 0.01`$ Å for the clean Ge(001)-(2$`\times `$1) reconstruction . The simple picture of an unperturbed germanium surface layer is supported by this fact. However, a bond length of $`2.68\pm 0.07`$ Å was determined for the germanium dimer 11-16. This is an increase of 9 % compared to the substrate bond length of 2.45 Å. Since the germanium atoms 11 and 16 have the smallest observed Ge-Pb distance of $`2.71\pm 0.06`$ Å to their next neighbor atoms 7 and 8 (drawn dark grey in figure 3) a stronger, partially covalent Pb-Ge interaction weakening the Ge-Ge dimer bonds may cause the observed increase in the dimer bond length.

To summarize the results for the Ge(001)-$`\left(\genfrac{}{}{0pt}{}{\text{2\hspace{0.33em}1}}{\text{0\hspace{0.33em}3}}\right)`$-Pb room-temperature phase, the following model is consistent with the observed structural parameters: A distorted, metallic Pb(111) overlayer is only weakly bonded within the layer and to the substrate and has an increased in-plane motion or disorder compared to the out-of-plane direction. A partially covalent interaction between one germanium dimer and the two next neighbor lead atoms per unit-cell stabilizes the overlayer. This interaction tends to make the overlayer commensurate with the substrate and to weaken the corresponding dimer bond under the lead layer. The remaining two substrate dimers per unit-cell continue the dynamic flipping movement characteristic of the clean surface. The delicate balance between covalent and metallic behavior of the surface atoms indicates that this system is likely to respond to small changes of the total energy, e.g., by simply changing the temperature. The structural parameters derived from the SXRD measurements suggest that this system is predisposed to structural changes as a function of the temperature.

### 3.2 Low-temperature phase

Apart from the additional reflections characteristic of the $`\left(\genfrac{}{}{0pt}{}{\text{2\hspace{0.33em}1}}{\text{0\hspace{0.33em}6}}\right)`$-periodicity the in-plane data measured at room and at low temperature look quite similar (see figure 1(a) and figure 2(a)) – mainly due to the characteristic ring pattern formed by the strongest reflections. This indicates that the fundamental framework of the reconstruction remains intact upon cooling. But despite the similarities there are definite differences as can be seen by the comparison of the fractional-order rod scans shown in figure 5. The differences are not artifacts caused by choosing incorrect scaling factors since they are not uniform along all rods and they are not just caused by changes in thermal motion since they occur in some cases only in a certain range of $`l`$-values (see, e.g., the (2,3/2,$`l`$)-rod for $`l>2`$) and in other cases along the whole rod (see, e.g., the (8/3, 5/3, $`l`$)-rod). The parameters determined in the detailed data analysis are listed in table 2 and table 3 and the model is shown in figure 4. The weakly bonded lead atom (number 10 in figure 3) remains in a high-symmetry site upon cooling (numbers 10 and 20 in figure 4). The differences to a similar structural model for the low-temperature $`\left(\genfrac{}{}{0pt}{}{\text{2\hspace{0.33em}1}}{\text{0\hspace{0.33em}6}}\right)`$-reconstruction with lead atom 2 on a high-symmetry trench site amounts only up to 0.09 Å for most atoms. But the height (i.e. $`z`$-position) difference of 1.48 Å for the high-symmetry site atoms 10 and 20 is necessary to describe the data adequately and therefore the chosen symmetry is unambiguously the one occurring in the low-temperature phase. For the lead layer the lowering of the $`z`$-position of atom 20, causing the described height difference, and the outward displacement of the next neighbor lead atoms 16 and 19 by 0.77 Å are the most pronounced changes in comparison to the room-temperature structure. All lead atoms are adequately described by one in-plane and one out-of-plane ADP. The ADP amplitudes of $`0.21\pm 0.03`$ Å in the surface plane are smaller than the $`0.26\pm 0.04`$ Å and $`0.34\pm 0.10`$ Å at room temperature, but not as low as one would expect solely from a thermal effect. Static disorder would explain the larger in-plane ADP-amplitudes. The counterpart of the “floating” lead overlayer at room temperature, which is only weakly pinned to its commensurate position, is probably a network made up of a variety of very similar low-temperature structures, which forms due to a broad minimum of the total energy as a function of the in-plane position.

The germanium atoms in the low-temperature phase are well described by an isotropic ADP with a value close to zero. This is clear evidence of a well-ordered substrate structure at low temperatures. The precise value of this ADP does not strongly influence the reduced $`\chi ^2`$ and therefore it was fixed at a constant small value.

The germanium dimers close to the lead atoms 10 and 20 are statically tilted in alternate directions and probably the dynamical motion is frozen out at low temperatures. The germanium dimer atoms next to lead atom 10 have a $`z`$-position difference of $`0.29\pm 0.06`$ Å. The germanium dimer atoms next to lead atom 20, which has a significantly lower $`z`$-position than at room temperature, have a larger $`z`$-position difference of $`0.54\pm 0.05`$ Å. Both values are smaller than the corresponding values for the clean surface: $`0.74\pm 0.15`$ Å and $`0.69\pm 0.04`$ Å  for the room-temperature phase and $`0.84\pm 0.02`$ Å  for the low-temperature phase. This indicates a stronger interaction of the adsorbate with the substrate at low temperature than at room temperature. In addition the dimers 21-27 and 30-32, which probably have a partially covalent interaction with the lead layer at room temperature are (within the uncertainty of the measurements) symmetric at low-temperature. A symmetrical configuration would be energetically unfavourable without germanium-lead interactions arising from the adsorbate layer. It is interesting that not only the distance from these germanium atoms in a symmetric configuration to their next neighbor lead atoms is close to the covalent bond distance, but also the approximately sp3 hybridized germanium atoms 26 and 31 of the statically tilted dimers have a distance similar to the covalent bond length to the lead atoms 1 and 3. Therefore, we conclude that the lead layer in the low-temperature phase is more firmly bonded to the substrate with a more covalent character than it is the case at room temperature.

The covalent bond length for lead is 2.94 Å. The lead atom pairs (1,6), (3,9), (11,16) and (13,19) have a distance of $`3.0\pm 0.05`$ Å between the atoms of each pair. This may indicate the presence of partially covalent bonds even within the overlayer. This assumption is supported by figure 6. In this figure the number of lead atoms with next neighbor distances (within the lead layer) below a certain threshold value versus the threshold value distance is shown. The only significant difference between the curves for the room and low-temperature phase (dashed and solid line respectively) is a plateau from $`3.0`$ Å to $`3.2`$ Å for the low-temperature phase curve. That means the small distances within the lead layer are a significant feature of the low-temperature phase. It is interesting to note that the lead atom pairs with small interatomic distances surround the lead atoms 10 and 20, which exhibit a pronounced height difference compared to their position at room temperature.

Assuming that there are partially covalent Ge-Pb and Pb-Pb bonds immediately raises the question of the hybridization of the lead atoms and the electron density distribution within the lead layer. As a first approximation one may assume that one germanium atom of each buckled dimer has approximately two electrons for a Pb-Ge bond and that the electron density within the lead layer is not strongly influenced by these bonds. Each germanium atom of the symmetric dimers has approximately one electron and the lead overlayer has to donate one electron for each covalent Pb-Ge bond. A completely covalent Pb-Pb bond would localize two electrons. Therefore, covalent bonds from symmetric germanium dimers to lead atoms and within the lead layer will lower the electron density in other parts of the overlayer which means that the number of partially covalent bonds in the energetically most favorable configuration is limited. The formation of stronger covalent bonds is inherently inhibited by the basic atomic structure of the system; the pathway for lowering the total energy of the system by the formation of covalent bonds is inherently frustrated. The only way out of this frustration would be the formation of a completely new structure with a lower total energy. However no such structure with a lead coverage of 1.5 ML has been observed and it is unlikely that the activation energy for this phase transition involving considerable mass transport would be available at low temperatures.

Considering the classification of the phase transition the results of the structure determination presented here suggest that it is of order-disorder type. The low-temperature phase with static germanium dimers and stronger bonding is the ordered form of the room-temperature phase with dynamically flipping dimers and a “floating” lead overlayer. We propose that within the Pb/Ge(001) system there is a delicate balance of different pathways for lowering the total energy of the system and that the true nature of the system is probably more complex than our simplistic bonding arguments would suggest.

It is instructive to compare the phase transition described here with the transition from the hexagonal incommensurate phase of lead on Si(111) to the Si(111)-$`\left(\genfrac{}{}{0pt}{}{\text{ 3\hspace{0.33em}2}}{\text{-1\hspace{0.33em}1}}\right)`$-Pb low-temperature structure . In both systems the adsorbate forms a Pb(111) overlayer. The overlayer is incommensurate (Pb/Si(111)) or weakly pinned to a commensurate position (Pb/Ge(001)). In both low-temperature phases interatomic distances indicate the presence of Pb-substrate bonds with predominantly covalent character. The major differences arise from the substrate; in the case of Si(111) the lowest energy structure, the (7$`\times `$7) reconstruction, is removed by lead deposition and annealing – the silicon substrate of the hexagonal incommensurate phase has essentially the undisturbed bulk-truncated structure which is not stable without the lead overlayer. In the case of Ge(001) the dimers in the (2$`\times `$1)-reconstruction of the clean surface are partially removed by lead deposition and reestablished themselves at larger lead coverages. Therefore, the lead overlayer of the Ge(001)-$`\left(\genfrac{}{}{0pt}{}{\text{2\hspace{0.33em}1}}{\text{0\hspace{0.33em}3}}\right)`$-Pb reconstruction stabilizes the reconstructed substrate and the lead overlayer of the hexagonal incommensurate phase on Si(111) stabilizes the bulk-truncated structure of the substrate. Both systems lower their total energy upon cooling by the formation of bonds with covalent character and the capability of forming such bonds plays a key role in determining the substrate structure.

## 4 Summary

At room temperature the topmost germanium layer of the Ge(001)-$`\left(\genfrac{}{}{0pt}{}{\text{2\hspace{0.33em}1}}{\text{0\hspace{0.33em}3}}\right)`$-Pb reconstruction is similar to the clean (2$`\times `$1) reconstructed surface which is described by dynamically flipping germanium dimers. The mainly metallic lead overlayer is probably pinned to a commensurate registry with the substrate by weak covalent interactions to one germanium dimer per unit-cell.

Upon cooling, the dimers with a weak covalent bond to the adsorbate at room temperature remain in a symmetrical position. The remaining dimers are statically tilted in alternate directions as in the Ge(001)-(2$`\times `$1) surface reconstruction at low-temperatures. Within the lead overlayer the main difference to the positions at room temperature is a height difference of the two lead atoms that are weakly bonded on high-symmetry sites at room temperature. Additional partially covalent interactions between the substrate and adsorbate and within the adsorbate layer stabilize the system.

The structural changes are consistent with the occurrence of an order-disorder phase transition in this system. The phase transition is partially frustrated by the mutual exclusiveness of different pathways for lowering the total energy.

Total energy calculations for this complex system with such a large unit-cell are definitely a challenge but, with regard to the interesting properties and the expected insight in the delicate energy balance, probably worth the effort. We hope that the determination of the detailed atomic structure reported here will pave the way for theoretical investigations of the Ge(001)-$`\left(\genfrac{}{}{0pt}{}{\text{2\hspace{0.33em}1}}{\text{0\hspace{0.33em}6}}\right)`$-Pb reconstruction at $`T=0`$ K and at elevated temperatures.

We thank the staff of HASYLAB for their technical assistance. Financial support from the Danish Research Council through Dansync, the Bundesministerium für Bildung, Wissenschaft, Forschung und Technologie (BMBF) under project no. 05 KS1GUC/3, the “Graduiertenkolleg Physik nanostrukturierter Festkörper” and the Volkswagen Stiftung is gratefully acknowledged.

## References

## Tables and table captions

## Figure captions


text:
# Characterising anomalous transport in accretion disks from X-ray observations

## 1 Introduction

Deeper understanding of the transport mechanisms that operate within accretion disks is important for a broad range of X-ray emitting astrophysical objects. In the absence of local measurements, the key questions are (1) to quantify the ways in which the unseen transport processes are anomalous as distinct from diffusive; and (2) to establish how this may be determined remotely from observations of global quantities such as the input and outflow of energy. This information can then be used to provide a constraint for turbulence/instability models of astrophysical accretion disks.

Here we analyse the total X-ray flux over several years from two accreting astrophysical objects – Cygnus X-1 and GRS1915+105, the stellar-mass analogues of disk-jet active galaxies powered by a massive black hole – and, for comparison, the nonaccreting Crab which is powered by a neutron star. The reasons for selecting these sources for statistical analysis are twofold. First, since 1996 February 20 they have been observed continually for several years by the All-Sky Monitor (ASM) on board the RXTE satellite (Swank et al. rxte (2001)), providing large data sets of around thirty thousand points. This enables us to seek correlations over several orders of magnitude up to the longest accessible timescales. Second, the source luminosities are sufficiently high to neglect instrument thresholds, uncertainties, and other sources in the field of view. The raw data are held at the Goddard Space Flight Center (GSFC) and can be accessed via their website111http://heasarc.gsfc.nasa.gov/docs/xte/asm\_products.html. Each point represents the total X-ray flux (measured by the number of counts during periods that last 90 seconds) in the range 1.3–12.2 $`keV`$, and the breakdown into three energy bands (1.3–3, 3–5 and 5–12.2 $`keV`$) is also available. We have analysed the different channels and found them to be similar to the total flux. Sampling intervals between the 90-second X-ray counting periods are distributed with means of 93 minutes for the Crab, 77 minutes for Cygnus X-1, and 96 minutes for GRS1915+105. To give an indication of the spread of sampling intervals, 90% of the intervals are below 186 minutes for the Crab, 187 minutes for Cygnus X-1, and 194 minutes for GRS1915+105. The implications for our techniques are explained where appropriate below (and see Appendix). Calibration is undertaken by the ASM/RXTE team and the processed data are freely accessible on their website222http://xte.mit.edu/XTE/asmlc/ASM.html (Bradt et al. asm (2001)). The RXTE counts are not directly proportional to luminosity, so the exact relationship between these two quantities would have to be accounted for in any model. However, this will affect only the nature of the PDFs and not the temporal correlation in the data. The three raw X-ray time-series are plotted in Figs. 1, 2 and 3.

Previous studies have focused on spatial and spectral structure, quasi-periodic oscillations, and modelling of temporal variability: see Weisskopf et al. (cr (2000)) and references therein for the Crab; Maccarone et al. (cy (2000)) and references therein for Cygnus X-1; and Belloni et al. (belloni (2000)) and Rao et al. (rao (2000)) for summaries of the spectral and temporal analyses of GRS1915+105. Dhawan et al. (dhawan (2000)) have studied the jet of this source, and the importance of microquasars in general is discussed by Mirabel and Rodríguez (mirabel (1999)). Nayakshin et al. (nay (2000)) model the gross variability of GRS1915+105 on all but the shortest timescales, and recent observations by Chaty et al. (chaty (2001)) have been used to search for interactions with the surrounding interstellar medium. As a complement to these techniques, we here examine three key statistical measures (described further in Sect. 2) for each time-series. These are the probability density function (PDF), the growth of range (extent of statistical self-similarity), and the differenced form, which we then compare with the well-known signatures of Gaussian noise and random walks. In particular, we test for signatures of long-range correlations that are the hallmarks of turbulent transport. Insofar as these techniques can be applied successfully to the observational signatures of anomalous transport in the present context, they may also be transferable to related questions in space and laboratory plasma physics.

There is already much evidence that accretion disks are locales for turbulent transport and instabilities. In order to explain typical accretion rates observed in a range of disk types, the standard disk model (Shakura & Sunyaev shak (1973); Pringle & Rees pringle (1972)) uses turbulent viscosity to produce an appropriate outward transport of angular momentum. Following extensive numerical simulations (reviewed by, for example, Gammie gamm (1998)), the source of this turbulence is now believed to be a magnetic shear instability (Balbus & Hawley balb (1991)). In addition, $`1/f`$ power spectra (inverse power-law frequency dependence with unspecified index) – suggestive of highly-correlated, possibly self-organised critical, behaviour – have been observed in diverse accretion systems (reviewed by, for example, Dendy et al. dendy (1998)). It is also well established that instabilities occuring in accretion disks (Dubus et al. dubus (2001) and references therein) can give rise to anomalous statistics. However, these statistics will differ from those of fully-developed turbulence, and hence we seek to constrain such models using the differencing and rescaling technique described below.

Self-similarity, non-Gaussianity and non-trivial temporal scaling together are strong indications of highly-correlated processes such as turbulence (Bohr et al. bohr (1998)). We will show how trivial scaling of near-Gaussian fluctuations in the Crab X-ray signal – evidence of diffusive transport – contrasts with non-trivial scaling of non-Gaussian fluctuations in the X-ray signals from Cygnus X-1 and GRS1915+105. Whilst there may be other methods we could apply, we are confident that the methods described here are sufficiently insensitive to the timing and counting errors in the data; many techniques are unreliable even when random errors are small (Sornette sorn (2000)).

## 2 Techniques

### 2.1 The Probability Density Function (PDF)

The first step in our analysis of each data set is the construction of its PDF. The PDF $`P(y)`$ of a variable $`Y`$ is defined such that the probability that $`Y`$ lies within a small interval $`dy`$ centred on $`Y=y`$, is equal to $`P(y)dy`$. $`P(y)`$ is normalised so that

$$_{y_{min}}^{y_{max}}P(y)𝑑y=1.$$

Here, we use $`dy=(y_{max}y_{min})/100`$ and further normalise each PDF to its mean $`y`$ and standard deviation $`\sigma `$ to enable comparison with theoretical distributions. We note that the PDF was used by Antar et al. (antar (2001)) to characterise turbulent fluctuations in tokamak edge plasmas, and by Bhavsar & Barrow (bhav (2001)) to model the magnitude distribution of the brightest cluster galaxies. Bramwell et al. (bram1 (1998, 2000, 2001)) discuss the PDFs of fluctuations in highly-correlated systems (see Sect. 4), and Burlaga (burlaga (2001)) presents a review of log-normal distributions in the turbulent solar wind.

### 2.2 Growth of range

Consider a self-similar function $`y\left(t\right)`$. The difference between the maximum and minimum values of $`y`$ during a time interval $`\mathrm{\Delta }t`$ defines its range for that interval, $`R\left(\mathrm{\Delta }t\right)`$. Since $`y`$ is self-similar, the ensemble-averaged value of $`R`$ will scale with $`\mathrm{\Delta }t`$. We may write

$$R\left(\mathrm{\Delta }t\right)=c\mathrm{\Delta }t^H$$
(1)
with $`c`$ and $`H`$ constants; $`H`$ here defines the Hurst exponent (Hastings & Sugihara frac (1993)). For data that are only approximately self-similar, we use the above relation to check their closeness to self-similarity and also to obtain an effective value for $`H`$ as follows. By moving the window $`\mathrm{\Delta }t`$ one point at a time through the raw data, an array of $`R\left(\mathrm{\Delta }t\right)`$ values is created from which the mean $`R`$ is found (reducing the effects of uneven sampling). This is repeated for a range of $`\mathrm{\Delta }t`$ within the length of the data set. A plot of $`\mathrm{log}R\left(\mathrm{\Delta }t\right)`$ against $`\mathrm{log}\mathrm{\Delta }t`$ will reveal any deviations from self-similarity whilst the slope will give the best estimate of $`H`$. We use linear regression to calculate the $`95\%`$ confidence interval for $`H`$ (see Figs. 6, 9 and 13).

Trivially, a function that is exactly constant over time has $`H=0`$. At the other extreme, $`H=1`$ indicates a function whose range increases linearly with time (for positive $`c`$ in Eq. 1). Intermediate values of $`H`$ are generated by fractal functions, random Gaussian noise ($`H0.2`$), and Gaussian random walks (whose next value in time is the sum of the previous value and a random Gaussian increment; $`H0.5`$). The value of $`H`$ does not uniquely establish correlation, however; uncorrelated series may present significant probabilities of observing greater values as the timescale increases. Consequently, the growth of range can be rather insensitive as a measure of correlation. We can in principle define a measure of correlation $`\beta `$ in terms of fractal exponents such as $`H`$ ($`\beta =2H+1`$). Following Malamud & Turcotte (1999), $`\beta =0`$ for uncorrelated noise and $`\beta =2`$ for a Gaussian random walk. However, the use of only one method to estimate an unknown exponent (and hence $`\beta `$) is to be avoided (Schmittbuhl et al. schmitt (1995)). For example, comparing values of $`H`$ obtained from the growth of range of synthesized series with known $`\beta `$, Malamud & Turcotte (mala (1999)) find that $`H0.2`$ for $`\beta <0.5`$, whereas the slopes of the power spectra are equal to $`\beta `$ for $`0<\beta <2`$. We therefore obtain, for comparison, estimates of the strengths of correlation from the slopes of power spectra with Thompson multi-tapering (Percival & Walden percival (1993)) to reduce the variance.

The Hurst exponent has been used to quantify solar magnetic complexity (Adams et al. adams (1997)), correlation in the dynamics of the upper photosphere (Hanslmeier et al. hans (2000)), and persistence in solar activity (Lepreti et al. lepreti (2000) and references therein). In conjunction with rescaling techniques, the Hurst exponent has also been used to quantify self-similarity and long-range correlations in turbulent fluctuations in magnetic fusion plasmas (Carreras et al. carr (1998)).

### 2.3 Differencing and rescaling

Starting from the raw data $`y\left(t\right)`$, we first form a set of differenced series $`Z(t,\tau )`$ for a range of values of the time-lag $`\tau `$:

$$Z(t,\tau )=y\left(t\right)y\left(t\tau \right).$$
(2)
From these we calculate a set of PDFs for $`Z(t)`$, one for each value of $`\tau `$, which we denote by $`P(Z,\tau )`$ (not normalised) – see Fig. 4a. It is often found to be the case that the peaks $`P(0,\tau )`$ scale approximately as $`\tau ^\alpha `$, as in Fig. 4b. Such scaling is, of course, sought using the values of $`P(0,\tau )`$ because they are the most accurate points of the PDFs, having the most data points in their bins. The question then arises whether the differenced series $`Z(t,\tau )`$ can all be derived from a single PDF that scales. That is, whether a common functional form emerges if we use the exponent $`\alpha `$ to rescale the axes such that

$$ZZ\tau ^\alpha =Z_s\mathrm{and}PP\tau ^\alpha =P_s,$$
(3)
causing the separate PDFs to collapse onto the same ($`\tau =1`$) curve – Fig. 4c (compare Mantegna & Stanley diff (1995)).

Clearly, a reasonable straight line in Fig. 4b (at least up to some maximum timescale $`\tau _{max}`$) is necessary but not sufficient for rescaling to succeed, since only the peaks of $`P(Z,\tau )`$ are used to determine $`\tau `$. Physically, successful rescaling of the differenced X-ray data would imply, on all timescales up to $`\tau _{max}`$, that (i) the sizes of X-ray fluctuations (that is, the *differences* between the observed values) are governed by a single type of process, and (ii) the total X-ray output is correlated, not random, in time. In the trivial case where $`P(Z,\tau )`$ is independent of $`\tau `$, we would infer the absence of temporal correlation in the X-ray output.

The value of $`\alpha `$ is given by $`\alpha =1/m`$ where $`m`$ is the slope of $`logP(0)`$ v. $`log\tau `$ as shown in Fig. 4(b). $`\alpha `$ characterises the common functional form of the distributions $`P(Z)`$ viz:

1. $`0<\alpha <2`$ Lévy (power-law with $`\sigma \mathrm{}`$)
2. $`\alpha =2`$ Gaussian (finite $`\sigma `$)
3. $`\alpha >2`$ power-law with finite $`\sigma `$

where $`P(Z)|Z|^{(1+\alpha )}`$ for $`Z\pm \mathrm{}`$ (Sornette sorn (2000)).

Thus the differencing and rescaling procedure not only reveals the timescales over which physical processes occur, but can also confirm any correlation suggested by the growth of range and inverse power-law form of the power spectra. Moreover, it quantifies the asymptotic behaviour of the distribution of fluctuations, which is essential for constraining turbulence/instability models (Bohr et al. bohr (1998)).

Differencing was used by Mantegna & Stanley (diff (1995)) to investigate fluctuations in the value of a financial index, but neither this technique nor the growth of range has yet (to our knowledge) been applied to astrophysical X-ray sources.

## 3 Results

### 3.1 The Crab

Fig. 5 shows the PDF of the total X-ray count-rate $`y\left(t\right)`$ that form the time-series for the Crab (Fig. 1).

The curve is close to Gaussian but with longer tails. Growth of the range is as low as that of Gaussian noise ($`H0.2`$, Fig. 6); the slope of the power spectrum is better at detecting weak correlation (see Sect. 2.2) and gives $`0.26\beta 0.28`$ at the $`95\%`$ confidence level.

Interestingly, the PDFs of the differenced data (Fig. 7) require no rescaling.

That is, not only is the type of distribution independent of $`\tau `$, but the differences themselves show no spread over time. We infer that the raw time-series is uncorrelated; if the correlation suggested by the power spectrum exists, it is too weak to be detected by this method.

### 3.2 Cygnus X-1

Fig. 8 shows the PDF of the raw count-rates from Cygnus X-1 (Fig. 2) over a continuous, relatively quiescent three-year period.

In contrast to Fig. 5, this PDF is clearly non-Gaussian. It appears possible to characterise some aspects of its non-Gaussian statistical properties in terms of, for example, the Gumbel and Fréchet distributions whose properties we outline in Sect. 4. Meanwhile, we note from Fig. 8 that the distribution of small amplitude events appears to fit the left-hand tail of Gumbel distributions with $`a=[1,2]`$ while the distribution of large amplitude events appears to fit the right-hand tail of a Fréchet distribution with $`a=1.1`$. This suggests that the total flux has contributions from different physical processes, perhaps arising in different parts of the accretion disk and its surroundings, whose individual PDFs could differ. Fitting to each of the many component distributions is problematic; here we simply note that the PDF in Fig. 8 could not be reproduced by summing Gaussian PDFs.

Cygnus X-1 has a substantially higher growth of range ($`H0.35`$, Fig. 9) than the Crab, but still below that of a Gaussian random walk ($`H0.5`$). This is confirmed by the power spectrum whose slope gives $`1.00\beta 1.02`$ at the $`95\%`$ confidence level (see Sect. 2.2).

This higher growth is evident in the unscaled PDFs of the differenced Cygnus X-1 data in Fig. 10.

By deriving a scaling exponent $`\alpha `$ from the peaks of Fig. 10, and using it to rescale as in Fig. 11, it is clear that the increments scale remarkably well over the full range of the data (four decades in time). This establishes both the existence of correlation in the X-ray output, and variations controlled by one type of process, on timescales up to three years.

However, with its slightly higher peak and broader tails, the PDF is distinctly non-Gaussian. In fact, the raw X-ray time-series of Cygnus X-1 is correctly described as a weakly leptokurtic random walk; that is, the PDF of the increments is long-tailed (Bouchard & Potters fin (2000)). The asymptotic form of the increments $`P(Z)`$ is quantified by $`\alpha `$; in this case, $`5.4\alpha 9.9`$ with $`95\%`$ confidence.

### 3.3 GRS1915+105

Fig. 12 shows the PDF of the raw data from GRS1915+105 (Fig. 3).

Like Cygnus X-1, GRS1915+105 is better described in terms of non-Gaussian statistics than Gaussian. For example, the extreme left- and right-hand tails fit Fréchet and Gumbel curves respectively, the reverse situation of Cygnus X-1. Note also that there are several peaks, suggesting a multi-component source involving at least two physical processes with different statistical properties. With $`H0.3`$ (Fig. 13), midway between Gaussian noise ($`H0.2`$) and a Gaussian random walk ($`H0.5`$), the growth of range of GRS1915’s count-rate lies between those of the Crab and Cygnus X-1. This is confirmed by the slope of the power spectrum: $`0.91\beta 0.93`$ at the $`95\%`$ confidence level (see Sect. 2.2).

The PDFs of the increments do not rescale over the full four decades, see Fig. 14.

The separate PDFs rescale approximately onto the same curve for $`\tau `$ up to only 1.5 decades (see Fig. 15) with stronger leptokurtosis than for Cygnus X-1.

However, we cannot unambiguously identify scaling in this régime given the variation in sampling interval; see Appendix.

## 4 Non-Gaussian statistical properties

The tails of the PDFs of Cygnus X-1 and GRS1915+105 are clearly non-Gaussian. Here, we compare these PDFs with those arising from extremal statistics. The two limiting distributions of interest are “Gumbel’s asymptote” and Fréchet (Fisher & Tippett fisher (1928); Gumbel gum (1958); Sornette sorn (2000)). In outline, the limiting distributions result from selecting the maximum value $`y_{max}`$ from each of a large number of large samples whose individual members are drawn from a distribution $`P(y)`$. When $`P(y)`$ decreases more rapidly than any power-law (as $`y\mathrm{}`$), “Gumbel’s asymptote” has the form

$$P_G(y_{max})=K(e^{ue^u})^a\mathrm{with}u=b(ys)$$
(4)
where in the limit of an infinite number of measurements $`a1`$; the constants $`K`$, $`b`$, and $`s`$ are fixed by normalisation as in Sect. 2.1.

Fréchet distributions $`P_F(y_{max})`$ arise in the same manner when the underlying PDF $`P(y)`$ is power-law. Mathematically, $`P_F(y_{max})`$ is defined by Eq. 4 but with $`u=\alpha +\beta \mathrm{ln}(1+y/G)`$, where $`K`$, $`b`$, $`s`$, $`\alpha `$, and $`G`$ are again fixed by normalisation as in Sect. 2.1, and $`\beta =(1a)^1`$. These curves exist for $`1<a<3/2`$.

Physically, the fact that the tails of the Cygnus X-1 and GRS1915+105 data can be fitted to Gumbel and Fréchet distributions may suggest that the observed signals have the character of maximal events. In this case, they would be the brightest among multiple events (whose PDF could be power-law) occuring within each observational time window. Interestingly, extremal statistics in a global measure are found in turbulent fluids and other highly-correlated systems (Bramwell et al. bram1 (1998, 2000, 2001)). In these cases, anomalous values of $`a`$ are found for Gumbel’s asymptote, and we plot $`a=1`$ and $`2`$ along with the Fréchet curve $`a=1.1`$ for comparison with the X-ray data. Also, since our global quantity is emitted flux rather than absorbed power, our curves show the opposite handedness to the results of Bramwell et al. (bram1 (1998)).

We also plot the log-normal PDFs having the same $`y`$ and $`\sigma `$ as the data, and these curves fit as closely as extremal distributions (see Figs. 8 and 12). The significance and origin of extremal and log-normal distributions is currently of considerable interest in statistical physics and turbulence studies (Bramwell et al. bram3 (2001); Burlaga burlaga (2001)).

## 5 Conclusions

Using three key statistical methods – the PDF, the growth of range, and differencing and rescaling – we have identified and quantified some fundamental contrasts in the character of the total X-ray output of the Crab, Cygnus X-1, and GRS1915+105 over multi-year time intervals. The Crab shows near-Gaussian behaviour with its low Hurst exponent and lack of temporal scaling. This is to be expected since the flux emanates from a very wide region over which correlation is unlikely, there being no observational evidence of an accretion disk in this source. In contrast, Cygnus X-1 has a time-series with scaling over three years, and resembles a random walk with leptokurtic increments whose PDF satisfies $`P(Z)|Z|^{(1+\alpha )}`$ for $`Z\pm \mathrm{}`$ and $`\alpha 8`$. Its PDF closely follows multi-component Gumbel and Fréchet curves, suggestive of the dominance of maximal events. Similarly, the output from GRS1915+105 lies close to such curves and its increments are also leptokurtic, although the roles of Fréchet and Gumbel components in the tails are reversed. However, as suggested by a lower Hurst exponent, results from this source are consistent with only short-range scaling with $`\alpha 6`$, and more data are required to establish this finding. Thus we have evidence that the two accreting objects display a degree of correlation in their X-ray time-series, which is absent from the nonaccreting Crab. This is a quantitative, observational, and model-independent measure of anomalous (non-diffusive) transport in accretion disks.

###### Acknowledgements.

We are grateful to John Kirk, Michel Tagger, and Nick Watkins for helpful suggestions. JG acknowledges a Research Studentship and SCC a Lecturer Fellowship from the UK Particle Physics and Astronomy Research Council. This work was also supported in part by the UK DTI and Euratom. S.C. acknowledges support from grant F/00-180/A from the Leverhulme Trust. Data provided by the ASM/RXTE teams at MIT and at the RXTE SOF and GOF at NASA’s GSFC.

## Appendix: Confidence limits for $`\tau `$