text:
# The Frequency Stability of Millisecond Oscillations in Thermonuclear X-Ray Bursts

## 1 Introduction

Millisecond oscillations have been observed during thermonuclear X-ray bursts from nine neutron stars in low mass X-ray binaries (LMXBs; see Strohmayer, 2001, for a review). The bursts are triggered by the unstable nuclear burning of accreted material on the neutron star’s surface (see Lewin, van Paradijs & Taam, 1993, for a review). It therefore has long been expected that anisotropies in the burning should produce pulsations at the stellar spin frequency (e.g Bildsten, 1995; Strohmayer et al., 1996). Indeed, the general characteristics of the observed oscillations suggest that they originate from a brightness anisotropy on the stellar surface. First, the oscillations are highly coherent towards the end of the burst (Strohmayer et al., 1996; Strohmayer & Markwardt, 1999). Second, the oscillations are frequently observed in the rise of bursts, when there is spectral evidence for a growing burning region (Strohmayer, Zhang, & Swank, 1997). Third, the frequencies of these oscillations are remarkably similar in bursts separated by several years (Strohmayer et al., 1998; Muno et al., 2000), which suggests that they originate from a stable clock.

If the burst oscillation frequency, $`\nu _{\mathrm{burst}}`$, is indeed the spin frequency of the neutron star, the oscillations are key to understanding several aspects of these systems. Pairs of kilohertz quasi-periodic oscillations (kHz QPOs) are observed in the persistent emission from many neutron star LMXBs. Their frequency difference $`\mathrm{\Delta }\nu _{\mathrm{kHz}}`$ remains comparable to $`0.5`$ or $`1\times \nu _{\mathrm{burst}}`$, even though the frequencies of the kHz QPOs vary by more than a factor of 2 (see van der Klis, 2000, for a review). It has been suggested that these QPOs result from signals at a Keplerian frequency in the disk and a beat between that frequency and the spin of the neutron star (Strohmayer et al., 1996; Miller, Lamb, & Psaltis, 1998).111Alternative models for the kHz QPOs cannot account naturally for the similarity between $`\nu _{\mathrm{burst}}`$ and $`\nu _{\mathrm{kHz}}`$ (e.g Psaltis, 2001). If this is the case, we can infer the rotation rate of the neutron stars in more than 20 systems (van der Klis, 2000). The inferred spin frequencies cluster around $`300`$ Hz. This confirms the suggestion that these neutron star LMXBs are the progenitors of recycled millisecond radio pulsars (Alpar et al., 1982; Radhakrishnan & Srinivasan, 1982). However, it also implies that some mechanism limits the frequencies to which these stars can be spun-up by accretion, such as a correlation between the surface magnetic field strength and X-ray luminosity (White & Zhang, 1997) or gravitational radiation from a quadrupole in the mass distribution of the star (Bildsten, 1998).

A closer examination of these burst oscillations reveals an even more complex picture (see also Psaltis, 2001, for a detailed discussion). First, the coherence of the oscillations is cast into some doubt by the detection of a sudden 0.25 cycle phase shift in one oscillation out of seven studied by Strohmayer (2001). Second, oscillations are not only observed in the rise of thermonuclear bursts, but also far into the decay (e.g. Smith, Morgan, & Bradt, 1997), even though by this time the fuel over the entire surface of the neutron star is expected to have been consumed by the burning front leaving no observable brightness anisotropies (e.g. Fryxell & Woosley, 1982). Third, while the spin of the neutron star must be constant during the course of a burst, the frequency of the oscillations is usually observed to increase and saturate at some maximum frequency (Strohmayer et al., 1997). It is this “asymptotic” frequency that appears stable in bursts separated by several years. The frequency drift is thought to originate from one of several mechanisms that might carry anisotropies in a retrograde motion about the neutron star, resulting in an observed frequency slightly lower than that of the star’s spin (Strohmayer et al. 1997; Cumming & Bildsten 2000; Spitkovsky, Levin, and Ushomirsky 2002; Heyl 2002). Finally, in a handful of bursts there are episodes of drifts to lower frequency (“spin-down”) in the tails of the bursts (Strohmayer, 1999; Miller, 2000; Muno et al., 2000) and strong detections of signals at second frequencies within a few Hz of the primary signal (Miller, 2000; Galloway et al., 2001).

If anything is clear, it is that more observational clues are required to understand burst oscillations. In this paper, we address the how the frequency of these oscillations evolves as a function of time, in order to establish the coherence of the underlying mechanism. In particular, we employ standard pulse phase analysis techniques adapted from pulsar timing studies (Manchester & Taylor, 1977).

## 2 Observations and Data Analysis

Our analysis used observations with the Proportional Counter Array (PCA; Jahoda et al. 1996) on the Rossi X-Ray Timing Explorer (RXTE). PCA consists of five identical gas-filled proportional counter units with a total effective area of 6000 cm2 and sensitivity in the 2.5–60 keV range. The detector is capable of recording photons with microsecond time resolution and 256-channel energy resolution. The data were recorded in a wide variety of data modes with different time and energy resolutions, depending upon the details of the original proposed programs and the available telemetry bandwidth. For all of the analysis presented here, we converted the TT photon arrival times at the spacecraft to Barycentric Dynamical Time (TDB)

at the solar system barycenter, using the Jet Propulsion Laboratory DE-200 solar system ephemeris (Standish et al. 1992).

We searched the entire RXTE public data archive for X-ray bursts from 8 neutron stars222Oscillations in bursts associated with MXB 1743$``$29 have been observed during observations of the bursting pulsar GRO J1744$``$28. A search for these bursts was not part of the analysis presented here. that are known to exhibit burst oscillations (see Table 1 and Muno et al. 2001). As of September 2001, we have identified a total of 159 X-ray bursts from these 8 sources. Each of these bursts was then searched for millisecond oscillations as described below.

For our analysis, we used data containing a single energy channel (2.5–60 keV) and $`2^{13}`$ s (122 $`\mu `$s) resolution. We note that the use of more restricted energy bands can sometimes allow oscillation trains to be detected for a second or two longer, since the oscillation amplitudes are often stronger at higher energies (e.g., Muno et al. 2000; Giles et al. 2002). However, in many cases, the data modes with higher energy resolution contained data gaps due to the limited size of the satellite memory buffers, making it difficult to produce a coherent phase model spanning the entire oscillation train (see Section 2.2). Data modes specifically designed to start recording during a burst (“burst catcher” modes) were analyzed in combination with the other modes in order to fill data gaps, but the burst catcher modes with high time resolution only recorded in a single 2.5–60 keV energy channel.

### 2.1 Fourier Timing

We produced Fourier power spectra of 1 s intervals of data for the first 16 s of each burst, and searched for signals within a frequency range of $`\pm 5`$ Hz from the frequencies listed in Table 1. Any oscillation with a single trial probability of less than $`2\times 10^5`$ that it is due to noise (a power of $`P>21.6`$ according to the normalization of Leahy et al. 1983) was considered a detection. This corresponds to a 1% probability that a noise signal will be that strong given a search of an entire burst, so that 1–2 noise signals may appear significant from the search of our entire sample of bursts. We detected 68 oscillation trains in this manner (Table 1). We also recorded the frequencies of each signal detected. The frequencies in Table 1 were determined recursively, by computing the median value (rounded to the nearest Hz) from the signals detected until that value did not change. The range of detected frequencies is reported in Section 3.1.

To qualitatively examine the frequency evolution, we computed power spectra of 2 s intervals of data every 0.25 s, and plotted contours of power as a function of both time and frequency. We tapered the data with a Welch function (Press et al., 1992) before computing each power spectrum, in order to suppress the power in sidebands due to the window function. A typical example of such a dynamic power spectrum is illustrated in Figure 2. We use dynamic power spectra to characterize the oscillations in Sections 3.2–3.3.

### 2.2 Phase Connection

In order to examine in detail whether the oscillations are coherent and drift to stable asymptotic frequencies, we developed frequency models using the phase connection method described in Muno et al. (2000). We modeled all oscillations that lasted continuously for more than 2 s using the $`2^{13}`$ s data (Table 1). A summary of the procedure is as follows. We fold the data in 0.25–0.5 s intervals according to a trial frequency model (which is the derivative of a phase model), and compute a Fourier transform of the resulting profile. If the power $`P`$ in the fundamental frequency of the pulsation333We have confirmed that all of the profiles were consistent with sinusoidal signals (Muno et al., 2000, 2001), justifying the use of only the power at the fundamental frequency. has less than a 10% chance of occurring randomly due to noise ($`P>4.6`$ in the normalization of Leahy et al. 1983), we measure its phase by a linear least squares fit of a sinusoid to the profile. The measured phase residuals $`\varphi _{\mathrm{obs}}(t_i)`$ are then related to the predicted phases $`\varphi _{\mathrm{model}}(t_i)`$ by

$$\mathrm{\Delta }\varphi (t_i)\varphi _{\mathrm{obs}}(t_i)\varphi _{\mathrm{model}}(t_i).$$
(1)
In the above convention, an increasing phase residual indicates that the instantaneous trial frequency ($`d\varphi _{\mathrm{model}}/dt_i`$) is too low.

The phase residuals $`\mathrm{\Delta }\varphi (t_i)`$ are then modeled via a least squares minimization, using a $`\chi ^2`$ statistic calculated from the phase residuals:

$$\chi ^2=\frac{(\mathrm{\Delta }\varphi (t_i))^2}{(\sigma _{\mathrm{\Delta }\varphi })^2}.$$
(2)
Here, $`\sigma _{\mathrm{\Delta }\varphi }`$ is the uncertainty on the phase measured from the folded profile, which can easily be determined from the diagonal values of the covariance matrix in the linear least squares fit (Press et al., 1992). We have confirmed that these uncertainties are accurate via Monte Carlo simulations, in which we (i) perturbed the counts in each bin of the profile within the range expected from Poisson counting noise and (ii) recorded the uncertainty as the standard deviation of phase measurements of 100 perturbed profiles. We also verified the values for $`\sigma _{\mathrm{\Delta }\varphi }`$ by visual inspection.

The phase connection technique provides a correction $`d(\mathrm{\Delta }\varphi )/dt`$ to the original frequency model $`\nu _{\mathrm{model}}`$, so that the best-fit frequency is

$$\nu (t)=\nu _{\mathrm{model}}(t)+\frac{d}{dt}(\mathrm{\Delta }\varphi )$$
(3)
This procedure can be iterated using the new frequency model, if it is composed of linear functions.

The advantages of this method are that (i) the form of the observed phase residuals naturally suggests the best model for the frequency evolution and provides a quantitative test of how well the trial phase evolution matches the data, (ii) there exist standard techniques to calculate the uncertainties on the model from the least squares fit to the phases, and (iii) the method utilizes the phase information that the power spectrum discards, which allows much finer frequency resolution on short time scales. In contrast, the $`Z^2`$ method of Strohmayer & Markwardt (1999) distinguishes between models based upon the amount of power in the folded profile, but provides little information about how the frequency of an oscillation deviates from the model assumed. The $`Z^2`$ statistic is still very useful for examining the amplitude and phase evolution of the oscillations on short (less than 0.1 s) time scales after the frequency evolution is modeled. We apply this technique and use the resulting models in Sections 3.4–3.7.

### 2.3 Energy Spectra

Finally, we produced energy spectra in 0.25 s intervals from each burst using available combinations of data modes that provide at least 32 energy channels. The detector response was estimated using PCARSP in FTOOLS version 5.1 (see http://heasarc.gsfc.nasa.gov/lheasoft/). We subtracted spectra from 15 s of emission prior to the burst to account for background, and fit each spectrum between 2.5–20 keV with a model consisting of a blackbody multiplied by a constant interstellar absorption factor. The absorbing column for each burst was taken to be the mean value from fits to each 0.25 s interval assuming a variable absorption. The model provides a color temperature ($`T_{\mathrm{col}}`$) and a normalization proportional to the square of the apparent radius of the burst emission surface ($`R_{\mathrm{col}}`$), and allows us to estimate the bolometric flux as a function of time. We compare these spectra to the properties of the oscillations in Sections 3.3 and 3.6.

## 3 Results

### 3.1 Frequency Range of Oscillations

In order to establish the range of frequencies over which the oscillations are observed, and in particular whether there are strict upper limits to their frequencies, we produced a histogram of the distribution of the frequencies at which signals were detected in non-overlapping 1 s intervals during the bursts from each source (Figure 2). For the majority of sources, the histograms are narrowly distributed with a width of 2-4 Hz, significantly smaller than the 10 Hz search window. This corresponds to frequency drifts ranging from 0.4% (KS 1731$``$260) to 1.2% (4U 1702$``$429).

The distribution from 4U 1636$``$536 appears to be double-peaked, possibly owing to two effects. First, oscillations are usually observed from this source only in the rise and tail, and the frequency gap between the two peaks could result from the fact that we do not observe the intermediate frequencies as the oscillations evolve (see Section 3.2). However, two oscillations are present simultaneously in several bursts from this source (see

Section 3.3 and Miller 2000), which could indicate there are indeed two fundamental frequencies.

One detected oscillation is unusual given the distributions in Figure 2, and has already been noted by Wijnands, Strohmayer, & Franco (2000). It is detected in MXB 1659$``$298 about 4 Hz above the next highest frequency, with a probability of $`4\times 10^4`$ that it is due to noise given all of the trial frequencies searched. The oscillation occurs at 571 Hz in a single interval about 15 s into the tail of a burst, while the oscillation detected earlier in the burst seemed to disappear about 5 s into the burst after reaching an asymptotic frequency of 567 Hz. If this

signal is a burst oscillation, it calls into question whether we are observing frequency saturation at all.

### 3.2 Instances of Unusual Frequency Evolution

We also searched the dynamic power spectra to confirm that the oscillations in our sample drift upward in frequency during the course of the burst, saturating at an approximately constant frequency late in the tail. This trend occurs in most bursts, and is illustrated in Figure 2. However, there are a few exceptions that have been noted in the literature, such as (i) clear evidence for a frequency decrease several seconds after the burst is first detected (spin-down) and (ii) signals present at two frequencies simultaneously. Spin-down episodes and simultaneous signals are not common.

We find only 3 clear examples of spin-down episodes during bursts, which we display in Figure 3.2. Two of these have been reported previously, from 4U 1636$``$536 (Strohmayer, 1999) and KS 1731$``$260 (Muno et al., 2000), and we find one additional example from 4U 1728$``$34. In all cases, the frequency decrease is on order 1 Hz. There is no clear pattern as to when in the burst the spin-down occurs. In KS 1731$``$260, the event occurs in the peak of the burst, while in 4U 1728$``$34 and 4U 1636$``$536 the spin-down occurs as the burst decays.

We also find two clear examples of bursts in which two signals are detected simultaneously at frequencies separated by about 1 Hz. The most significant example has already been reported by Miller (2000) from 4U 1636$``$536, and we find a second example in 4U 1702$``$429. Both are displayed in Figure 3.3. The second oscillation in 4U 1636$``$536 occurs between 2–3 s into the burst, and has a chance probability $`<3\times 10^{10}`$ of occurring randomly due to noise (see also Miller, 2000). The second signal from 4U 1702$``$429 occurs 11 s into the burst, and is only marginally significant, with a single-trial chance probability of $`5\times 10^5`$. Less significant examples, with chance probabilities of $`0.01`$, have also been reported in

two more bursts from 4U 1636$``$536 by Miller (2000), and in one burst from 4U 1916$``$053 by Galloway et al. (2001). In all cases, the secondary frequencies are 1–3 Hz lower than the main signal.

### 3.3 How Radius Expansion Affects Oscillations

In many bursts, photospheric radius expansion is evident at the start of the burst, during which $`R_{\mathrm{col}}`$ increases and $`T_{\mathrm{col}}`$ decreases such that the bolometric flux remains constant, presumably at the Eddington limit (see Lewin et al. 1993). We plot examples of bursts with and without radius expansion in Figure 5. We examined whether radius expansion affects the times in the burst when the oscillations are observed. We find that all 16 of the bursts in which oscillation trains are observed continuously throughout the rises, peaks, and tails fail to exhibit radius expansion (e.g. Figure 5a). Of the 51 bursts where the oscillation train is not observed during the peak, 40 exhibit radius expansion (e.g. Figure 5b). One burst from 4U 1608$``$522 (1998 March 27 at 14:08:30 TDB) could not

be examined, because there were gaps in the high time resolution data during the peak. Since it also has been noticed previously that oscillations in KS 1731$``$260 (Smith, Morgan, & Bradt, 1997) and MXB 1743$``$29 (Strohmayer et al., 1997) appear only after radius expansion ends, we conclude that radius expansion either obscures or interrupts the oscillations.

### 3.4 Mathematical Form of the Frequency (Phase) Evolution

We modeled the frequency evolution of the oscillations using the phase connection technique described in Section 2.2. We tried both polynomial models (up to fifth order) and a saturating exponential such that the frequency evolution is of the form used by Strohmayer & Markwardt (1999):

$$\nu (t)=\nu _0+\mathrm{\Delta }\nu \mathrm{exp}(t/\tau ).$$
(4)
Here, $`\nu _0`$ is the asymptotic frequency of the oscillation, $`\mathrm{\Delta }\nu `$ is the frequency drift referenced to the start of the burst, and $`\tau `$ is the time scale for the frequency drift. When modeling the data with polynomials, we iterated the fitting procedure until no improvement in the fits to the phase residuals was found. The fits using the exponential model were executed on phase residuals found assuming a constant trial frequency, as the corrections must add linearly to the trial phase model in order to use the iterative process. Table 2 lists the number of oscillations best fit with each model, as determined from the value of reduced $`\chi ^2`$. If the reduced $`\chi ^2`$ was greater than one, and if both the number of degrees of freedom and the value of $`\chi ^2`$ decreased in comparing one model to another, an $`F`$-test (Bevington, 1969) was used to ensure that the decrease in $`\chi ^2`$ had less than a 5% chance of being random.

Figure 6a illustrates the importance of considering the phases in modeling the frequency evolution. Judging from power spectra alone, the frequency of this oscillation from 4U 1636$``$536 appears roughly constant (top panel). However, when one examines the phase residuals from a constant frequency model it is clear that additional evolution has occurred (center panel). An acceptable fit to the frequency evolution only can be achieved using a fifth order polynomial (bottom panel). Although some of this frequency evolution can be inferred from the power spectra alone (see also Miller, 2000), by tracking the phase we are able to achieve high frequency resolution and to determine how coherent the oscillations are.

We modeled 59 oscillations from 6 different sources (Table 1). The results are summarized in Table 2. Thirty-seven oscillations required either an exponential model or a polynomial of third degree or higher. Exponential models are only favored over polynomials in only 6 out of 37 cases. The exponential models are formally inconsistent (greater than 99% confidence) in 22 of the remaining 31 cases, while all such frequency models are inconsistent with the data in 12 out of 37 cases.

The polynomial models generally reproduce the frequency evolution better than exponential models because there is slower frequency evolution at the start of the oscillation train, inflection in the evolution during the train, or a decrease in the oscillation frequency as the trains disappear (see also Miller, 2000). Figure 6 exhibits all of these variations in the frequency evolution. In addition to the upward drifts in frequency that saturate after a few seconds, the frequencies also tend to wander by about 0.1 Hz (see, e.g., the evolution in the top panel of Figure 6a).

We see no evidence for a rapid decrease in frequency at the start of the oscillation that is comparable in magnitude to the later upward drift in frequency. Such a frequency decrease could occur as the burning front spreads around the surface of the neutron star (e.g., Strohmayer et al., 1997; Spitkovsky et al., 2002), but may be too rapid to observe using our analyses.

### 3.5 Short-term Stability of the Burst Oscillations

In general, we find many more large values of $`\chi ^2`$ than would be expected if the phase evolution were described by smooth, low-order polynomial and exponential functions. In particular, 30% of the oscillation trains are inconsistent with smooth models with 99% confidence. This could indicate that the underlying mechanism is not perfectly stable. We tried higher-order polynomials in instances where our best-fit model is statistically unacceptable, but we find that we achieve only marginal improvements in $`\chi ^2`$. Higher-order polynomials can smooth residuals in the middle of the oscillation train, but tend to create larger residuals at the endpoints. Figure 7 illustrates two examples of oscillation trains that we were unable to fit with

exponential or polynomial models. Piecewise smooth phase functions would provide better fits in some instances, such as between 2–4 s in Figure 7a, and 3–6 s in Figure 7b. However, the resulting discontinuities in the derivatives of the phases would imply that the frequencies of the oscillations shift in less than 0.25 s.

Moreover, some of the residuals appear to be discrete changes in the phases of the oscillations. One possible explanation for the phase residuals in the first 1–2 s of Figure 7a is that two signals are present simultaneously within about 1 Hz, the resolution of the folding technique. Simultaneous signals with only slightly larger frequency separations have been seen in other bursts from 4U 1636$``$536 (Miller, 2000). In Figure 7b, there appear to be discrete jumps of about 0.1 cycle in the phase of the oscillation from 7–10 s into the burst, similar to that observed by Strohmayer (2001) during an oscillation in the rise of a burst from 4U 1636$``$536. We confirmed that the sudden changes in the phases of the oscillations are not due to variations in the pulse profile, which is always sinusoidal. The presence of phase jumps suggests that the oscillations are not strictly coherent.

### 3.6 Magnitude and Time Scale of Frequency Drift

Our models also allowed us to compare the frequency evolution to other properties of the bursts. We first measured the observed change in frequency $`\mathrm{\Delta }\nu `$ for the continuous portion of each oscillation. We then measured the time scale $`\tau `$ that represents how long the oscillation took to evolve in frequency by $`(1e^1)\mathrm{\Delta }\nu =0.632\times \mathrm{\Delta }\nu `$ from the minimum observed frequency (compare Equation 2). As measures of the burst time scales, we have fit the decay in flux of each of these bursts with either one or two exponential functions ($`\mathrm{exp}[t/t_d]`$). We also computed the peak bolometric flux and the fluence of each burst, using the exponential fit to determine the flux contribution at late times. We find that neither $`\mathrm{\Delta }\nu `$ nor $`\tau `$ are correlated with the above burst properties.

We have also compared $`\mathrm{\Delta }\nu `$ and $`\tau `$ with the duration of the oscillations and the time at which they began. We find that $`\tau `$ is not correlated with the start time or duration of the oscillations. Instead, $`\tau `$ is distributed uniformly between 0.75–5.5 s. We plot the frequency drift $`\mathrm{\Delta }\nu `$ as a function of these parameters in Figure 3.6. Oscillations that are first observed early in the burst evolve over a larger frequency range (top panel), but their frequency drift is not correlated with how long the oscillations last (bottom panel). This suggests that the evolution is initiated by the start of the burst.

Cummings & Bildsten (2000) previously noticed that the fractional frequency drifts $`\mathrm{\Delta }\nu /\nu `$ in sources with $``$600 Hz oscillations are a factor of two smaller than those in $``$300 Hz oscillations. We have indicated the data from the faster ($``$600 Hz) oscillations in Figure 3.6 with filled circles. The continuous portions of the fast oscillations tend to start later in the burst (with only two exceptions) and to drift by a smaller absolute frequency range. We believe that this is due to the fact that the faster oscillations tend to occur in bursts with radius expansion (Muno et

al. 2001),which prevents the oscillations from being observed earlier in the bursts (see Section 3.1). Indeed, the fast oscillators in Figure 2 do not in general have a narrower distribution of observed frequencies when considering the non-continuous portions of the oscillations as well. However, we do confirm that the fractional frequency drift is a factor of two smaller in the $``$600 Hz oscillators than in the $``$300 Hz ones, but only because the frequencies $`\nu `$ are larger.

If the frequencies of the oscillations are associated with underlying stable clocks, then the fact that they evolve in frequency implies that the phases of the oscillations drift by several cycles with respect to those clocks. If we take the maximum observed frequency ($`\nu _{\mathrm{max}}`$) to represent the stable clock (see Section 3.7), then we can define a phase loss by

$$\mathrm{\Delta }\varphi _{\mathrm{loss}}=\nu _{\mathrm{max}}t\nu (t)𝑑t,$$
(5)
where $`\nu (t)`$ is the derivative of our phase model. We find that oscillations typically lose between 1–5 cycles, depending on the magnitude of the frequency drift, the duration of the oscillations, and the time scale for the frequency evolution. One exceptional oscillation, illustrated in Figure 6b, drifts by more than 17 cycles with respect to $`\nu _{\mathrm{max}}`$. This oscillation is unusual for its combination of a long duration (10 s), a large frequency change ($`\mathrm{\Delta }\nu =4.1`$ Hz), and a slow drift time scale ($`\tau =5.5`$ s).

### 3.7 “Asymptotic” Frequencies of Burst Oscillations

Aside from the exceptions noted above, the fundamental features that led to the adoption of an exponential model for the frequency evolution by Strohmayer & Markwardt (1999) are present in about $`70\%`$ of burst oscillations: the oscillations increase in frequency when they are first observed, but appear to reach stable frequencies on the time scale of seconds. Several authors have noted that the asymptotic frequencies of burst oscillations detected several years apart are remarkably similar (Strohmayer et al., 1998; Strohmayer & Markwardt, 1999; Muno et al., 2000). Therefore, we measured the maximum frequencies that the oscillations reached according to our best fit models, in order to determine whether the apparent saturation frequencies of the oscillations are stable. The 1-$`\sigma `$ uncertainties on the maximum frequency were found by a search through $`\chi ^2`$ space (e.g Lampton, Margon, & Bowyer, 1976). In those instances where the best available fit was statistically unacceptable, we inflated the uncertainties on the phase measurements to for the reduced $`\chi ^2`$ to equal 1 before the search in $`\chi ^2`$ space. The maximum frequency was constrained to occur after a minimum, as in several cases the absolute highest frequency was observed at the start of the oscillation train (e.g Miller, 2000, and Figure 3a). Our measurements are consistent with those found using the $`Z^2`$ statistic by Strohmayer et al. (1998), Strohmayer & Markwardt (1999), and Giles et al. (2002). We compare our results for 4U 1636$``$536 with those of Giles et al. (2002) in Section 4.2

In Figure 9, we plot how the maximum frequencies vary with time. To indicate those oscillations where the frequency has apparently saturated, we have indicated which ones required at least an exponential or a 3rd order polynomial phase model by diamonds. The measurements taken from first or second order polynomial models are plotted only with error bars. We have calculated the fractional standard deviation in the maximum frequencies, $`\sigma _\nu /\nu _{\mathrm{max}}`$, for those sources with more than two values, and listed them in Table 3. For Aql X-1 and KS 1731$``$260, the frequency evolution appeared to saturate during two bursts each, so we report the fractional difference between the measurements. The deviations in these frequencies $`\sigma _\nu /\nu _{\mathrm{max}}`$ are quite small, less than one part in 1000.

## 4 Discussion

The frequency evolution of burst oscillations can be described by smooth phase models in all but 30% of the pulsations that we have studied (Section 3.5; see also Strohmayer 2001). The discrepant cases appear as phase jumps in Figure 7. If the oscillations result from brightness asymmetries on the surface of the neutron star, it appears that some of the following occur: (i) there are multiple anisotropies or modes that propagate with slightly different velocities relative to the surface, which results in two simultaneous signals that interfere with each other (see also Section 3.2), (ii) the velocity of the pattern changes suddenly on time scales of 0.25 s, or (iii) the anisotropies change their longitude in less than a tenth of a second, resulting in sudden phase jumps of about 0.1 cycles.

On time scales of years, the oscillations are stable to a few parts in 1000. This is consistent with the Doppler

shifts that are expected from the binary orbit. Since this raises the possibility of measuring a mass function for the companion star, we now examine this frequency dispersion in detail.

### 4.1 The Stability of the Maximum Frequencies

The change in the apparent spin frequency due to the orbital motion of the neutron star is given by

$$\frac{\mathrm{\Delta }\nu _{\mathrm{max}}}{\nu _{\mathrm{max}}}=2.0\times 10^3\frac{M_c\mathrm{sin}i}{P_{hr}^{1/3}(M_x+M_c)^{2/3}},$$
(6)
where $`M_x`$ and $`M_c`$ are the mass of the neutron star and companion respectively in solar units, $`P_{hr}`$ is the binary orbital period in hours, and $`i`$ is the inclination of the orbital plane normal to the line of sight. Assuming a star with solar density that over-fills it Roche lobe, the mass and orbital period are related by $`M_c0.11P_{\mathrm{hr}}`$ (Frank, King, & Raine 1995). For $`P_{hr}=5`$, $`\mathrm{sin}i=0.87`$, $`M_x=1.4`$, and $`M_c=0.5`$, the half-amplitude of the Doppler variation would be $`\mathrm{\Delta }\nu _{\mathrm{max}}/\nu _{\mathrm{max}}=4\times 10^4`$. This is a factor of 2–5 smaller than the dispersions from 4U 1636$``$536, 4U 1702$``$429, and 4U 1728$``$34 (Figure 9), but it is of the right order of magnitude.

Only three oscillation sources have known orbital ephemerides: MXB 1659$``$298 (Wachter, Smale, & Bailyn, 2000), Aql X-1 (Welsh, Robinson, & Young, 2000), and 4U 1636$``$536 (Giles et al., 2002). We did not measure enough oscillations from MXB 1659$``$298 and Aql X-1 to constrain the parameters of the binary orbit. We have measured eleven oscillations that appear to saturate from 4U 1636$``$536 (the diamonds in Figure 9), so we have attempted to fit a sinusoid to the maximum frequencies as a function of orbital phase (Figure 4.1). The best-fit sinusoid, allowing the reference phase to vary, has a $`\chi ^2`$ of 221 for 8 degrees of freedom. Therefore, sinusoidal orbital modulations cannot cause all of the observed dispersion. In Figure 4.1 we plot a sinusoid with an amplitude equal to the standard deviation of the points fit, $`\sigma _\nu /\nu _{\mathrm{max}}=7.7\times 10^4`$. The disagreement between the data and the sinusoid is primarily due to three oscillations with low maximum frequencies. Two of these lasted for more than 5 s, one of which is displayed in the right panel of Figure 6 (1999 June 10). We will discuss this further in the next section.

There are no other effects that ordinarily change the apparent spin frequencies of neutron stars (i.e. pulsars) that can produce the dispersion in Figures 9 and 4.1. Accretion torques can alternately quicken or slow the rotation of the neutron star, and would provide the next largest effect (Frank et al., 1995). The torque on a neutron star can be estimated by assuming that its magnetic field disrupts the disk at some fiducial radius $`R_M`$ and removes all of its angular momentum. At $`R_M`$, we have $`I\dot{\mathrm{\Omega }}=\dot{M}(GMR_M)^{1/2}`$. We take $`R_M`$ to be the co-rotation radius, where

$`2\pi \nu _{\mathrm{kHz}}=(GM/R_M^3)^{1/2}`$, so that

$$\dot{\nu }=7\times 10^5\dot{M}R_6M_x^{2/3}I_{45}^1\nu _{\mathrm{kHz}}^{1/3}\mathrm{Hz}\mathrm{yr}^1$$
(7)
where $`\dot{M}`$ is the accretion rate in Eddington units, $`R_6`$ is the radius of the neutron star in units of 10 km, $`M_x`$ is its mass in solar units, $`I_{45}`$ is its moment of inertia in units of $`10^{45}`$ g cm2, and $`\nu _{\mathrm{kHz}}`$ is the stellar spin frequency in kHz. The values are scaled to those expected for a burster. For a spin frequency of 200 Hz, Equation 7 corresponds to a $`\mathrm{\Delta }\nu _{\mathrm{max}}/\nu _{\mathrm{max}}6\times 10^7`$ in one year. Clearly, the observed dispersion in the maximum frequencies on time scales of a month (Figure 9) is much larger than can be explained by accretion torques (Equation 7).

### 4.2 Frequency Dispersion in 4U 1636$``$536

Giles et al. (2002) recently examined a larger sample of oscillations from 4U 1636$``$536, and found that there is a subset of the burst oscillations whose maximum frequencies are distributed as a Gaussian about 581.64 Hz with a width of 0.08 Hz. We cannot independently test this, as the uncertainties that we report on the maximum frequencies are up to a factor of 5 larger than those in Giles et al. (2002). This is because the phase-connection technique we used to measure the maximum frequencies accounts for the frequency evolution of the oscillations, while the dynamic $`Z^2`$ technique Giles et al. (2002) used assumes that the frequency is constant during each 2 s interval that they test. In our analysis, the maximum frequency is usually observed at the end of the oscillation train where the data can tolerate relatively large variations in the values of the polynomial models. 444On the other hand, we may be over-estimating our uncertainties when we model the oscillations with higher-order polynomials, which tend to have large derivatives at their endpoints. A more physical model might assume that the frequency change in the oscillations is much smaller at the endpoints. However, the oscillation that appears 4 Hz from the main signal at the end of the burst on 1999 Apr 14 from MXB 1659$``$298 (Wijnands et al., 2001) makes us reluctant to apply this assumption in our analysis, and leads us to assume larger uncertainties. In the analysis of Giles et al. (2002), if significant frequency evolution occurs during a 2 s interval, they measure the mean frequency in that interval, but do not consider how much larger the frequency could be at the end of the interval. As a result the magnitude of our uncertainties, 0.05–0.4 Hz, are often larger than the width of the Gaussian of Giles et al. (2002).

If we only include those oscillations that are part of the Gaussian in Giles et al. (2002), we find that the data are consistent with no modulation. Fitting a sinusoid to the data, we derive a fractional amplitude of $`(5.0\pm 2.1)\times 10^4`$ ($`\chi ^2=6.7`$ for 5 degrees of freedom). This is equivalent to a Doppler shift of 150$`\pm `$63 km s-1, which is significantly larger than the 90% upper limit to the velocity from Giles et al. (2002; 55 km s-1). The larger uncertainties that we assumed result in larger upper limits on the velocity. These systematic uncertainties also present an additional hurdle in measuring the orbital motion of the neutron star.

### 4.3 Models for the Frequency Evolution

Three models have been proposed to explain the frequency evolution of the burst oscillations, all of which assume that anisotropies develop in the surface brightness of the neutron star (Strohmayer et al., 1997; Cumming & Bildsten, 2000; Heyl, 2002; Spitkovsky et al., 2002). The simplest model was proposed by Strohmayer et al. (1997), who suggested that the oscillations originate from hot regions on a burning layer that expands and decouples from the neutron star at the start of the burst. The oscillations are observed when the burning layer begins to contract, causing them to increase in frequency as the layer re-couples to the neutron star. However, recent calculations assuming that the burning layer is rigidly rotating predict a frequency drift a factor of three smaller than that observed (Cumming et al., 2002). Moreover, Cumming & Bildsten (2000) predicted that there should be a fractional dispersion of at most $`\mathrm{\Delta }\nu /\nu 10^4`$ in the observed asymptotic frequencies, if the height of burning layer is smaller when it re-couples to the star than when it decouples. This is a factor of 10 smaller than the residual dispersion from 4U 1636$``$536 in Figure 4.1.

As a result, it is important to consider other models for the frequency evolution of burst oscillations. Spitkovsky et al. (2002) have suggested that a brightness contrast could result from a hydrodynamical instability in a geostrophic flow, similar to Jupiter’s Great Red Spot, while Heyl (2002) has proposed that Rossby waves propagating retrograde on the neutron star could produce traveling brightness patterns on the stellar surface. The frequency evolution in both models is caused by changes in the velocity at which the pattern propagates around the star. The velocity in turn is sensitive to the temperature, density, and vertical structure of the surface layers. Thus, variations in the properties of the surface layers near the end of a burst can naturally explain the observed dispersion in the asymptotic frequencies of the oscillations (Figures 9 and 4.1). On short time scales, instabilities in either mechanism conceivably could cause the apparent phase jumps in Figure 7. Moreover, Heyl (2002) also has predicted that Rossby waves with different numbers of radial nodes could produce simultaneous signals at multiple frequencies, as are seen in rare cases (Section 3.2; compare Section 3.5). These models are promising, although how these mechanisms produce brightness variations on the stellar surface remains to be studied in detail.

## 5 Conclusions

We have examined 68 burst oscillations observed in 159 bursts from 8 different sources, and have modeled the frequency evolution of 59 oscillations from 6 sources (Table 2). The most notable result is that there appears to be slight instability in the mechanism generating the oscillations on time scales of both seconds and years. About 30% of the oscillation trains do not exhibit smooth frequency evolution, as both low order polynomials and exponential models are statistically inconsistent with the data (Figure 7). The possible explanations are that (i) two signals may be present simultaneously, invalidating our assumption that there is a single signal in our analysis, (ii) there are discrete phase jumps that occur on time scales less than a second, and (iii) the frequencies of the oscillations shift dramatically on time scales of 0.25 s. On longer time scales, the maximum frequencies observed during these oscillation trains exhibit a fractional dispersion of $`\mathrm{\Delta }\nu _{\mathrm{max}}/\nu _{\mathrm{max}}4\times 10^3`$ (Figure 9). For 5 of the sources studied, this dispersion could be consistent with Doppler shifts of signals originating from rotating neutron stars as they orbit the centers of mass of their binaries. However, the dispersion in the maximum frequencies in 4U 1636$``$536 is uncorrelated with its known orbital period (Figure 4.1).

We now are able to accurately characterize the general properties of the frequency evolution in our large sample of bursts:

1) The frequency of the oscillations drifts by up to 1.2% (Figure 2 and 2), but generally reaches a stable value on the time scales of seconds.

2) Both spin-down in the oscillation train and obvious examples of simultaneous signals separated by around 1 Hz are relatively rare, occurring in 3 and 2 bursts out of 65 respectively (Figure 3.2 and 3.3).

3) The amounts of the frequency drifts appear larger when oscillations are observed earlier in the bursts (Figure 3.6), but are uncorrelated with the durations of the oscillations. This implies that the starts of the bursts set in motion the eventual changes in the oscillation frequencies, as opposed to some mechanism that operates when the oscillations themselves appear.

4) Photospheric radius expansion appears to temporarily interrupt oscillations. All of the 16 oscillation trains that appear continuously from the rise to the tail of the burst fail to exhibit radius expansion. When radius expansion does occur, oscillations are often seen in the rise and/or the tail of the burst, but rarely during the radius expansion (Figure 5).

5) Aside from the above points, the time scale and magnitude of the frequency evolution does not appear correlated with other properties of the burst, such as its decay time scale, peak flux, or fluence.

6) In about 60% of the oscillation trains, the frequency evolution displays inflections that cannot be described by an exponential model. Two illustrative examples are shown in Figure 6. Thus, in addition to the upward frequency drift, the frequency tends to wander by on order 0.1 Hz in a random manner.

Several models have been proposed in which the oscillations originate from an anisotropy in the emission from the neutron star’s surface. Each models the frequency evolution as originating from different mechanisms (Strohmayer et al., 1997; Cumming et al., 2002; Heyl, 2002; Spitkovsky et al., 2002). The apparent lack of stability in the underlying clock producing these oscillations favors the models in which the brightness asymmetry originates from hydrodynamic instabilities (Spitkovsky et al., 2002) or modes excited in the neutron star ocean (Heyl, 2002). Future studies of the amplitude, harmonic structure, energy spectrum, and phase lags of the oscillations will further constrain the location of the hot spot and probe the atmosphere of the neutron star.

We thank Lars Bildsten and Fred Lamb for useful comments, and the referee for careful reading of this manuscript. We also thank Pavlin Savov for developing the burst detection software, and Derek Fox for developing much of the phase connection software. This work was supported by NASA, under contract NAS 5-30612 and grant NAG 5-9184.


format: text

text:
# Generation of Flows in the Solar Chromosphere Due to Magnetofluid Coupling

## 1 Introduction

Recent observations, strongly fortified by immensely improved measuring and interpretive capabilities, have convincingly demonstrated that the solar corona is a highly dynamic arena replete with multiple–scale spatiotemporal structures (Aschwanden et al. 2001a). A major new advance is the discovery that strong flows are found everywhere — in the subcoronal (chromosphere) as well as in the coronal regions (see e.g. (Schrijver et al. 1999; Winebarger, LeLuca and Golub 2001; Wilhelm 2001; Aschwanden et al. 2001a; Aschwanden et al. 2001b; Seaton et al. 2001; Winebarger et al. 2002) and references therein). Equally important is the growing belief and realization that the plasma flows may complement the abilities of the magnetic field in the creation of the amazing richness observed in the coronal structures.

Interestingly enough, even before the observational mandate, theoretical efforts in harnessing the plasma flows to solve some of the riddles of solar physics had already begun. In particular the dynamics of flow–based structure creation and heating was the subject of Mahajan et al. (1999, 2001); in this model the flows provided the basic material as well as energy for the primary heating of the coronal loops. A systematic treatment of loop models that include flows was also developed by Orlando, Peres and Serio (1995a, 1995b), and by Mahajan et al. (1999, 2001).

If flows are to play an important and essential role in determining the dynamics and structure of the solar corona, we must immediately face the problem of finding sources and mechanisms for the creation of these flows. Catastrophic models of flow production in which the magnetic energy is suddenly converted into bulk kinetic energy (and thermal energy) are rather well–known; various forms of magnetic reconnection (flares, micro and nanoflares) schemes permeate the literature (see e.g. (Wilhelm 2001; Christopoilou, Georgakilas and Koutchmy 2001) for chromosphere up–flow generations). A few other mechanism of this genre also exist: Uchida et al. (2001) proposed that the major part of the supply of energy and mass to the active regions of the corona may come from a dynamical leakage of magnetic twists produced in the subphotospheric convection layer; Ohsaki, Shatashvili, Yoshida and Mahajan (2001,2002) have shown how a slowly evolving closed structure (modelled as a double–Beltrami two–fluid equilibrium) may experience, under appropriate conditions, a sudden loss of equilibrium with the initial magnetic energy appearing as the mass flow energy. Another mechanism, based on loop interactions and fragmentations and explaining the formation of loop threads, was given in Sakai and Furusawa (2002). These mechanisms, though extremely interesting, are an unlikely source to account for the observed ubiquity of plasma flows. One should, perhaps, look for relatively gentler, more widespread, and steadier mechanisms.

Before exploring new avenues for flow generation, and deciding which region of the solar atmosphere should be the first target for investigation, we seek some guidance from phenomenology. Based on estimates of energy fluxes required to heat the chromosphere and the corona, Goodman (2001) has shown that the mechanism which transports mechanical energy from the convection zone to the chromosphere (to sustain its heating rate) could also supply the energy needed to heat the corona, and accelerate the solar wind (SW). The coronal heating problem, in this context, is shifted to the problem of the dynamic energization of the chromosphere. In the latter process the role of flows is found to be critical as warranted by the following observations made in soft X–rays and extreme ultraviolet (EUV) wavelengths, and recent findings from the Transition Region and Coronal Explorer (TRACE): the overdensity of coronal loops, the chromospheric up–flows of heated plasma, and the localization of the heating function in the lower corona (e.g. Schrijver, et al. 1999; Aschwanden et al. 2001a; Aschwanden 2001b) and references therein). The main message then, is that the coronal heating problem may only be solved by including processes (including the flow dynamics) in the chromosphere and the transition region. The challenge, therefore, is to develop a semi–steady state theory of flow generation in the chromosphere.

There are only two obvious energy sources that could power flow generation in the chromosphere: the magnetic field, and the thermal pressure of the plasma. We have already mentioned a few examples of the magnetically driven transient but sudden flow–generation. Looking for a more quiescent pathway, we shall now concentrate on the latter source. To convert thermal energy into kinetic energy, some variant of the Bernoulli mechanism (the existence of larger kinetic energy in regions of lower pressure) must be invoked. We will soon show that the double–Beltrami–Bernoulli states accessible to a two–fluid system in which the velocity field is formally treated at par with the magnetic field (Mahajan and Yoshida 1998), can readily provide the necessary framework for this conversion.

## 2 Model

In astrophysics (particularly in the physics of the solar atmosphere), it is useful to assign at least two distinct connotations to the plasma “flow”: 1) the flow is a primary object whose dynamics bears critically on the phenomena under investigation. The problems of the formation and the original heating of the coronal structure, the creation of channels for particle escape, for instance, fall in this category, 2) The flow is a secondary feature of the system, possibly created as a by-product (e.g. see Ohsaki et al. (2002)) and/or used to drive or suppress an instability. Since the generation of flows which will eventually create the coronal loops is the theme of this effort, the flows here are fundamental.

As pointed out earlier, our theoretical model is based on a rather simple application of the magnetofluid theory developed over the last few years. We plan to restrict ourselves to almost steady state considerations (for a steady and continuous supply of plasma flows emerging from the chromosphere). Very near the photospheric surface, the influence of neutrals and ionization processes (and processes of flux emergence etc.) would not permit a quasi–equilibrium approach. A little farther distance ($`\mathrm{\Delta }r`$) from the surface, however, we expect that the quasi–equilibrium two–fluid model will capture the essential physics of flow generation.

¿From recent observational data (see e.g. Goodman (2000); Aschwanden et al. (2001a); Socas-Navarro and Almeida (2002) and references therein) we could obtain the following average plasma density and temperature at $`\mathrm{\Delta }r\left(5002000\right)`$km: $`n\left(10^{15}10^{11}\right)\mathrm{cm}^3;T1`$eV (for simplicity we will assume equal electron and ion temperatures). The information about the magnetic field is a little harder to extract because of the low sensitivity and lack of high spatial resolution of the measurements coupled with the inhomogeneity and co–existence of small– and large–scale structures with different temperatures in nearby regions. At these distances we have different values for the network and for the internetwork fields. The network plasmas have typically short–scale fields in the range $`B_0\left(7001500\right)`$G, have more or less uniform density and will be prone to explosive/eruptive analysis of the kind carried out in Ohsaki et al. (2002). The internetwork fields, on the other hand, are generally smaller (with some exceptions (Socas-Navarro and Almeida (2002)) — $`B_o500`$G, and are embedded in larger–scale plasma structures with inhomogeneous densities. The theory of steady creation of flows in the lower chromosphere will be based on these latter objects.

## 3 Generation of flows in lower chromosphere

To illustrate the basic physics of flow creation, we deal with the simplest two–fluid equilibria with $`T=\mathrm{const}`$ leading to $`n^1pT\mathrm{ln}n`$. The generalization to a homentropic fluid with $`p=\mathrm{const}\mathrm{n}^\gamma `$ is straightforward and is done in the numerical work.

The dimensionless equations describing the model equilibrium can be read off from Mahajan et al. (2001)

$$\frac{1}{n}\times 𝐛\times 𝐛+\left(\frac{r_{A0}}{r}\beta _0\mathrm{ln}n\frac{V^2}{2}\right)+𝐕\times \left(\times 𝐕\right)=0,$$
(1)
$$\times \left(𝐕\frac{\alpha _0}{n}\times 𝐛\right)\times 𝐛=0,$$
(2)
$$\left(n𝐕\right)=0,$$
(3)
$$𝐛=0,$$
(4)
where the notation is standard with the following normalizations: the density $`n`$ to $`n_0`$ prevalent at about $`500`$km (and farther) from the solar surface (this is where the plasma is created), the magnetic field to the ambient field strength at the same distance, and velocities to the Alfvén velocity $`V_{A0}`$. The parameters $`r_{A0}=GM_{}/V_{A0}^2R_{}=2\beta _0/r_{c0},\alpha _0=\lambda _{i0}/R_{},\beta _0=c_{s0}^2/V_{A0}^2`$ are defined with $`n_0,T_0,B_0`$. Here $`c_{s0}`$ is a sound speed, $`R_{}`$ is the solar radius, $`r_c=GM_{}/2c_{s0}^2R_{}`$, and $`\lambda _{i0}=c/\omega _{i0}`$ is the collisionless skin depth.

The double–Beltrami solutions, expressing the simple physics that the electrons follow the field lines, while the ions, due to their inertia, follow the field lines modified by the fluid vorticity, are contained in the pair,

$$𝐛+\alpha _0\times 𝐕=dn𝐕,𝐛=an\left[𝐕\frac{\alpha _0}{n}\times 𝐛\right],$$
(5)
where $`a`$ and $`d`$ are dimensionless constants related to the two ideal invariants, the magnetic and the generalized helicities (Mahajan and Yoshida 1998; Mahajan et al. 2001),

$$h_1=\left(𝐀𝐛\right)d^3x,$$
(6)
$$h_2=\left(𝐀+𝐕\right)\left(𝐛+\times 𝐕\right)d^3x.$$
(7)
On substituting (6)–(7) into (1)–(2), one obtains the Bernoulli condition

$$\left(\frac{2\beta _0r_{c0}}{r}\beta _0\mathrm{ln}n\frac{V^2}{2}\right)=0,$$
(8)
relating the density with the flow kinetic energy, and solar gravity. Equations (1), (5), and (8) represent a close system, and may be easily manipulated to yield ($`g\left(r\right)=r_{c0}/r`$)

$$\frac{\alpha _0^2}{n}\times \times 𝐕+\alpha _0\times \left[\left(\frac{1}{an}d\right)n𝐕\right]+\left(1\frac{d}{a}\right)𝐕=0,$$
(9)
$$\alpha _0^2\times \left(\frac{1}{n}\times 𝐛\right)+\alpha _0\times \left[\left(\frac{1}{an}d\right)𝐛\right]+\left(1\frac{d}{a}\right)𝐛=0,$$
(10)
$$n=\mathrm{exp}\left(\left[2g_0\frac{V_0^2}{2\beta _0}2g+\frac{V^2}{2\beta _0}\right]\right),$$
(11)
a set ready to be solved for the density, the velocity and the magnetic field. We must point out that this time–independent set is not suitable for studying chromospheric heating processes (primary heating at lower heights). The main thrust of this paper is to uncover mechanisms which create flows in the chromosphere — flows that will supply matter and energy needed to create the coronal structures, and provide their primary heating. The creation and heating problem, of course, requires a fully time–dependent treatment (Mahajan et al. 2001).

We have to resort to numerical methods to obtain detailed solutions for the coupled nonlinear system (9), (10), and (11). We have carried out a 1D simulation (the relevant dimension being the height “Z” from the center of the Sun; $`Z_0=R_{}+\mathrm{\Delta }r`$ is the surface at which the boundary conditions are applied) for a variety of boundary conditions. The boundary surface is so chosen that at this height $`Z_0`$ the ionization is expected to be complete. For estimates of $`\mathrm{\Delta }r`$ given earlier, the relevant heights lie within $`\left[\left(1+0.710^3\right)\left(1+2.810^3\right)\right]R_{}`$.

The simulation results are presented in Figs. 1–2. These are the plots of various physical quantities as functions of the height. The first figure consists of three frames (a–b, c–d, and e–f) each consisting of two pictures — one for the density and the magnetic field and the other for the velocity field. The parameters defining different frames are: $`n_010^{14}\mathrm{cm}^3`$, $`B_0500`$G, $`V_{A0}110`$km/s implying $`\beta _00.011`$ and $`r_{c0}=900`$ for a–b; $`n_010^{16}\mathrm{cm}^3`$, $`B_01000`$G, $`V_{A0}22`$km/s implying $`\beta _00.2<1`$ and $`r_{c0}=950`$ for c–d; and $`n_010^{17}\mathrm{cm}^3`$ and $`B_01500`$G, $`V_{A0}11`$km/s, $`\beta _01`$ and $`r_{c0}=1000`$ for the frame e–f. In each frame there are three sets of curves labelled by $`\alpha _0`$ (1–2–3 corresponding respectively to $`\alpha _0=10^5,\mathrm{\hspace{0.17em}10}^3,\mathrm{\hspace{0.17em}10}^1`$), the measure of the strength of the two–fluid Hall currents.

For all our runs the temperature at the boundary was taken to be $`1`$eV ($`c_{s0}10`$km/s), and the boundary conditions, $`\left|b_0\right|=1,V_0=1`$km/s (with $`V_{x0}=V_{y0}=V_{z0}`$) were imposed. Notice that we begin with just a small residual flow speed. The choice, $`da100`$ and $`\left(ad\right)/a^210^6`$ for the parameters characterizing the double Beltrami state, reflects the physical constraint that we are dealing with a sub–Alfvénic flow with a very small $`\alpha _0`$ (Mahajan et al. 1999). We must admit that the values of $`\alpha _0`$ chosen for the simulation are much larger than their actual values ($`10^8`$ for corona and $`10^{11}`$ for chromosphere); our present code cannot resolve the equivalent short lengths, though, we hope to do better in future. We believe, however, that the nature of the final results is properly captured by these artificial values of $`\alpha _0`$.

The most remarkable result of the simulation is that for small and realistic values of $`\alpha _0`$ (curves labelled 1), there exists some height where the density begins to drop precipitously with a corresponding sharp rise in the flow speed. The effect is even stronger for the low beta (a–b are the lowest beta frames) plasmas. It is also obvious that at very short distances, the stratification is practically due to gravity, but as we approach the velocity “blow–up” height, the self–consistent magneto–Bernoulli processes take over and control the density (and hence the velocity) stratification.

An examination of the Bernoulli condition (11) readily yields an indirect estimate for the height at which the observed shock formation may take place. For a low–beta plasma, the sharp fall in density is expected to occur where the local flow kinetic energy exceeds the kinetic energy specified at the boundary (this is true for all $`\alpha _0`$), i.e,

$$\left|𝐕\right|^2V_0^2>2\beta _0.$$
(12)
For the current simulation, at $`\beta _0=0.01`$, it occurs approximately at $`\left|𝐕\right|^2>0.02`$ or at $`\left|V\right|0.14`$. This analytically–predicted value is very close to the simulation result (see Fig. 2(b)). Simulation results also confirm that the velocity blow–up distance depends mainly on $`\beta _0`$, and that the final velocity is greater for greater $`\beta _0`$ (Fig. 2(a)). The data presented in Fig. 1 and Fig. 2 corresponds to a uniform temperature plasma. For this case, the variations in plasma pressure are entirely due to the variations in density. Since the magnetic energy remains practically uniform over the distance, sharp decrease in density with a corresponding sharp rise in the flow–speed (flow energy changes are of the order of $`n^{1/2}`$) is nothing but the expression of the commonly understood Bernoulli effect. We must emphasize that the general results remain unchanged in our extensive simulations in which the temperature is allowed to vary (but we have to use a homentropic equation of state to analytically derive the Beltrami states. The final parameters, naturally, depend upon the adiabaticity index $`\gamma `$).

Taking into account the fact that the fiducial height $`Z_0`$ is different for different cases (larger for smaller $`\beta _0`$), one expects that the plots for the “blow–up” distance will approximately match one other, and there exists some minimum distance (possibly of the order $`500`$km or so) from the photosphere below which steady up–flows can not be generated. This expectation is simply supported by the numerical results given in Fig. 2(a). Observations clearly indicate that there exists a narrow layer above the photosphere where no spicules, mottles etc. are seen.

To check whether the generated flows are predominantly radial or somewhat more isotropic (to explain the observational constraints) we studied in detail the relatively large $`\beta _0`$ case (fixing $`\beta _0`$ is quite difficult due to complications like ionization) and found that the flows tend to be mostly radial only for large $`\alpha _0`$ (see, for example, plots labeled 2 and 3 in Fig. 1(b,d,f)). The situation could change considerably when we deal with a more inclusive time–dependent dynamical model with dissipation. Plasma heating, then, could result from the dissipation of the perpendicular energy so that at larger distances, the flows would have larger radial components. Heating would also keep $`\beta (𝐫,t)`$ large at upper heights shifting the velocity blow–up distance further or eliminating it all together; we know from Fig. 1 that as $`\beta _0`$ goes up, the density fall(velocity amplification) becomes smoother. These issues will be dealt with later in a more detailed work.

Note, that if one were to ignore the flow term in (8) (a totally wrong assumption commonly used in many studies), we will end up finding essentially radial flows. The magnitude of these flows, however, remains small; there is no region of sharp rise (12), and the generated flows achieve reasonable energies at heights typically 10 times greater than the heights at which the correct Bernoulli condition would do the trick.

## 4 Conclusions and Summary

We have shown a possible pathway for a steady generation of flows in the quasi–equilibrium stage established through the interaction of the emerging magnetic fluxes with the existing cold solar chromosphere (when an ionized $`1`$eV plasma is trapped in internetwork structures). The suggested mechanism is a straightforward application of the recently–developed magnetofluid model (Mahajan and Yoshida 1998; Mahajan et al. 1999, 2001); a generalized Bernoulli mechanism (a necessary condition for the double–Beltrami magnetofluid equilibrium) allows the pressure energy to be very effectively transformed to flow kinetic energy as the plasma moves away from the sun. We find that at reasonable heights and for realistic plasma parameters, there is a precipitous pressure fall with a sharp amplification of the flow speed. In the presence of dissipation, these flows are likely to play a fundamental role in the heating of the inner and upper chromosphere, although our explicit purpose in this paper was to create a steady source of matter and energy for the formation and primary heating of the corona. Our preliminary results agree with the observation data, and lend promise to attempts, based on the existence of subcoronal flows, to tackle unresolved problems like the coronal heating and origin of the solar wind.

## ACKNOWLEDGEMENTS

Authors acknowledge the help of Dr. K.I. Sigua for preparing the simulation pictures and thank Dr. E. Marsch for interesting and useful discussions. All the authors thank Abdus Salam International Centre for Theoretical Physics, Trieste, Italy. The work of SMM was supported by USDOE Contract No. DE–FG03–96ER–54366. Work of KIN was supported by a Grant No. 02–02–16199a of Russia Fund of Fundamental Research (RFFR). The work of NLS was supported by ISTC Project G–633.


text:
# Infrared spectroscopy of the largest known trans-neptunian object 2001 KX76

## 1 Introduction

The trans-neptunian object (TNO) 28976 = 2001 KX76 was discovered in May 2001, and it is the largest known TNO (Millis et al. Millis (2001)). Assuming a geometric albedo of 4% its diameter is $``$1100 km. TNOs are remnant planetesimals from the early solar system formation stages, populating a region just beyond the orbit of Neptune (Edgeworth Edgeworth49 (1949); Kuiper Kuiper51 (1951)), called the Edgeworth-Kuiper belt (EKb). They comprises three dynamical classes: objects in the 3:2 mean motion resonance with Neptune have been described as “Plutinos” (Jewitt & Luu JewLuu96 (1996)), those beyond about $`41AU`$ as “Classical EKb Objects”, and those with a much larger semi-major axis and higher eccentricity than the previous classes are known as “Scattered Disk Objects” (Luu et al. Luu97 (1997)). 2001 KX76 has orbital elements that fit with those of the Classical EKb Objects.

The low temperature in the EKb and the low density of TNOs imply they are probably the most pristine objects in the Solar System. So they can provide key information on the composition and early conditions of the pre-solar nebula. Though the study of their surface properties is very important from a cosmogonical point of view, as it could provide important clues to understand the conditions existing at the beginning of the solar system.

A powerful method for remote determination of the composition of volatile surface component of the outer solar system objects is the Near-infrared spectroscopy (Brown & Cruikshank BrownCrui97 (1997)). Due to the small size and large distance of this objects, this technique is actually limited to the largest members of the EKb. Very few near-infrared spectra of TNOs have been already published: Brown et al. (Brownetal97 (1997)) present the spectrum of 1993 SC; Luu & Jewitt (LuuJew98 (1998)) the spectrum of 1996 TL66; Brown et al. (BrownCrPe99 (1999)) the spectrum of 1996 TO66; Brown et al. (Brownetal2000 (2000)) the spectrum of 2000 EB173; Licandro et al. (Lic (2001)) the spectra of (20000) Varuna (= 2000 WR106) and 2000 EB173; and Jewitt & Luu (JewLuu2001 (2001)) the spectra of 2000 EB173, 1999 DE9, 1996 TS66, and 1993 SC. Even the number of published spectra of TNOs is very low, large different surface properties have been inferred among them. Strong color variations, from neutral to very red objects, have been also reported photometrically by several authors (e.g. Jewitt & Luu JewLuu2001 (2001), Gil-Hutton & Licandro GilLic (2001), Davies et al. Davies00 (2000), Barucci et al. Baruccietal00 (2000)). This color diversity is confirmed in the reported infrared spectra. But also some absorption bands possibly due to hydrocarbons and/or water ice were observed in the infrared spectrum of some objects (1993 SC, 1999 DE9, 1996 TO66, 20000 Varuna), while the spectra of other TNOs (1996 TL66, 1996 TS66, 2000 EB173) are featureless in the observed range. This diversity of surface compositions has been observed also in similar objects like Centaurs and irregular satellites (Brown Brown2000 (2000)). The observed surface diversity can be atributed to intrinsically different compositions among TNOs and/or to some resurfacing processes (Jewitt & Luu JewLuu2001 (2001); Gil-Hutton Gil (2002)).

## 2 Observations

We have obtained low resolution spectra of 2001 KX76 on July 6.0 UT and July 12.0 UT, 2001, with the 3.6m Telescopio Nazionale Galileo (TNG), using NICS, the near-infrared camera and spectrometer (see Baffa et al. Baf (2001)). Among the many imaging and spectroscopic observing modes, NICS offers a unique, high throughput, low resolution spectroscopic mode with an Amici prism disperser (Oliva oliva01 (2001)), which yields a complete 0.9-2.4 $`\mu `$m spectrum. A $`1.5\mathrm{"}`$ width slit corresponding to a spectral resolving power $`R34`$ and quasi-constant along the spectrum, has been used. The low resolution together with the high efficiency of the Amici prism (about 90% across the useful wavelength range) allowed us to obtain spectra of faint objects like TNOs with a four meter class telescope for the first time (Licandro et al. Lic (2001)), and with the advantage of having the whole infrared range measured simultaneously.

The identification of the TNO was done by taking series of images through the Js filter separated by one our, and by comparing them. The object was identified as a moving object at the predicted position and with the predicted proper motion. The slit was oriented in the paralactic angle, and the tracking was at the TNO proper motion. The acquisition consisted of a series of 5 images of 60 seconds exposure time in one position (position A) of the slit and then offsetting the telescope by $`30\mathrm{"}`$ in the direction of the slit (position B). This process was repeated and a number of ABBA cycles were acquired. The total exposure time was 3600 seconds each night. The reduction and calibration of the spectra was done as in Licandro et al. (Lic (2001)).

To correct for telluric absorption and to obtain the relative reflectance, the G star Land (SA) 110-361 (Landolt Landolt (1992)), which has visible colors very similar to that of the Sun, was observed during the same night just after 2001 KX76, and at a similar airmass. Land (SA) 110-361 was observed also in previous nights together with the solar analogue star P330E (Colina & Bohlin Col (1997)). Both stars present similar spectra in the infrared region, so we used Land (SA) 110-361 as a solar analogue. The spectrum of the TNO was divided by the spectrum of the solar analogue star, and then normalized to unity around 1.6 $`\mu `$m, thus obtaining the relative reflectance plotted in Fig. 1. Around the telluric water band absorptions the S/N of the spectrum is very low, additionally, the telluric absorption varies between the TNO spectra and the standard stars spectra introducing false spectral features. Therefore, these parts are not included in the final spectra.

## 3 Discussion

The spectra of 2001 KX76 are shown in Fig. 1. Both are identical withing the noise. This means that the object was observed at a very similar rotational phase or that the composition of its surface is very homogeneous. The spectrum of 2001 KX76 is featureless and, compared with the very red slope spectra of Varuna and 2000 EB173 (Licandro et al. Lic (2001)), is almost neutral (see Fig. 2).

The spectrum of 2001 KX76 do not show the water ice absorption features at 1.5 and 2.0 $`\mu `$min observed in Varuna. Similar objects with neutral colors present a large variation in the detected amount of water ice in the surface. Brown et al. (BrownCrPe99 (1999)) detected water ice in 1996 TO66 while Luu & Jewitt (LuuJew98 (1998)) did not in 1996 TL66. Also red TNOs presents this variation in the detected amount of water ice in their surfaces (see Licandro et al. Lic (2001)). Though it seems that the diversity between red and neutral TNOs has no relation with the presence of water ice in the surface.

The observed very red color of some TNOs is probably due to an evolved surface which is the result of long term irradiation by solar radiation, solar wind, and galactic cosmic-rays. This results in the selective loss of hydrogen and the formation of an “irradiation mantle” of carbon residues (Moore et al. Mooreetal83 (1983); Johnson et al. Johnetal84 (1984); Strazzulla & Johnson StrazzullaJo91 (1991)). This process makes that an initially neutral color and high albedo ice becomes redish. But further irradiation gradually reduces the albedo at all wavelengths, and the material becomes very dark, neutral in color, and spectrally featureless (Andronico et al. Androetal87 (1987); Thompson et al. Thompsonetal87 (1987)). Gil-Hutton (Gil (2002)) shows that the total amount of radiation received by TNOs during the age of the solar system can produce very dark, neutral surfaces. The neutral color and the lack of ice features in the spectrum of 2001 KX76 suggest that the surface of this object has a very evolved irradiation mantle, that has not suffered from a strong collisional resurfacing process. In contrary, the presence of water ice in the surface of Varuna suggest that it has undergone a strong resurfacing due to collisions, capable of retrieve “fresh” material from below the irradiation mantle. Even if the collisional resurfacing process were not the primary cause of the observed color dispersion of the TNOs according to Jewitt & Luu (JewLuu2001 (2001)), it could be the reason to explain why water ice is observed in the very evolved surface of these objects. Thus, more infrared spectroscopy of a large number of TNOs is needed.

## 4 Conclusions

Two spectra of 2001 KX76, obtained in different nights, are presented. Both are almost identical withing the noise, suggesting that the TNO was observed at a very similar rotational phase or that the composition of its surface is very homogeneous. The spectra are featureless and correspond to a neutral colored object. This suggest that the surface of 2001 KX76 is probably highly evolved by long term radiation and that any possible collisional resurfacing process did not erode significantly its irradiation mantle.

###### Acknowledgements.

JL thanks Bobby Bus and Humberto Campins for their usefull suggestions on faint solar analogue stars to be observed for calibrations purposes with large telescopes. This paper is based on observations made with the Italian Telescopio Nazionale Galileo (TNG) operated on the island of La Palma by the Centro Galileo Galilei of the CNAA (Consorzio Nazionale per l’Astronomia e l’Astrofisica) at the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofisica de Canarias. We are grateful to all the technical staff and telescope operators for their assistance during the commissioning phase of NICS.


text:
# Spatial solitons in an optically pumped semiconductor microresonator

\[

## Abstract

We show experimentally and numerically the existence of stable spatial solitons in an optically pumped semiconductor microresonator. We demonstrate that the pump substantially reduces the light intensity necessary to sustain the solitons and thereby reduces destabilizing thermal effects. We demonstrate coherent switching on and off of bright solitons. We discuss differences between pumped and unpumped below-bandgap-solitons.

\] Spatial solitons in semiconductor microresonators  are self-trapped light beams which can form due to both transverse (e.g., self-focusing) and longitudinal (nonlinear resonance ) nonlinear effects. These nonlinear effects can act in the same sense or oppositely, with the consequence of reduced soliton stability in the later case. Thus choice of nonlinear resonator parameters suited best for sustaining stable solitons is important. The best-suited schemes are unpumped microresonator excited above bandgap  and optically pumped microresonator excited below bandgap .

Here we pump a semiconductor microresonator optically in a range up to the lasing threshold and study formation of switched structures. We show experimentally and numerically the existence of the resonator solitons in the pumped microresonators. We demonstrate that the pumping substantially reduces the light intensity necessary to sustain and switch the solitons (diminishing thermal problems as a side effect) so that semiconductor laser diodes are sufficient for sustaining solitons. We discuss differences between pumped and unpumped solitons below bandgap.

We pump a bistable quantum-well (GaAs/GaAlAs) microresonator using either a tunable Ti:sapphire laser or a high power multimode laser diode. The pumped area of the resonator is illuminated additionally by a focused beam from a single-mode laser diode that provides quite large Fresnel number ($``$100) and near resonant illumination of the resonator. The main control parameters in the experiment are the resonator detuning and intensities of the illumination and pump.

Observations are done in reflection since the substrate of the microresonator structure is opaque at the working wavelengths: dark switched structures in reflection correspond to bright ones in transmission and vice versa. Switched structures formed in the illuminated beam cross section were monitored in the plane of the microresonator in two ways: (i) A CCD camera with electro-optical shutter recorded 2D snapshots of switched structures. (ii) A fast and small aperture photo detector monitored local dynamics.

Switched structures observed (Fig. 1 a,b) manifest themselves as resonator solitons: 1) they are of the size ($``$10 $`\mu `$m) expected for such solitons ; 2) they are round spots whose size and shape are independent on the intensity and shape of the illuminating beam (e.g., elliptical beam shape as shown in Fig. 1c); 3) they are robust against perturbations of the illuminating light intensity; 4) they are bistable, i.e. they can be switched on (Fig. 2a) and off (Fig. 2b) by sharply focused address pulses.

This soliton nature of the observed switched structures is supported by numerical simulations of the intracavity field structures (in 2D) using the model equations for the pumped semiconductor microresonator driven by a plane wave . Fig. 3 shows typical examples of calculated resonator solitons below bandgap. Bright solitons have a large existence range in the pumped case (Fig. 3b), dark solitons exist, though with smaller range of stability, in the unpumped case (Fig. 3a). Parameter domains of existence of resonator solitons are related to those of optical bistability for plane waves, are shown in Fig. 3.

Analysis shows that increase of the pump intensity leads to shrinking of the resonator solitons’ existence domain and shifting towards low intensity of the light sustaining the solitons (Such reduction of the sustaining light intensity was observed experimentally in ).

When the pump intensity approaches the transparency point of the semiconductor material, the resonator solitons’ domain of existence disappears. It reappears above the transparency point. In the experiment we have quite strong contribution of the imaginary part (absorption/gain) of the complex nonlinearity at the working wavelength (854 nm). Therefore the transparency point is very close to the lasing threshold so that inversion without lasing is difficult to realize.

Slightly above threshold we observe in presence of illumination structures (Fig. 4) reminiscent of the solitons in electrically pumped resonators .

In summary, optically pumped semiconductor resonators are well suited for sustaining solitons below bandgap:(i) background light intensity necessary to sustain and switch resonator solitons is substantially reduced by the pumping and therefore destabilizing thermal effects are minimized, (ii) above the transparency point only the dispersive part of semiconductor nonlinearity stabilizes a soliton, therefore the domain of existence of ”below bandgap” (purely dispersive) bright solitons can be quite large. Moreover optical as opposed to electrical pumping allows more homogeneous pumping conditions . This suggests that optically pumped resonators lend themselves more readily for localization and motion control of solitons then electrically pumped ones.

Acknowledgment
This work was supported by Deutsche Forschungsgemeinschaft under grant We743/12-1.