text:
# Long-term optical and X-ray observations of the old novae DI Lacertae and V841 OphiuchiBased on observations with the Apache Point Observatory 3.5-m telescope which is owned and operated by the Astrophysical Research Consortium, and on observations with the WIYN Observatory 3.5-m telescope which is jointly operated by the University of Wisconsin, Indiana University, Yale University, and the National Optical Astronomy Observatories.

## 1 Introduction

Cataclysmic variables (CVs) are semi-detached interacting binary stars composed of a white dwarf (WD) primary star and a low mass ($`0.5`$M) main sequence secondary star, with typical orbital periods of $`1`$ d. The Roche-lobe-filling secondary star loses mass through the inner Lagrangian point into an accretion disk formed around the WD (for non-magnetic CVs). Classical novae are a subclass of CV in which a thermonuclear runaway is triggered in a reservoir of matter that has been gradually accreted onto the WD. The resultant outburst produces a peak brightness increase of $``$ 6–15 mag, and releases $`10^{44}10^{46}`$ erg (see review in Warner, 1995, ch. 5).

The outbursts of novae are often well-observed, and the long term behavior of many dwarf novae and novalike CVs is monitored by the American Association of Variable Star Observers. However, there have been few programs for monitoring the variability of novae at quiescence. One such program, using an automated telescope called RoboScope (see §2.1), has monitored 22 old novae and 42 novalikes, most of them for over 9 years now. RoboScope has found several kinds of unusual photometric behavior. About one-third of the old novae have shown quasi-periodic variability for a year or two, interspersed with stable light curves. These variations do not appear stochastic as they repeat at similar periods when they reappear. Preliminary analysis of RoboScope light curves spanning several years for the old nova DI Lacertae suggested brightness oscillations on a time scale of $`35`$ d (Honeycutt et al., 1995), while the old nova V841 Ophiuchi also displayed prominent variability (Honeycutt et al., 1994). In order to determine if these oscillations could be related to a disk instability mechanism (as operates in dwarf novae) or to mass transfer or magnetic effects on the secondary star, we attempted a detailed study of DI Lac and V841 Oph, which are among the most active and brightest of the systems showing the oscillations.

DI Lac (= Nova Lac 1910) was a moderately fast nova that reached a maximum brightness of $`m_V(\mathrm{max})=4.3`$ mag during its outburst. V841 Oph (= Nova Oph 1848) was a slow nova that reached a similar brightness, $`m_V(\mathrm{max})=4.2`$ mag (Kukarkin et al., 1971). The reddening for both objects has been measured from International Ultraviolet Explorer ($`IUE`$) ultraviolet and ground-based optical spectra. The range for DI Lac is $`E(BV)=0.15`$–0.41 (Bruch, 1984; Cassatella & Gonzalez-Riestra, 1990), while for V841 Oph it is $`E(BV)=0.30`$–0.58 (Bruch, 1984; Cassatella & Gonzalez-Riestra, 1990; Weight et al., 1994; Verbunt et al., 1997). Nova shells were not detected in H$`\alpha `$ for either CV (Cohen, 1985). In X-rays, DI Lac was undetected in the $`ROSAT`$ PSPC All Sky Survey, with a $`2\sigma `$ upper limit of 0.013 cts s-1 in the 11–201 channel (0.1–2.0 keV) energy range during 299 s total exposure time, while V841 Oph was a marginal detection with $`0.017\pm 0.007`$ cts s-1 in the “hard” 0.9–2.0 keV ROSAT energy range ($`0\pm 0.025`$ cts s-1 in the full ROSAT energy range) during 402 s total exposure time (Verbunt et al., 1997).

We report here on the results of ground-based optical photometry and spectroscopy, and Rossi X-ray Timing Explorer ($`RXTE`$) X-ray observations of DI Lac and V841 Oph. The optical photometry comprises an almost decade-long light curve for each star, while the contemporaneous spectroscopy and X-ray observation repeatedly sampled each nova during separate intervals of $`45`$–55 d in length.

## 2 Observations and Analysis

Our long-term optical photometry of DI Lac and V841 Oph was accomplished from 1990–1998, while the spectroscopic amd X-ray observations took place in 1997. These observations are summarized in Tables 1 and 2, and are discussed in detail below.

### 2.1 Optical Photometry

Our optical photometry data were acquired by RoboScope, a 41-cm telescope in Indiana equipped for automated differential CCD stellar photometry (Honeycutt & Turner, 1992). All observatory operations (including data reductions) are accomplished as fully unattended and unsupervised tasks, which makes practical the acquisition of long homogeneous data streams. Typically, RoboScope obtains one or two 4-min exposures per clear night for each of $`140`$ program stars. The data are reduced using the method of incomplete ensemble photometry (Honeycutt, 1992). For DI Lac, 85 ensemble stars were used in 987 exposures over 9 observing seasons from 1990 November to 1998 November. For V841 Oph, 24 ensemble stars were used in 611 exposures over 8 observing seasons from 1991 May to 1998 September. The zero-points of the differential light curves were established using secondary standard stars from Henden & Honeycutt (1995, 1997). Six such secondary standards were used for V841 Oph, while 13 were employed for DI Lac. Typical $`1\sigma `$ uncertainties in the calibrated magnitudes are on the order of 0.01–0.05 mag. Over the entire multiyear range of the RoboScope light curves, the mean magnitudes were $`V14.55`$ for DI Lac and $`V13.6`$ for V841 Oph, with full ranges of variability of $`\mathrm{\Delta }V0.9`$ mag and $`\mathrm{\Delta }V1`$ mag, respectively. The RoboScope light curves of DI Lac and V841 Oph are shown in Figures 1 and 2.

With only limited excursions to anomalously fainter magnitudes, the typical minimum magnitude of $`14.7`$ for DI Lac has a sharp boundary at any given epoch. However, the minimum brightness displays a slight, apparently linear, trend towards increasing brightness over the 8 years of RoboScope observation. If we ignore for the moment the small number of data points that fall below the well-defined low brightness limit seen in the top panel of Figure 1, then the typical minimum brightness of DI Lac changes from $`V14.75`$ mag at the beginning of the RoboScope coverage, to $`V14.67`$ mag at the end. This corresponds to a mean increase of the minimum brightness by $`\mathrm{\Delta }V0.01`$ mag yr-1. The maximum brightness of DI Lac is much less uniform, and it is difficult to determine if the maximum brightness behaves in the same manner as the minimum brightness. A linear fit to the complete DI Lac light curve (with each data point weighted by the inverse square of its $`1\sigma `$ uncertainty) gives the rate of change of the mean magnitude as $`<+0.001`$ mag yr-1; that is, the mean brightness of DI Lac is essentially constant. However, if we exclude the 2 excursions to fainter magnitudes at HJD 2450600 and HJD 2451000 (as well as the excursion to brighter magnitudes at HJD 2450300), then the slope of the linear fit changes to $`\mathrm{\Delta }V=0.004`$ mag yr-1; that is, a trend towards increasing mean brightness.

The faint magnitude excursions in the light curve of DI Lac are very similar to the “dips” seen by Honeycutt et al. (1998a) in long-term RoboScope light curves of five old novae and novalike CVs. In the latter systems, the dips were found to often be paired with a preceding or following outburst, although the dips also sometimes occurred as isolated events. Honeycutt et al. (1998a) do not find a clear mechanism (e.g. disk instabilities, truncated disks, mass transfer modulation) that is responsible for the dips, but conclude that the overall photometric behavior (dips and outbursts) of the novae and novalikes that they studied is likely to be governed by either a combination of disk and mass transfer events, or by mass transfer events alone.

In contrast to the linear trend in the minimum brightness of DI Lac, the minimum brightness of V841 Oph appears to follow an almost sinusoidal trend that completes approximately 1.5–2 cycles during the 7.5 year RoboScope coverage. As seen in the top panel of Figure 2, the minimum brightness of V841 Oph varies from $`V14.5`$ mag to $`V13.8`$ mag. The maximum brightness of V841 Oph, which ranges from $`V13.1`$ mag to $`V13.4`$ mag, displays the same behavior; that is, when the minimum brightness at a given epoch is at its faint (bright) level, the maximum brightness is also at its faint (bright) level. We note that the magnitude range spanned by the minimum brightness in V841 Oph is about twice as large as that spanned by the maximum brightness. Although a linear fit to this overall trend is perhaps not the optimum choice, the corresponding rate of mean magnitude change is $`\mathrm{\Delta }V=+0.010`$ mag yr-1. It is possible (but not required by the data!) that the apparently linear trend in the light curve of DI Lac may also be a sinusoid, but with a much longer cycle length than in V841 Oph.

We searched for periodicities in the range 10–100 days in the RoboScope light curve data in two ways: by applying the phase dispersion minimization (PDM) algorithm (Stellingwerf, 1978), as well as performing an independent power spectrum analysis using the clean algorithm (Roberts et al., 1987). We performed the period search on both the entire data set for each CV, as well as on the data subsets from individual observing seasons (indicated in Figures 1 and 2). The most prominent period detected in each data set is listed in Table 3 (in several cases, most notably for the combined data sets for each nova, two significant periods are listed). In general, the two period search methods gave equivalent results to within better than 1–4 days. All of the detected time scales appear to correspond to quasiperiodic behavior rather than truly periodic variability such as that resulting from, for example, an eclipsing orbit. The time scales have significant “flexibility” (on the order of 1–5 days for periods $`50`$ d, 5–20 days for periods $`50`$ d) over the multi-year length of the light curves (as shown by substantially broadened dips and peaks in the PDM and clean analyses, respectively). Phase-binned light curves of both CVs folded on the shorter period found for each combined data set are shown in Figure 3. We note that the true amplitude of variation in a given cycle or observing season will be larger than the 0.1–0.2 mag suggested by Figure 3, in which data from many cycles with slightly different amplitudes and periods have been averaged together. Inspection of Figures 1 and 2 shows amplitudes of 0.4–0.8 mag for the full range of variability within each observing season. The long-term sinusoidal trend in the light curve of V841 Oph is also well-fit by a periodicity in the range 1800–1900 d or, with slightly less agreement, 1250–1300 d. (We note that these two period ranges are approximately related by the ratio 3:2, so we are likely seeing an aliasing effect due to the fact that only 1.5–2 cycles of these long periods are contained in the light curves.)

### 2.2 Optical Spectroscopy

We obtained optical spectra of both CVs using the Double Imaging Spectrograph on the Apache Point Observatory (APO) 3.5-m telescope (Gillespie et al., 1995)111also see http://www.apo.nmsu.edu/ during University of Washington share time. Spectra of DI Lac were obtained in July and August of 1997; spectra of V841 Oph were obtained in May of 1997. The APO spectra have a resolution of $`2`$Å and cover simultaneous wavelength ranges of $`4200`$–5000Å and $`5800`$–6800Å. The raw spectrum images were reduced in the standard fashion using IRAF. The instrumental response was removed via spectra of standard stars (Massey et al., 1988) obtained on the same nights; however, slit losses due to guiding errors rendered the absolute flux calibration unreliable. Consequently, we have normalized the spectra to a constant continuum level of 1.0.

A representative spectrum of DI Lac is shown in Figure 4.

Note the prominent Balmer absorption troughs containing narrow, central emission cores. Our other two APO spectra of DI Lac (not shown) are essentially identical to this one, with the exception that the Balmer emission cores in the August spectrum are somewhat stronger (relative to He II) than in the July spectra. The He I $`\lambda 4471`$ feature (see inset) is present with the same profile shape in all three of our spectra. It is suggestive of the absorption trough with central emission as seen in the Balmer lines. The feature labeled with “?” is likely an instrumental or reduction artifact since its wavelength does not match any identifiable line and it occurs in all of our APO spectra (those of V841 Oph included).

Additional optical spectra of both novae were obtained using the Hydra Multi-fiber Positioner + Bench Spectrograph (Barden et al., 1993, 1994) on the WIYN 3.5-m telescope222see http://www.noao.edu/wiyn/. Spectra of DI Lac were obtained in 1997 September, and spectra of V841 Oph were obtained in 1997 May. The WIYN spectra have a resolution of $`1`$Å and cover (usable) wavelength ranges of $`6200`$–6800Å for DI Lac and $`4400`$–5000Å for V841 Oph. For the WIYN spectra of DI Lac, no standard star observations were available, so it was not possible to correct for the instrumental response or perform a flux calibration. We do not show these spectra here, but note that they are qualitatively similar to the red APO spectrum shown in Figure 4 – no differences corresponding to the different brightness levels in the RoboScope light curve are apparent. The V841 Oph WIYN data were reduced and calibrated as with the APO spectra. Again, slit losses and, additionally, the presence of cirrus made the absolute flux calibration unreliable, so we have normalized the spectra to a constant continuum level of 1.0. All of the spectra of V841 Oph are shown in Figure 5 (the three WIYN spectra have been averaged together to improve the S/N).

There is no indication of the absorption troughs seen in the spectra of DI Lac. Additionally, the C III/N III + He II emission complex appears to be stronger, with broader lines, in V841 Oph than in DI Lac.

Flux-calibrated optical spectra of both DI Lac and V841 Oph are shown in Williams (1983), with resolutions of 5.1Å and 6.7Å, respectively, and in Ringwald et al. (1996), with resolution of 15–20Å for both. Williams (1983) obtained his spectra in 1981, while Ringwald et al. (1996) obtained theirs in 1990. In both papers, the spectra of both novae show steep blue continua; the observed continuum of DI Lac is bluer than that of V841 Oph (but different reddening values could eliminate or reverse this characteristic). Ringwald et al. (1996) fit a power law of the form $`F_\lambda \lambda ^\alpha `$ to the continuum regions 5600–6500Å and 6650–8000Å in their spectra. These fits yielded indices of $`\alpha =1.81(2)`$ – dereddened to $`\alpha _0=3.10(3)`$ using $`E(BV)=0.41`$ – for DI Lac and $`\alpha =1.45(1)`$ – dereddened to $`\alpha _0=2.79(1)`$ using $`E(BV)=0.58`$ – for V841 Oph. Cassatella et al. (1989) performed this power law fit to the wavelength region 1200–3200Å in $`IUE`$ UV spectra of DI Lac and V841 Oph. Using reddenings at the low ends of the ranges given in §1, $`E(BV)=0.15`$ for DI Lac and $`E(BV)=0.30`$ for V841 Oph, they calculated dereddened slopes of $`\alpha _0=1.5`$ and $`\alpha _0=2.0`$, respectively. Comparison with the slope expected for the continuum flux distribution of a steady-state accretion disk composed of parcels radiating as blackbodies, $`\alpha =2.3`$ (Pringle, 1981), suggests that the true reddening values for DI Lac and V841 Oph are near the middle of the ranges given in §1.

In 1936, both of these novae were observed spectroscopically by Humason (1938), who noted their strong blue continua, but also noted the absence of any lines, emission or absorption. Since that time, the available published spectra of DI Lac and V841 Oph in Greenstein (1960), Williams (1983), and Ringwald et al. (1996) are consistent with our current spectra (see spectroscopic histories summarized in Ringwald et al., 1996), with the exception that Williams (1983) did not report the detection of any He II $`\lambda 4686`$ emission in either star.

### 2.3 X-ray Observations

DI Lac and V841 Oph were each observed ten times during July–September and April–June of 1997, respectively, with the Proportional Counter Array (PCA; Jahoda et al., 1996) on the $`RXTE`$ satellite (e.g. Bradt et al., 1993)333also see http://heasarc.gsfc.nasa.gov/docs/xte/. The PCA consists of five xenon-methane proportional counters effective over the range 2–60 keV (with 18% energy resolution at 6 keV). Each $`RXTE`$ visit constituted a Good-Time-Interval on-source lasting $`1000`$–2500 s. We extracted the net average X-ray count rate from each visit and the X-ray spectrum from the ten combined visits to each star using the FTOOLS/LHEASOFT (v5.0) software package. We utilized the faint source background models (pca\_bkgd\_faintl7\_e3v19990824.mdl and pca\_bkgd\_faint240\_e3v19990909.mdl) and an up-to-date South Atlantic Anomaly history file (pca\_saa\_history\_20000223). Extraction of the count rates and spectra was carried out as described in “The RXTE Cookbook”444 http://heasarc.gsfc.nasa.gov/docs/xte/recipes/cook\_book.html.

The X-ray light curves constructed from the mean count rates at each visit for DI Lac and V841 Oph are shown in Figures 6 and 7, respectively, along with the corresponding sections of their RoboScope $`V`$ band light curves.

The light curve counts were summed in the energy range 0–15 keV. Inspection of the X-ray spectra (see below) revealed an essentially zero count rate dominated by noise at energies higher than 15 keV. The mean X-ray count rates averaged over all ten visits to each target are $`1.30\pm 0.32`$ ct s-1 for DI Lac and $`0.56\pm 0.73`$ ct s-1 for V841 Oph. Considering these very low count rates, it is not clear whether or not the variability seen in the X-ray light curves is real.

Additional analysis and model-fitting (DI Lac only) of the X-ray spectra of these novae was performed using the routine XSPEC (v11). Although we extracted the full range of available energy channels in the X-ray spectra, only the range 0–15 keV was used for fitting the models to the DI Lac data. The X-ray spectra for DI Lac and V841 Oph are shown in Figures 8 and 9.

The X-ray spectrum of V841 Oph is shown on the same scale as that of DI Lac; the former is clearly much weaker than the latter. We attempted to extract a better spectrum of V841 Oph by using only the data from the four visits with the largest count rates (visits 3, 5, 9, and 11 – see Table 2). Unfortunately, this spectrum was virtually identical to the total combined spectrum (i.e. flat and featureless), with the exception of a small upward offset of the mean count rate, from about 0.02 cts s-1 keV-1 in the total combined spectrum to about 0.04 cts s-1 keV-1 in the “high” count rate spectrum. Consequently, no attempt was made to fit a model to the spectrum of V841 Oph.

For DI Lac, we fit three simple models to the X-ray spectrum: (1) blackbody, (2) Raymond-Smith thermal plasma, and (3) bremsstrahlung. Each model was modified by a multiplicative component representing photoelectric (H column) absorption. More complex models were not warranted owing to the small number of counts per energy channel in the data. In all cases, the data were weighted using the XSPEC-recommended scheme appropriate for small count numbers, $`W_\mathrm{i}=1+(N_\mathrm{i}+0.75)^{0.5}`$ (Gehrels, 1986). The results of these model fits are summarized in the first three rows of Table 4. All of the models produce nearly indistinguishable fits to the observed spectrum. The bremsstrahlung model is formally the best-fitting (with lowest reduced $`\chi ^2`$ and highest null hypothesis probability) and is plotted over the observed spectrum in Figure 8.

All of the models were only able to constrain the hydrogen column density to an upper limit; the values of $`N_\mathrm{H}`$ quoted in the first three rows of Table 4 are the upper ends of the 1-parameter 90% confidence intervals. We used the relation

$$N_\mathrm{H}=E(BV)\times 5.8\times 10^{21}\mathrm{cm}^2\mathrm{mag}^1$$

(Bohlin et al., 1978) to estimate values of $`N_\mathrm{H}`$ spanning the reddening range given in §1. We then re-fit the models to the X-ray spectrum of DI Lac with only two free parameters: temperature and normalization. The results are listed in the bottom six rows of Table 4. While slightly different from those obtained from the original model fits, they are completely consistent with them. The bremsstrahlung model yields a slightly lower nominal temperature of $`kT4.1`$–4.2 keV (vs. $`kT=4.5`$ keV when $`N_\mathrm{H}`$ is a free parameter) and is still the best-fitting model for both values of $`N_\mathrm{H}`$.

We note that the uncertainties quoted in Table 4 for $`kT`$ and $`C_{\mathrm{norm}}`$ are 1-parameter 90% confidence intervals; the interdependent parameter uncertainties will be slightly different. For example, for the bremsstrahlung model (with free $`N_\mathrm{H}`$), the XSPEC routine steppar555 This routine iteratively refits the model to the data while stepping the values of selected parameters through a given range and produces a plot in parameter space of the 2-dimensional $`\chi ^2`$ contours. gives a 2-parameter 90% confidence interval of

$$kT=4.5_{1.5}^{+1.9}\mathrm{keV}$$

for $`kT`$ with respect to $`N_\mathrm{H}`$, and

$$kT=4.5_{1.7}^{+1.8}\mathrm{keV}$$

for $`kT`$ with respect to $`C_{\mathrm{norm}}`$.

The X-ray count rates were too low to reliably fit spectrum models to any of the ten individual $`RXTE`$ visits to DI Lac. We constructed “low” and “high” state X-ray spectra by combining the data from visits with the lowest mean count rates (1, 2, 3, 4, 7) and from visits with the highest mean count rates (5, 6, 9) to see if there were any differences between the low and high X-ray brightness states666 Visits 8 and 10 were not used. The former was deemed to be an “intermediate” brightness state, and the latter an anomalously low state.. Other than an $`1`$ cts s-1 keV-1 increase at the peak energy of the high state spectrum, the two spectra are essentially indistinguishable from each other (and from the combined spectrum from all ten visits). We conclude from this that either the spectral energy distribution of the X-ray emitting source in DI Lac does not change as the overall brightness of the X-ray source changes, or that the variability seen in the X-ray light curve of DI Lac is not real (i.e. it is a side-effect of the very low count rates).

## 3 Discussion

### 3.1 Optical Photometric Variability

Inspection of the RoboScope light curves of DI Lac and V841 Oph reveals variability on three distinct time scales. First, in the range 10–100 d, quasiperiods with typical lengths of $`30`$–50 d are present in each of these CVs. Yet, the presence (or lack) of periodic and/or quasiperiodic variability on time scales up to a few hundred days in V841 Oph has been contested in the literature. (To the best of our knowledge, no extensive database of photometric observations of DI Lac has been published prior to this work.) Using 420 archival visual observations of V841 Oph spanning 28 yr, Della Valle & Rosino (1987) quote an average period of 51.5 d with individual cycles ranging from 45 d to 57 d. Shara et al. (1989) observed V841 Oph once per night for 31 nights during a 45 night interval in 1986 using a CCD camera + $`B`$ filter, and found a period of $`40`$ d. All of these values are comparable to the $`35`$–50 d we found for our complete RoboScope light curve of V841 Oph, as well as to the range of periods found for individual observing seasons in the RoboScope data (see Table 3). On the other hand, Duerbeck (1992) notes the presence of considerable variability in archival light curves of V841 Oph (including some of the same data used by Della Valle & Rosino, 1987), but finds no evidence for periodicities up to 200 d.

It is quite interesting that the variability time scales in these two novae are similar to those that characterize dwarf nova outbursts (Warner, 1995, ch. 3). Old novae are presumed to have mass transfer rates, $`\dot{M}`$, above the threshold level ($`\dot{M}_{\mathrm{crit}}`$) at which the accretion disk thermal instability mechanism operates in dwarf novae to produce their outbursts (Osaki, 1996). However, the hibernation theory of cyclic evolution between CV variability types (Shara, 1989, and references therein), predicts that when $`\dot{M}`$ eventually decreases back into the instability regime some decades after the nova outburst, then dwarf nova outbursts should resume. Both GK Persei (= Nova Per 1901; Bianchini et al., 1986) and V446 Herculis (= Nova Her 1960; Honeycutt et al., 1998b) have displayed dwarf nova eruptions in their post-nova stages, and similar behavior has been suggested for a number of other old novae, commencing $``$ 50–200 yr after their outbursts (Livio, 1989; Warner, 1995, ch. 4). The amplitudes of variability displayed by DI Lac and V841 Oph (see Figures 1 and 2) are smaller than those typically observed for dwarf novae ($`\mathrm{\Delta }V`$ 2–4 mag), but we cannot rule out the possibility of a disk instability mechanism being in operation – after all, it has been 90 years since the nova outburst of DI Lac and over 150 years since that of V841 Oph. Further, the classic thermal instability mechanism is not well-explored for high-$`\dot{M}`$ disks, and other disk instabilities might operate as well (e.g. Whitehurst, 1988; Mineshige, 1993; Godon, 1998).

The amplitudes of the $`30`$–50 d oscillations in DI Lac and V841 Oph are similar to those of the spaced stunted outbursts seen in some novalike CVs (Honeycutt et al., 1998a). One of the suggestions for stunted outbursts made in Honeycutt et al. (1998a) is that they are dwarf nova eruptions seen against a background of brighter light in the system. If that suggestion is correct, then the oscillations reported here (as well as similar 25 d oscillations reported for RW Tri, Honeycutt et al., 1994) may be examples of dwarf nova type outburst behavior at relatively short recurrence times, $`T_\mathrm{n}50`$ d ($`T_\mathrm{n}`$ is typically $`10`$ d to hundreds of days – even tens of years in extreme cases – for dwarf novae; Warner, 1995, ch. 3). Among the true dwarf novae, the Z Camelopardalis stars have outburst recurrence times of $`10`$–30 d, somewhat shorter than the time scales for the oscillations observed in these novae. The Z Cam stars are thought to have high mass transfer rates (near $`\dot{M}_{\mathrm{crit}}`$). This suggests that if the oscillations observed in novae are analogous to dwarf nova outbursts, then they may be linked to the presence of relatively high mass transfer rates that are, nonetheless, smaller than both $`\dot{M}_{\mathrm{crit}}`$ and $`\dot{M}`$ in the Z Cam stars.

Second, V841 Oph displays a distinct sinusoidal variation with a best period of 1800–1900 d (4.9–5.2 yr) or a slightly less preferred period of 1250–1300 d (3.4–3.6 yr). DI Lac does not appear to display similar sinusoidal variability, unless it occurs with a much longer period than in V841 Oph. Bianchini (1990) reports a long-term period in V841 Oph of 3.4 yr with mean amplitude of $`0.3`$ mag, which was determined from 420 visual observations spanning 30 yr (essentially the same early- to mid-20th century data utilized by Duerbeck, 1992, see below). Richman et al. (1994) re-analyzed the light curves presented by Bianchini (1990) and consider the case for a multi-year periodicity in V841 Oph to be very weak. However, with our RoboScope light curve, sinusoidal variability on comparable multi-year times scales has now been observed in two separate data sets for V841 Oph spanning almost 80 years of observation. While this still does not provide firm evidence of strictly periodic behavior, it certainly points to the presence of a mechanism operating in this CV that modulates its brightness on a characteristic time scale of several years.

Photometric variability in old novae on time scales of several years has generally been attributed to solar-type magnetic cycles on the secondary star that might be able to control the rate of mass loss through the inner Lagrangian point (Bianchini, 1990; Applegate, 1992). Although Richman et al. (1994) find that no CV studied to date displays strictly periodic behavior over multi-year intervals, they also note that the observed amplitudes ($`0.2`$ mag) and “apparent” time scales of variability of 5–40 yr are “plausible consequences from solar-type magnetic cycles.” We note that the shorter time scale ($``$30–50 d) variability discussed above might also be linked to modulation of $`\dot{M}`$ due to secondary star magnetic activity. If solar-type magnetic cycles can affect $`\dot{M}`$ on time scales of years, then starspots induced by the magnetic cycle that migrate under the $`L_1`$ point might also affect $`\dot{M}`$ on shorter time scales, in a manner suggested by Livio & Pringle (1994) to explain the very low brightness states seen in some CVs.

Third, the mean magnitudes of both CVs are changing slowly, at rates measured in a few millimagnitudes (mmag) per year. The mean brightness of DI Lac is increasing by 4 mmag yr-1 since 1990. Duerbeck (1992) analyzed a large set of visual, photoelectric, and CCD observations of post-novae available in the literature (see references in Duerbeck, 1992). These data cover a large part of the 20th century, from the 1920’s through the 1980’s. Duerbeck (1992) summarizes the long-term behavior of DI Lac as exhibiting a steady decline of $`13\pm 1`$ mmag yr-1 in visual observations obtained between 1921 and 1952. Sparse photoelectric and CCD data from 1953 to 1981 suggest a possible small brightness increase, which is supported by our recent RoboScope observations. For V841 Oph, visual observations from 1919–52 and 1978–91 yield a brightness decline of $`7\pm 1`$ mmag yr-1, while “no definitive brightness decline” (Duerbeck, 1992) was found in photoelectric data from 1954 to 1988. Our post-1991 observations of V841 Oph suggest a trend towards decreasing mean brightness comparable to that noted by Duerbeck (1992) in the archival visual observations, although the obvious sinusoidal variation that dominates our light curve at long time scales makes this conclusion somewhat suspect.

In the “classical” hibernation scenario (Shara, 1989), novae stay bright for $``$ 50–300 yr following outburst, gradually declining in brightness during this time to a hibernating state characterized by low brightness ($`M_V10`$ or fainter) and low mass accretion rate ($`\dot{M}10^{12}`$ M yr-1). However, the photometric record of DI Lac from the 1920’s (a decade after its outburst) to the 1990’s does not show a simple decline in brightness. Instead, the brightness of DI Lac decreased for about four decades after outburst, but then apparently leveled off for several decades. In the last decade (covered in our RoboScope light curve), the brightness of DI Lac has been increasing. It is possible that we are seeing the influence of a mechanism (such as magnetic activity of the secondary star – see above – and/or accretion-induced irradiation of the secondary star) that modulates $`\dot{M}`$ (and, hence, the brightness) in DI Lac over a long-term cycle (with period on the order of many decades) that obscures the general decline predicted by the hibernation scenario.

### 3.2 Optical Spectrum Variability

The optical spectra of DI Lac and V841 Oph are quite different, and this can possibly be ascribed to their difference in post-outburst ages. The spectrum of V841 Oph is similar to those of novalike CVs (e.g. Warner, 1995, ch. 4) – this suggests that V841 Oph still has a high accretion rate, even $`150`$ yr after outburst. The narrow emission components in DI Lac are suggestive of emission originating from the irradiated inner face of the secondary star, while the absorption troughs imply the presence of optically thick material in the system. The latter feature could be material ejected during the outburst of this younger post-nova; the lack of any detected H$`\alpha `$ emission shell (Cohen, 1985) does not preclude the existence of denser, non-emitting circumstellar material. The former feature offers a possible explanation for the gradual brightness increase seen in our RoboScope light curve of DI Lac (discussed in §3.1) if the secondary star is being slowly heated via irradiation and is, in turn, increasing the rate of mass transfer through the $`L_1`$ point.

Although the RoboScope coverage is incomplete, we infer from the adjacent data that DI Lac was faint during our July spectra and bright during our August spectrum (see Figure 6). As mentioned in §2.2, the only difference between our July and August spectra is that the emission cores were somewhat stronger in August (bright state) than July (faint state). This is consistent with the hypothesis that these narrow emission features originate on the irradiated face of the secondary star if we make the logical assumption that the irradiation increases when DI Lac is bright. It does not, however, illuminate the exact mechanism producing the irradiation (i.e. whether the “excess” flux in the high state originates near the disk center – presumably due to a disk instability producing increased accretion onto the WD – or in the outer disk – presumably due to an increase in mass transfer through the $`L_1`$ point.) If we assume that the spectrum shown in Figure 4 (from July 27 09:01:47 UT) was obtained at an arbitrary orbital phase of 0.0, then the August 17 spectrum was obtained at a relative orbital phase of 0.9 (using $`P_{\mathrm{orb}}=0.543773`$ d; Ritter & Kolb, 1998). Thus, these two spectra were obtained at similar orbital phases, and we do not expect the difference in the narrow emission component strength to be only due to system orientation.

The RoboScope coverage of V841 Oph during our spectroscopic (and X-ray) observations is more sparse than for DI Lac, but we can infer that the CV was returning to the faint state during our May 07 spectrum, was in its faint state during our May 14 and 17 spectra, and was likely near the bright state during our May 27 spectrum (see Figure 7). (The two elevated brightness points at HJD 2450610 and HJD 2450625 suggest that V841 Oph may have returned to a bright state sometime between the low states bracketing HJD 2450590 and HJD 2450630.) If this is the case, then the somewhat stronger H$`\alpha `$ emission of the May 07 and May 27 spectra (see Figure 5) may be linked to the bright state of the CV. As with DI Lac, we calculated relative orbital phases for each of these spectra (using $`P_{\mathrm{orb}}=0.60423`$ d; Ritter & Kolb, 1998), and obtained $`\varphi =0.0,0.6,0.75,0.1`$ for May 07, 14, 17, 27, respectively. Unfortunately, this casts some doubt on the link between H$`\alpha `$ emission strength and brightness state, since both of the bright state spectra (phases 0.0, 0.1) were obtained at different orbital phases than the faint state spectra (phases 0.6, 0.75). So, we cannot rule out the influence of system orientation effects in these spectra.

### 3.3 X-ray Variability

Unfortunately, because of the weak X-ray emission from these CVs, little can be firmly stated about their X-ray variability. The count rates during visits 5 and 6 to DI Lac, which took place when we infer from the RoboScope light curve that the CV was in a bright state, are slightly elevated compared to the preceding visits. However, visit 7 also occurred during this optical bright state and does not show an elevated count rate. The mean X-ray count rate during each visit to V841 Oph also does not display any strong correlation with the corresponding optical state. The lack of large changes in X-ray flux or spectrum during the optical variations (as usually evident in dwarf novae; e.g. Szkody et al., 1999), argues against a disk instability scenario. Alternatively, this could indicate that the optical brightness is determined by activity in the outer disk only and, therefore, is not reflected in the X-ray behavior; however, the X-ray count rates are too low (and their corresponding error bars too large) to make any firm conclusions.

## 4 Conclusions

Our long-term optical light curves of the novae DI Lac and V841 Oph obtained with RoboScope reveal quasiperiodic variability with a characteristic time scale of $`30`$–50 d in both CVs. In addition, the light curve of V841 Oph displays evidence for sinusoidal variability with a period of 3.5–5 yr. The latter cannot be said to be strictly periodic since our data set covers only $`1.5`$ cycles; however, when this detection is added to past reports of multi-year periodicities in V841 Oph that have been reported in the literature, a strong case can be made for the presence of repeating multi-year variability with a preferred time scale in this CV. The most likely origin of such behavior is a solar-type magnetic cycle of the secondary star that modulates the rate of mass transfer through the $`L_1`$ point. If this is the case, then the shorter $`30`$–50 d quasiperiodic variability in this system might also be related to the magnetic activity on the secondary star; for example, due to starspots that migrate under the $`L_1`$ point and temporarily throttle mass transfer. DI Lac does not show evidence for cyclic multi-year variability, unless it occurs with a much longer period than the length of our RoboScope coverage ($`10`$ yr). This casts some doubt on the origin of the $`30`$–50 d quasiperiodic variability as a facet of the secondary star’s magnetic activity, since both V841 Oph and DI Lac show the $`30`$–50 d quasiperiodic behavior, but only V841 Oph displays evidence for a multi-year solar-type magnetic cycle on the secondary star. On the other hand, the $`30`$–50 d time scale is very reminiscent of that expected for dwarf nova type behavior (although the amplitude of variability in these novae is smaller than in typical dwarf nova outbursts). This raises the possibility that this variability is caused by a disk instability (either the thermal disk instability that leads to dwarf nova outbursts operating at a low level or some other form of disk instability).

The X-ray spectrum of DI Lac is fit almost equally well by the three models we tried: a simple blackbody, a Raymond-Smith thermal plasma, and bremsstrahlung emission. More complicated models (e.g. involving multiple components, lines, etc.) are unwarranted due to the low X-ray flux. In addition to being the most physically plausible X-ray emission emission mechanism in a (non-magnetic) CV, the bremsstrahlung model is formally the best-fitting. Our X-ray spectrum model fit parameters ($`kT4`$ keV) are consistent with those obtained from $`ROSAT`$ X-ray spectra of a sample of 37 disk-accreting CVs (Richman, 1996). Unfortunately, the X-ray count rates in both of these systems were too low for any conclusions to be made about their time-resolved X-ray behavior. X-ray observations sample the innermost disk region, and we would expect to see different time-resolved behavior as the novae go into their optically bright state if the transition is triggered by a disk instability vs. a change in $`\dot{M}`$ from the secondary star. Additional time-resolved X-ray observations using the more sensitive Chandra and/or XMM X-ray satellites may be necessary to illuminate the inner workings of these old novae.

PS, DWH, and VD acknowledge support from NASA grant NAG5-4791. DWH thanks the NOAO librarian Mary Guerrieri for her valuable assistance locating several papers cited herein. This research made use of NASA’s Astrophysics Data System Abstract Service and the SIMBAD database operated by CDS, Strasbourg, France.


format: text

text:
# Control and Tracking of Unstable Orbits in Hamiltonian Flows: The Diamagnetic Kepler Problem

## Abstract

We report on the application of chaos control to the irregular motion of an electron under the combined influence of a Coulomb and a magnetic field, the so-called “diamagnetic Kepler problem” (DKP). We show how to stabilize the classical chaotic orbits into regular motion and present a new method to follow the unstable orbits as the energy is increased from the regular to the chaotic regime.

The dynamics of an hydrogen atom in a strong magnetic field is known as the “diamagnetic Kepler problem” (DKP). After more than two decades of research, the DKP is currently the most intensively studied autonomous, classically chaotic, atomic physics system. As first realized by Edmonds , the DKP is important in such diverse areas as solid state physics, astrophysics, and Rydberg atoms and occupies central stage in classical and quantum chaos research . The classical flow covers a wide range of Hamiltonian dynamics reaching from bound, nearly integrable behavior to completely chaotic and unbound motion as the scaled energy is varied . This letter reports the numerical implementation of a control scheme to stabilize this chaotic behavior. The methodology is sufficiently general to apply to all Hamiltonian flows and , with slight modifications, is also applicable to ballistic dynamics (e.g. billiards) and area-preserving mappings.

Ever since Ott, Grebogi and York (OGY)  introduced the notion of chaos control a decade ago, the number of experimental verifications has continuously and rapidly expanded to include a wide variety of applications ranging from magnetoelastic ribbons , electronic circuits  , lasers  , chemical reactions  to thermal convection loops . The OGY algorithm was written explicitly for dissipative systems and most studies (experimental and theoretical) have focused their attention on this class of systems. Non-trivial modifications are needed to adapt the OGY technique to conservative flows (or maps) and Lai, Ding and Grebogi  have provided the necessary changes in the context of area-preserving mappings.

Our system is a continuous, 2 degrees of freedom (4D phase space) Hamiltonian flow representing the motion of an electron under the combined influence of a Coulomb and an external magnetic field. We use scaled semi-parabolic coordinates and write the resulting scaled (pseudo-) Hamiltonian as (for angular momentum $`L=0`$)

$$\widehat{h}_{DK}=\frac{1}{2}(p_\nu ^2+p_\mu ^2)ϵ(\nu ^2+\mu ^2)+\frac{1}{8}\nu ^2\mu ^2(\nu ^2+\mu ^2)2.$$
(1)
The scaled energy $`ϵ`$ is related to the physical energy $`E`$ by $`ϵ=\gamma _0^{2/3}E`$ where the parameter $`\gamma _0=B/B_0`$ denotes the strength of the magnetic field relative to the unit $`B_0\mathrm{2.35\; 10}^5T`$. The structure of the dynamics of the Hamiltonian (1) depends solely on the value of the scaled energy $`ϵ`$. For instance, it is known that for $`0.5<ϵ<0.13`$ the system exhibits bounded motion with mixed chaotic and regular motion and for $`0.13<ϵ<0.0`$ the last large stable island disappears and the dynamics is mostly chaotic .

The dimension reduction (from 4D to 2D) and discretization is performed by observing the dynamics on the Poincaré section defined by $`\mu =0,\dot{\mu }>0`$. The energy shell is then mapped to an area bounded by the condition $`p_\nu ^22ϵ\nu ^2=4`$ which represents an ellipse in the $`(\nu ,p_\nu )`$ plane. Figure (1) shows the collection of points $`\{\nu _n,p_{\nu ,n}\}`$ obtained by numerical integration of the equations of motion for $`ϵ=0.5`$. One notices, for this energy, that the mixed phase space is still mostly regular with a small stochastic area filled by the successive piercings of one chaotic trajectory. We recall that a periodic trajectory would show up on the Poincaré section as a finite number of crossings.

In attempting to bring order to the DKP dynamics, we had to overcome a number of difficulties not encountered in previous studies. We had in mind the adaptive control of chaotic DKP trajectories under changing conditions (e.g. drift in the magnetic field) through the stabilization of preselected unstable periodic orbits (UPOs) (monitored on the Poincaré section) via small programmed perturbations to the scaled energy. The difficulties are five-fold. First, a typical trajectory spends a lot of time away from the Poincaré section where the flow is discretized and because of the sensitivity of the chaotic dynamics we had to device an efficient variable step symplectic integrator thereby preserving the geometrical structure of the Hamiltonian . Second, the locations of a number of UPOs (our target states) and the corresponding Jacobian matrices for each member of the periodic orbits must be obtained numerically. It was found necessary to intervene at every crossing of the Poincaré section to assure a robust stabilization. Third, because of anomalous transport in Hamiltonian phase space , the waiting time for a trajectory to enter the neighborhood of a predefined target UPO is unduly long and a targeting strategy is in general necessary to steer the trajectory towards the appropriate region of phase space. This problem has been addressed in  and is assumed to be solved. Fourth, the eigenvalues of area-preserving Jacobians are often complex and the stable and unstable manifolds, at the core of the OGY procedure, are no longer along the directions of their eigenvectors. A new method had to be implemented along the lines described in . Fifth, a tracking algorithm for updating the control parameters (positions of UPOs, Jacobian matrices …) was designed to insure that control is kept while external conditions are changing. Part of our solutions to these problems are discussed below and the details will be reported elsewhere.

To locate the UPOs on the Poincaré section, we have used two methods: a modified recurrence method (MRM) and an extension to conservative mappings of a stability transform algorithm (STA) introduced in . The MRM uses the basic idea proposed by Auerbach et al.  for extracting periodic orbits from a chaotic time series $`\{𝐗_i\}_{i=1,N}`$. However, instead of scanning a very long $`N1`$ time series for pairs of points separated in time by $`m`$ iterates (for a period-m UPO) and within a preassigned spatial distance $`r`$ from each other, $`M`$ points are launched on the Poincaré section and integrated forward until the $`m`$-th return on the section where the pair of points are compared. In both methods, the accepted pairs, a distance $`r`$ apart, are distributed among the possibly different UPO candidates, the center of mass of each subset is calculated and the coordinates are stored. For a given relative accuracy, $`10^410^6`$, we find $`MN`$ by several orders of magnitude. In Hamiltonian cases, $`N`$ must be at least $`10^6`$ in order to reach an acceptable precision. In the MRM, one can also easily take advantage of the symmetries on the Poincaré section for the flow. For example, figure (1) shows that the mapping is symmetric with respect to both $`\nu `$ and $`p_\nu `$ axes. Therefore one only needs to spread the initial conditions over one quarter of phase space for a smaller set of initial conditions or a denser distribution. This optimization is essential if one considers that between every two consecutive points on the Poincaré section, a segment of trajectory in phase space must be integrated numerically. Furthermore, for the DKP, this segment grows rapidly as the scaled energy approaches zero where the trajectory spends an increasingly long time in the branches (they grow as $`|ϵ|^{1/2}`$) of the pseudo-potential.

The STA amounts to the following. Consider a discrete chaotic dynamical system $`S`$ defined in a $`d`$-dimensional space by

$$𝐗_{n+1}=𝐅(𝐗_n,p)$$
(2)
where $`p`$ is some external parameter. The goal is to construct from equation (2) other dynamical systems $`S_k`$ with the same periodic orbits (in number and positions), but whose stability has changed, i.e unstable orbits become stable, and stable orbits remain stable. The new systems are defined by the linear transformations (written for a period $`m`$ point, $`𝐅^{(m)}(𝐗^{},p)=𝐗^{}`$)

$$𝐗_{n+1}=𝐗_n+\lambda 𝐂_k[𝐅^{(m)}(𝐗_n,p)𝐗_n]$$
(3)
where $`0<\lambda 1`$ is adjusted to improve convergence and $`𝐂_k`$ are non-singular constant $`d\times d`$ matrices (see  for their explicit forms). Because of the stability of the periodic orbits of the constructed system $`S_k`$, every trajectory of $`S_k`$ converges after a finite number of iterations to a periodic $`m`$ point $`𝐗^{}`$. Per construction, the $`𝐗^{}`$ are also periodic points of the original system $`S`$. We have tested this procedure for a number of area preserving mappings and flows and the global convergence of the method is excellent. For flows, where the mapping in (3) is replaced by a numerical integration for $`m`$ returns to the Poincaré section, the STA has proved equally reliable. In this respect, our adaptive symplectic integrator  has insured that the (long) trajectories away from the Poincaré section were kept sufficiently accurate. The position of the fixed point in figure (2) (bottom panel) has been reached after only 250 iterations starting from an arbitrary initial point: this is representative of the efficiency of the STA. For the same relative accuracy ($`10^6`$), the MRM needs approximately $`10^4`$ initial points.

In order to stabilize a chaotic trajectory around one of these UPOs, we have implemented a numerical version of the OGY method as modified in  for area-preserving mappings. This method is believed to be dynamically optimal in that it explicitly uses the local geometry of the underlying system. Given a dynamical system of the type (2), and a target UPO of period $`m`$, $`\{𝐗(i,p_0)\}_{i=1,m}`$ , at some nominal parameter value $`p_0`$, one characterizes the local stable and unstable manifolds by the vectors $`𝐞_{s,i}`$ and $`𝐞_{u,i}`$ respectively as well as their contravariant counterparts $`𝐟_{s,i}`$ and $`𝐟_{u,i}`$ satisfying $`𝐟_{u,i}𝐞_{u,i}=𝐟_{s,i}𝐞_{s,i}=1`$ and $`𝐟_{u,i}𝐞_{s,i}=𝐟_{s,i}𝐞_{u,i}=0`$. The stabilizing perturbations $`\delta p_np_np_0`$ are then obtained by firstly linearizing the dynamics in a $`\delta `$-neighborhood of a member of the periodic orbit, say $`𝐗(k,p_0)`$ , and around $`p_0`$, namely

$$𝐗_{n+1}𝐗(k+1,p_n)𝐔_k\left[𝐗_n𝐗(k,p_n)\right]$$
(4)
where the $`d\times d`$ Jacobian matrix $`𝐔D_𝐗𝐅(𝐗,p)`$ is evaluated at $`[𝐗=𝐗(k,p_0),p=p_0]`$ and it is understood that $`𝐗_n𝐗(k,p_0)\delta 1`$ and $`|\delta p_n||p_0|`$.

Secondly, the control criterion is imposed that $`𝐗_{n+1}=𝐅(𝐗_n,p_n)`$ should lie along the stable direction at $`𝐗(k+1,p_0)`$, i.e. $`𝐟_{u,k+1}[𝐗_{n+1}𝐗(k+1,p_0)]=0`$ which, together with the parametric variation of the periodic points, $`𝐗(k,p_0+\delta p)𝐗(k,p_0)+𝐠_k\delta p`$ , leads to the following expression for the parameter perturbation at the $`n`$-th iteration:

$$\delta p_n=\frac{𝐟_{u,k+1}\left\{𝐔_k\left[𝐗_n𝐗(k,p_0)\right]\right\}}{𝐟_{u,k+1}(𝐠_{k+1}𝐔_k𝐠_k)}.$$
(5)
The individual Jacobian matrices $`𝐔_i`$ were obtained by least-square minimization  of the linear relationship between subsequent pairs of iterates on the Poincaré section for over 1000 points drawn from the close neighborhood of each $`𝐗(i,p_0)`$. For the dissipative cases, the eigenvectors of the Jacobian matrix $`𝐔_k^{(m)}=D_𝐗𝐅^{(m)}(𝐗_k,p_0)=_{j=0}^{m1}𝐔_{k+j}`$ would determine the manifold directions, but as first realized in , the conservative cases lead generically to complex eigenvalues and eigenvectors and a different strategy is required. The solution adopted to generate the unstable (stable) direction $`𝐞_u`$ ($`𝐞_s`$) is to rotate an arbitrary unit vector $`𝐮_0`$ ($`𝐬_0`$) according to $`𝐔_k^{(qm)}𝐮_0`$ and $`𝐔_k^{(qm)}𝐬_0`$ which converge respectively to $`𝐞_{u,k}`$ and $`𝐞_{s,k}`$ for $`q1`$. The vectors are periodically renormalized to prevent overflow and in general, $`q10`$ is sufficient for convergence. This procedure is very reliable and robust.

Putting all the ingredients together, we have achieved the stabilization of the chaotic behavior of the diamagnetic-Kepler Hamiltonian (1) for different values of $`ϵ`$ and different UPOs. Figure (2) presents a sample of results for the control of period-2, period-3 and period-1 at $`ϵ`$ = -0.3, -0.2 and -0.1 respectively. The control mechanism was switched on when an iterate entered a neighborhood of size $`\delta =10^3`$ ( in the 3 panels, one has removed the initial transient) and the perturbations on the scaled energy (alternatively the magnetic field) were never larger than $`\delta ϵ=10^2`$. After 3000 iterates, we release the control, $`\delta ϵ_n=0`$ , and the trajectory returns to its natural chaotic behavior. One also notes (left panels) that the dynamics is increasingly chaotic as $`ϵ`$ approaches zero and by $`ϵ=0.1`$ no obvious trace of regularity is present. We emphasize that in performing the present control task no a priori knowledge of the governing map on the Poincaré section is known and the complete procedure is purely numerical from beginning to end.

Under parameter variation, the location of the UPOs changes as well as their associated Jacobian matrices, and unless action is taken the control over the chaotic dynamics may be lost. The strategy of detecting parameter drift and updating accordingly the control parameters is known as tracking . We have developed a new tracking algorithm which is simple and efficient. For concreteness, we will concentrate the presentation on the updating localization of the UPOs.

The method is based on a judicious construction of a family of extrapolating rational functions . One starts from a number $`s`$ of UPO coordinates (the seeds), e.g. $`\{\nu _i,ϵ_i\}_{i=1,s}`$, from which a rational function, $`R_0(s,ϵ)`$ is constructed. $`R_0`$ is then evaluate at $`ϵ_{s+1}=ϵ_s+\mathrm{\Delta }ϵ`$ with $`\mathrm{\Delta }ϵ1`$ to give an extrapolated value $`\nu _{s+1}`$. A new set of seeds is formed $`\{\nu _{i+1},ϵ_{i+1}\}_{i=1,s}`$ , removing $`\{\nu _1,ϵ_1\}`$ and keeping the number of seeds fixed , and a new extrapolating function $`R_1(s,ϵ)`$ is calculated. The process is repeated until a final value $`ϵ_f`$ is reached. The same procedure may also be applied to the pairs $`\{p_{\nu ,i},ϵ_i\}_{i=1,s}`$ to provide initial conditions $`\{\nu _i,p_{\nu ,i}\}`$ for the integration of the complete phase-space trajectory. Two remarks are in order. First, this approach is much more accurate at $`ϵ_f`$ than the direct extrapolation $`R_0(s,ϵ_f)`$, although all the original seeds are replaced by extrapolated values after $`s`$ iterations. Second, the algorithm permits to cover a wide range of parameter variations beyond the initial seeds (3-5 in practice) while keeping high accuracy. Since the individual extrapolation steps are chosen small, we call the method adiabatic localization (AL) or more generally adiabatic tracking (AT).

We have applied AT to the dissipative Hénon map for the same range chosen in . The results are in complete agreement with those of  and with the precise values of the position of UPOs obtained by STA. We have also applied this algorithm to the standard map, the Hénon-Heiles Hamiltonian and the DKP and more complete results will be reported at a later time. Figure (3) shows the tracking of a period-1 orbit for the DKP. We started from five seeds for $`ϵ`$: - 0.3, - 0.29, - 0.28, -0.27, -0.26 and set $`\mathrm{\Delta }ϵ=0.001`$. The tracking process is extended to $`ϵ=0.075`$. The results from the STA confirms the accuracy of the adiabatic tracking. The right panel clearly indicates that the result from AT for $`ϵ=0.075`$ is indeed a period 1 orbit.

We have presented for the first time the application of a numerical implementation of the OGY control to a realistic conservative flow, the DKP Hamiltonian. On our way to a complete solution, we have i. developed a new variable step symplectic integrator, ii. modified the standard recurrence method and extended a powerful technique (STA) to locate the positions of the unstable periodic orbits, iii. obtained accurate numerical Jacobian matrices and the local manifolds, and iv. provided an efficient approach (AT) for tracking the UPOs under parameter changes. We are presently considering two important mesoscopic applications: nonlinear dynamics in electromagnetic traps , and chaotic ray dynamics in lasing microcavities . It still remains an open question however if manipulations of an external parameter (e.g. a magnetic field) to induce stabilization of a classical unstable orbit can be extended to the semi-classical regime, for example, for the control of Rydberg wave packet dynamics. Research along these lines is actively being pursued.


text:
# Conservation laws for the nonlinear Schrödinger equation in Miwa variables.

## 1 Introduction.

In the present note we would like to discuss once more the conservation laws for the nonlinear Schrödinger equation (NLSE),

$`i{\displaystyle \frac{q}{t}}+{\displaystyle \frac{^2q}{x^2}}+2q^2r`$ $`=`$ $`0`$ (1)
$`i{\displaystyle \frac{r}{t}}+{\displaystyle \frac{^2r}{x^2}}+2qr^2`$ $`=`$ $`0`$ (2)
Existence of an infinite number of conserved quantities is an characteristic feature of integrable partial differential equations (PDE), and this question in context of the NLSE was discussed in the very first works devoted to this model, where authors established integrability of the NLSE and elaborated the corresponding inverse scattering transform (IST). The IST is based on a representation of the equation in question as a compatibility condition for a overdetermined linear system (the so-called zero curvature representation) which in the case of the NLSE can be written as

$`{\displaystyle \frac{}{x}}\mathrm{\Psi }(x,t;\lambda )`$ $`=`$ $`U(x,t;\lambda )\mathrm{\Psi }(x,t;\lambda )`$ (3)
$`{\displaystyle \frac{}{t}}\mathrm{\Psi }(x,t;\lambda )`$ $`=`$ $`V(x,t;\lambda )\mathrm{\Psi }(x,t;\lambda )`$ (4)
Here $`\mathrm{\Psi }`$ is a $`2`$-column, the matrix $`U`$ is given by

$$U=i\left(\begin{array}{cc}\lambda & r\\ q& \lambda \end{array}\right)$$
(5)
and $`V`$ is some $`2\times 2`$ matrix which is a second order polynomial in $`\lambda `$ (in what follows we will not need its explicit form).

The $`U`$-$`V`$ representation (3), (4) of the NLSE is enough to obtain answers to a wide range of questions related to this equation. In particular, one can derive from (3), (4) an infinite number of integrals of motion of the NLSE. This can be done as follows (here we will only outline some key moments, more elaborated description of the IST one can find in various textbooks on this topic, as, e.g., ). Introducing the so-called Jost functions $`\mathrm{\Psi }_\pm `$ of the scattering problem (3) (i.e. solutions of (3) satisfying different boundary conditions) and the scattering matrix $`T(x,t;\lambda )`$, $`\mathrm{\Psi }_+=\mathrm{\Psi }_{}T`$, one can obtain from the zero-curvature representation that the diagonal elements of the $`2\times 2`$ matrix $`T`$ do not depend on time. Hence they (or their logarithms) can be used as generating functions of constants of motion (the later are the coefficients of the power series in $`\lambda `$ ($`\lambda ^1`$) of the former).

Another approach stems from the viewpoint when one considers an integrable equation (the NLSE in our case) as a member of an integrable hierarchy (the AKNS in our case). The $`U`$-$`V`$ representation of the equations of the AKNS hierarchy can be written as

$`{\displaystyle \frac{}{x}}\mathrm{\Psi }`$ $`=`$ $`U\mathrm{\Psi }`$ (6)
$`{\displaystyle \frac{}{t_n}}\mathrm{\Psi }`$ $`=`$ $`V_n\mathrm{\Psi }`$ (7)
where $`U`$ is given again by (5) and $`V_n`$’s are different matrices ($`V_n`$ is some $`n`$th order polynomial in $`\lambda `$). All equations of the hierarchy are compatible and one can solve them simultaneously, i.e. one can think of $`q`$ and $`r`$ as functions of an infinite set of variables $`q,r=q,r(t_1,t_2,t_3,\mathrm{})`$ with the evolution with respect to $`t_k`$ being given by the $`k`$th NLSE ($`k`$th member of the AKNS hierarchy). In some situations such standpoint leads to more transparent results, and the main aim of this note is to apply it to the question of the description of the constants of motion of the NLSE (read constants of motion of the AKNS hierarchy).

## 2 Functional representation of the AKNS hierarchy.

Our starting point is the so-called functional representation of the AKNS hierarchy which can be written as

$`i\zeta _1q`$ $`=`$ $`qq^{}+\zeta ^2q^2r^{}`$ (8)
$`i\zeta _1r`$ $`=`$ $`rr^++\zeta ^2q^+r^2`$ (9)
Here $`q`$ and $`r`$ are functions of an infinite set of times, $`q=q(\mathrm{t})`$, $`r=r(\mathrm{t})`$,

$$f(\mathrm{t})=f(t_1,t_2,t_3,\mathrm{})$$
(10)
$`_1=/t_1`$ and the designation $`f^\pm `$ stands for the function with shifted arguments (Miwa shifts),

$`f^\pm `$ $`=`$ $`f(\mathrm{t}\pm i[\zeta ])`$ (11)
$`=`$ $`f(t_1\pm i\zeta ,t_2\pm i\zeta ^2/2,t_3\pm i\zeta ^3/3,\mathrm{})`$ (12)
This equation, which can be termed as ’AKNS hierarchy in Miwa variables’, may be derived in different ways. First this can be done by careful analysis of the linear problems (6), (7). One can find the shifts $`t_kt_k\pm i\zeta ^k/k`$ in some textbooks on the IST (see, e.g., chapters 3,4 of the book ). Another way is to use the (generalized) Hirota’s bilinear identities which are one of the most important formulae of the Kyoto school approach to the integrable systems. Explicitly equations (8) and (9) have been written down in the paper . One should also mention paper  where similar representation has been derived for the Davey-Stewartson system, from which one can easily obtain one for the AKNS hierarchy. We will not repeat here derivation of (8) and (9) and only demonstrate that the first equations of the AKNS hierarchy (the NLSE in particular) can be easily obtained from them. Indeed, using the multidimensional Taylor series for $`f^\pm =f(\mathrm{t}\pm i[\zeta ])`$,

$$f^\pm =f\pm i\zeta _1f+\frac{\zeta ^2}{2}\left(\pm i_2f_{11}f\right)+\frac{\zeta ^3}{6}\left(\pm \mathrm{\hspace{0.17em}2}i_3f3_{21}fi_{111}f\right)+\mathrm{}$$
(13)
(here $`_k`$ stands for $`/t_k`$, $`_{jk}`$ for $`^2/t_jt_k`$, etc) and expanding equations (8), (9) in power series in $`\zeta `$ one will obtain that the functions $`q`$ and $`r`$ satisfy an infinite number of PDEs. The first non-trivial equations (the $`\zeta ^2`$ terms),

$`i_2q+_{11}q+2q^2r`$ $`=`$ $`0`$ (14)
$`i_2r+_{11}r+2qr^2`$ $`=`$ $`0`$ (15)
are nothing else than the NLSE. Hereafter we will identify variables $`x`$, $`t`$ and $`t_1`$, $`t_2`$

$$x=t_1,t=t_2$$
(16)
The next equations

$`2_3q3i_{21}q_{111}q6q^2_1r`$ $`=`$ $`0`$ (17)
$`2_3r+3i_{21}r_{111}r6r^2_1q`$ $`=`$ $`0`$ (18)
can be rewritten using (14), (15) as

$`_3q+_{111}q+6qr_1q`$ $`=`$ $`0`$ (19)
$`_3r+_{111}r+6qr_1r`$ $`=`$ $`0`$ (20)
These are the third-order NLSE. Thus one can view (8), (9) as a ’condensed’ form of the AKNS hierarchy.

The key moment is that, if we deal not only with solutions of the NLSE, $`q,r(x,t)`$ but consider them as solutions of all equations of the hierarchy, $`q=q(t_1,t_2,\mathrm{})`$, $`r=r(t_1,t_2,\mathrm{})`$, then we can formally solve the auxiliary linear problem. Indeed, it follows from (8), (9) that the $`2\times 2`$ matrix

$$\mathrm{\Psi }=\left(\begin{array}{cc}1& \zeta r^{}\\ \zeta q^+& 1\end{array}\right)\left(\begin{array}{cc}\mathrm{exp}(iu_1)& 0\\ 0& \mathrm{exp}(iu_2)\end{array}\right)$$
(21)
where

$`u_1`$ $`=`$ $`{\displaystyle \frac{x}{2\zeta }}+\zeta {\displaystyle dxq^+r}`$ (22)
$`u_2`$ $`=`$ $`{\displaystyle \frac{x}{2\zeta }}+\zeta {\displaystyle dxqr^{}}`$ (23)
solves (6) with $`\lambda =(2\zeta )^1`$. Hence we can now rewrite the results which were presented in terms of solutions of (6) (i.e. in terms of the Jost or Baker-Akhiezer functions) in terms of $`q`$, $`r`$ themselves. This is also valid for the generating function of the integrals of motion.

Of course, matrix (21) is a formal solution of (6) and one must be ready to face some problems when, e.g., one will try to construct the Jost functions (i.e. to satisfy some boundary conditions). However, for our purposes this is not an obstacle. Moreover, we will not repeat the ’classical’ algorithm, Jost functions $`\mathrm{\Psi }_\pm (\lambda )`$ $``$ scattering matrix $`T(\lambda )`$ $``$ generating function $`\mathrm{ln}T_{11}(\lambda )`$. Knowing the answer, we will first present the final result, and then, using the functional representation (8), (9) of the AKNS hierarchy, will prove it.

## 3 Conservation laws.

The main result of this work can be presented as follows: the function

$$J(\mathrm{t},\zeta )=q(\mathrm{t}+i[\zeta ])r(\mathrm{t})$$
(24)
is the generating function for the constants of motion.

Indeed, it follows from (14), (15) that

$$i\frac{}{t}q^+r=\frac{}{x}\left(q^+\frac{r}{x}\frac{q^+}{x}r\right)2q^+r\left(q^+r^+qr\right)$$
(25)
Using again (14), (15), this time with shifted arguments, one can easily get

$$i\zeta \frac{}{x}q^+r=q^+r^+qr$$
(26)
which leads to

$$\frac{}{t}J(\zeta )=\frac{}{x}F(\zeta )$$
(27)
where

$$F=i\frac{q^+}{x}rq^+\frac{r}{x}\zeta \left(q^+r\right)^2$$
(28)
(recall that $`x=t_1`$ and $`t=t_2`$).

Thus we have obtained an infinite number of divergent-like conservation laws

$$\frac{}{t}J_m=\frac{}{x}F_m$$
(29)
where $`J_m`$’s and $`F_m`$’s are coefficients of the Taylor series for $`J(\zeta )`$ and $`F(\zeta )`$

$`J(\zeta )`$ $`=`$ $`{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}J_m\zeta ^m`$ (30)
$`F(\zeta )`$ $`=`$ $`{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}F_m\zeta ^m`$ (31)
Some first of the conserved densities are given by

$`J_0`$ $`=`$ $`qr`$ (32)
$`J_1`$ $`=`$ $`i{\displaystyle \frac{q}{x}}r`$ (33)
$`J_2`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(i_2q_{11}q\right)r=\left({\displaystyle \frac{^2q}{x^2}}+q^2r\right)r`$ (34)
$`J_3`$ $`=`$ $`\left({\displaystyle \frac{i}{3}}_3q{\displaystyle \frac{1}{2}}_{21}q{\displaystyle \frac{i}{6}}_{111}q\right)r=i\left({\displaystyle \frac{^3q}{x^3}}+4qr{\displaystyle \frac{q}{x}}+q^2{\displaystyle \frac{r}{x}}\right)r`$ (35)
Note that to present $`J_m`$ for $`m=2,\mathrm{}`$ in a standard way, i.e. in terms of $`q`$, $`r`$ and their derivatives with respect to $`x`$, one has to use evolution equations of the hierarchy (19), (20) and higher. However, it is possible not to use these equations but instead to ’iterate’ the identity (8), which can be rewritten as

$`q^+`$ $`=`$ $`q+i\zeta q_x^+\zeta ^2\left(q^+\right)^2r`$ (36)
$`=`$ $`q+i\zeta q_x\zeta ^2\left[q_{xx}^++\left(q^+\right)^2r\right]i\zeta ^3\left[\left(q^+\right)^2r\right]_x`$ (37)
$`=`$ $`\mathrm{}`$ (38)
or to return to the traditional inverse scattering scheme: it follows from (36) that $`J(\mathrm{t},\zeta )`$ satisfies

$$J=qr+i\zeta r\frac{}{x}\frac{J}{r}\zeta ^2J^2$$
(39)
which leads to the recurrence relation

$`J_0=qr,J_1=q_xr`$ (40)
$`J_{m+1}=r{\displaystyle \frac{}{x}}{\displaystyle \frac{J_m}{r}}{\displaystyle \underset{l=0}{\overset{m1}{}}}J_lJ_{m1l},m1`$ (41)
One can easily identify these equations with the standard for the inverse scattering approach equations for the generating function. Thus, the main result (24) of this paper can be interpreted as follows. Equation (39), if considered as an ordinal differential equation, is the famous Riccati equation which cannot be solved explicitly for arbitrary functions $`q`$ and $`r`$. However in our case $`q`$ and $`r`$ are not arbitrary but related by an infinite number of PDEs of the hierarchy, and it turns out that in this situation, though these restrictions do not determine $`q`$ and $`r`$ uniquely, equation (39) can be solved formally and this solution is given by (24).

Equations (29) can be rewritten as

$$\frac{}{t}I_m=0$$
(42)
where $`I_m`$’s are integrals of the densities $`J_m`$. In the case when $`q`$ and $`r`$ vanish (sufficiently rapidly) as $`x\pm \mathrm{}`$ $`I_m`$’s are given by

$$I_m=_{\mathrm{}}^{\mathrm{}}dxJ_m$$
(43)
In the periodical case, $`q,r(x+L)=q,r(x)`$,

$$I_m=_0^LdxJ_m$$
(44)
while in the case of non-trivial boundary conditions, say the finite-density ones, the integral in the right-hand side should be in some way regularized.

## 4 Conclusion.

To conclude we would like to note the following. The aim of this note is twofold. First we want to demonstrate that the Kyoto school approach can be used not only to reveal and study some mathematical structures behind integrable equations, but also to solve some ’practical’ problems, as one discussed above. Another purpose of this paper is to attract attention to the functional representation of the AKNS hierarchy. This system is one of the first studied integrable models for which the IST has been developed, and since 70’s till now the IST is the main tool to study the NLSE and related problems. As to the methods which were developed later for such equations as, e.g., KP and 2D Toda equations, for our knowledge, their application to the AKNS hierarchy is rather limited. At the same time the formulation of a problem in terms of the functional equations seems be perspective for, e.g., developing some perturbation schemes for various NLSE-related problems.


text:
# Recent X-ray Observations of Disk Galaxies: Tracing the Dynamic Interstellar Medium

## 1 Introduction

X-ray observations of nearby disk galaxies are important to the understanding of the hot interstellar medium (ISM) and its role in galaxy evolution. We are conducting an extensive study of the face-on spiral galaxy M101 based partly on an ultra-deep (230 ks) ROSAT HRI image, and we have also obtained a Chandra observation on the edge-on galaxy NGC4631. These observations, together with data in other wavelength bands, provide us with complementary perspectives of hot gas and its interaction with other galactic components in late-type disk galaxies.

## 2 M101

This face-on galaxy allows us to study the production of hot gas not only under normal circumstances but also under localized starburst conditions. Diffuse $`0.25`$ keV emission is present in the central ($`R10`$ kpc) region of the galaxy (Snowden & Pietsch 1995). Enhanced X-ray emission is also detected in several giant HII complexes (GHCs; Williams & Chu 1995). Wang et al. (1999) have reported the results on discrete X-ray sources detected in the M101 field. But some of these sources apparently represent peaks of extended X-ray-emitting features.

Hot gas (a few times $`10^6`$ K) as traced by the RHRI correlates well with recent star forming regions in the galaxy (e.g., Fig. 1). In particular, enhanced diffuse X-ray emission features line up well with the active southeast arm, which harbors such GHCs as NGC 5461 and NGC 5462, which are more luminous than the 30 Dor complex in the LMC. We further find that the radial distribution of “diffuse” X-ray emission is nearly identical to that of the far-UV radiation and is substantially flatter than that of optical light, which is produced chiefly by old stellar populations. Thus the hot gas originates primarily in recent star formation regions. In comparison, the cooler gas ($`10^6`$ K) as traced by the ROSAT PSPC 0.25 keV image (Snowden & Pietsch 1995) is seen only in the inner 4 radius, consistent with an origin in the galaxy’s corona or a combination of blown-out bubbles.

Fig. 2 provides a close-up of a region at the western end of the two H$`\alpha `$ branches of the southeast spiral arm (Fig. 1). This region happens to be covered by archival HST WFPC2 observations. An apparent ring of HII regions coincides spatially with an HI supershell, as proposed by Kamphuis (1993). Whereas the RHRI emission is primarily associated with the HII ring, the softer PSPC 0.25 keV emission peaks in the interior of the ring. Thus, the X-ray-emitting gas within the ring appears to be cooler than that in the HII regions. If the HI supershell is indeed powered by the mechanical energy input from massive stars, its dynamic age is several times $`10^7`$ yrs, about an order of magnitude greater than that of the HII regions. This can naturally be explained by the propagation of star formation, which is likely triggered by the expansion of the supershell.

The HST image also presents morphological evidence for outflows from HII regions into the interior of ring. These frothy HII regions appear being torn apart, being stripped away, and mass-loading the interior of the supershell. Similar feedback activities are also evident in the GHC NGC 5462.

## 3 NGC4631

ROSAT observations have shown that X-ray emission, particularly in the 0.25 keV band, extends beyond both the cold and warm gas disks of the galaxy (Wang et al. 1995). The Chandra ACIS-S observation of the galaxy provides us with a new opportunity to further the study. While the analysis of the data is still on-going, Fig. 3 shows unambiguously the presence of extraplanar diffuse X-ray emission. The X-ray morphology resembles the radio halo of the galaxy (e.g., Hummel et al. 1991), indicating a close connection among outflows of hot gas, cosmic rays, and magnetic field. A spatially-resolved spectroscopy will further enable us to study the detailed chemical and thermal properties of the hot gas, constraining the dynamics of the corona. The results will have strong implications for our understanding of the feedback from galaxies as well as their structure and evolution.