text:
# CO on Titan: More Evidence for a Well-Mixed Vertical Profile

## 1. Introduction

The atmosphere of Titan exhibits a complex photochemistry, and many nitriles and hydrocarbons have been detected by Voyager spacecraft and from Earth. Until recently, however, only two oxygen-bearing species had been detected on Titan: $`CO_2`$ (observed by Voyager 1; Samuelson et al.  (1983)) and CO (observed from Earth; Lutz et al.  (1983)).

The presence of oxygenated molecules is interesting because the atmosphere of Titan is strongly reducing. The cold temperatures of the lower stratosphere and the troposphere imply that $`CO_2`$ condenses out of the lower atmosphere and is continuously deposited on the surface. To sustain the carbon dioxide abundance a source of oxygen is needed, and it is generally assumed to be supplied in water from bombardment of the upper atmosphere by icy grains. In this model vaporized water is quickly photolyzed to produce OH, and OH reacts with hydrocarbon radicals such as $`CH_3`$ to produce CO. CO in turn reacts with OH to produce $`CO_2`$ (Samuelson et al.  (1983), Yung et al.  (1984), Toublanc et al.  (1995), Lara et al.  (1996)). While $`CO_2`$ has a short lifetime (order 103–104 years), the photochemical lifetime of $`CO`$ in the atmosphere of Titan is estimated to be very long ($`10^9`$ years; Yung et al.  (1984), Chassefière and Cabane  (1991)).

Observationally the missing piece of the oxygen chemistry has been the source, water. Recently, water vapor was detected in the upper atmosphere of Titan by the Short Wavelength Spectrometer (SWS) aboard the Infrared Satellite Observatory (Coustenis et al.  (1998)). With observations of the three major components of oxygen chemistry, it is now possible to check the internal consistency of photochemical models, and to compare the oxygen chemistry and water infall rate of Titan with the other giant planets, particularly Saturn (Feuchtgruber et al.  (1997), Coustenis et al.  (1998)).

Understanding the oxygen chemistry relies on accurate knowledge of the abundance and distribution of each species. A longstanding discussion regarding the CO distribution in Titan’s atmosphere, spanning more than a decade, has been primarily directed toward determining if CO is well-mixed (Marten et al.  (1988), Gurwell and Muhleman (1995), Hidayat et al.  (1998)). Since the residence lifetime of CO is long compared to transport timescales, the molecular weight of $`CO`$ is the same as for the dominant $`N_2`$ gas, and the atmosphere is never cold enough for CO to condense, carbon monoxide should be uniformly mixed in the Titan atmosphere to high altitudes.

Observational data, however, give conflicting results. Table 1 provides data on the CO abundance as measured by ground-based observers over the past 17 years. These observations have been sensitive to either the troposphere (near- and mid-IR) or the stratosphere (millimeter). The data in Table 1 show that no clear consensus has emerged regarding the CO abundance, either in the troposphere or the stratosphere.

In this Note, we present an analysis of new interferometric observations of the $`CO(21)`$ line on Titan. The results of this study have important implications for our understanding of the oxygen budget and photochemistry of the stratosphere of Titan.

## 2. Observations and Data Reduction

Observations of the $`CO(21)`$ rotational transition (rest frequency $`\nu _0`$=230.5380 GHz) on Titan were made on November 11 and 12, 1999 with the Owens Valley Radio Observatory Millimeter Array, located near Big Pine, California. The Titan ephemeris data was generated using the Jet Propulsion Laboratory’s Horizons on-line system (Giorgini et al.  (1996)). Titan was approaching eastern elongation with respect to Saturn, with a separation increasing from $``$145′′ to more than 200′′ over the two day period.

The interferometer was aligned in a fairly compact configuration, providing a synthesized beam of roughly 2″$`\times `$2.5″ at the observing frequency, while Titan’s apparent surface diameter was 0.86′′ (at a distance of 8.2165 AU). Titan was observed on each night over a period of about 7 hours, when it was above 30 elevation. A single measurement on each baseline consisted of a three minute integration during which the complex visibility of the source was recorded. Amplitude and phase gain variations were monitored through observations of 0235+164 approximately every 20 minutes, and antenna pointing was checked about every two hours using 3C84. The total integration time spent on Titan equaled 238 minutes on each night.

The signal was detected for each antenna pair in two correlator systems: a wide-band analog cross correlator ($``$1 GHz bandwidth) and a digital spectrometer. The CO line is significantly wider than this system bandwidth, and we utilized two local oscillator tunings to provide better coverage of the line: on November 11 the digital spectrometer measured (in two secondary LO tunings) the line in the upper sideband from $`\mathrm{\Delta }\nu =656`$ MHz $`+32`$ MHz, and on November 12 from $`\mathrm{\Delta }\nu =32`$ MHz $`+656`$ MHz. Spectra in the image sideband, approximately 3 GHz lower in frequency, were also recorded. The sideband signals were isolated to better than 20 dB using a phase-switching cycle. The combination of two first and second LO tunings allowed us to ultimately measure $`\pm `$650 MHz of the center frequency of each sideband at 4 MHz resolution, and $`\pm `$16 MHz at 0.5 MHz resolution in the line core.

Calibration of the digital correlator passband was done through observations of 3C273 and an internal correlated noise source. The relative calibration of the sidebands was accurately measured from observations of Uranus, 3C273, 3C84, and 0235+164 (to $`1`$%), since nearly all weather and instrumental effects impact each sideband in a similar manner. We note that the ability to isolate the astronomical signal sidebands and record independent spectra in each sideband represents a significant advantage for interferometric relative to single-dish observations of CO on Titan because the emission line is significantly broader than the spectrometer bandwidth. In this case the most precise measure of the line-to-continuum (LTC) ratio is provided by separating the sidebands. This relative sideband calibration then allows for the production of an accurate LTC spectrum, with one sideband sensing the continuum, and the other the line.

The Titan signal strength was sufficient for the application of phase self-calibration (see Thompson, Moran, and Swenson (1986)) to remove atmospheric phase variations, which cause decorrelation of the signal, on timescales shorter than the standard calibration cycle. After calibration, a complete spatially unresolved (e.g. “zero-spacing”) spectrum for each day was obtained by fitting the observed complex visibilities for each channel with a model of the Titan visibility function, correcting for the spatial sampling of the interferometer. The absolute flux scale was provided by scaling the continuum sideband intensity to equal the radiative transfer model flux of Titan at 227.5 GHz, corrected for the date and time of the observations (see below). The same scaling factor was applied to the emission line sideband, preserving the relative calibration.

The resulting combined spectrum of the $`CO(21)`$ line on Titan is shown in Figure 1. The data clearly shows the $`CO(21)`$ line is a strong emission feature in the spectrum of Titan. The image sideband spectrum is essentially flat except for a weak emission line due to the $`HC_3N(2524)`$ rotational transition at 227.419 GHz. This is particularly important because it shows that the sideband isolation procedure was effective to well below the noise level of the spectrum.

## 3. Modeling & Analysis

The radiative transfer model used to analyze the new CO data is nearly identical to the one discussed in Gurwell and Muhleman (1995), and we only highlight important aspects.

The basic parameters of the Titan atmosphere were derived from revised Voyager 1 radio occultation results (Lindal et al.  (1983), Lellouch (1990), Coustenis and Bézard (1995)), including an atmospheric base at 2575 km from the center of Titan, with a surface pressure and temperature of 1460 millibar and 96.7 K. For the thermal profile of the atmosphere we used an equatorial profile determined by Coustenis and Bézard (1995, their profile A) based upon the occultation results and Voyager 1/IRIS spectra (Fig. 2) combined with model J of Yelle (1991) for the upper atmosphere; this same model was used by Hidayat et al. (1998) to analyze their results. This model is appropriate since the observations reported here are unresolved (whole-disk) spectra, which are heavily weighted by emission from equatorial and low latitudes.

The millimeter continuum opacity on Titan is due to collision induced dipole absorption by $`N_2N_2`$, $`N_2`$Ar, and $`N_2`$CH4, and was modeled according to the results of Courtin (1988). The spectroscopic parameters for the $`CO(21)`$ line were taken from the JPL catalog (Pickett et al.  (1992); see also http://spec.jpl.nasa.gov). The full Voigt lineshape profile calculation using a fast computational method (Hui et al.  (1978)), integrated in pressure over atmospheric layers of constant temperature, was done using a collisional line-broadening coefficient for $`CO(21)`$ in $`N_2`$ of $`\gamma =2.21(T/300K)^{0.74}`$ MHz mbar-1 (Semmoud-Monnanteuil and Colmont (1987)). Radiative transfer calculations at appropriate frequencies were performed for a variety of radial steps, including limb-sounding geometries, and integrated over the apparent disk to provide the model whole-disk spectrum.

The contribution functions ($`W(z)=e^\tau d\tau /dz`$) for several frequency offsets from the $`CO(21)`$ line center are shown in Fig. 2, for a single raypath at the disk center. This function describes the relative contribution of different regions of the atmosphere to the emitted radiation at each frequency. The plotted functions assume a CO abundance of 50 ppm, constant with altitude. The $`\mathrm{\Delta }\nu =3000`$ MHz contribution function corresponds to the middle of the continuum (lower) sideband, and is dominated by the collision induced opacity of $`N_2`$. The other functions correspond to the emission (upper) sideband, and are dominated by $`CO(21)`$ opacity. The full line senses the atmosphere from 40 km (the tropopause) to 400 km. However, the range from 200 to 400 km is sounded mostly in the inner 4 MHz of the line core. At 4 MHz spectral resolution we are limited to sensing the CO abundance from the tropopause to $``$200 km. The 0.5 MHz spectrum of the line core pushes this upper bound to near 350 km in the absence of noise. Thermal noise on the spectral measurements in practice limit our sensitivity to $`300`$ km.

### 3.1. Best-fit Uniform CO Distribution

The radiative transfer model was run for a series of uniform CO distributions from $`q`$(CO)=10 to 90 ppm, in steps of 10 ppm. Resulting spectra are shown in Fig. 1 (in steps of 20 ppm for clarity). The model spectra have been convolved to the measurement spectral resolution of the data in each panel. The model calculations show that the continuum (lower) sideband emission is essentially unaffected by the CO distributions considered, and is an excellent continuum measurement. The model gives a flux at 227.5 GHz of 1.565 Jy for the geometry of the observations, equal to a disk-average Rayleigh-Jeans brightness temperature of 71.4 K.

Even by eye, the 50 ppm uniform model provides an exceptionally good fit to the data at both resolutions. The 50 ppm model gives an RMS residual of 86.8 mJy for the 4 MHz data, with models of 40 and 60 ppm giving RMS residuals that are factors of 1.4 and 1.1 times larger, respectively. Given that a large number of channels are involved (324), even an 10% increase in RMS residuals is quite significant. The 0.5 MHz spectrum is also consistent with this model. A rigorous least-squares analysis for the best-fit uniform profile gives a formal solution of 52$`\pm `$2 ppm from 40 to 300 km.

### 3.2. Best-fit Non-Uniform CO Distribution

An iterative least-squares inversion algorithm (following Gurwell and Muhleman (1995)) based on the radiative transfer model was utilized to solve for a best-fit non-uniform CO distribution. The logarithm of the CO distribution was constrained to be a linear function of altitude. This constrained solution tests whether a gradient in the CO distribution is consistent with the observed spectrum.

We find that the best-fit non-uniform profile, with formal error, is 48$`\pm `$4 ppm at 40 km, rising to 60$`\pm `$10 ppm at 300 km. The RMS residual is 86.1 mJy, representing less than 1% improvement in the residual over the best-fit uniform profile.

### 3.3. Error Estimates and the Best-fit CO Distribution

The formal errors quoted in the above sections are the direct results of the least-squares analyses, and therefore do not take into account errors in the radiative transfer modeling or the calibration of the spectrum.

To test whether the continuum emission model is a serious source of error (since the spectrum is calibrated by referencing to the continuum sideband data), we recomputed the continuum emission at 227.5 GHz, scaling the collision induced dipole absorption calculated from the data of Courtin (1988) by factors of 0.5 and 2. The calculated emission results were indistinguishable from the nominal model, which can be explained as the result of two factors. First, the collision induced continuum absorption scales as the square of pressure, and is therefore a very steep function of altitude. Therefore, increasing (or decreasing) the absorption coefficient even by factors of two will only increase the peak of the contribution function by a small fraction of a scale height. Second, the peak of the contribution function is right at the tropopause, where the temperature gradient is near zero. The result is that the emission change is very small, and we estimate that this error is about 1%.

The spectrum sidebands were calibrated assuming the QSO calibration sources had a spectral index of -0.5 (e.g. flux $`\nu ^{0.5}`$). However, the spectral index of these types of sources vary over the range of 0 to -1, and could lead to a calibration error of approximately 1% in the relative calibration over the 3 GHz difference in the sidebands.

Adding the calibration errors in quadrature, we find an error in the relatively calibrated spectrum of about 1.4%. Using the uniform distribution models discussed in section 3.1, we find that a 3% error in the relative calibration could lead to an error of roughly 10 ppm in a worst-case situation (we note that this does not include a refitting of the lineshape, which would tend to reduce this error; hence this is a worst-case estimate). The calibration error is then about 6 ppm using this scaling. For the uniform model solution, the formal error is significantly smaller than this calibration error estimate, and we believe that the error of our measurement is therefore about 6 ppm.

Noting that the non-uniform solution only improves the RMS residuals by 1% at best, and that the formal errors on the non-uniform solution encompass our uniform solution, we favor the uniform model for the CO distribution, which is in agreement with the current understanding of the chemistry of CO in the atmosphere of Titan. The high resolution data provides the information on altitudes above 200 km, and as can be seen in Fig. 1 this data has a higher RMS noise (by a factor of $``$2); this increases our error estimate by a factor of about two over this altitude range. We therefore find that the $`CO(21)`$ spectrum is best fit by a uniform profile of 52 ppm, with estimated errors of 6 ppm (40 to 200 km) and 12 ppm (200 to 300 km).

## 4. Discussion

The results presented here are nearly identical to our previous estimate of the CO distribution based on observations of the $`CO(10)`$ transition (Gurwell and Muhleman (1995)) and consistent with the original measurement of tropospheric CO (Lutz et al.  (1983)). Taken together, these measurements suggest a vertical profile of CO that is constant with altitude, at about 52 ppm, from the surface to at least 300 km.

These results are at odds with the recent measurements of Noll (1996), who found a tropospheric abundance of 10 ppm, and Hidayat et al. (1998), who found a stratospheric CO abundance of around 27 ppm (Table 1). Noll (1996) explored the possibility that their simple reflecting layer was not the surface, but a higher altitude ’haze’ layer. If the reflecting layer was at 0.9 bar (14 km) the spectrum was best fit with a CO abundance of 60 ppm. However, based on other evidence they found this model less satisfactory than a surface reflecting layer. The results of Hidayat et al. come from an analysis of several lines of CO, including the $`CO(10)`$ and $`CO(21)`$ lines; the discrepancy between their results and ours does not appear to be due to differences in modeling the atmosphere of Titan, but derives from differences in the measurement techniques and the resulting calibrated spectra (A. Marten, personal communication). However, we do point out that the interferometric method does offer advantages over single-dish observations for measuring the very broad lines of CO from the atmosphere of Titan.

We find the model of a uniform distribution of CO in the atmosphere of Titan provides a good fit to our data, but we cannot rule out a difference between the tropospheric and stratospheric CO abundance, since our data is insensitive to the lower atmosphere. A final confirmation of the abundance of CO and its vertical distribution requires further near- and mid-IR measurements of CO in the troposphere.

Acknowledgements

This work was supported in part by NASA grant NAG5-7946.

TABLE 1. Observations of CO in Titan’s Atmosphere

| Altitude | Mixing ratio (ppm)a | Wavelength | Reference |
| --- | --- | --- | --- |
| Troposphere | 48$`{}_{32}{}^{}{}_{}{}^{+100}`$ | 1.57$`\mu m`$ | Lutz et al.  (1983) |
| Stratosphere | 60$`\pm `$40b | 2.6 mm | Muhleman et al.  (1984) |
| Stratosphere | 2$`{}_{1}{}^{}{}_{}{}^{+2}`$ | 2.6 mm | Marten et al.  (1988) |
| Stratosphere | 50$`\pm `$10 | 2.6 mm | Gurwell and Muhleman (1995) |
| Troposphere | 10$`{}_{5}{}^{}{}_{}{}^{+10}`$ | 4.8$`\mu m`$ | Noll et al.  (1996) |
| Stratosphere | 27$`\pm `$5c | 2.6, 1.3, 0.9 mm | Hidayat et al.  (1998) |
| Stratosphere | 52$`\pm `$6 | 1.3 mm | this work |
| aMixing ratio defined as N(CO)/N(Total), i.e. not referenced to $`N_2`$. | | | |
| bReanalyzed by Paubert et al. (1984): 75$`{}_{45}{}^{}{}_{}{}^{+105}`$ ppm | | | |
| cNon-uniform model: 29$`\pm `$5 ppm (60 km), 24$`\pm `$5 ppm (175 km), 4.8$`\pm `$2 ppm (350 km) | | | |

Figure Captions


format: text

text:
# Growth of sumsets in abelian semigroups11footnote 1Supported in part by grants from the PSC–CUNY Research Award Program and the NSA Mathematical Sciences Program.

Let $`S`$ be an abelian semigroup, written additively, that contains the identity element 0. Let $`A`$ be a nonempty subset of $`S.`$ The cardinality of $`A`$ is denoted $`|A|.`$ For any positive integer $`h`$, the sumset $`hA`$ is the set of all sums of $`h`$ not necessarily distinct elements of $`A`$. We define $`hA=\{0\}`$ if $`h=0`$. Let $`A_1,\mathrm{},A_r,`$ and $`B`$ be nonempty subsets of $`S`$, and let $`h_1,\mathrm{},h_r`$ be nonnegative integers. We denote by

$$B+h_1A_1+\mathrm{}+h_rA_r$$
(1)
the set of all elements of $`S`$ that can be represented in the form $`b+u_1+\mathrm{}+u_r,`$ where $`bB`$ and $`u_ih_iA_i`$ for all $`i=1,\mathrm{},r.`$ If the sets $`A_1,\mathrm{},A_r,`$ and $`B`$ are finite, then the sumset (1) is finite for all $`h_1,\mathrm{},h_r.`$ The growth function of this sumset is

$$\gamma (h_1,\mathrm{},h_r)=|B+h_1A_1+\mathrm{}+h_rA_r|.$$

For example, let $`S`$ be the additive semigroup of nonnegative integers $`𝐍_0`$, and let $`A_1,\mathrm{},A_r,`$ and $`B`$ be nonempty, finite subsets of $`𝐍_0`$, normalized so that $`0BA_1\mathrm{}A_r`$ and $`\mathrm{gcd}(A_1\mathrm{}A_r)=1`$. Let $`b^{}=\mathrm{max}(B)`$ and $`a_i^{}=\mathrm{max}A_i`$ for $`i=1,\mathrm{},r`$. Han, Kirfel, and Nathanson  determined the asymptotic structure of the sumset $`B+h_1A_1+\mathrm{}+h_rA_r`$. They proved that there exist integers $`c`$ and $`d`$ and finite sets $`C[0,c2]`$ and $`D[0,d2]`$ such that

$$B+h_1A_1+\mathrm{}+h_rA_r=C[c,b^{}+\underset{i=1}{\overset{r}{}}a_i^{}h_id]\left(b^{}+\underset{i=1}{\overset{r}{}}a_i^{}h_iD\right).$$

for $`\mathrm{min}(h_1,\mathrm{},h_r)`$ sufficiently large. This implies that the growth function is eventually a multilinear function of $`h_1,\mathrm{},h_r,`$ that is, there exists an integer $`\mathrm{\Delta }`$ such that

$$|B+h_1A_1+\mathrm{}+h_rA_r|=a_1^{}h_1+\mathrm{}+a_r^{}h_r+b^{}+1\mathrm{\Delta }$$

for $`\mathrm{min}(h_1,\mathrm{},h_r)`$ sufficiently large. The explicit determination of the sets $`C`$ and $`D`$ is a difficult unsolved problem in additive number theory. In the case $`r=1`$, it is called the linear diophantine problem of Frobenius. For a survey of finite sumsets in additive number theory, see Nathanson .

The theorem about sums of finite sets of integers generalizes to sums in an arbitrary abelian semigroup $`S.`$ We shall prove that if $`A_1,\mathrm{},A_r,`$ and $`B`$ are finite, nonempty subsets of $`S`$, then the growth function $`\gamma (h_1,\mathrm{},h_r)`$ is eventually polynomial, that is, there exists a polynomial $`p(z_1,\mathrm{},z_r)`$ such that

$$\gamma (h_1,\mathrm{},h_r)=|B+h_1A_1+\mathrm{}+h_rA_r|=p(h_1,\mathrm{},h_r)$$

for $`\mathrm{min}(h_1,\mathrm{},h_r)`$ sufficiently large. The case $`r=1`$ is due to Khovanskii . We use his method to extend the result to the case $`r2.`$ The idea of the proof is to show that the growth function is the Hilbert function of a suitably constructed module graded by the additive semigroup $`𝐍_0^r`$ of $`r`$–tuples of nonnegative integers.

We need the following result about Hilbert functions. Let $`R`$ be a finitely generated $`𝐍_0^r`$–graded connected commutative algebra over a field $`E`$. Then $`R=_{h𝐍_0^r}R_h`$. Suppose that $`R`$ is generated by $`s`$ homogeneous elements $`y_1,\mathrm{},y_s`$ with $`y_iR_{\delta _i}`$, that is, the degree of $`y_i`$ is $`\mathrm{deg}y_i=\delta _i𝐍_0^r`$. Let $`M`$ be a finitely generated $`𝐍_0^r`$–graded $`R`$–module. For $`h=(h_1,\mathrm{},h_r)𝐍_0^r`$, we define the Hilbert function

$$H(M,h)=dim_E\left(M_{(h_1,\mathrm{},h_r)}\right).$$

For $`z=(z_1,\mathrm{},z_r)`$, we define

$$z^h=z_1^{h_1}\mathrm{}z_r^{h_r}.$$

Consider the formal power series

$$F(M,z)=\underset{h𝐍_0^r}{}H(M,h)z^h.$$

Then there exists a vector $`\beta `$ with integer coordinates and a polynomial $`P(M,z)=P(M,z_1,\mathrm{},z_h)`$ with integer coefficients such that

$$F(M,z)=\frac{z^\beta P(M,z)}{_{i=1}^s(1z^{\delta _i})}.$$

(This is Theorem 2.3 in Stanley \[7, p. 33\]).

###### Theorem 1

Let $`A_1,\mathrm{},A_r,`$ and $`B`$ be finite, nonempty subsets of an abelian semigroup $`S`$. There exists a polynomial $`p(z_1,\mathrm{},z_r)`$ such that

$$|B+h_1A_1+\mathrm{}+h_rA_r|=p(h_1,\mathrm{},h_r)$$

for all sufficiently large integers $`h_1,\mathrm{},h_r.`$

Proof. For $`i=1,\mathrm{},r,`$ let

$$A_i=\{a_{i,1},\mathrm{},a_{i,k_i}\},$$

where

$$|A_i|=k_i1.$$

We introduce a variable $`x_{i,j}`$ for each $`i=1,\mathrm{},r`$ and $`j=1,\mathrm{},k_i`$. Fix a field $`E`$. We begin with the polynomial ring

$$R=E[x_{1,1},\mathrm{},x_{r,k_r}]$$

in the $`s=k_1+\mathrm{}+k_r`$ variables $`x_{i,j}`$. The algebra $`R`$ is connected since it is an integral domain (cf. Hartshorne \[2, Exercise 2.19, p. 82\]). For each $`r`$–tuple $`(h_1,\mathrm{},h_r)𝐍_0^r`$ we let

$$R_{(h_1,\mathrm{},h_r)}$$

be the vector subspace of $`R`$ consisting of all polynomials that are homogeneous of degree $`h_i`$ in the variables $`x_{i,1},\mathrm{},x_{i,k_i}`$. In particular, $`E=R_{(0,\mathrm{},0)}`$. Then

$$R=\underset{(h_1,\mathrm{},h_r)𝐍_0^r}{}R_{(h_1,\mathrm{},h_r)}.$$

The multiplication in the algebra $`R`$ is consistent with this direct sum decomposition in the sense that

$$R_{(h_1,\mathrm{},h_r)}R_{(h_1^{},\mathrm{},h_r^{})}R_{(h_1+h_1^{},\mathrm{},h_r+h_r^{})},$$

and so $`R`$ is graded by the semigroup $`𝐍_0^r`$.

Next we construct an $`𝐍_0^r`$–graded $`R`$–module $`M.`$ To each $`r`$–tuple $`(h_1,\mathrm{},h_r)𝐍_0^r`$ we associate a finite-dimensional vector space $`M_{(h_1,\mathrm{},h_r)}`$ over the field $`E`$ in the following way. To each element

$$uB+h_1A_1+\mathrm{}+h_rA_r$$

we assign the symbol

$$[u,h_1,\mathrm{},h_r].$$

Let $`M_{(h_1,\mathrm{},h_r)}`$ be the vector space consisting of all $`E`$–linear combinations of these symbols. Then

$$dim_EM_{(h_1,\mathrm{},h_r)}=|B+h_1A_1+\mathrm{}+h_rA_r|.$$
(2)
Let

$$M=\underset{(h_1,\mathrm{},h_r)𝐍_0^r}{}M_{(h_1,\mathrm{},h_r)}.$$

This is an $`𝐍_0^r`$–graded vector space over $`E.`$

To make $`M`$ a module over the algebra $`R,`$ we must construct a bilinear multiplication $`R\times MM.`$ We define the product of the variable $`x_{i,j}R`$ and the basis element $`[u,h_1,\mathrm{},h_r]M`$ as follows:

$$x_{i,j}[u,h_1,\mathrm{},h_r]=[u+a_{i,j},h_1,\mathrm{},h_{i1},h_i+1,h_{i+1},\mathrm{},h_r].$$

This makes sense since

$$uB+h_1A_1+\mathrm{}+h_iA_i+\mathrm{}+h_rA_r$$

and so

$$u+a_{i,j}B+h_1A_1+\mathrm{}+(h_i+1)A_i+\mathrm{}+h_rA_r.$$

This induces a well-defined multiplication of elements of $`M`$ by polynomials in $`R`$ since, if $`i<i^{}`$,

$`x_{i^{},j^{}}\left(x_{i,j}[u,h_1,\mathrm{},h_r]\right)`$
$`=`$ $`x_{i^{},j^{}}[u+a_{i,j},h_1,\mathrm{},h_i+1,\mathrm{},h_r]`$
$`=`$ $`[u+a_{i,j}+a_{i^{},j^{}},h_1,\mathrm{},h_i+1,\mathrm{},h_i^{}+1,\mathrm{},h_r]`$
$`=`$ $`[u+a_{i^{},j^{}}+a_{i,j},h_1,\mathrm{},h_i+1,\mathrm{},h_i^{}+1,\mathrm{},h_r]`$
$`=`$ $`x_{i,j}[u+a_{i^{},j^{}},h_1,\mathrm{},h_i^{}+1,\mathrm{},h_r]`$
$`=`$ $`x_{i,j}\left(x_{i^{},j^{}}[u,h_1,\mathrm{},h_r]\right).`$
The case $`ii^{}`$ is similar. Note that this is the only place where we use the commutativity of the semigroup $`S.`$ It follows that $`M`$ is an $`R`$–module. Moreover,

$$R_{(h_1,\mathrm{},h_r)}M_{(h_1^{},\mathrm{},h_r^{})}M_{(h_1+h_1^{},\mathrm{},h_r+h_r^{})},$$

and so $`M`$ is a graded $`R`$–module. Furthermore, the finite set

$$\{[b,0,\mathrm{},0]:bB\}M$$

generates $`M`$ as an $`R`$–module.

Since $`x_{i,j}R_{\delta _{i,j}}`$, where $`\mathrm{deg}(x_{i,j})=\delta _{i,j}`$ is the $`r`$–tuple whose $`i`$–th coordinate is 1 and whose other coordinates are 0, and since

$$\frac{1}{(1z_i)^{k_i}}=\underset{h_i=0}{\overset{\mathrm{}}{}}\left(\genfrac{}{}{0pt}{}{h_i+k_i1}{k_i1}\right)z_i^{h_i},$$

we have

$`F(M,z)`$ $`=`$ $`{\displaystyle \underset{h𝐍_0^r}{}}H(M,h)z^h`$
$`=`$ $`{\displaystyle \frac{z^\beta P(M,z)}{_{i=1}^r_{j=1}^{k_i}(1z^{\delta _{i,j}})}}`$
$`=`$ $`{\displaystyle \frac{z^\beta P(M,z)}{_{i=1}^r(1z_i)^{k_i}}}`$
$`=`$ $`z^\beta P(M,z){\displaystyle \underset{i=1}{\overset{r}{}}}{\displaystyle \underset{h_i=0}{\overset{\mathrm{}}{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{h_i+k_i1}{k_i1}}\right)z_i^{h_1}`$
$`=`$ $`z^\beta P(M,z){\displaystyle \underset{h=(h_1,\mathrm{},h_r)𝐍_0^r}{}}{\displaystyle \underset{i=1}{\overset{r}{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{h_i+k_i1}{k_i1}}\right)z^h.`$
This implies that the Hilbert function $`H(M,h)`$ is a polynomial in $`h_1,\mathrm{},h_r`$ for $`\mathrm{min}(h_1,\mathrm{},h_r)`$ sufficiently large. By (2), the growth function is the Hilbert function of $`M`$. This completes the proof.


text:
# Black hole X-ray binaries: A new view on soft-hard spectral transitions

## 1 Introduction

For a decade it has been known that the spectra of X-ray novae show changes from a soft state at high luminosity to a hard state when the luminosity has declined during the outburst (Tanaka 1989). The persistent canonical black hole system Cyg X-1 also undergoes occasional transitions between its standard low luminosity (hard) state and a soft state (see Fig. 1). Such changes between the two spectral states have been observed for several systems, regardless of whether the compact object is a neutron star (Aql X-1, 1608-522) or a black hole (GS/GRS 1124-684, GX 339-4) (Tanaka & Shibazaki 1996). Here we concentrate on black hole sources. Observations show that the phenomenon always occurs at a luminosity around $`10^{37}\mathrm{erg}/\mathrm{s}`$, which corresponds to a mass accretion rate of about $`10^{17}\mathrm{g}/\mathrm{s}`$ (Tanaka 1999).

The two spectral states are thought to be related to different states of accretion: (1) the soft spectrum originates from a thin disk which extends down to the last stable orbit plus a corona above the disk, (2) the hard spectrum originates from a thin disk outside a transition radius $`r_{tr}`$ and a coronal flow/ ADAF inside. The spectral transitions of Nova Muscae 1991 and Cygnus X-1 were modelled based on this picture by Esin et al. (1997, 1998). The value of $`r_{tr}`$ was taken as the maximal distance $`r`$ for which an ADAF with that accretion rate can exist (“strong ADAF proposal”, Narayan & Yi 1995). We determine the location of the inner edge of the thin disk from the equilibrium between it and the corona above.

## 2 Generation of the coronal flow

### 2.1 Evaporation

The equilibrium between the cool accretion disk and the corona above (Meyer & Meyer-Hofmeister 1994) is established in the following way. Frictional heat released in the corona flows down into cooler and denser transition layers. There it is radiated away if the density is sufficiently high. If the density is too low, cool matter is heated up and evaporated into the corona until an equilibrium density is established (Meyer 1999).

Mass drained from the corona by an inward drift is replaced by mass evaporating from the thin disk as the system establishes a stationary state. When the evaporation rate exceeds the mass flow rate in the cool disk the disk terminates. Inside only a hot coronal flow exists.

### 2.2 Physics of the corona

Mass flow in the corona is similar to that in the thin disk. Differential (Kepler-like) rotation causes transfer of angular momentum outwards and mass flow inwards. The corona is geometrically much thicker than the disk underneath. Therefore sidewise energy transport is not negligible. Sidewise advection, heat conduction downward, radiation from the hot optically thin gas flow and wind loss are all important for the equilibrium between corona and thin disk. A detailed description would demand the solution of a set of partial differential equations in radial distance $`r`$ and vertical height $`z`$. In particular a sonic transition requires treatment of a free boundary condition on an extended surface.

From simplified modelling and analysis we find the following pattern of coronal flow. When a hole in the thin disk exists there are three regimes with increasing distance from the black hole. (1) Near the inner edge of the thin disk the gas flows towards the black hole. (2) At larger $`r`$ wind loss is important taking away about 20% of the total matter inflow. (3) At even larger distances some matter flows outward in the corona as a consequence of conservation of angular momentum. One might compare this with the flow in a “free” thin disk without the tidal forces acting in a binary. In such a disk matter flows inward in the inner region and outward in the outer region, with conservation of the total mass and angular momentum (Pringle 1981).

### 2.3 Model

We model the equilibrium between corona and thin disk in a simplified way. This is possible since the evaporation process is concentrated near the inner edge of the thin disk. Thus the corona above the innermost zone of the disk dominates the global structure. Further inward there is no thin disk anymore. The representative dominant region from $`r`$ to $`r`$+$`\mathrm{\Delta }r`$ has to be chosen such that evaporation further outward is not important. One incorporates the effects of frictional heat generation, conduction, radiation, sidewise loss of energy and wind loss at large height into this one zone ( “one-zone-model”). A set of ordinary differential equations for mass, motion, and energy with boundary conditions at the bottom (downward thermal flux - pressure relation) and at the top (sonic transition) uniquely determine mass accretion rate, wind loss and temperature in the corona as a function of radius. We restrict the analysis to a stationary corona.

The evaporation process was first investigated for disks in dwarf nova systems (Meyer & Meyer-Hofmeister 1994, Liu et al 1995). The situation is similar for disks around black holes (Meyer 1999). The coronal gas flowing into the hole and replaced by evaporation from the disk is understood as the supply for an ADAF which was used successfully to model the spectra of several black hole sources. A recent review by Narayan et al. (1998) gives a detailed description of accretion in the vicinity of a black hole.

## 3 Computational results

### 3.1 The critical mass flow rate $`\dot{M}_{\mathrm{crit}}`$

We use the same equations as Liu et al. (1995). The efficiency of evaporation at given distance $`r`$ from the compact star determines the location of the inner edge of the thin disk $`r_{\mathrm{tr}}`$. The relation between the mass flow rate $`\dot{M}`$ in the disk and $`r_{\mathrm{tr}}`$ was now computed also for black hole systems. In Fig. 2 we show this relation for a 6 $`M_{}`$ black hole (viscosity parameter $`\alpha `$=0.3).

Up to now only the decreasing branch was known and investigated. The interesting new feature is that the efficiency of evaporation reaches a maximum. This means that as the mass accretion rate in the disk is increased the inner edge moves inward, but if the rate exceeds a critical value $`\dot{M}_{\mathrm{crit}}`$ the thin disk can no longer be fully depleted by evaporation (for this accretion rate the inner disk edge is at about 340 Schwarzschild radii). The thin disk then extends inward towards the last stable orbit.

The temperature in the corona increases with decreasing radius, but reaches a saturation value where the coronal mass flow reaches maximum. The value h in Fig. 2 is the height where the pressure has decreased by 1/e.

### 3.2 What causes the maximum of the coronal mass flow rate?

A change in the physical process that removes the heat released by friction is the cause for the maximum of the coronal mass flow rate seen in Fig. 2. A dimensional analysis of the equations yields the following result. For large inner radii coronal heating is balanced by inward advection and wind loss. This fixes the coronal temperature at about 1/8 of the virial temperature $`T_\mathrm{v}`$ ($`\mathrm{}T_\mathrm{v}/\mu =GM/r`$, $`\mathrm{}`$ gas constant, $`\mu `$ molecular weight, $`G`$ gravitational constant) (see Fig. 2). Downward heat conduction and subsequent radiation in the denser lower region play a minor role for the energy loss though they always establish the equilibrium density in the corona above the disk.

With rising temperature, thermal heat conduction removes an increasing part of the energy released and finally becomes dominant. For optically thin bremsstrahlung the temperature saturates at a universal value defined by a combination of conductivity and radiation coefficients, the Boltzmann and the gas constant, and the non-dimensional $`\alpha `$-parameter of friction, (see Fig. 2). Dimensional analysis of the equations yields the rate of mass accretion through the corona as a function of temperature divided by the Kepler frequency $`(GM/r^3)^{1/2}`$. For small radii this gives the $`r^{3/2}`$ law in Fig. 2.

The maximum accretion rate occurs where the sub-virial temperature for large radii reaches the saturation temperature for small radii. Since the virial temperature is proportional to $`M/r`$, this radius $`r_{\mathrm{crit}}`$ is proportional to $`M`$. Then the accretion rate, proportional to the inverse of the Kepler frequency, also becomes proportional to $`M`$.

### 3.3 Approximations used for our model

Synchrotron and Compton cooling have been left out. Synchrotron cooling is non-dominant as long as the magnetic energy density stays below roughly 1/3 of the pressure. Compton cooling and heating by photons from the disk surface and from the accretion centre are non-dominant at all distances larger than that of the peak of the coronal mass flow rate, $`rr_{\mathrm{crit}}`$ ($`340R_s`$, Fig. 2). They become important for smaller radii.

The conductive flux remains small compared to the upper limit, the transport by free streaming electrons, so that classical thermal heat conduction is a good approximation. We have neglected lateral heat inflow by thermal conduction. This term is small compared to the dominant advective and wind losses at large radii, and vanishes when the temperature becomes constant at small radii.

Temperature equilibrium between electrons and ions requires that the collision times between them remains shorter than the heating timescale. This holds for $`rr_{\mathrm{crit}}`$, but the condition fails for $`r<r_{\mathrm{crit}}`$ where a two temperature corona can develop.

Tangled coronal magnetic fields could reduce electron thermal conductivity. We note however that reconnection and evaporation tend to establish a rather direct magnetic path between disk and corona.

## 4 Spectral transitions

### 4.1 Predictions from the evaporation model

At maximum luminosity of an X-ray nova outburst the mass accretion rate is high and the thin disk extends inward to the last stable orbit. A corona exists above the thin disk, but the mass flow in the thin disk is so high that no hole appears. In decline from outburst the mass accretion rate decreases. When $`\dot{M}_{\mathrm{crit}}`$ is reached a hole forms at $``$ 340 $`R_s`$ and the transition soft/hard occurs. If the mass accretion rate varies up and down as in high-mass X-ray binaries we expect hard/soft and soft/hard transitions. In Fig. 3 we show the expected behaviour schematically.

The descending branch for smaller $`r`$ indicates the possibility that an interior disk could form. We note that a gap exists between the exterior standard thin disk and an interior disk. In this gap the flow assumes the character of an ADAF with different temperature of ions and electrons, due to its high temperature and poor collisional coupling. This provides the possibility that the interior disk fed by this flow has a two temperature corona on top, different from a standard thin disk plus corona in the high state. We will discuss this in a further investigation.

### 4.2 Comparison with observations

The three persistent (high-mass) black hole X-ray sources LMC X-1, LMC X-3 and Cyg X-1 show a different behaviour. LMC X-1 is always in the soft state (Schmidtke et al. 1999). For LMC X-3, most of the time in the soft state, recently recurrent hard states have been detected (Wilms et al. 1999). Cyg X-1 spends most of its time in the hard state with occasional transitions to the soft state (see e.g. Fig. 1). This can be interpreted as caused by different long-term mean mass transfer rates: the highest rate (scaled to Eddington luminosity) in LMC X-1, the lowest in Cyg X-1, and in between in LMC X-3. Transient sources show a soft/hard transition during the decay from outburst. The best studied source is the X-ray Nova Muscae 1991 (Cui et al. 1997).

The transition always occurs around $`L_X10^{37}`$ erg/s (Tanaka 1999). Our value for the critical mass accretion rate for a 6 $`M_{}`$ black hole, $`10^{17.2}`$g/s, corresponds to a standard accretion disk luminosity of about $`10^{37.2}`$ erg/s. This is very close agreement.

For accretion rates below $`\dot{M}_{\mathrm{crit}}`$ the location of the inner edge of the standard thin disk derived from the evaporation model also agrees with observations (Liu et al. 1999).

At the moment of spectral transition our model predicts the inner edge near 340 Schwarzschild radii. The observed timescale for the spectral transition of a few days (Zhang et al. 1997) agrees with the time one obtains for the formation of a disk at 340 $`R_s`$ with an accretion rate $`\dot{M}_{\mathrm{crit}}`$.

But even in the low state X-ray observations of a reflecting component indicates the existence of a disk further inward, at 10 to 25 $`R_s`$ (Gilfanov et al. (1998), Zycki et al. (1999)). This might point to a non-standard interior disk as discussed above and explain why the spectral transitions in Cygnus X-1 could be well fitted by Esin et al. (1998) with a disk reaching inward to $`100R_s`$.

## 5 Conclusions

We understand the spectral transition as related to a critical mass accretion rate. For rates $`\dot{M}\dot{M}_{\mathrm{crit}}`$ (the peak coronal mass flow rate) the standard disk reaches inward to the last stable orbit and the spectrum is soft. Otherwise the ADAF in the inner accretion region provides a hard spectrum. At $`\dot{M}_{\mathrm{crit}}`$ the transition between dominant advective losses further out and dominant radiative losses further in occurs. Except for the difference between the sub-virial temperature of the corona and the closer-to-virial temperature of an ADAF of the same mass flow rate, this same critical radius is predicted by the “strong ADAF proposal” (Narayan & Yi (1995). In general however, the strong ADAF proposal results in an ADAF region larger than that which the evaporation model yields.

The transition between the two spectral states has been observed for black hole and neutron star systems, in persistent and transient sources (Tanaka & Shibazaki 1996, Campana et al.1998). This points to similar physical accretion processes. Menou et al. (1999) already discussed the accretion via an ADAF in neutron star transient sources. Our results should also be applicable.

The relations for a 6 $`M_{}`$ black hole plotted in Fig. 2 can be scaled to other masses: in units of Schwarzschild radii and Eddington accretion rates the plot is universal. The application to disks around supermassive black holes implies interesting conclusions for AGN.

###### Acknowledgements.

We thank Marat Gilfanov, Eugene Churazov and Michael Revnivtsev for the spectral data of Cygnus X-1.


text:
# RUNAWAY OF LINE-DRIVEN WINDS TOWARDS CRITICAL AND OVERLOADED SOLUTIONS

## 1 Introduction

Atmospheres of hot luminous stars and accretion disks in active galactic nuclei and cataclysmic variables form extensive outflows due to super-Eddington radiation fluxes in UV resonance and subordinate lines. An understanding of these winds is hampered by the pathological dependence of the driving force on the flow velocity gradient. Castor, Abbott, & Klein (1975; CAK hereafter) found that line-driven winds (hereafter LDWs) from O stars should adopt a unique, critical state which corresponds to maximum mass loss rate. The equation of motion for a 1-D, spherically symmetric, polytropic outflow subject to a Sobolev line force allows for two infinite families of so-called shallow and steep solutions. However, none of these families can provide for a global solution alone. Shallow solutions do not reach infinity, while steep solutions do not extend into the subsonic regime including the photosphere. The critical wind starts then as the fastest shallow solution and switches at the critical point in a continuous and differentiable manner to the slowest steep solution. Hence the critical point and not the sonic point determines the bottleneck in the wind. This description in principle applies equally to winds from stars and accretion disks.

A physical interpretation of the CAK critical point was given by Abbott (1980), who derived a new type of radiative-acoustic waves (hereafter Abbott waves). These waves can propagate inward, in the stellar rest frame, only from below the CAK critical point. Above the critical point, they are advected outwards. Hence, the CAK critical point serves as an information barrier, much as the sonic or Alfvén points in thermal and hydromagnetic winds. Abbott’s analysis was challenged by Owocki & Rybicki (1986) who found for a pure absorption LDW the signal speed to be the sound speed and not the much faster Abbott speed. As noted already by these authors, this should be a consequence of assuming pure line absorption, which does not allow for any radiatively modified, inward wave mode. Meanwhile there is ample evidence for Abbott waves in time-dependent wind simulations (Owocki & Puls 1999).

Shallow solutions fail to reach infinity because they cannot perform the required spherical expansion work, implying that the flow starts to decelerate. Since this usually occurs very far out in the wind, the local wind speed is much larger than the local escape speed, and the wind escapes to infinity. Thus, a simple generalization of the CAK model allowing for flow deceleration renders shallow solutions globally admissible. This raises a fundamental question of why the wind would adopt the critical solution at all, and attain the critical mass loss rate and velocity law, as proposed by CAK.

In this Letter we analyze a physical mechanism which drives shallow solutions towards the critical one, and discuss under what conditions this evolutions does not terminate at the CAK solution, but continues into the realm of overloaded solutions. We find that simulations so far were affected by numerical runaway towards the critical solution, by not accounting for Abbott waves in the Courant time step.

## 2 Abbott waves

Abbott waves are readily derived by bringing the wind equations into characteristic form. We consider a 1-D planar wind of velocity $`v(z,t)`$ and density $`\rho (z,t)`$, assuming zero sound speed. The continuity and Euler equation are,

$$\frac{\rho }{t}+v\frac{\rho }{z}+\rho \frac{v}{z}=0,$$
(1)
$$E\frac{v}{t}+v\frac{v}{z}+g(z)CF(z)\left(\frac{v/z}{\rho }\right)^\alpha =0.$$
(2)
Here, $`g(z)`$ and $`F(z)`$ are gravity and radiative flux, respectively. The CAK line force is given by $`g_\mathrm{l}CF(z)(v^{}/\rho )^\alpha `$ (with $`v^{}v/z`$), with constant $`C`$ and exponent $`0<\alpha <1`$. The unique, stationary CAK wind, $`v_\mathrm{c}(z),\rho _\mathrm{c}(z)`$, is found by requiring a critical point at some $`z_\mathrm{c}`$. The number of solutions for $`vv^{}(z)`$ changes from 2 to 1 at $`z_\mathrm{c}`$ (which is a saddle point), hence $`E/(vv^{})_\mathrm{c}=0`$ holds. Writing $`C`$ in terms of critical point quantities, the Euler equation becomes,

$$\frac{v}{t}+v\frac{v}{z}+g(z)$$

$$\alpha ^\alpha (1\alpha )^{(1\alpha )}\frac{F(z)}{F(z_\mathrm{c})}g(z_\mathrm{c})^{1\alpha }(\rho _\mathrm{c}v_\mathrm{c})^\alpha \left(\frac{v/z}{\rho }\right)^\alpha =0.$$
(3)
Note that for stationary planar winds, $`\rho v`$ is constant. If, in addition, $`g`$ and $`F`$ are taken constant with height, and $`\rho _\mathrm{c}v_\mathrm{c}v^{}/\rho `$ is replaced by $`vv^{}/\dot{m}`$, with normalized mass loss rate $`\dot{m}\rho v/\rho _\mathrm{c}v_\mathrm{c}`$, one finds that $`E`$ does no longer depend explicitly on $`z`$ for stationary solutions. Hence, $`vv^{}`$ is independent of $`z`$, too. This implies that $`z_\mathrm{c}`$ is ill-defined, and every point of the CAK solution is a critical point. CAK removed this degeneracy by introducing gas pressure terms. Here we take a different approach and assume $`g=z/(1+z^2)`$. A situation with roughly constant radiative flux and gravity showing a maximum at finite height could be encountered above isothermal disks around compact objects (cf. Feldmeier & Shlosman 1999). The critical point is determined by the regularity condition, $`dE/dz_\mathrm{c}=0`$, hence $`z_\mathrm{c}=1`$ and the critical point coincides with the gravity maximum. For simplicity also we chose $`\alpha =1/2`$ from now on, which is reasonably close to realistic values $`\alpha 2/3`$ (Puls et al. 1999). None of our results should depend qualitatively on the assumptions made so far. The Euler equation is

$$\frac{v}{t}+v\frac{v}{z}+g(z)2\sqrt{g_\mathrm{c}\rho _\mathrm{c}v_\mathrm{c}}\sqrt{\frac{v/z}{\rho }}=0,$$
(4)
where $`g_\mathrm{c}g(z_\mathrm{c})`$. The stationary solutions for wind acceleration are given by

$$vv^{}(z)=\frac{g_\mathrm{c}}{\dot{m}}\left(1\pm \sqrt{1\frac{\dot{m}g(z)}{g_\mathrm{c}}}\right)^2,$$
(5)
where plus and minus signs refer to steep and shallow solutions, respectively. For $`\dot{m}1`$, shallow and steep solutions are globally, i.e., everywhere, defined. For $`\dot{m}>1`$, solutions are called overloaded, and become imaginary in a neighborhood of the gravity maximum. These winds carry too large mass loss rates and eventually stagnate.

Next we put the Euler equation into quasi-linear form, which does not mean to linearize it. Differentiating $`E`$ with respect to $`z`$ (Courant & Hilbert 1962; Abbott 1980) and introducing $`fv/z`$, eqs. (1, 4) become,

$`\left[{\displaystyle \frac{}{t}}+v{\displaystyle \frac{}{z}}\right]\rho +\rho f=0,`$ (6)
$`\left[{\displaystyle \frac{}{t}}+(v+v_\mathrm{A}){\displaystyle \frac{}{z}}\right]{\displaystyle \frac{f}{\rho }}+{\displaystyle \frac{1}{\rho }}{\displaystyle \frac{g}{z}}=0,`$ (7)
with inward Abbott speed in the rest frame, $`v_\mathrm{A}\sqrt{g_\mathrm{c}v/\dot{m}v^{}}`$. In the WKB approximation, individual spatial and temporal variations are much larger than the inhomogeneous term $`g^{}/\rho `$ in eq. (7), and the latter can be neglected. Consequently, $`v^{}/\rho `$ is a Riemann invariant propagating at characteristic speed $`v+v_\mathrm{A}`$. Perturbations of $`v^{}/\rho `$ correspond to the amplitude of a wave propagating at phase speed $`v+v_\mathrm{A}`$. Note that $`v^{}/\rho `$ is proportional to the Sobolev line optical depth, indicating that this wave is a true radiative mode.

The second characteristic is determined by the continuity equation (6). In the advection operator in square brackets, $`v`$ has to be read as $`v+0`$ in the zero-sound speed limit. This outwards propagating invariant corresponds to a sound wave, with amplitude $`\rho `$ scaling with gas pressure.

At the critical point, $`\dot{m}=1`$ and $`vv^{}(z_\mathrm{c})=g(z_\mathrm{c})`$ after eq. (5), hence $`v_{\mathrm{Ac}}=v_\mathrm{c}`$ (where we introduced $`v_{\mathrm{Ac}}v_\mathrm{A}(z_\mathrm{c})`$). Abbott waves stagnate at the critical point, in analogy with sound waves at the sonic point. For shallow solutions, $`\dot{m}<1`$ and $`vv^{}<v_\mathrm{c}v_\mathrm{c}^{}`$ from eq. (5), hence $`v+v_\mathrm{A}<0`$. Shallow LDW solutions are therefore the subcritical analog to solar wind breezes.

Because in the rest frame, the inward Abbott mode can propagate at larger absolute speeds than the outward sound mode, Abbott waves can determine the Courant time step in time-explicit hydrodynamic simulations. Violating the Courant step results in numerical instability. Despite this fact, Abbott waves along shallow solutions were never considered in the literature.

## 3 Wind convergence towards the critical solution

We turn our attention to a physical mechanism which can drive LDWs away from shallow solutions, and towards the critical one. Starting from an arbitrary shallow solution as initial condition, we explicitly introduce perturbations at some fixed location in the wind and study their evolution. In order to keep unperturbed shallow solutions stable in numerical simulations, we fix one outer boundary condition, according to inward propagating Abbott waves. Either a constant mass loss rate at the outer boundary or non-reflecting boundary conditions (Hedstrom 1977) serve this aim. At the inner, subcritical boundary, we also fix one boundary condition, according to incoming sound waves. Non-reflecting boundary conditions and $`\rho =const`$ give similar results.

Wind convergence towards the critical solution is then triggered by negative flow velocity gradients. Allowing for $`v^{}<0`$ turns the inward Abbott mode of phase speed $`v+v_\mathrm{A}<0`$ in the rest frame into an outwards propagating mode. This is readily seen for a line force which is zero for negative $`v^{}`$, i.e., when all photons are absorbed at a resonance location between the photosphere and the wind point. The Euler equation simplifies to that for an ordinary gas, with characteristic speed $`v0>0`$ in the zero sound speed limit. At the other extreme, for a purely local line force where the unattenuated stellar or disk radiation field reaches the wind point, $`g_\mathrm{l}\sqrt{|v^{}|}`$. Here the Abbott phase speed is found to be $`v+v_\mathrm{A}`$, with $`v_\mathrm{A}=+\sqrt{g_\mathrm{c}v/\dot{m}v^{}}`$ for $`v^{}<0`$.

Consider then a sawtooth-like velocity perturbation (sinusoidal perturbation lead to similar results). Slopes $`v^{}>0`$ propagate inwards, slopes $`v^{}<0`$ propagate outwards. Hence, as a kinematical consequence, a sawtooth which is initially symmetric with respect to the underlying stationary velocity law evolves towards larger velocities. This is demonstrated in Figure 1 where, in course of time, a periodic sawtooth perturbation is introduced at $`z=2`$. The line force is assumed to be $`\sqrt{|v^{}|}`$, and the initial shallow solution has $`\dot{m}=0.8`$. The figure shows 2$`\frac{1}{2}`$ perturbation cycles. For upward pointing kinks the slopes propagate apart and a flat velocity law develops between them. At each time step $`dt`$, a new increment $`dv=4dt\delta v/T`$ ($`\delta v`$ and $`T`$ being the amplitude and period of the sawtooth) is added at $`z=2`$, hence the flattening velocity law does not show up in region A of Figure 1. Overall, the wind speed at the perturbation site evolves towards larger values during these phases. On the other hand, for downward pointing kinks of the sawtooth, $`\delta v`$, the two approaching slopes merge, and the wind speed evolves back towards its unperturbed value after each decrement $`dv=4dt\delta v/T`$. The wind velocity hardly evolves during these phases, cf. region B of Figure 1. Over a full perturbation cycle, the wind speed clearly increases.

Essentially any perturbation which introduces negative $`v^{}`$ will accelerate the wind. The amplitude of the perturbation is rather irrelevant since, with decreasing perturbation wavelength, negative $`v^{}`$ occur at ever smaller amplitudes. However, in more realistic winds, dissipative effects may smear out short-scale perturbations before they can grow. Details of the physical mechanism will be discussed elsewhere.

If the perturbation lies downstream from the critical point, the wind converges to the critical solution. Namely, as soon as the perturbation site comes to lie on the supercritical part of the CAK solution during its evolution, positive velocity slopes propagate outwards, and combine with negative slopes to a full wave train. No information is propagated upstream. This unconditional stability of the outer CAK solution is shown in the lower panel of Fig. 1.

## 4 Wind convergence towards overloaded solutions

Wind runaway towards larger speeds as caused by perturbations introduced upstream from the critical point does not terminate at the critical CAK solution. For low-lying perturbations, communication with the wind base is still possible once the subcritical branch of the CAK solution is reached. The wind gets further accelerated into the domain of mass-overloaded solutions (where $`vv^{}>v_\mathrm{c}v_\mathrm{c}^{}`$ and hence $`v>v_\mathrm{c}`$ for $`z<z_\mathrm{c}`$ according to eq. 5), until a generalized critical point develops, which prevents inward propagation of Abbott waves and adjustment of the mass loss rate. Such generalized critical points are given by ‘termination’ points, $`z_\mathrm{t}`$, of overloaded solutions, where the velocity becomes imaginary. At $`z_\mathrm{t}`$, the number of real solutions $`vv^{}(r)`$ changes from 2 (shallow and steep) to 0. Hence, termination points are defined by the same condition as the CAK critical point (at which the number of solutions changes from 2 via 1 to 2), $`E/(vv^{})_\mathrm{t}=0`$. From the stationary version of eq. (4), $`v_{\mathrm{At}}=v_\mathrm{t}`$, hence Abbott waves stagnate at termination points, and the latter are generalized critical points.

The fact that perturbations with negative $`v^{}`$ accelerate the wind either to the critical or an overloaded state can be casted into black-hole conjecture (Penrose 1965): a LDW avoids a ‘naked’ base, and encloses it with a critical surface.

Since to each $`z_\mathrm{t}`$ there corresponds a unique, supercritical mass loss rate, the latter is determined by the perturbation location alone. Using $`v_{\mathrm{At}}=v_\mathrm{t}`$, one finds $`\dot{m}_\mathrm{t}=g_\mathrm{c}/g_\mathrm{t}>1`$ for a planar wind with constant radiative flux.

At a termination point, $`vv^{}`$ jumps to the decelerating branch, $`vv^{}<0`$. Beyond a well-defined location above the gravity maximum, the super-CAK mass loss rate can again be lifted by the line force, and $`vv^{}`$ jumps back to the accelerating branch. Hence, two stationary kinks occur in the velocity law. Figure 2 shows a hydrodynamic simulation of the evolution towards an overloaded solution. Sawtooth-type velocity perturbations were introduced at $`z=0.8`$. Correspondingly, $`\dot{m}=1.025`$ for the overloaded solution, using $`g=z/(1+z^2)`$.

Future work has to clarify whether LDWs show deep-seated perturbations. It seems however unlikely that they would occur at a unique location. Hence, overloaded winds should be non-stationary and show a range of supercritical mass loss rates.

More fundamentally, time-dependent overloaded solutions occur already for single, unique perturbation sites, once the latter lie below a certain height. For the present wind model, this is at $`z0.66`$. The overloading is then so severe and the decelerating region so broad that negative wind speeds result (cf. Poe, Owocki, & Castor 1990). The corresponding mass loss rates are still only a few percent larger than the CAK value. The gas which falls back towards the photosphere collides with outflowing gas, and a time-dependent situation develops. Within each perturbation period, a shock forms in the velocity law, supplemented by a dense shell. These shocks and shells propagate outwards (Feldmeier & Shlosman 2000).

Although strong perturbations introducing negative velocity gradients can appear already in O star winds, accretion disk winds are the prime suspects. The reason for this is that accretion processes and their radiation fields in cataclysmic variables and galactic nuclei are intrinsically variable on a range of timescales (Frank, King, & Raine 1992), and that disk LDWs are driven by a combination of uncorrelated, local and central radiation fluxes.

## 5 Summary

We find that shallow solutions to line-driven winds are subcritical with respect to Abbott waves (sub-abbottic). These waves cause shallow solutions to evolve towards larger speeds and mass loss rates because of the asymmetry of the line force with regard to positive and negative velocity gradients and because perturbations with opposite signs of $`dv/dz`$ propagate in opposite directions. Steep velocity slopes propagate towards the wind base, steepen the inner wind and lift it to higher mass loss rates. In the presence of enduring wind perturbations, this proceeds until a critical point forms and Abbott waves can no longer penetrate inwards.

The resulting solution does not necessarily correspond to the CAK wind. For perturbations which originate below the critical point, the developing Abbott wave barrier is found to be the termination point of a mass-overloaded solution. The velocity law acquires a kink at the termination point, where the wind starts to decelerate. Whether the wind converges to a critical or overloaded solution depends entirely on the location of perturbations, and not, e.g., on boundary conditions at the wind base.

If Abbott waves are not accounted for in the Courant time step of hydrodynamic simulations, we find that numerical runaway can drive the solution towards the critical CAK wind. A detailed discussion of this will be given elsewhere.

Future work has to clarify whether and where perturbations causing local flow deceleration, $`dv/dz<0`$, can occur in LDWs. Overloaded winds may be detected observationally. While their mass loss rates should still be close to CAK values, broad regions of decelerating flow could be identified in P Cygni line profiles. Furthermore, shocks occurring in overloaded solutions with infalling gas may contribute to the X-ray emission from LDWs, besides shocks from the line-driven instability (Lucy 1982; Owocki et al. 1988). Note that the present wind runaway occurs already in the lowest order Sobolev approximation, and is therefore, unrelated to the line-driven instability which depends on velocity curvature terms (Feldmeier 1998).

###### Acknowledgements.

We thank R. Buchler, J. Drew, R. Kudritzki, C. Norman, S. Owocki, and J. Puls for intense blackboard discussions, and the referee, Stan Owocki, for suggestions improving the manuscript. This work was supported in part by PPA/G/S/1997/00285, NAG 5-3841, WKU-522762-98-6 and HST GO-08123.01-97A.


text:
# Estimation of initial conditions from a scalar time series

## Abstract

We introduce a method to estimate the initial conditions of a mutivariable dynamical system from a scalar signal. The method is based on a modified multidimensional Newton-Raphson method which includes the time evolution of the system. The method can estimate initial conditions of periodic and chaotic systems and the required length of scalar signal is very small. Also, the method works even when the conditional Lyapunov exponent is positive. An important application of our method is that synchronization of two chaotic signals using a scalar signal becomes trivial and instantaneous.

A trajectory that a given dynamical system traverses in its state space depends on the particular set of initial conditions with which it starts. In particular, the state of a chaotic system at a latter time is exponentially sensitive to changes in its initial state . This defining feature of a chaotic system leads to a complex behaviour in state space that appears random yet is deterministic, which means that an initial state uniquely fixes the future course of its evolution. Though their are several invariant measures of a chaotic system which are not sensitive to the intial conditions, the exact trajectory crucially depends on the intial state and hence is difficult to reproduce due to sensitivity to initial conditions.

In light of these facts, it is interesting and important to ask whether the complete set of initial conditions of a given multivariable dynamical system can be estimated from a given scalar time series for a single state space variable. We show that this question can be answered in the affirmative and present a novel and simple method to estimate the initial conditions. Our method is based on a modified multidimensional Newton-Raphson method  where we include the time evolution of the system. The length of the time series required for the calculations is typically very small.

Our results raise some interesting issues regarding the information content of a time series. In standard embedding techniques, a vector space is constructed from successive iterates of a single variable and a trajectory is reconstructed in this space . While embedding, it is crucial to choose an appropriate time delay so that the successive iterates are well resolved and contain qualitatively different information. In our method we use a very small time series and the total duration is typically much less than the standard delay time in embedding techniques. It is interesting that we can recover the initial conditions and hence the trajectory from such a small duration of time series.

An important application of our method is in the problem of synchronization of two chaotic systems . This problem itself has attracted wide attention in recent times due to its potential application to secured communication  and parameter estimation . Our estimation of the initial conditions makes the problem of synchronization almost trivial. We also find that our method works for most of the cases where other methods fail .

Let us consider an autonomous dynamical system given by,

$$\dot{𝐱}=𝐅(𝐱),$$
(1)
where $`𝐱=(x_1,x_2,\mathrm{},x_d)`$ is a $`d`$-dimensional state vector whose evolution is governed by the function $`𝐅=(F_1,F_2,\mathrm{},F_d)`$. Given an initial state vector $`𝐱(0)`$ at time $`t=0`$, the time evolution $`𝐱(t)`$ is uniquely determined by Eq. (1). Now let us assume that only one component of the state vector is known to us and we take it to be $`x_1(t)`$ without loss of generality. The problem that we address is to obtain the initial state vector $`𝐱(0)`$ from the knowledge of the scalar signal $`x_1(t)`$.

Let $`𝐲(0)`$ denote a random initial state vector and $`𝐲(t)`$ its time evolution obtained from Eq. (1). Let $`𝐰(t)`$ denote the difference

$$𝐰(t)=𝐲(t)𝐱(t).$$
(2)
We look for the solution of the equation

$$𝐰(t)=0.$$
(3)
Noting that the initial state vectors $`𝐲(0)`$ and $`𝐱(0)`$ uniquely determine the difference $`𝐰(t)`$, one of the solutions of Eq. (3) is $`𝐲(\mathrm{𝟎})𝐱(0)=0`$ and this is the solution that we are searching for.

We now introduce the notation $`𝐰^n=𝐰^n(𝐲^0,𝐱^0)=𝐰(n\mathrm{\Delta }t)`$, where $`\mathrm{\Delta }t`$ is a small time interval. Similarly, $`𝐲^n=𝐲(n\mathrm{\Delta }t)`$ and $`𝐱^n=𝐱(n\mathrm{\Delta }t)`$. With this notation condition (3) can be written as $`𝐰^n=0`$.

Our approach to the solution of Eq. (3) is a modified Newton-Raphson method  which includes the time evolution of the system.

Let us first consider $`𝐰^1`$. We have

$`0`$ $`=`$ $`𝐰^1(𝐱^0,𝐱^0),`$ (4)
$`=`$ $`𝐰^1(𝐲^0+\delta 𝐲^0,𝐱^0),`$ (5)
$`=`$ $`𝐰^1(𝐲^0,𝐱^0)+(\delta 𝐲^0\mathbf{}\mathbf{}_{𝐲^\mathrm{𝟎}})𝐰^1(𝐲^0,𝐱^0)+𝒪((\delta 𝐲^0)^2),`$ (6)
where $`\delta 𝐲^0=𝐱^0𝐲^0=𝐰^0`$ and the last step is a Taylor series expansion in $`\delta 𝐲^0`$. For small $`\mathrm{\Delta }t`$, we can write

$$𝐰^1(𝐲^0,𝐱^0)=𝐰^0+\mathrm{\Delta }t[𝐅(𝐲^0)𝐅(𝐱^0)]+𝒪((\mathrm{\Delta }t)^2).$$
(7)
Substituting Eq. (7) in Eq. (6) and neglecting higher order terms, we get,

$$𝐰^1(𝐲^0,𝐱^0)=𝐰^0\mathrm{\Delta }t(𝐰^0\mathbf{}\mathbf{}_{𝐲^\mathrm{𝟎}})𝐅(𝐲^0).$$
(8)
It is convenient to write the above equation in a matrix form as

$`W^1`$ $`=`$ $`(I+\mathrm{\Delta }tJ^0)W^0,`$ (9)
$`=`$ $`A^0W^0,`$ (10)
where $`W^n`$ is the column matrix corresponding to the vector $`𝐰^n`$, $`I`$ is the identity matrix, $`A^n=I+\mathrm{\Delta }tJ^n`$, and the elements of the Jacobian matrix $`J^n`$ are given by

$$J_{ij}^n=\frac{F_i(𝐲^n)}{y_j^n}.$$
(11)
Next we consider $`𝐰^2`$ or $`W^2`$. Proceeding as above, we get (see Eq. (10)),

$`W^2`$ $`=`$ $`(I+\mathrm{\Delta }tJ^1)W^1`$ (12)
$`=`$ $`(I+\mathrm{\Delta }tJ^1)(I+\mathrm{\Delta }tJ^0)W^0`$ (13)
$`=`$ $`A^1A^0W^0.`$ (14)
Similarly, the equation for $`W^n`$ is

$`W^n`$ $`=`$ $`(I+\mathrm{\Delta }tJ^{n1})(I+\mathrm{\Delta }tJ^{n2})\mathrm{}(I+\mathrm{\Delta }tJ^0)W^0`$ (15)
$`=`$ $`A^{n1}A^{n2}\mathrm{}A^0W^0.`$ (16)
We now concentrate on the first component of the signal whose time series is assumed to be known. For a $`d`$-dimensional system we need $`d1`$ equations to determine the initial state vector $`𝐱^0`$. Eqs. (10), (14) and (16) give us the required relations.

$`W_1^1`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{d}{}}}A_{1i}^0W_i^0,`$ (17)
$`W_1^2`$ $`=`$ $`{\displaystyle \underset{i,j=1}{\overset{d}{}}}A_{1i}^1A_{ij}^0W_j^0,`$ (18)
$`\mathrm{}`$ (19)
$`W_1^{d1}`$ $`=`$ $`{\displaystyle \underset{i,\mathrm{},l,m=1}{\overset{d}{}}}A_{1i}^{d2}\mathrm{}A_{lm}^0W_m^0,`$ (20)
These are $`d1`$ simultaneous equations for $`W^0`$.

The numerical procedure is as follows. We set the intial state of system (1) to a random initial guess vector $`\left(𝐲^0\right)_{old}`$ with $`\left(y_1^0\right)_{old}=x_1^0`$ and evolve it using Eq. (1). Using this vector $`𝐲(t)`$ we write down $`d1`$ simultaneous equations (Eqs. (20)) which can be solved for $`d1`$ unknown components of $`𝐰^0=\delta 𝐲^0`$. Also, $`\delta y_1^0=0`$. Thus the initial guess vector can be improved by

$$\left(𝐲^0\right)_{new}=\left(𝐲^0\right)_{old}+\delta 𝐲^0.$$
(21)
This sets up an iterative scheme giving us better and better estimates of the initial vector which converge to $`𝐱^0`$.

We note that as in Newton-Raphson method, the choice of the initial guess vector can be very important . In some cases, the iterative procedure of Eq. (21) may not converge or converge to a wrong root. In such cases, a different choice of initial guess vector can be useful.

We further note the similarity of our method with the so called method of variational equations in analytical dynamics . The method of variational equations can be appied to a known Hamiltonian system to determine an unknown neighbouring trajectory to an already known one. There, the method requires a complete particular solution of a known set of Hamiltonian equations of motion. In contrast, we have used our method for dissipative chaotic systems. In such systems an analytical solution of the equations of motion cannot be known. Further our method requires only one component of a complete trajectory to be sampled. This has important consequences in the problem of synchronization using a scalar signal.

We now demonstrate our method of estimating the initial state. As our first example we discuss the Rössler system given by ,

$`\dot{x}_1`$ $`=`$ $`x_2x_3,`$ (22)
$`\dot{x}_2`$ $`=`$ $`x_1+ax_2,`$ (23)
$`\dot{x}_3`$ $`=`$ $`b+x_3(x_1c).`$ (24)
First, we consider a case when the time series for $`x_1`$ is given and we want to estimate $`(x_2^0,x_3^0)`$. We chose the parameters $`(a,b,c)`$ such that the system is in the chaotic regime and the initial state $`𝐱^0`$ is on the chaotic attractor. We start with an arbitrary initial state $`𝐲^0=(y_1^0,y_2^0,y_3^0)`$ with $`y_1^0=x_1^0`$. From Eqs. (20) we get a pair of simultaneous equations as,

$`w_1^1`$ $`=y_1^1x_1^1=`$ $`\mathrm{\Delta }t\delta y_2^0+\mathrm{\Delta }t\delta y_3^0`$ (25)
$`w_1^2`$ $`=y_1^2x_1^2=`$ $`\left(2\mathrm{\Delta }t+a(\mathrm{\Delta }t)^2\right)\delta y_2^0`$ (27)
$`+\left(2\mathrm{\Delta }t+(y_1^0c)(\mathrm{\Delta }t)^2\right)\delta y_3^0`$
which can be solved for $`(\delta y_2^0,\delta y_3^0)`$. With $`\delta y_1^0=0`$ we use these in an iterative manner (Eq. (21)) to obtain the correct intial conditions.

Table 1 shows successively corrected $`(y_2^0,y_3^0)`$ obtained using the iterative process as discussed. These are the successive estimates for $`(x_2^0,x_3^0)`$. Let $`e_i=|y_i^0x_i^0|`$ denote the absolute error in the estimation of $`x_i^0`$. In Fig. 1(a) we plot a graph corresponding to Table 1 showing errors $`e_2`$ and $`e_3`$ (on logarithmic scale) plotted against the number of iterations of our method (Eq. (21)). From Table 1 and Fig. 1(a) we see that the successive estimates converge to the correct values of $`(x_2^0,x_3^0)`$. Using only two data points in the given time series $`x_1(t)`$, we can thus readily estimate the full initial state $`𝐱^0`$. We also note that the rate of convergence is very good. In about 8 to 10 iterates we obtain the initial values $`(x_2^0,x_3^0)`$ to within computer accuracy. If we write the deviations of the succesive iterates from the correct values in the form

$$(e_i)_n=\left|\left(y_i^0\right)_nx_i^0\right|e^{\alpha n},$$
(28)
where $`n`$ is the number of iterations, then the value of the parameter $`\alpha `$ is found to be $`2.01`$ for $`e_2`$ and $`2.02`$ for $`e_3`$. This is consistant with the fact that Newton-Raphson method has a quadratic convergence .

We note that the largest Lyapunov exponent for the subsystem $`(y_2^0,y_3^0)`$ (conditional or subsystem Lyapunov exponent) is positive . The success of our method does not depend on whether this Lyapunov exponent is positive or negative. This is important for synchronization of chaotic signals as will be discussed afterwards.

We next present cases where time series for the variables $`x_2`$ and $`x_3`$ of the Rössler system are given. The procedure is similar to the case of time series for $`x_1`$ as discussed above. Fig. 1 (b) shows the errors $`e_1`$ and $`e_3`$, when time series for $`x_2`$ is given, plotted against the number of iterations. The parameter $`\alpha `$ (Eq. (28)) is $`1.97`$ for $`e_1`$ and $`1.95`$ for $`e_3`$. This again indicates a quadratic convergence. Similarly, Fig. 1 (c) shows the quantities $`e_1`$ and $`e_2`$ when time series for $`x_3`$ is given, as a function of the number of iterations. The parameter $`\alpha `$ (Eq. (28)) is $`1.28`$ for $`e_1`$ and 1.30 for $`e_2`$ which shows a convergence slower than quadratic. We note that the largest subsystem Lyapunov exponent is negative when time series for $`x_2`$ is given and is positive when time series for $`x_3`$ is given .

As our next example we consider the Chua’s circuit which in its dimensionless form is given by ,

$`\dot{x}_1`$ $`=`$ $`\alpha (x_2x_1f(x_1)),`$ (29)
$`\dot{x}_2`$ $`=`$ $`x_1x_2+x_3,`$ (30)
$`\dot{x}_3`$ $`=`$ $`\beta x_2,`$ (31)
$`f(x_1)`$ $`=`$ $`bx_1+{\displaystyle \frac{1}{2}}(ab)\left[|x_11||x_1+1|\right].`$ (32)
We chose the parameters $`a,b,\alpha `$ and $`\beta `$ such that the attractor is a limit cycle. The initial state $`𝐱^0`$ is chosen in the basin of attraction of this limit cycle. Fig. 2 shows the errors $`e_2`$ and $`e_3`$ when a time series for the variable $`x_1`$ is given as a function of the number of iterations of our method. The parameter $`\alpha `$ (Eq. (28)) in this case takes values $`1.29`$ for $`e_2`$ and $`1.25`$ for $`e_3`$ showing a slower than quadratic convergence.

We have also applied our method for the cases when time series for the variables $`x_2`$ and $`x_3`$ are given and also for cases when parameters are such that the attractor is chaotic. In all the cases we are able to estimate the full initial state vector.

We have successfully applied our method to estimate the initial state vector using a given scalar time series for many other dynamical systems as well. These include Lorenz system  in its periodic, chaotic or intermittent regimes, the disk dynamo system modelling a periodic reversal of earth’s magnetic field , a 3-d plasma system formed by a three wave resonant coupling equations  and a four dimensional phase converter circuit .

Now we will discuss an important application of our estimation method in the problem of synchronization of two identical chaotic systems coupled unidirectionally by a scalar signal. Let us suppose that Eq. (1) describes a chaotic system and let us consider a replica of it given by, $`\dot{𝐲}=𝐅(𝐲)`$. Without losing generality let’s further assume that a scalar output signal $`x_1(t)`$ is given. The aim is to synchronize vector $`𝐲(t)`$ with $`𝐱(t)`$ using this scalar signal. Our method of estimating initial conditions makes this procedure trivial. Using the scalar signal we estimate the initial vector $`𝐱^0`$ and set $`𝐲^0=𝐱^0`$. This clearly leads to an instantaneous synchronization of the two trajectories. As we have demonstrated in the case of Rössler system, our method works even when the largest conditional Lyapunov exponent is positive which is the case when other methods of synchronization are known to fail .

To summarise, we have introduced a novel yet simple method to estimate initial conditions of a multivariable dynamical system from a given scalar signal. Our method is based on a multidimensional Newton-Raphson method where we include the time evolution of the system. The method gives a reasonably fast convergence to the correct initial state. The required length of the time series is very small. The method works even when the largest conditional Lyapunov exponent is positive. An important consequence of the method is that the problem of synchronization of identical chaotic systems using scalar signal becomes trivial since evolution of two such systems can then be started from identical initial states.