text:
# High Resolution 4.7 𝜇m Keck/NIRSPEC Spectra of Protostars. I: Ices and Infalling Gas in the Disk of L1489 IRS1footnote 11footnote 1The data presented herein were obtained at the W.M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W.M. Keck Foundation.

## 1 Introduction

In the process of low mass star formation, a mixture of gas, dust, and ices accumulates in protostellar envelopes and disks. The fate of this molecular material is diverse. Most of it will fall toward the protostar and dissociates in the inner disk region or stellar photosphere. Some material will be blown away and destroyed by the stellar wind. Some may also survive and be the building material for comets and planets. Major aspects of this complicated process are not well understood, and poorly observationally constrained. For example, do the ices that form comets still resemble ices of the original pristine molecular clouds or are new ices of different composition being formed in the envelope or disk? The type of ices being formed depends on the composition of the gas that accretes onto grains. Reducing environments produce H2O-rich (‘polar’) ices, while in cold inert environments ‘apolar’ ices rich in CO, N2 and O2 can be formed (Tielens & Hagen, 1982). Depending on the composition, ices evaporate between temperatures of 18 and 90 K. Also, heat can change the solid state structure of ices by for example crystallization. Energetic particles (e.g. cosmic rays) and ultraviolet (UV) radiation are able to initiate reactions in ices and form new species. Dynamics and shocks within disks may be able to destroy ices as well.

Clearly, to determine the relative importance of these ice formation and destruction processes, knowledge of the physical conditions and structure of envelopes and disks is crucial. Much theoretical and observational work on this topic has been done over the last $``$10 years. Molecular gas was detected in a suite of protostellar disks by millimeter wave observations sensitive to emission over radii of several hundred AU (Dutrey, Guilloteau, & Guelin 1997; Thi et al. 2001). Gas phase abundances were found to be reduced by factors of 5 to several 100, depending on the source and the sublimation temperature of the molecules. Models of disk mid-planes indeed show high depletions because of the formation of icy mantles on grains (Aikawa et al., 1997; Willacy et al., 1998). The predicted depletions were in fact higher than observed and thus desorption mechanisms are needed to explain the millimeter wave observations (Goldsmith, Langer, & Velusamy, 1999). It was realized that the outer parts of disks are heated more efficiently when they are flared (Kenyon & Hartmann, 1987). Thus, by the influence of the stellar radiation a layer with ‘super-heated’ dust is formed in which molecules have been dissociated (Chiang & Goldreich, 1997). The layer below that is warm enough to evaporate the ices, but not dissociate the released molecules. The importance of this warm layer, and relative gas phase molecular abundances, depends strongly on how effective ice desorption mechanisms are (Willacy et al., 1998). Recent studies indicate that desorption by UV and X-ray photons may be strong enough to explain observations of molecular gas by millimeter wave telescopes (Willacy & Langer, 2000; Najita, Bergin, & Ullom, 2001). This idea is confirmed by multi transition molecular line observations which indicate temperatures ($`>2040`$ K) and densities that typically occur in this warm layer (van Zadelhoff et al., 2001). In a third layer, the disk mid-plane, cold, dense conditions prevail, resulting in extreme depletions of gas phase molecules many orders of magnitude larger than in quiescent dense molecular clouds. Indeed, recent absorption line observations failed to detect gas phase CO in the edge-on disk around the protostar Elias 18, thus indicating an enormous depletion in the mid-plane (Shuping et al., 2001).

In this Paper we report high spectral resolution ($`R=`$25,000) 4.7 $`\mu \mathrm{m}`$ M band observations of the obscured protostar L1489 IRS (IRAS 04016+2610) in the Taurus Molecular Cloud. L1489 IRS is a low luminosity object (3.7 $`L_{}`$), with a spectral energy distribution resembling that of an embedded class I protostar (Kenyon, Calvet, & Hartmann 1987). Detailed millimeter wave line and continuum studies show that L1489 IRS is surrounded by a large, 2000 AU radius rotating thick disk-like structure (Hogerheijde et al., 1998; Saito et al., 2001), rather than an inside-out collapsing envelope (Hogerheijde & Sandell, 2000). The rotation is sub-Keplerian, and the disk as a whole is contracting. Thus it was suggested that L1489 IRS represents a short-lasting (2$`\times 10^4`$ yr) transitional phase between embedded YSOs that have large envelopes and small (few hundred AU) rotationally supported disks, and T Tauri stars which have no envelopes and fully rotationally supported 500-800 AU size disks (Hogerheijde, 2001). This circumstellar (or circumbinary: Wood et al. 2001) disk is seen close to edge-on (60 to $`<`$90o) in scattered light images (Whitney et al., 1997; Padgett et al., 1999). A CO outflow emanates from the object (Myers et al., 1988) with Herbig Haro objects lying along it (Gómez, Whitney, & Kenyon, 1997).

Low spectral resolution infrared observations of L1489 IRS show that deep H2O (Sato et al., 1990) and CO (Chiar et al. 1998; Teixeira, Emerson, & Palumbo 1998) ice bands are present along the line of sight. In this Paper we will use the newly available spectrometer NIRSPEC at the Keck II telescope to obtain high resolution M band spectra ($`\mathrm{\Delta }v=`$12 $`\mathrm{km}\mathrm{s}^1`$ at 4.7 $`\mu \mathrm{m}`$) of this source. The large array of NIRSPEC allows both the vibrational band of solid CO and the surrounding ro-vibrational transitions of gas phase CO to be observed in the same high resolution spectrum. This offers a new view on this system, both on the origin and evolution of ices and the interrelationship of gas and ices as well as on the kinematics and structure of the young, contracting, close to edge-on disk. It is a unique view, because infrared absorption line studies trace all gas and solid state material at all radii from the star, while present day millimeter wave observations are limited by their relatively low spatial resolution ($``$100 AU). Thus one of the questions that will be answered in this Paper is whether the large scale inward motions seen in millimeter wave emission lines continue to smaller radii (1 AU or less) from the star.

Previous studies have already shown that rich astrophysical information can be obtained from high spectral resolution observations in the atmospheric M band. Mainly massive, luminous protostars were observed, however (Mitchell et al., 1990), and the few observations of low mass protostars cover a small wavelength range containing only a few gas phase lines and not the solid CO band (Shuping et al., 2001; Carr, Mathieu, & Najita, 2001). This is Paper I in a series on high resolution M band spectroscopy of protostars, initiated by the availability of the NIRSPEC spectrometer at Keck II with which weak, low mass protostars can be routinely observed at high spectral resolution over a large wavelength range covering both the solid and gas phase CO features.

The reduction of the long slit spectra is discussed in §2. In §3.1, we analyze the observed gas lines, which at this resolution even give dynamical information. We use standard curve of growth and rotation diagram techniques to get a first idea of gas column and temperatures. In order to analyze the solid CO band profile in this line of sight, a detailed discussion of available laboratory experiments of solid CO is given in §3.2. In §4.1 we apply the infalling disk model of Hogerheijde (2001) to explain the observed gas phase $`{}_{}{}^{12}\mathrm{CO}`$ and $`{}_{}{}^{13}\mathrm{CO}`$ line profiles, and constrain the physical conditions and structure of the disk. The possibility of binarity is briefly discussed in §4.2. The gas phase analysis is linked to the the solid CO results to determine the origin and thermal history of solid CO in §4.3. We conclude with suggestions for future work in §5.

## 2 Observations

The infrared source L1489 IRS was observed with the NIRSPEC spectrometer (McLean et al., 1998) at the Keck II telescope atop Mauna Kea on the nights of 28/29 and 29/30 January 2001. The sky was constantly clear and dry, and the seeing was reasonable ($`0.50.7^{\prime \prime }`$ at 2.2 $`\mu \mathrm{m}`$). NIRSPEC was used in the echelle mode with the 0.43$`\times 24^{\prime \prime }`$ slit, providing a resolving power of $`R=\lambda /\mathrm{\Delta }\lambda =25,000`$ ($``$12 $`\mathrm{km}\mathrm{s}^1`$)222This nominal spectral resolution was verified by measuring the Gaussian width of absorption lines in a variety of astrophysical sources. with three Nyquist sampled settings covering the wavelength range 4.615–4.819 $`\mu \mathrm{m}`$ in the atmospheric M transmission band.

The data was reduced in a standard way, using IDL routines. The thermal background emission was removed by differencing the nodding pair. Each nodding position was integrated on for 1 minute, before pointing the telescope to the other nodding position. The positioning had to be done by hand because condensations on the dewar window prevented us from using the image rotator, and the instrument had to be used in the non-standard ‘stationary guiding mode’. In this mode, the condensations were well removed by subtracting nodding pairs, although the sensitivity of Keck/NIRSPEC was reduced by several magnitudes.

The most critical step in the reduction of this data is correction for atmospheric absorption features, and separate telluric and interstellar CO absorption features. For our observations of L1489 IRS this proved to be relatively easy, because at the time of our observations the deep interstellar CO lines are shifted by as much as 43 $`\mathrm{km}\mathrm{s}^1`$ with respect to telluric CO lines. The low airmass of the source, 1.00–1.05, also made this step easier. The standard stars HR 1380 (A7V) and HR 1497 (B3V) were observed at similar airmass, and their spectral shape and hydrogen absorption features were divided out with Kurucz model atmospheres. A good telluric correction was achieved, although residuals are seen near deep atmospheric lines. To be safe, we therefore removed the parts of the spectrum which have less than 50% of the maximum transmission in each setting. This does not affect most of the features in the L1489 IRS spectrum, because the velocity shift of 43 $`\mathrm{km}\mathrm{s}^1`$ fully separates them from telluric features. The final signal-to-noise value on the unsmoothed spectrum is $`45`$, after integration times of 11, 14, and 20 minutes on each of the three settings of this $`M=4.8`$ magnitude (Kenyon et al., 1993) source. The different settings were wavelength calibrated on the atmospheric CO emission lines, and subsequently the three settings were combined by applying relative multiplication factors. We have not attempted to flux calibrate the spectrum, since we are interested in absorption features only.

## 3 Results

The fully reduced echelle spectrum of L1489 IRS shows, in great detail, many deep narrow absorption lines of gas phase $`{}_{}{}^{12}\mathrm{CO}`$ and $`{}_{}{}^{13}\mathrm{CO}`$, and a few weak lines of $`\mathrm{C}^{18}\mathrm{O}`$ (Fig. 1). These lines were identified, using the line frequencies in the HITRAN catalogue (Rothman et al., 1992). The broad absorption feature between 2122-2149 $`\mathrm{cm}^1`$ can be attributed to the stretching vibration mode of $`{}_{}{}^{12}\mathrm{CO}`$ in ices along this line of sight.

In order to analyze these gas and solid state absorption features, a shallow, second order polynomial continuum was applied to derive the optical depth spectrum. The solid CO band was then analyzed using available laboratory experiments (§3.2). The derivation of physical parameters from the gas phase lines is highly model dependent. First, we will derive temperatures and column densities using the standard curve of growth and rotation diagram techniques (§3.1). Then, we will independently test an astrophysically relevant power law model in §4.1. This information is combined in §4.3 to discuss the origin and thermal history of the solid CO seen in this line of sight.

### 3.1 Gas Phase CO

The $`{}_{}{}^{12}\mathrm{CO}`$ lines have a complicated profile. Deep lines are present at a velocity of $`+43`$ $`\mathrm{km}\mathrm{s}^1`$ with respect to earth, which is the systemic velocity of L1489 IRS at the date these observations were done (taking $`v_{\mathrm{lsr}}=`$ 5 $`\mathrm{km}\mathrm{s}^1`$ from millimeter emission lines; e.g. Hogerheijde & Sandell 2000). Each of the $`{}_{}{}^{12}\mathrm{CO}`$ lines is accompanied by a fairly prominent wing at the red shifted side (Fig. 3.1). As an exploring step in the analysis, we decomposed the main $`{}_{}{}^{12}\mathrm{CO}`$ component and its wing by fitting two Gaussians. They are separated by on average 23$`\pm `$6 $`\mathrm{km}\mathrm{s}^1`$, and the main feature and its wing have widths of FWHM=20$`\pm `$3 $`\mathrm{km}\mathrm{s}^1`$ and 53$`\pm `$17 $`\mathrm{km}\mathrm{s}^1`$ respectively. The main $`{}_{}{}^{12}\mathrm{CO}`$ feature at the systemic velocity is resolved. Excess absorption is visible at the blue and red sides with respect to $`{}_{}{}^{13}\mathrm{CO}`$ (Fig. 3.1). The $`{}_{}{}^{13}\mathrm{CO}`$ lines have Gaussian shapes with a width equal to the instrumental resolution (FWHM=12 $`\mathrm{km}\mathrm{s}^1`$). The main absorption feature of the $`{}_{}{}^{12}\mathrm{CO}`$ lines is therefore only in part responsible for the same gas seen in $`{}_{}{}^{13}\mathrm{CO}`$.

We derived equivalent widths for the isotopes and $`{}_{}{}^{12}\mathrm{CO}`$ components (Table 1), and applied a standard curve of growth technique to calculate column densities for each $`J`$ level. The comparison of equal $`J`$ levels of $`\mathrm{C}^{18}\mathrm{O}`$, $`{}_{}{}^{13}\mathrm{CO}`$, and main $`{}_{}{}^{12}\mathrm{CO}`$ component provides a handle on the intrinsic line width $`b_\mathrm{D}`$ (=FWHM/$`2\sqrt{\mathrm{ln2}}`$), in the assumption that all material absorbs at the same velocity (however, see §4.1). One also has to assume that the isotope ratios are constant along the line of sight ($`{}_{}{}^{12}\mathrm{CO}`$/$`{}_{}{}^{13}\mathrm{CO}`$=80 and $`{}_{}{}^{12}\mathrm{CO}`$/$`\mathrm{C}^{18}\mathrm{O}`$=560; Wilson & Rood 1994). The equivalent widths of $`\mathrm{C}^{18}\mathrm{O}`$ and $`{}_{}{}^{13}\mathrm{CO}`$ are then simultaneously fit at $`b_\mathrm{D}`$$`>`$0.8 $`\mathrm{km}\mathrm{s}^1`$. Lower $`b_\mathrm{D}`$ significantly ($`>3\sigma `$) underestimates the $`{}_{}{}^{13}\mathrm{CO}`$ lines, with respect to $`\mathrm{C}^{18}\mathrm{O}`$. As mentioned above, the main $`{}_{}{}^{12}\mathrm{CO}`$ component is clearly contaminated by gas not seen in $`{}_{}{}^{13}\mathrm{CO}`$, both on the blue and red shifted sides. As a first order correction we lowered the $`{}_{}{}^{12}\mathrm{CO}`$ equivalent widths in Table 1 with 40%, which corresponds to unresolved lines at the observed peak optical depth. Then, we find that all isotopes are best fit simultaneously at $`b_\mathrm{D}`$=1.4$`\pm 0.1`$ $`\mathrm{km}\mathrm{s}^1`$. This would be an upper limit if the correction for contamination of the $`{}_{}{}^{12}\mathrm{CO}`$ lines were too small. Hence, this curve of growth analysis of the CO isotopes shows that, in the assumption that all gas absorbs at the same velocity, $`b_\mathrm{D}`$ is limited to 0.8$`<`$ $`b_\mathrm{D}`$$`<`$1.5 $`\mathrm{km}\mathrm{s}^1`$.

With the column densities per $`J`$ level at hand, a rotation diagram was constructed for the $`{}_{}{}^{13}\mathrm{CO}`$ lines (Fig. 3.1) to derive the total column density and temperature at a number of allowed $`b_\mathrm{D}`$ values (Table 2). Clearly, the rotation diagram shows a double temperature structure, much resembling that of high mass objects (Mitchell et al., 1990): cold ($`T15`$ K), and warm gas ($`T250`$ K) are present along the same line of sight. The column density of the cold component toward L1489 IRS is a particularly strong function of $`b_\mathrm{D}`$, increasing with an order of magnitude from 0.7 to 1.3 $`\mathrm{km}\mathrm{s}^1`$. An independent CO column of 1.4$`\times 10^{19}`$ $`\mathrm{cm}^2`$ can be estimated from $`A_\mathrm{V}`$=29 (Myers et al., 1987), assuming the dust and gas are co-spatial. This would suggest $`b_\mathrm{D}`$ is in the 1.0-1.3 $`\mathrm{km}\mathrm{s}^1`$ range.

For the highly red shifted $`{}_{}{}^{12}\mathrm{CO}`$ wing a rotation diagram is constructed as well. Here, we assume that the absorption is optically thin, because the wings are absent in $`{}_{}{}^{13}\mathrm{CO}`$, as evidenced by the high signal-to-noise average line profile (Fig. 3.1). The rotational temperature of this highly red shifted gas is similar to that of the warm $`{}_{}{}^{13}\mathrm{CO}`$ component (250$`{}_{100}{}^{}{}_{}{}^{+500}`$ K), however the column is an order of magnitude lower (Table 2). Finally, assuming that the absorption in the wings of the resolved $`{}_{}{}^{12}\mathrm{CO}`$ main component is optically thin, the column of this blue and red shifted low velocity gas (within 10 $`\mathrm{km}\mathrm{s}^1`$ of the systemic velocity) is about 25% of the column of the high velocity red shifted gas at a similar, although poorly determined, temperature.

Summarizing, this basic analysis indicates that, from $`{}_{}{}^{13}\mathrm{CO}`$, both a large column of cold 15 K gas and a significant amount of warm gas ($`T250`$ K) are present within $`3`$ $`\mathrm{km}\mathrm{s}^1`$ of the systemic velocity (Table 2). The $`{}_{}{}^{12}\mathrm{CO}`$ lines show that warm gas at $`T250`$ K is also present at highly red shifted velocities ($`20100`$ $`\mathrm{km}\mathrm{s}^1`$), but at a factor 10 lower column. A small amount of warm gas is present at low red and blue shifted velocities as well (within 10 $`\mathrm{km}\mathrm{s}^1`$). As further described in §4.1, the gas components at the red shifted and systemic velocities can be fitted within the same physical model of a contracting disk, but the origin of the warm gas at low blue shifted velocities is more difficult to explain.

### 3.2 Solid CO

The broad absorption feature between 2122-2149 $`\mathrm{cm}^1`$ (Fig. 1) can be entirely attributed to the stretching vibration mode of $`{}_{}{}^{12}\mathrm{CO}`$ in circumstellar ices (§4.3). The high spectral resolution allows, for the first time, to unambiguously separate the gas phase CO lines from the solid state absorption and study the solid CO band profile in great detail. In accordance with previous, low resolution studies (Chiar et al., 1998; Teixeira et al., 1998), a distinct narrow feature is seen at 2140 $`\mathrm{cm}^1`$, and a significantly broader component at longer wavelengths. Our data however, indicates the presence of a new, third component on the blue side, separate from the narrow 2140 $`\mathrm{cm}^1`$ feature, most notable by a change of the blue slope at 2142 $`\mathrm{cm}^1`$ (Fig. 1).

In order to explain the shape of this CO absorption profile, we have taken laboratory experiments of solid CO from the literature (Sandford et al. 1988; Schmitt, Greenberg, & Grim 1989; Palumbo & Strazzulla 1993; Trotta 1996; Ehrenfreund et al. 1997; Elsila, Allamandola, & Sandford 1997; Baratta & Palumbo 1998; Teixeira et al. 1998). We determined peak position and width of the laboratory profiles as a function of ice composition, temperature, and cosmic ray bombardment intensity. Additionally, for CO–rich ices we used optical constants to calculate the absorption profile as a function of particle shape in the small particle, Rayleigh limit (for details, see Ehrenfreund et al. 1997). The absorption profile of CO–poor ices (concentration $`<`$30%) is not affected by these particle shape effects. Changing these various parameters gives a wide variety of peak positions (2135–2144 $`\mathrm{cm}^1`$) and widths (FWHM=1.5–14 $`\mathrm{cm}^1`$) in the laboratory. We limit ourselves here by identifying general trends, as summarized in Fig. 3.2. These trends are then related to the Gaussian peak position and width of the three aforementioned components observed toward L1489 IRS (also given in Fig. 3.2).

The absorption band of a thin film of pure, solid CO at $`T=`$10 K peaks at $``$2137 $`\mathrm{cm}^1`$, and is very narrow ($``$2 $`\mathrm{cm}^1`$). This clearly does not fit any of the components observed toward L1489 IRS. However, for ellipsoidally shaped grains (more precisely, a distribution of ellipsoidal shapes, ‘CDE’; Bohren & Huffman 1983), this sensitive strong band shifts to shorter wavelengths and becomes broader. Now it exactly fits the 2140 $`\mathrm{cm}^1`$component observed toward L1489 IRS, both in peak position and width (Fig. 3.2). In view of a recent controversy on optical constants (Ehrenfreund et al., 1997), it is worth noting that good fits are obtained with optical constants from the works of Ehrenfreund et al. (1997) and Baratta & Palumbo (1998). The optical constants of Trotta (1996) do not induce strong particle shape effects and therefore pure CO does not provide a good fit to L1489 IRS for any particle shape.

Broadening of the laboratory profile, in order to fit the 2140 $`\mathrm{cm}^1`$ feature in L1489 IRS, is also achieved by adding a small amount of CO2, O2, or H2O molecules. To avoid a too large broadening, and minimize aforementioned particle shape effects, this mixture needs be diluted in N2. This astrophysically relevant molecule does not broaden the feature, and gives a small blue shift (Ehrenfreund et al., 1997; Elsila et al., 1997), required to fit the 2140 $`\mathrm{cm}^1`$ feature in L1489 IRS. Thus, both this mixture, as well as ellipsoidally shaped pure CO ice grains, provide good fits to the central 2140 $`\mathrm{cm}^1`$ feature.

Now, the width significantly increases, and the peak shifts to longer wavelengths by diluting CO in a mixture of molecules with large dipole moments such as H2O or CH3OH (Sandford et al., 1988; Tielens et al., 1991). This particular behavior is needed to fit the broad long wavelength wing seen toward L1489 IRS. Solid CH3OH has an abundance less than a few percent of solid H2O toward low mass objects (Chiar, Adamson, & Whittet, 1996). H2O seems the best dilutant, because of its large interstellar abundance. The use of H2O:CO mixtures requires interesting additional constraints. Although a low temperature, unprocessed H2O–rich ice does provide a good fit to the red wing, it can be excluded based on the presence of a prominent second absorption at $``$2150 $`\mathrm{cm}^1`$ in the laboratory, which is clearly absent toward L1489 IRS. This second peak is caused by CO molecules located in pockets in an amorphous ice. These CO molecules are weakly bound, and the $``$2150 $`\mathrm{cm}^1`$ peak disappears rapidly at higher $`T`$ or as a result of cosmic ray hits (Sandford et al., 1988). Thus, the H2O ice responsible for the long wavelength wing toward L1489 IRS must be thermally ($`T>`$ 50 K) or energetically processed.

The blue wing seen at $``$2143 $`\mathrm{cm}^1`$ toward L1489 IRS can only be explained by an apolar ice. Adding a significant amount of CO2 to a CO ice (CO2/CO$`>1`$) results in the blue shift and broadening required to fit the observed wing. Somewhat less CO2 is needed (CO2/CO$`0.5`$) when a large amount of O2 is present. N2 may be added as well, but is not essential except as a dilutant to reduce the effects of particle shape. A good fit is obtained by the mixture N2:O2:CO2:CO=1:5:0.5:1, as proposed in Elsila et al. (1997). If CO2/CO$`<0.5`$ the band peaks at too high wavelength. In CO–rich ices this problem can however be overcome by particle shape effects (Fig. 3.2). In view of this effect it is not possible to constrain the relative molecular abundances of this interstellar component in more detail, but it is clear that an apolar CO2 or O2 ice is needed, different from the distinct 2140 $`\mathrm{cm}^1`$ feature.

A three component fit to the entire CO ice band of L1489 IRS is shown in Fig. 3.2. Although this is not a unique fit, it does obey the global trends that we identified in the laboratory experiments.

Finally, the solid CO column density is derived by dividing the integrated optical depth over the band strength $`A`$. We take $`A=1.1\times 10^{17}`$ cm molecule-1 independent of ice composition (Gerakines et al., 1995), and thus find $`N`$(solid CO)=6.5$`\times 10^{17}`$ $`\mathrm{cm}^2`$. The main source of uncertainty here is in $`A`$, which is about 10%. CO in polar ices contributes 3.5$`\times 10^{17}`$ $`\mathrm{cm}^2`$ to the total column and the apolar components at 2140 and 2142 $`\mathrm{cm}^1`$ contribute each 1.5$`\times 10^{17}`$ $`\mathrm{cm}^2`$.

## 4 Discussion

### 4.1 An Infalling Disk

The astrophysical meaning of the apparent two component temperature structure seen in the $`{}_{}{}^{13}\mathrm{CO}`$ rotation diagram (Fig. 3.1) requires further investigation. For high mass protostars it was found that similar rotation diagrams can be ‘mimicked’ by power law models of spherical envelopes (van der Tak et al., 2000).

For L1489 IRS, the detection of molecular gas at a range of temperatures and red shifted velocities could indicate the presence of infalling gas at a range of radii from the protostar. Indeed, a 2000 AU radius contracting, disk-like structure was found in millimeter wave interferometer data (Hogerheijde & Sandell, 2000). In a detailed follow-up study, Hogerheijde (2001) adopts a flared-disk model based on Chiang & Goldreich (1997) with a radial power-law distribution for the temperature

$$T=34(R/1000\mathrm{AU})^{0.4}\mathrm{K},$$
(1)
and a density distribution that has a power-law drop-off with radius and a vertical exponential drop-off with scale height $`h`$

$$\rho (R,z)=\rho _0(R/1000\mathrm{A}\mathrm{U})^{1.5}\mathrm{e}^{z^2/h^2}\mathrm{kg}\mathrm{cm}^3.$$
(2)
The scale-height $`h`$ is assumed to be a simple function of R, $`h=R/2`$. An inward-directed radial velocity field described as

$$V_{\mathrm{in}}=1.3(R/100\mathrm{AU})^{0.5}\mathrm{km}\mathrm{s}^1$$
(3)
is inferred, in addition to Keplerian rotation around a 0.65 M central star.

Can this contracting disk model, based on (sub-) millimeter emission observations with angular resolution of $`4`$$`8^{\prime \prime }`$, reproduce the observed infrared CO absorption line profiles measured along a pencil beam? The absorption lines are modeled with the radiative transfer code of Hogerheijde & van der Tak (2000); the high densities in the disk ensure LTE excitation for the lines involved, and line trapping is neglected in the excitation calculation. The model spectra include dust opacity at a standard gas/dust ratio, as well as a $`N(\mathrm{CO})=1\times 10^{18}`$ cm-2 column of cold foreground material (15 K; Hogerheijde 2001); both factors do not affect the spectra in any significant way. The calculated spectrum is convolved with a Gaussian of FWHM=12 $`\mathrm{km}\mathrm{s}^1`$, which is the NIRSPEC instrumental resolution.

We find that, while keeping all other parameters the same as in Hogerheijde (2001), the assumed density profile sensitively influences the wings of the 12CO lines. This is enhanced by the fact that we are observing the flared disk of L1489 IRS at an inclination between $`60^{}`$ and $`<90^{}`$ (cf., Padgett et al. 1999), and the pencil beam crosses the disk at a few scale heights. Small changes in the density profile, for example induced by the thermal structure, have a large effect on the absorption line profile. In the model of Hogerheijde (2001) the scale height increases linearly with distance from the star, and thus the density $`\rho (l)`$ along the line of sight $`l`$ follows the density in the mid-plane (Eq. 2) reduced by a factor $`e^{4/\mathrm{tan}^2(\alpha )}`$, with $`\alpha `$ the inclination. Here, we include the effect of density variations, or deviations from the adopted scale height $`h=R/2`$, by relaxing the values of the density along the line of sight $`l`$ by fitting $`\rho (l)=\rho _0(l/1000\mathrm{AU})^p`$ to the data.

This initial model successfully fits the peak velocity and depth of both high and low $`J`$ 13CO lines (Fig. 2). Its rotation diagram is quite different from that of the curve of growth analysis (Fig. 3.1), showing that rotation diagrams must be interpreted with great care. Our model also matches the range of velocities observed in the red wings of the 12CO lines, when taking $`p=0.55\pm 0.15`$. This is a much shallower density profile compared to that derived from millimeter wave data ($`p=1.5`$; Eq. 2) and indicates that the scale height increases more than linear, i.e. the disk flares more than assumed in Hogerheijde (2001). With this result, it is possible to determine the important relation of disk scale height $`h(R)=a.R^b`$ as a function of $`R`$, but only if the disk inclination is a priori known. Unfortunately the inclination is not better constrained than within the range of $`60^{}`$ and $`<90^{}`$ imposed by near-infrared data (Padgett et al. 1999). We can therefore not distinguish between low and high values of $`a`$ and corresponding high and low inclinations respectively. In either case, the total 12CO column along the pencil beam is $`1.2\times 10^{19}`$ $`\mathrm{cm}^2`$, with 58% of the CO mass at a temperature of $`T=2060`$ K, 15% at $`6090`$ K, and 27% at $`6090`$ K. This result is of importance in §4.3 in the interpretation of the solid CO observations, and in particular in assessing the thermal history of ices. The total column of our model is in good agreement with the column derived from the visual extinction ($`1.4\times 10^{19}`$ $`\mathrm{cm}^2`$; §3.1). It is also of the same order of magnitude as the total column through the mid-plane, calculated from dust and line emission ($`N[\mathrm{CO}]=6\times 10^{18}`$ $`\mathrm{cm}^2`$; Hogerheijde 2001), and confirms the relatively edge-on orientation of the disk.

However, apart from these successes, the $`{}_{}{}^{12}\mathrm{CO}`$ lines show that our infalling disk model produces too much warm gas at high velocities (Fig. 2). The 12CO lines are a factor of 2.5 deeper, and, in contrast to 13CO, they peak at a too high velocity (+10 $`\mathrm{km}\mathrm{s}^1`$) with respect to the observations. In principle, one could make the $`{}_{}{}^{12}\mathrm{CO}`$ lines less deep by assuming that $`1\%`$ of the original, unextincted continuum flux (corresponding to 30% of the extincted continuum) reaches the slit without passing through the disk, by scattering on large grains. The shift in peak velocity however requires a solution of a more fundamental origin. Perhaps the infall velocity function is shallower, and the disk is more rotationally supported at lower radii. The amount of warm gas at high velocities can also be lowered by assuming that only part of the disk participates in the high velocity inflow, such as a thin hot surface layer, or gas accelerated in magnetic field tubes directed from the inner disk to the stellar photosphere. Such a two component model is consistent with the rotation diagram derived from the curve of growth (Fig. 3.1), and also with the low observed mass accretion rate. If we take the inflow at face value, and assume that the entire disk participates, the mass accretion rate would be $`10^6`$ M, generating 7 L in accretion luminosity. The star’s Lbol is estimated at 3.7 L which also contains the stellar luminosity. It is therefore indeed likely that the mass accretion onto the star is significantly lower, as is also traced through the lack of the hydrogen Pf$`\beta `$ emission line in our spectrum (2148.8 $`\mathrm{cm}^1`$; Fig. 1) and the weakness of Br$`\gamma `$ emission (Muzerolle, Hartmann, & Calvet, 1998).

### 4.2 Binarity?

An entirely different explanation for the line profiles may lie in the possibility that L1489 IRS is a protobinary system. A protobinary nature of L1489 IRS is suggested by various pieces of evidence (Lucas, Blundell, & Roche 2000; Wood et al. 2001 and references therein). The presence of a quadrupolar outflow system is inferred from K band polarization images, C18O emission line profiles, Herbig-Haro knots that are scattered throughout the L1489 IRS environment, and a very complex near infrared scattered light pattern. Three dimensional models, in which the axisymmetry of the infalling circumstellar envelope is broken by multiple outflow cavities that are perpendicular to each other, are able to account for the observed morphology.

The putative binary itself, however, has not been resolved so far. An upper limit on the projected separation from near infrared images has been set on $`<`$ 20 AU (Padgett et al., 1999). If the CO absorption line profile is in any way related to a binary system, then the large observed velocities ($``$23 $`\mathrm{km}\mathrm{s}^1`$; §3.1) may indeed favor a close binary system. The line profile is then expected to vary on a time scale of a few months, which can easily be tested. In this case, much of the observed warm gas might be present in two small circumstellar disks, which are in a close orbit around each other. Some of the warm gas may also be present at low density in the central cavity created by the binary. The large column of cold gas may originate in the circumbinary disk. The binary extracts momentum from the 2000 AU circumbinary disk, setting up the inward motion seen in millimeter wave emission lines. We leave further investigation of this topic for future studies.

### 4.3 The Origin and Evolution of Ices

In order to establish if the solid CO observed toward L1489 IRS originates in foreground clouds or in a circumstellar (or binary) disk, it is worth to compare with ices observed in lines of sight not affected by star formation. Observations of field stars obscured by intervening quiescent material of the Taurus Molecular Cloud have revealed that solid CO is not present when the extinction $`A_\mathrm{V}5`$ (e.g. Teixeira & Emerson 1999). The solid CO toward L1489 IRS can therefore not be associated with foreground clouds, which have a gas column of $`N`$(12CO)= 1$`\times 10^{18}`$ $`\mathrm{cm}^2`$ (Hogerheijde, 2001), corresponding to $`A_\mathrm{V}2`$.

Thus, the solid CO must be present in the disk of L1489 IRS. The absorption profile is intriguingly different from that seen in quiescent clouds. The broad red wing has a depth of $``$30% with respect to the narrow 2140 $`\mathrm{cm}^1`$ peak, which is significantly more than toward all measured background stars (10%; Chiar et al. 1995). This may well be an effect of thermal processing along the L1489 IRS line of sight, because the sublimation temperatures of polar and apolar ices, causing the broad and narrow features respectively, are very different (90 versus 18 K). However, a chemical origin of an increased abundance of polar ices in disks cannot be excluded, because the apparently edge-on system Elias 18 in the Taurus Molecular Cloud has an extremely large CO depletion factor (solid/\[gas+solid\]$``$100% versus 7% for L1489 IRS), but a deep red ‘polar’ CO wing is present as well (Shuping et al., 2001; Chiar et al., 1998). On the other hand, energetic processing may take place even in the cold disk of Elias 18 (Whittet et al., 2001). Clearly, it is necessary to observationally characterize the ices in circumstellar disks in much more detail.

If for now we assume the sublimation scenario, we can do some general extrapolations which can be compared with the results of our gas phase study (§4.1). By scaling the long wavelength wing of solid CO of background field stars to that of L1489 IRS, we find that a column of 6$`\times 10^{17}`$ $`\mathrm{cm}^2`$ of CO has evaporated from the apolar ice component in the part of the L1489 IRS disk along the pencil absorption beam where $`T<90`$ K (the sublimation temperature of polar ices). Then the column of solid CO that went from the quiescent cloud into building this part of the disk is 12.5$`\times 10^{17}`$ $`\mathrm{cm}^2`$. Extrapolating this further, we use the observed CO depletion factor of 30% toward field stars behind the Taurus Molecular Cloud (Chiar et al., 1995) to calculate that the original quiescent gas column must have been of the order of 3$`\times 10^{18}`$ $`\mathrm{cm}^2`$. Adding the evaporated column, the expected present day gas column at $`T<90`$ K is 3.6$`\times 10^{18}`$ $`\mathrm{cm}^2`$. This is of the same order of magnitude as the CO column below 90 K in our collapsing disk model ($`N`$\[CO\]=8.7$`\times 10^{18}`$ $`\mathrm{cm}^2`$), which may indicate that no chemical change in the apolar/polar CO ice ratio and no significant additional depletion has occurred in the evolution from quiescent Taurus Molecular Cloud material to the formation of the L1489 IRS disk. This contrasts strongly with the very large depletions found in the (older) disks of T Tauri stars (Dutrey et al., 1997). The low CO depletion along the pencil beam toward L1489 IRS (7%) and the supposed signs of thermal processing (see below) may be due to the fact that our line of sight does not cross the disk mid-plane, i.e. the system is not exactly edge-on. The ice processing we see takes place higher in the disk atmosphere, perhaps in the warm layer below the super-heated dust layer responsible for millimeter wave line emission (van Zadelhoff et al., 2001). It must be noted that in the model of Hogerheijde (2001) the gas temperatures are larger than 25 K, prohibiting the formation of apolar ices and large CO depletions anywhere in the disk. The observed presence of apolar CO ices thus indicates that, as already suggested in §4.1, the line of sight may cross the cold, rotationally supported disk interior not traced in the observations and infall model of Hogerheijde (2001).

Apart from evaporation of apolar ices, other hints of thermal processing include the aforementioned absence of the 2150 $`\mathrm{cm}^1`$ absorption (§3.2), which occurs in cold unprocessed polar CO ices but disappears at temperatures $`T>50`$ K. Also, the blue apolar wing may be a consequence of thermal processing. If the central 2140 $`\mathrm{cm}^1`$ peak is due to a mixture of O2, N2 and CO instead of pure CO (spectroscopically these cannot be distinguished), thermal or energetic processing (UV radiation from the ISRF, UV induced by H2 cosmic ray collisions, or direct hits of cosmic rays) could efficiently produce CO2. This could cause the band to broaden and shift to the position of the observed blue wing. Chemical models indicate that energetic processing of molecules in disks takes place on a time scale of 106 yrs (Aikawa et al., 1999), which is somewhat longer than the age of the disk of L1489 IRS ($`5\times 10^5`$ yrs). This however applies to the disk mid-plane, and the time scale may well be shorter in the lower density higher disk layers that our observations of L1489 IRS trace. A possible problem with the energetic processing interpretation is the absence of a feature adjacent to the short wavelength side of the CO ice band, usually attributed to energetically produced C$``$N bondings (Whittet et al., 2001). Another spectroscopic tracer of thermal processing is the signature of crystallization in the band profiles of H2O and CO2 ices. Our infalling disk model predicts that only 15% of the gas is within the temperature range at which ices crystallize (60–90 K), and thus crystallization is not expected to play a significant role in the disk of L1489 IRS. This model prediction can be tested with future high quality H2O and CO2 spectra of L1489 IRS.

In summary, several pieces of evidence indicate that the CO ices in the disk of L1489 IRS have experienced thermal or energetic processing. The strongest arguments are the low depletion factor and the low ratio of apolar to polar ices with respect to the quiescent Taurus Molecular Cloud material. This may be explained by the fact that the disk of L1489 IRS is seen under an angle, and our pencil absorption beam traces the warm upper disk layers.

## 5 Summary and Future Work

We have shown that valuable and unique information is obtained from high resolution spectroscopy of the CO fundamental at 4.7 $`\mu \mathrm{m}`$ toward the low mass class I protostar L1489 IRS in the Taurus Molecular Cloud. At a resolution of $`R=25,000`$ (12 $`\mathrm{km}\mathrm{s}^1`$) this object shows a multitude of deep ro-vibrational absorption lines of $`{}_{}{}^{12}\mathrm{CO}`$, as well as $`{}_{}{}^{13}\mathrm{CO}`$ and $`\mathrm{C}^{18}\mathrm{O}`$. The isotopes trace large columns of warm and cold gas in the circumstellar disk at or within 3 $`\mathrm{km}\mathrm{s}^1`$ of the systemic velocity, while the $`{}_{}{}^{12}\mathrm{CO}`$ line profiles show warm gas that is red shifted at a range of velocities of up to 100 $`\mathrm{km}\mathrm{s}^1`$. Both the line depth of the isotopes and the extent of the red shifted warm gas seen in $`{}_{}{}^{12}\mathrm{CO}`$ are well explained by an infalling flared disk model with power laws for the temperature, infall velocity, and density (a small column of blue shifted gas seen in the $`{}_{}{}^{12}\mathrm{CO}`$ line wing however remains unexplained). These observations show that the inward motions inferred on scales of several hundred AU through millimeter wave interferometry, continue to within 0.1 AU of the star, where the velocity model of Hogerheijde (2001) predicts inward velocities exceeding several tens of $`\mathrm{km}\mathrm{s}^1`$. A detailed comparison of our power law infall model however overestimates the amount of warm, high velocity infalling gas. Much of this gas must therefore be rotationally supported, and only a thin disk surface layer is infalling, or gas is accelerated along magnetic fields in the inner parts of the disk. High spatial resolution millimeter wave observations (with ALMA) are needed to test our model, e.g. to refine the determination of the velocity field in the inner disk parts and the dependence of the disk scale height on radius $`h(R)`$. Finally, high spatial resolution infrared interferometer observations would be able to see if L1489 IRS is a close binary system ($`<`$ 20 AU), which is essential to assess the importance of this aspect to the observed line profiles and disk evolution.

In the same spectrum, a deep CO ice band is seen toward L1489 IRS, which we conclude to originate in the disk. The high signal-to-noise and spectral resolution allowed us to separate it from circumstellar CO absorption lines, and to do a detailed comparison with laboratory mixtures. It is, for the first time, found that besides the well known long wavelength wing of polar H2O:CO ices and the central narrow absorption peak at 2140 $`\mathrm{cm}^1`$, an additional, separate component is present on the blue side of the ice band. Both the central peak and the blue shifted component are due to CO in apolar ices, one due to pure CO and one due to a mixture of CO2 and perhaps O2 and N2. The relatively large depth of the long wavelength, polar wing relative to the apolar components may suggest that the ices in the disk of L1489 IRS are thermally processed with respect to the quiescent Taurus Molecular Cloud material, which is supported by the low depletion factor ($``$7%) in this line of sight. This is likely due to the fact that we see the disk under an angle, and our pencil absorption beam traces the upper disk layers.

The present work shows that high spectral resolution 4.7 $`\mu \mathrm{m}`$ observations are a great tool to better understand protostellar disks, which define the initial conditions of planet and comet formation both in the solid and gas phase. This is an exploratory study and needs to be followed up by observing a larger sample of protostellar disks at 4.7 $`\mu \mathrm{m}`$ to investigate the influence of parameters such as foreground contribution, disk inclination, age, and binarity on gas and ice band profiles.

The UCLA NIRSPEC instrument team is acknowledged for constructing an excellent instrument. We thank M. Brown (Caltech, USA) for providing a procedure to rectify 2D spectra in IDL, T. Teixeira (University of Århus, Denmark) for generously providing her compilation of laboratory CO spectra to us, and the anonymous referee for careful reading of the manuscript which helped to clarify the description of several aspects of the analysis. The research of ACAB at the Caltech Submillimeter Observatory is funded by the NSF through contract AST-9980846. The research of MRH is supported by the Miller Institute for Basic Research in Science. NASA support to GAB is gratefully acknowledged. The authors wish to extend special thanks to those of Hawaiian ancestry on whose sacred mountain we are privileged to be guests. Without their generous hospitality, none of the observations presented herein would have been possible.


format: text

text:
# Darboux Transformations, Infinitesimal Symmetries and Conservation Laws for Nonlocal Two–Dimensional Toda Lattice

## 1 Introduction

It is well known that some of nonlocal differential equations of great physical significance (e.g., the Benjamin–Ono (BO) equation , the intermediate long wave (ILW) equation ) possess the same mathematical properties as local nonlinear equations integrable in the frameworks of the inverse scattering transformation (IST) method . Such the nonlocal equations were found to have the multi–soliton solutions ,, infinitely many conservation laws ,, the Bäcklund transformation , to pass the Painlev$`\stackrel{´}{\text{e}}`$ test , to be represented in the bilinear form , and as the compatibility condition of the overdetermined linear system (Lax pair) , to possess the algebro–geometric solutions . An effective technique of Darboux transformation (DT)  was also applied to nonlocal equations, such as the ILW equation, nonlocal analogue of the Kadomtsev–Petviashvili equation  and nonlocal Toda equation . This technique allows one to obtain the infinite hierarchies of solutions of Lax pairs and associated nonlinear equations to arise from the compatibility condition. Carrying out the proper sequence of the Darboux transformations adds the one–soliton component of solution to initial solution of nonlinear integrable equations.

The noncommutative (”quantum”) generalization of the spectral parameter of usual Lax pair was suggested in  to construct the hierarchies of nonlocal counterparts of the nonlinear equations admitting the compatibility condition representation. This way the hierarchies of the ILWn equations, modified ILWn equations and nonlocal Toda lattice were obtained . Nonlocal bilinear equations for the last system, which were proposed in work , should be modified for the multi–soliton solutions to be derived by applying the ordinary procedure of the perturbation theory to a constant solution . It is also inconvenient in the frameworks of the bilinear approach to determine a set of exponents, whose mutual products appear no in the expansion for the multi–soliton solutions, and, consequently, to describe the soliton solutions.

In the present paper we extend the DT technique to nonlocal partners of Lax pairs and associated nonlocal nonlinear equations. The dual Lax pair with quantized spectral parameter is introduced in Sec.II. The connections between the spaces of solutions of direct and dual Lax pairs are found. The reduction constraint on the Lax pairs coefficients that leads to nonlocal $`A_{n1}`$ Toda lattice $`(n𝐍)`$ is discussed in Sec.III. This nonlinear system is written here in the bilinear form, which is suitable for exploiting the usual perturbation theory and whose one–soliton expansion, as it is shown in the next section, has infinitely many exponential terms. In Sec.IV the theorem establishing the covariance of the Lax pairs with quantized spectral parameter with respect to the DT of direct pair is presented. For the case of the Toda lattice, we formulate sufficient conditions to keep the reduction constraint on the Lax pair coefficients under performing this transformation. Similarly we construct the DT of dual pair preserving the reduction, whose product with the DT of direct pair yields the formulas of binary DT. Iterations of these transformations on zero background give multi–soliton solutions that depend on infinitely many arbitrary parameters. The formulas of infinitesimal Darboux transformation and the expansions in series of the Lax pair solutions at the singular points neighborhood are presented in Sec.V. A simple way to produce infinite hierarchies of the infinitesimal symmetries and conservation laws of nonlocal $`A_{n1}`$ Toda lattice is given there.

## 2 Lax Pairs with Quantized Spectral Parameter

Let us consider the overdetermined linear system

$$\{\begin{array}{c}\mathrm{\Psi }_x=JT\mathrm{\Psi }\mathrm{\Lambda }+U\mathrm{\Psi }\hfill \\ \mathrm{\Psi }_t=AT^1\mathrm{\Psi }\mathrm{\Lambda }^1\hfill \end{array}$$
(1)
for matrix $`nn`$ function $`\mathrm{\Psi }\mathrm{\Psi }(x,t,\mathrm{\Lambda })`$. Here $`\mathrm{\Lambda }=\text{diag}(\lambda _j)`$ is the constant matrix, $`T`$ is the translation operator:

$$T=\mathrm{exp}(h_x)$$

($`h`$ is a constant), matrices $`J`$, $`U`$ and $`A`$ are independent of $`\mathrm{\Lambda }`$. This system coincides with usual Lax pair studied in the frameworks of the IST method  if $`h=0`$. Matrix $`\mathrm{\Lambda }`$ in this case is called the matrix spectral parameter. Joint action of operator $`T`$ and matrix $`\mathrm{\Lambda }`$ on solution $`\mathrm{\Psi }`$ can be regarded as a notion of the quantized spectral parameter . We will quote Eqs.(1) as direct Lax pair in that follows. Compatibility condition $`\mathrm{\Psi }_{xt}=\mathrm{\Psi }_{tx}`$ of direct Lax pair gives system of nonlocal matrix equations

$$\{\begin{array}{c}J_t=0\hfill \\ U_t=JTAAT^1J\hfill \\ A_x=UAAT^1U\hfill \end{array}.$$
(2)
It is remarkable that Eqs.(2) are also derived from the compatibility condition of the overdetermined linear system that is ”dual” to Lax pair (1). This system (dual Lax pair) is written in following manner

$$\{\begin{array}{c}\mathrm{\Xi }_x=K\left(T^1\mathrm{\Xi }J\right)\mathrm{\Xi }U\hfill \\ \mathrm{\Xi }_t=K^1\left(T\mathrm{\Xi }A\right)\hfill \end{array}$$
(3)
with $`\mathrm{\Xi }\mathrm{\Xi }(x,t,K)`$ and $`K=\text{diag}(\ae _j)`$ being respectively matrix solution and matrix spectral parameter of dual system. The spaces of solutions of direct and dual Lax pairs are connected in the local case. Such the connection turns out to exist for the Lax pairs considered here. Namely, matrix function

$$R(\mathrm{\Xi },\mathrm{\Psi })=K_x^{x+h}\left(T^1\mathrm{\Xi }J\right)\mathrm{\Psi }𝑑x+\mathrm{\Xi }\mathrm{\Psi }$$

is independent on variables $`x`$ and $`t`$ if $`K=\mathrm{\Lambda }=\lambda E`$ ($`\lambda `$ is scalar spectral parameter). Another relation between the spaces of solutions is provided by the closure of differential 1–form

$$d\omega (\mathrm{\Xi },\mathrm{\Psi })=\mathrm{\Xi }JT\mathrm{\Psi }dx+K^1\left(T\mathrm{\Xi }A\right)\mathrm{\Psi }\mathrm{\Lambda }^1dt.$$
(4)
## 3 Nonlocal Two–Dimensional Toda Lattice

Let matrix $`J`$ be defined as given

$$J=\{\delta _{j,k1}\}$$
(5)
(the indices are supposed hereafter to be equal on modulo $`n`$). System (2) is valid if we put

$$U=\sigma _x,A=\mathrm{exp}(\sigma )J^1\mathrm{exp}(T^1\sigma ),$$
(6)
where

$$\sigma =\text{diag}(\sigma _j)$$
(7)
and functions $`\sigma _j\sigma _j(x,t)`$ ($`j=1,\mathrm{},n`$) are the solutions of nonlocal generalization of the two–dimensional periodic $`A_{n1}`$ Toda lattice

$$\sigma _{j,xt}=\mathrm{exp}\left(T\sigma _{j+1}\sigma _j\right)\mathrm{exp}\left(\sigma _jT^1\sigma _{j1}\right).$$
(8)
Changing variables

$$\sigma _j=\mathrm{log}\frac{T\tau _{j+1}}{\tau _j}$$

results Eqs.(8) in bilinear form:

$$D_xD_t\tau _j\tau _j+2\tau _j^2=2\left(T\tau _{j+1}\right)T^1\tau _{j1}$$
(9)
($`D_x`$ and $`D_t`$ denote the Hirota derivatives ). In the local case these equations give the bilinear representation of the two–dimensional Toda lattice (compare with Eqs.(4.3) in ).

The important feature of nonlocal direct and dual Lax pairs is that matrix $`U`$ of general form cannot be led to the algebra $`SL(n)`$ case by means of the gauge transformation . However, the coefficients of linear systems (1,3) admit new reduction constraint that is in keeping with system (8):

$$\sigma =\rho E,$$
(10)
where single dependent variable $`\rho `$ solves equation ($`A_0`$ Toda lattice)

$$\rho _{xt}=\left(T1\right)\mathrm{exp}\left(\rho T^1\rho \right).$$
(11)
It should be stressed that this equation, whose local analogue is trivial, appears as the compatibility condition of the overdetermined linear systems of arbitrary matrix dimension. Making the dependent variable transformation in Eq.(11)

$$\rho =(T1)\mathrm{log}\tau ,$$

one obtains bilinear equation for $`\tau `$

$$D_xD_t\tau \tau +2\tau ^2=2\left(T\tau \right)T^1\tau .$$
(12)
Eq.(11) was derived in  as continious limit of the lattice equations describing the transfer of energy of plasma oscillations. The nonlocal equation that can be written in the form of Eq.(11) with operator $`T`$ defined in different manner was considered in work .

The soliton solutions of Eqs.(9) and Eq.(11) can be constructed applying the usual procedure of the Hirota method . However, the explicit form of the one–soliton solution is not obvious in this approach. To generate the hierarchy of solutions of the nonlocal Toda lattice we develop below another method.

## 4 Darboux Transformation Technique and Solitons

The underlying property of the DT technique is an existence of the kernel of transformations of the Lax pairs solutions for some value of the spectral parameter. This property is exploited in the present paper to generalize the technique considered for the nonlocal equations. For the sake of convenience we introduce notation

$$\mathrm{\Omega }(\mathrm{\Xi },\mathrm{\Psi })=_{(x_0,t_0)}^{(x,t)}𝑑\omega (\mathrm{\Xi },\mathrm{\Psi })+C(\mathrm{\Xi },\mathrm{\Psi }),$$

where it is supposed that constant matrix $`C(\mathrm{\Xi },\mathrm{\Psi })`$ can be determined from equation

$$KC(\mathrm{\Xi },\mathrm{\Psi })C(\mathrm{\Xi },\mathrm{\Psi })\mathrm{\Lambda }=R(\mathrm{\Xi },\mathrm{\Psi })|_{(x_0,t_0)}.$$

THEOREM Let $`\mathrm{\Phi }`$ be a solution of direct Lax pair (1) with matrix spectral parameter $`M`$. Then matrix

$$\stackrel{~}{\mathrm{\Psi }}=\mathrm{\Psi }_x\mathrm{\Phi }_x\mathrm{\Phi }^1\mathrm{\Psi }$$
(13)
and, if there exists appropriate matrix $`C(\mathrm{\Xi },\mathrm{\Phi })`$, matrix

$$\stackrel{~}{\mathrm{\Xi }}=\mathrm{\Omega }(\mathrm{\Xi },\mathrm{\Phi })\left(T\mathrm{\Phi }^1\right)J^1$$
(14)
are solutions respectively of direct and dual Lax pairs

$$\{\begin{array}{c}\stackrel{~}{\mathrm{\Psi }}_x=JT\stackrel{~}{\mathrm{\Psi }}\mathrm{\Lambda }+\stackrel{~}{U}\stackrel{~}{\mathrm{\Psi }}\hfill \\ \stackrel{~}{\mathrm{\Psi }}_t=\stackrel{~}{A}T^1\stackrel{~}{\mathrm{\Psi }}\mathrm{\Lambda }^1\hfill \end{array},\{\begin{array}{c}\stackrel{~}{\mathrm{\Xi }}_x=K\left(T^1\stackrel{~}{\mathrm{\Xi }}J\right)\stackrel{~}{\mathrm{\Xi }}\stackrel{~}{U}\hfill \\ \stackrel{~}{\mathrm{\Xi }}_t=K^1\left(T\stackrel{~}{\mathrm{\Xi }}\stackrel{~}{A}\right)\hfill \end{array},$$

whose coefficients are

$$\stackrel{~}{U}=U+J_xJ^1+J\left(T\mathrm{\Phi }_x\mathrm{\Phi }^1\right)J^1\mathrm{\Phi }_x\mathrm{\Phi }^1,$$
(15)
$$\stackrel{~}{A}=J(T\mathrm{\Phi })M\mathrm{\Phi }^1A(T^1\mathrm{\Phi })M^1\mathrm{\Phi }^1(T^1J^1).$$
(16)
The direct Lax pair is covariant with respect to transformation $`\{\mathrm{\Psi },U,A\}\{\stackrel{~}{\mathrm{\Psi }},\stackrel{~}{U},\stackrel{~}{A}\}`$ owing to the existence of the kernel: $`\stackrel{~}{\mathrm{\Psi }}0`$ if $`\mathrm{\Psi }=\mathrm{\Phi }`$. The covariance of the dual Lax pair is proven by applying identity

$$K\left(T^1\mathrm{\Omega }(\mathrm{\Xi },\mathrm{\Phi })\right)=\mathrm{\Omega }(\mathrm{\Xi },\mathrm{\Phi })M+\mathrm{\Xi }\mathrm{\Phi }.$$

It follows from the compatibility conditions of transformed direct or dual Lax pairs that their coefficients $`J`$, $`\stackrel{~}{U}`$ and $`\stackrel{~}{A}`$ are new solutions of system (2). Transformation (1316) carried out with solution $`\mathrm{\Phi }`$ of Eq.(1) is called the DT of direct pair. Reduction constraint (6) on the Lax pairs coefficients is kept under performing this DT by imposing the following conditions on solution $`\mathrm{\Phi }`$ and its spectral parameter $`M`$:

$$\mathrm{\Phi }=\text{diag}(\phi _j)B,$$
(17)
$$B_{jk}=\mathrm{exp}(2\pi i(j1)(k1)/n)$$

$$(j,k=1,\mathrm{},n),$$

$$M=\mu \text{diag}(\mathrm{exp}(2\pi i(j1)/n)).$$
(18)
Here $`\phi \phi (x,t,\mu )=(\phi _1,\mathrm{},\phi _n)^T`$ is a vector solution of Eq.(1) with scalar spectral parameter $`\mu `$. If $`\psi =(\psi _1,\mathrm{},\psi _n)^T`$ and $`\xi =(\xi _1,\mathrm{},\xi _n)`$ are the vector solutions of direct and dual Lax pairs with scalar spectral parameters $`\lambda `$ and $`\ae `$ respectively, then Eqs.(13,14) under conditions (17,18) give the next expressions for the components of transformed vector solutions

$$\stackrel{~}{\psi }_j=\psi _{j,x}\frac{\phi _{j,x}}{\phi _j}\psi _j,$$
(19)
$$\stackrel{~}{\xi }_j=\frac{\mathrm{\Omega }_{j+1}(\xi ,\phi )}{T\phi _{j+1}}$$
(20)
$`(j=1,\mathrm{},n)`$ with

$$\mathrm{\Omega }_j(\xi ,\phi )=_{(x_0,t_0)}^{(x,t)}\xi _{j1}T\phi _j𝑑x+\mathrm{exp}\left(T\sigma _j\sigma _{j1}\right)\frac{\left(T\xi _j\right)\phi _{j1}}{\ae \mu }dt+\frac{\rho _j(\xi ,\phi )}{\ae ^n\mu ^n}|_{(x_0,t_0)},$$

$$\rho _j(\xi ,\phi )=\underset{k=1}{\overset{n}{}}\ae ^{k1}\mu ^{nk}\left(\ae _x^{x+h}\left(T^1\xi _{jk1}\right)\phi _{jk}𝑑x+\xi _{jk}\phi _{jk}\right).$$

Transformation of solutions of nonlocal Toda lattice (8), which corresponds to transformations (19,20) of solutions of associated Lax pairs, has the form

$$\stackrel{~}{\sigma }_j=\sigma _j+\mathrm{log}\frac{T\phi _{j+1}}{\phi _j}(j=1,\mathrm{},n).$$
(21)
Similarly one can construct the formulas of DT of dual pair, using the matrix solution of Eq.(3). It can be shown that the transformations of direct and dual pairs commute. The conditions analogous to Eqs.(17,18) are imposed on the matrix solution and its spectral parameter for reduction constraint (6) to be inherited under carrying out DT of dual pair. Corresponding formulas for the transformations of the Lax pairs solutions and ones of nonlocal Toda lattice are

$$\stackrel{~}{\psi }_j=T^1\frac{\mathrm{\Omega }_j(\chi ,\psi )}{\chi _{j1}},$$
(22)
$$\stackrel{~}{\xi }_j=\xi _{j,x}\frac{\chi _{j,x}}{\chi _j}\xi _j,$$
(23)
$$\stackrel{~}{\sigma }_j=\sigma _j+\mathrm{log}\frac{\chi _j}{T^1\chi _{j1}}$$
(24)
$`(j=1,\mathrm{},n)`$, where $`\chi =(\chi _1,\mathrm{},\chi _n)`$ is the vector solution of dual Lax pair with scalar spectral parameter $`\nu `$. The product of transformations (1921) and (2224) yields the formulas of so–called binary DT

$$\stackrel{~}{\psi }_j=\psi _j\frac{T^1\mathrm{\Omega }_j(\chi ,\psi )}{T^1\mathrm{\Omega }_j(\chi ,\phi )}\phi _j,$$
(25)
$$\stackrel{~}{\xi }_j=\xi _j\frac{\mathrm{\Omega }_{j+1}(\xi ,\phi )}{\mathrm{\Omega }_{j+1}(\chi ,\phi )}\chi _j,$$
(26)
$$\stackrel{~}{\sigma }_j=\sigma _j+\mathrm{log}\frac{\mathrm{\Omega }_{j+1}(\chi ,\phi )}{T^1\mathrm{\Omega }_j(\chi ,\phi )}$$
(27)
$`(j=1,\mathrm{},n)`$.

To keep reduction constraint (10) we have to impose additional conditions on the solutions of Lax pairs, which are used in performing DT’s. For example, if we put $`\phi _j=\alpha ^j\vartheta `$ in transformation (1921), where $`\alpha =\text{exp}(2\pi ik/n)`$ $`(k=0,\mathrm{},n1)`$ and function $`\vartheta `$ solves system

$$\{\begin{array}{c}\vartheta _x=\alpha \mu T\vartheta +\rho _x\vartheta \hfill \\ \vartheta _t=\frac{\text{exp}(\rho T^1\rho )}{\alpha \mu }T^1\vartheta \hfill \end{array},$$

then transformed solution of Eq.(11) is

$$\stackrel{~}{\rho }=\rho +(T1)\mathrm{log}\vartheta .$$

Let us consider the zero background $`(\sigma =0)`$. Vectors $`\phi `$ and $`\chi `$, which were exploited in constructing the DT’s, satisfy nonlocal linear systems

$$\{\begin{array}{c}\phi _x=\mu JT\phi \hfill \\ \phi _t=\mu ^1J^1T^1\phi \hfill \end{array},$$
(28)
$$\{\begin{array}{c}\chi _x=\nu T^1\chi J\hfill \\ \chi _t=\nu ^1T\chi J^1\hfill \end{array}.$$
(29)
These systems have solutions depending on infinitely many parameters of the next form

$$\phi _j=\underset{k=0}{\overset{n1}{}}\underset{m_k}{}c_k^{(m_k)}\mathrm{exp}\left(2\pi i(j1)k/n\right)\mathrm{exp}\left(p_k^{(m_k)}xt/p_k^{(m_k)}\right),$$
(30)
$$\chi _j=\underset{k=0}{\overset{n1}{}}\underset{n_k}{}d_k^{(n_k)}\mathrm{exp}\left(2\pi i(j1)k/n\right)\mathrm{exp}\left(q_k^{(n_k)}xt/q_k^{(n_k)}\right),$$
(31)
$`(j=1,\mathrm{},n)`$. Here $`p_k^{(m_k)}`$ and $`q_k^{(n_k)}`$ satisfy equations

$$p_k^{(m_k)}+\mu \mathrm{exp}(2\pi ik/n)\mathrm{exp}(p_k^{(m_k)}h)=0,$$

$$q_k^{(n_k)}\nu \mathrm{exp}(2\pi ik/n)\mathrm{exp}(q_k^{(n_k)}h)=0$$

$`(k=0,\mathrm{},n1)`$, $`c_k^{(m_k)}`$ and $`d_k^{(n_k)}`$ are arbitrary constants. Carrying out of DT accordingly to Eqs.(1921) (or Eqs.(2224)) gives the one–soliton solutions, which depend on infinitely many free parameters. It means that in general the set $`\tau _j`$ $`(j=1,\mathrm{},n)`$ of the one–soliton $`\tau `$–functions of Eqs.(9) and the one–soliton $`\tau `$–functions of Eq.(12) contain infinitely many exponential terms. The products of these exponents appear no in the expression of the multi–soliton solutions of the bilinear equations. One–soliton solutions in nonlocal case form two families to be produced by DT’s of direct and dual pairs respectively, while in local case transformation (2224) can be obtained iterating transformation (1921).

In the case $`n=2`$, the components of solution of system (28) are represented in the following manner

$$\phi _j=\underset{p_+}{}c_{p_+}\mathrm{exp}\left(p_+xt/p_+\right)+(1)^{j1}\underset{p_{}}{}c_p_{}\mathrm{exp}\left(p_{}xt/p_{}\right)$$

$`(j=1,2)`$, where $`p_\pm `$ are the roots of equations

$$p\pm \mu \mathrm{exp}(ph)=0,$$

and the summations are supposed over all solutions of the last equations. Substitution of vector solution $`\phi `$ such, that $`\phi _1=\pm \phi _2`$, into Eq.(21) leads to the solutions of $`A_0`$ Toda lattice (11). For real $`h`$, simplest real nonsingular solution is

$$\rho =(T1)\mathrm{log}\vartheta ,$$

where

$$\vartheta =c_1\mathrm{exp}\left(p_1xt/p_1\right)+c_2\mathrm{exp}\left(p_2xt/p_2\right),$$

$`c_1`$, $`c_2`$ and $`\mu `$ are real constants, $`c_1c_2>0`$, $`p_1`$ and $`p_2`$ are real solutions of equation

$$p=\mu \mathrm{exp}(ph).$$

Iterations of transformations (1927) allow one to construct the infinite hierarchies of solutions of Eqs.(8) and corresponding solutions of Lax pairs (1,3). Final expressions for the transformed quantities are brought to the determinant form. Taking into account

$$\lambda \rho _{j+1}(\xi ,\psi )=\ae \rho _j(\xi ,\psi )+(\lambda ^n\ae ^n)\left(\ae _x^{x+h}\left(T^1\xi _{j1}\right)\psi _j𝑑x+\xi _j\psi _j\right),$$

$$\ae T^1\mathrm{\Omega }_j(\xi ,\psi )=\lambda \mathrm{\Omega }_{j+1}(\xi ,\psi )+\xi _j\psi _j,$$

one obtains

$$\rho _j(\stackrel{~}{\xi },\stackrel{~}{\psi })=(\ae ^n\lambda ^n)\frac{T\psi _j}{T\phi _j}\mathrm{\Omega }_j(\xi ,\phi )\rho _j(\xi ,\psi ),$$

$$\mathrm{\Omega }_j(\stackrel{~}{\xi },\stackrel{~}{\psi })=\frac{T\psi _j}{T\phi _j}\mathrm{\Omega }_j(\xi ,\phi )\mathrm{\Omega }_j(\xi ,\psi ),$$

where $`\stackrel{~}{\psi }`$ and $`\stackrel{~}{\xi }`$ are defined accordingly to Eqs.(19,20). The formulas of $`N`$-th iteration of the DT of direct pair are written by means of these equations as given

$$\stackrel{~}{\psi }_j=\frac{\mathrm{\Delta }^{(j)}[N+1]}{\mathrm{\Delta }^{(j)}[N]},$$

$$\stackrel{~}{\xi }_j=\left|\begin{array}{ccccc}T\phi _{j+1}^{(1)}& T\phi _{j+1,x}^{(1)}& \mathrm{}& T\phi _{j+1,(N2)x}^{(1)}& \mathrm{\Omega }_{j+1}(\xi ,\phi ^{(1)})\\ T\phi _{j+1}^{(2)}& T\phi _{j+1,x}^{(2)}& \mathrm{}& T\phi _{j+1,(N2)x}^{(2)}& \mathrm{\Omega }_{j+1}(\xi ,\phi ^{(2)})\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ T\phi _{j+1}^{(N)}& T\phi _{j+1,x}^{(N)}& \mathrm{}& T\phi _{j+1,(N2)x}^{(N)}& \mathrm{\Omega }_{j+1}(\xi ,\phi ^{(N)})\end{array}\right|/T\mathrm{\Delta }^{(j+1)}[N],$$

$$\stackrel{~}{\sigma }_j=\sigma _j+\mathrm{log}\frac{T\mathrm{\Delta }^{(j+1)}[N]}{\mathrm{\Delta }^{(j)}[N]},$$

where

$$\mathrm{\Delta }^{(j)}[N]=det[\phi _{j,(m1)x}^{(k)}],$$

$`\phi _j^{(k)}`$ are the components of the vector solutions of Eq.(1) with scalar spectral parameter $`\mu ^{(k)}`$, $`\phi ^{(N+1)}=\psi `$ $`(j=1,\mathrm{}n;k,m=1,\mathrm{},N)`$. The reduction constraint (10) is inherited under the iterations by supposing $`N=ln`$ $`(l=1,2,\mathrm{})`$, $`\phi ^{(ml+k)}=(\phi _{m+1}^{(k)}\mathrm{},\phi _{m+n}^{(k)})^T`$, $`\mu ^{(ml+k)}=\mu ^{(k)}`$ $`(m=0,\mathrm{},n1;k=1,\mathrm{},l)`$.

The iterations of the transformations presented in this section on zero background give the multi–soliton solutions of nonlocal $`A_{n1}`$ Toda lattice that depend on infinitely many arbitrary parameters. The interaction of the one–soliton components in the multi–soliton solution causes the shifts of the parameters of the one–soliton solutions.

## 5  Infinitesimal Symmetries and Conservation Laws

Taking limit $`\nu \mu `$ in Eqs.(2527) one finds solutions of the linearizations of direct and dual Lax pairs and nonlocal Toda lattice

$$\delta \psi _j=\mu \left(T^1\mathrm{\Omega }_j(\chi ,\psi )\right)\phi _j,$$
(32)
$$\delta \xi _j=\mu \mathrm{\Omega }_{j+1}(\xi ,\phi )\chi _j,$$
(33)
$$\delta \sigma _j=\chi _j\phi _j$$
(34)
($`j=1,\mathrm{},n`$). These formulas establish infinitesimal DT $`\psi \psi +\epsilon \delta \psi `$, $`\xi \xi +\epsilon \delta \xi `$, $`\sigma \sigma +\epsilon \delta \sigma `$ ($`\epsilon =\nu \mu `$). The closure of differential 1–form (4) immediately yields

$$(\mathrm{\Xi }JT\mathrm{\Psi })_t=(K^1\left(T\mathrm{\Xi }A\right)\mathrm{\Psi }\mathrm{\Lambda }^1)_x.$$

In terms of vector solutions $`\phi `$ and $`\chi `$ of Lax pairs with scalar spectral parameter $`\mu `$ this identity reads as

$$\left(\underset{j=1}{\overset{n}{}}\chi _{j1}T\phi _j\right)_t=\mu ^2\left(\underset{j=1}{\overset{n}{}}\left(T\chi _{j+1}\right)\phi _j\text{e}^{T\sigma _j\sigma _{j1}}\right)_x$$

or, equivalently,

$$T_t+X_x=0,$$
(35)
where we use notations

$$T=\underset{j=1}{\overset{n}{}}\chi _{j1}T\phi _j,$$
(36)
$$X=\mu ^2\underset{j=1}{\overset{n}{}}\left(T\chi _{j+1}\right)\phi _j\text{e}^{T\sigma _j\sigma _{j1}}.$$
(37)
The hierarchies of infinitesimal symmetries and conservation laws of nonlocal generalization of two–dimensional periodic $`A_{n1}`$ Toda lattice are obtained substituting in Eqs.(3437) the expansions of the Lax pairs solutions at the neighborhood of singular points on the spectral parameter plane. The components of vector solutions $`\phi `$ and $`\chi `$ of direct and dual Lax pairs (1,3), whose coefficients are defined by Eqs.(57), are represented in the neighborhood of point $`\mu =\mathrm{}`$ in the following manner

$$\phi _j=\left(1+\underset{k=1}{\overset{\mathrm{}}{}}\frac{A_j^{(k)}}{\mathrm{\Omega }^k}\right)\left(\frac{\mathrm{\Omega }}{\mu }\right)^{j1}\text{e}^{\mathrm{\Omega }\left(j1\right)h\mathrm{\Omega }x},$$
(38)
$$\chi _j=\left(1+\underset{k=1}{\overset{\mathrm{}}{}}\frac{B_j^{(k)}}{\mathrm{\Omega }^k}\right)\left(\frac{\mu }{\mathrm{\Omega }}\right)^{j1}\text{e}^{\mathrm{\Omega }\left(j1\right)h+\mathrm{\Omega }x}$$
(39)
($`j=1,\mathrm{},n`$), where $`\mathrm{\Omega }`$ solves equation

$$\mathrm{\Omega }^n\text{exp}(n\mathrm{\Omega }h)=\mu ^n$$

and coefficients $`A_j^{(k)}`$ and $`B_j^{(k)}`$ ($`k𝐍`$) satisfy equations

$$\{\begin{array}{c}A_j^{(k)}+A_{j,x_{}}^{(k1)}=TA_{j+1}^{(k)}+\sigma _{j,x}A_j^{(k1)}\hfill \\ A_{j,t}^{(k)}=\text{e}^{\sigma _jT^1\sigma _{j1}}T^1A_{j1}^{(k1)}\hfill \end{array},$$

$$\{\begin{array}{c}B_j^{(k)}+B_{j,x_{}}^{(k1)}=T^1B_{j1}^{(k)}\sigma _{j,x}B_j^{(k1)}\hfill \\ B_{j,t}^{(k)}=\text{e}^{T\sigma _{j+1}\sigma _j}TB_{j+1}^{(k1)}\hfill \end{array}.$$

For the neighborhood of singular point $`\mu =0`$ we have

$$\phi _j=\underset{k=0}{\overset{\mathrm{}}{}}C_j^{(k)}\alpha ^{1jk}\mu ^k\text{e}^{\alpha t/\mu },$$
(40)
$$\chi _j=\underset{k=0}{\overset{\mathrm{}}{}}D_j^{(k)}\alpha ^{jk1}\mu ^k\text{e}^{\alpha t/\mu }$$
(41)
($`j=1,\mathrm{},n`$). Here $`\alpha `$ is a root of equation

$$\alpha ^n=1,$$

$`C_j^{(k)}`$ and $`D_j^{(k)}`$ are solutions of the systems

$$\{\begin{array}{c}C_{j,x_{}}^{(k+1)}=TC_{j+1}^{(k)}+\sigma _{j,x}C_j^{(k+1)}\hfill \\ C_{j,t}^{(k)}+C_j^{(k+1)}=\text{e}^{\sigma _jT^1\sigma _{j1}}T^1C_{j1}^{(k+1)}\hfill \end{array},$$

$$\{\begin{array}{c}D_{j,x_{}}^{(k+1)}=T^1D_{j1}^{(k)}\sigma _{j,x}D_j^{(k+1)}\hfill \\ D_{j,t}^{(k)}D_j^{(k+1)}=\text{e}^{T\sigma _{j+1}\sigma _j}TD_{j+1}^{(k+1)}\hfill \end{array}.$$

One can check by straightforward calculations that the systems determining coefficients $`A_j^{(k)}`$, $`B_j^{(k)}`$, $`C_j^{(k)}`$ and $`D_j^{(k)}`$ are compatible. The formulas for the expansions in series at singular points neighborhood of the Lax pairs solutions in nonlocal case seem to be new. Note that in the case $`\mu =\mathrm{}`$ we have infinitely many types of asymptotic behavior. It follows after substituting expansions (38,39) and (40,41) in Eqs.(34,36,37) that the solution of the linearization of nonlocal Toda lattice, conserved densities and currents admit representation

$$\delta \sigma _j=\underset{k=0}{\overset{\mathrm{}}{}}\delta \sigma _j^{(k,\mathrm{})}\mathrm{\Omega }^k,\delta \sigma _j=\underset{k=0}{\overset{\mathrm{}}{}}\delta \sigma _j^{(k,0)}\mu ^k,$$

$$T=\underset{k=0}{\overset{\mathrm{}}{}}\underset{j=1}{\overset{n}{}}T_j^{(k,\mathrm{})}\mathrm{\Omega }^k,T=\underset{k=0}{\overset{\mathrm{}}{}}\underset{j=1}{\overset{n}{}}T_j^{(k,0)}\mu ^k,$$

$$X=\underset{k=0}{\overset{\mathrm{}}{}}\underset{j=1}{\overset{n}{}}X_j^{(k,\mathrm{})}\mathrm{\Omega }^k,X=\underset{k=0}{\overset{\mathrm{}}{}}\underset{j=1}{\overset{n}{}}X_j^{(k,0)}\mu ^k,$$

whose coefficients form the infinite hierarchies of infinitesimal symmetries and conservation laws. The first few nontrivial coefficients are

$$\delta \sigma _j^{(1,0)}=\sigma _{j,t},\delta \sigma _j^{(1,\mathrm{})}=\sigma _{j,x},\delta \sigma _j^{(2,\mathrm{})}=\sigma _{j,x}^2+\stackrel{t}{}\left(\text{e}^{T\sigma _{j+1}\sigma _j}+\text{e}^{\sigma _jT^1\sigma _{j1}}\right)_x𝑑t.$$

$$T_j^{(2,\mathrm{})}=\stackrel{t}{}\left(\text{e}^{T\sigma _j\sigma _{j1}}\stackrel{t}{}\left(\text{e}^{T^2\sigma _{j+1}T\sigma _j}+\text{e}^{\sigma _{j1}T^1\sigma _{j1}}\right)𝑑t\right)𝑑t\left(\stackrel{t}{}\text{e}^{T\sigma _j\sigma _{j1}}𝑑t\right)^2,$$

$$X_j^{(2,\mathrm{})}=\text{e}^{T\sigma _j\sigma _{j1}},$$

## 6 Conclusion

In this paper the multi–soliton solutions of nonlocal partner of two–dimensional Toda lattice were obtained. Unlike to local case, these solutions depend on infinitely many free parameters. This feature of multi–soliton solutions follows in the frameworks of the Darboux transformation technique from analogous property of the solutions of the Lax pairs of nonlocal Toda lattice. We presented also the formulas of the expansions in series on the spectral parameter powers of the solutions of the Lax pairs. These expansions were used to construct the hierarchies of infinitesimal symmetries and conservation laws of nonlocal two–dimensional Toda lattice.

## 7  Acknowledgements

I thank Dr. Heinz Steudel for useful discussions and hospitality. I thank also Gottlieb Daimler– und Karl Benz–Stiftung for financial support.


text:
# Discovery of Two Types of X-ray Dips in Cyg X-1

## 1 Introduction

Cyg X-1 was the first astronomical system discovered to show strong evidence for a stellar mass black hole (Bolton 1972; Webster $`\&`$ Murdin 1972). It is a binary system with an orbital period of 5.6 days. The results from radial velocity measurements imply that the mass of the compact object exceeds the upper limit on the mass of a neutron star, so the compact object is inferred to be a black hole. A more recent study, based on spectrum synthesis, has derived a mass of about 10 M for the black hole (Herrero et al. 1995). The companion star has been identified as an O9.7 Iab supergiant (Walborn 1973; Gies $`\&`$ Bolton 1986), with a mass of about 20 M (Herrero et al. 1995). Cyg X-1 is, therefore, intrinsically different from most known black hole candidates which contain only a low-mass companion star. This difference might be related to the fact that Cyg X-1 is a persistent X-ray source while those with a low-mass companion are exclusively transient sources.

The X-ray emission from Cyg X-1 is likely powered by accretion of material from the companion star by the black hole. In this case, the accretion flows are thought to follow a pattern that is intermediate between that of Roche-lobe overflow, as in low-mass systems, and that of wind accretion, which is common for high-mass systems (Gies $`\&`$ Bolton 1986). Such an accretion process is sometimes referred to as “focused wind accretion”, which can occur when the companion star is close to filling its Roche lobe (Friend & Castor 1982). The X-ray observations of Cyg X-1 have revealed that it has two distinct spectral states, hard and soft (see reviews by Oda 1977, Liang $`\&`$ Nolan 1984, Tanaka & Lewin 1995, Cui 1998, and Liang 1998). The source is usually found in the hard state, when its X-ray spectrum is relatively flat (with a typical power-law photon index of 1.5), but, occasionally, it makes a transition to the soft state, when the spectrum steepen significantly (to a typical photon index of 2.5).

In the hard state, the X-ray intensity of Cyg X-1 is strongly modulated by the binary motion (Wen et al. 1999; Brocksopp et al. 1999; Priedhorsky, Brandt, $`\&`$ Lund 1995; Holt et al. 1979), probably caused by varying amount of absorbing material along the line of sight through the stellar wind of the companion star. As expected, the intensity is found to reach a minimum at the times of superior conjunction of the black hole (when the companion star is in front of the black hole with respect to the line of sight), which is usually defined as orbital phase $`0`$. It is, however, a very broad minimum, spanning about 26$`\%`$ of the orbit. In contrast, no X-ray orbital modulation is measurable in the soft state (Wen et al. 1999), which might be due to a significant change in the physical conditions of the wind (e.g., much higher ionization level) and/or in the geometry of the X-ray emitting region.

Besides the global orbital modulation of the X-ray emission, X-ray intensity dips are often seen in the light curves of Cyg X-1 (e.g., Pravdo et al. 1980; Remillard & Canizares 1984; Kitamoto et al. 1984; Bałucińska & Hasinger 1991; Ebisawa et al. 1996). The dips vary in duration from minutes to hours. The spectrum of the source hardens during a dip, implying that the dips probably originate in the photoelectric absorption of X-rays by the ambient medium. This is further supported by the detection of Fe K absorption edge in the dip spectrum in some cases (Kitamoto et al. 1984). The absorbing column density can increase by more than one order of magnitude during a dip (Kitamoto et al. 1984; Bałucińska-Church et al. 1997). The distribution of dips against orbital phases shows a prominent peak around phase $`0`$ (Bałucińska-Church et al. 2000), implying an origin of the dips in the stellar wind from the companion star.

In this paper, we present results from a long observation of Cyg X-1 with the Rossi X-Ray Timing Explorer (RXTE) over one complete 5.6-day orbital cycle. Many intensity dips were detected. The large collecting area of the Proportional Counter Array (PCA) aboard RXTE made it possible to conduct a more detailed investigation of spectral evolution of the source during some of the strong dips, as well as to carry out orbital-phase-resolved spectroscopy to quantify the variation in the column density over an orbital cycle.

## 2 Data

For six consecutive days (2000 January 5–11) Cyg X-1 was observed with the large-area detectors abroad RXTE, with usual interruptions because of earth occultation and passage of the satellite through the South Atlantic Anomaly. The data were collected for an effective exposure time of about 240 ks. Judging from the long-term ASM/RXTE light curves of the source (both the flux and hardness ratios), Cyg X-1 was in its usual hard state when our pointed RXTE observations were made.

Besides the PCA, RXTE also carries another pointing instrument, the High-Energy X-ray Timing Experiment, which covers a nominal energy range of 15–250 keV. For this investigation, only the PCA data (2–60 keV) were used, because we are primarily interested in spectral changes at low energies. The PCA consists of five Proportional Counter Units (PCUs). Not all PCUs were turned on throughout the observation, because of technical difficulties with some of the PCUs. Consequently, the number of PCUs used varied between 3 and 5. Because of this, it was necessary to break up the observation into a number of segments and to analyze data for each segment separately.

## 3 Analyses and Results

Both spectral and timing analyses were carried out using the Standard 2 data (which has a timing resolution of 16 seconds). We used ftools v5.0, as well as the calibration files and background models that accompanied this release of the software, to extract X-ray spectra and light curves of Cyg X-1. The quality of the background models was examined by comparing an overall spectrum (source plus background) with the corresponding model background spectrum. The two should overlap each other at the highest energy channels, since the counts in those channels should be entirely due to the background. We found that background modeling was satisfactory, except for a couple of cases where renormalization of the background was necessary.

### 3.1 X-ray Dips

We made a light curve for each of the following energy bands: 3–6 keV, 6–15 keV, and 15–60 keV, using data from all PCUs (and all xenon layers) that were turned on for a segment of the observation. A corresponding background light curve was then constructed from the models and was subtracted from the overall light curve to obtain the light curve of the source for a particular energy band. The light curve of each segment was then normalized to count rates per PCU, and all segments were pieced together to obtain the final light curve for the entire observation. We also computed hardness ratios by taking the ratio of the light curves in two energy bands (higher energy band over lower energy band). The hardness ratios are useful in providing crude information about the spectral properties of a source.

The light curves of Cyg X-1 reveal frequent occurrences of a brief decrease in the intensity of the source. These are X-ray dips. In most cases, the hardness ratios indicate that the spectrum of Cyg X-1 hardens during a dip, implying that the dips are caused by additional absorption. As an example, Figure. 1 shows a segment of the light curves where two dips are clearly visible (on top of an X-ray flare). As a matter of fact, the dips manifest themselves much more prominently as peaks in the hardness ratio plots, which are also shown in the figure. Note that we have plotted orbital phases instead of times to emphasize the distribution of dips over the orbital cycle. For this purpose, we have used the most updated ephemeris of Cyg X-1 (Brocksopp et al. 1999). The range of orbital phases covered by our observations is between 0.50 and 1.60.

To obtain a complete sample of dips, we used the 3–6 keV light curve alone, because it is most sensitive to absorption dips and the data is of the highest quality (compared to light curves in higher energy bands). We adopted the following algorithm to detect dips. First, we smooth the light curve by a running boxcar with a width of 3 time bins (or 48 s), to minimize intrinsic variability of the source. Second, we scan over the entire light curve for all “local minima”. A local minimum is defined as a point with two higher adjacent neighbors. Third, we compare each local minimum with its next closest neighbor on each side. This process continues until two points on each side are found to be more than 20$`\sigma `$ higher than the local minimum itself (where $`\sigma `$ is the error bar of the local minimum). A dip is then identified. On the other hand, the process stops, if one neighbor on either side is lower than the local minimum or if the time interval examined exceeds 300 time bins (or 4.8 ks) or reaches the overall length of a data segment, before the 20$`\sigma `$ threshold is reached. We then move onto the next local minimum, until all minima in the list have been examined. The choice of a relatively high detection threshold is to minimize confusion due to the intrinsic variability of the source. Since the main objective of this investigation is to study the spectral properties of dips in a more detailed manner than what has been done previously (rather than focusing on the frequency of the dips; see (Bałucińska-Church et al. 2000), the reliability of detection is much more important than the number of dips that we can detect. We note that our algorithm is clearly biased against dips of very short durations (less than 80 s) or those of very long durations (greater than 4.8 ks or the duration of a data segment). However, it is very difficult to improve the situation, because the detection of short-duration dips is limited by the intrinsic variability of the source and the detection of long-duration ones is complicated by coverage gaps in the data.

We detected a total of 33 dips. To examine their distribution over the binary orbit, we divided the orbit into 20 phase bins and computed the normalized number of dips in each phase bin (which is the number of dips detected divided by the fractional coverage of the phase bin). The normalization is necessary because there are gaps in the light curve and thus the coverage is incomplete for all phase bins (note that each phase bin corresponds to about 0.28 days or 6.7 hours). The results are shown in Figure 2. The overall shape of the distribution, such as the apparent asymmetry about phase $`0`$, is in general agreement with that derived previously (Bałucińska-Church et al. 2000), although the statistics are very limited here. The asymmetry in the distribution might be a manifestation of enhanced “clump” formation in the leading hemisphere of the stellar wind, because of the compression of the wind by orbital motion.

The duration of the dips varies greatly. Figure 3 shows a dip that has the longest duration among all detected. It should be noted that this dip was not detected by the algorithm, because of its long duration as well as the data gaps present. It is interesting that it occurred almost exactly at phase $`0`$. Unfortunately, the presence of gaps in the data makes it impossible to derive its duration accurately. The profile of a dip can be quite complex, certainly not necessarily symmetric. Figure 4 shows an example of such a dip. Sub-structures within the dip are quite apparent (although it can also be argued that it is a superposition of two smaller dips).

Upon a closer examination of the dips detected, we found that they fell into two distinct categories, based on their spectral properties. While the majority of the dips show strong spectral variation over the duration of dipping activity (see the hardness ratio plots in Figures 1, 3, and 4), there are some dips that show little spectral variation. We refer to the former as type A and the latter as type B. As examples, Figures 5 and 6 show the light curve and hardness ratios of two type B dips of quite different durations. It is clear that type B dips can be comparable in depth to type A ones (e.g., comparing Figure 1 and Figure 6) but they show much less spectral variation (if any at all). Of the 33 dips detected, 4 dips are quite convincingly type B. They seem to distribute randomly over the orbit: one is around phase 0.637, one around phase 0.830, one around phase 0.995, and one around 1.148, but we do not have sufficient statistics to address this issue with confidence.

### 3.2 Dip Spectra

We proceeded to carry out more detailed spectral analyses and modeling for some of the strong dips. We divided the overall light curve into dip and non-dip time intervals (with gaps between the two to ensure the quality of the dip spectra). We then extracted a source spectrum for each time interval by combining data from all PCUs (and all xenon layers) that were on during the observation. To strengthen some of the conclusions drawn from the hardness ratios in a model-independent manner, we simply divided the spectrum of the source during a dip by the average non-dip spectrum. For illustration, Figure 7 shows the results for a type A dip and a type B dip, respectively.

We can see that the spectrum of the type A dip rolls over much faster than the non-dip spectrum below about 6 keV, indicating the need for additional absorbing material intrinsic to the binary system. It is interesting to notice that the dip spectrum never fully recovers to the non-dip spectrum even at high energies. This seems to signify the importance of Thompson scattering in attenuating X-ray flux (at a few percent level). The presence of sufficient number of electrons to participate in the scattering process would imply that the absorbing medium is perhaps highly ionized, which has been suggested previously (Pravdo et al. 1980; however, see Kitamoto et al. 1984; Bałucińska-Church et al. 1997).

To be more quantitative, we fit the spectra with an empirical model that was used in previous studies of Cyg X-1 (Ebisawa et al. 1996; Cui et al. 1997). First, we fit the non-dip spectrum with a broken power law (plus a Gaussian function, which is needed to minimize residuals between 6–7 keV). A detailed discussion of physical models is beyond the scope of this paper. Here we only focus our attention on how the column density varies in a dip. The model describes the spectrum quite well. The inferred column density is $`0.8\pm 0.4\times 10^{22}\text{ }`$cm-2, which is in general agreement with previous measurements (e.g., Bałucińska & Hasinger 1991; Ebisawa et al. 1986; Cui et al. 1998). We then applied the best-fit model for the non-dip spectrum to the dip spectrum by fixing all model parameters except for the column density. This procedure is necessary because the quality of the dip spectrum is not adequate for constraining so many parameters in the model, but there is no physical justification for doing so. This only represents an attempt to quantify the amount material in the binary system that caused an absorption dip. The model fits the dip spectrum quite well. The inferred column density for the dip is $`1.7\pm 0.1\times 10^{22}\text{ }`$cm-2, which roughly doubles the interstellar value. We note that this is only the average column density for the dip; the peak value can be much higher.

The situation for type B dips is quite different. For the type B dip shown in Figure 5, the ratio of the dip and non-dip spectra is nearly flat up to at least 20 keV, above which the statistics become quite poor. If we average all data points above 20 keV, the average level seems to lie above that for data points below 20 keV. The lack of a strong energy dependence of the ratio at low energies indicates that the observed attenuation cannot be due to photoelectric absorption. There are at least three possible scenarios that might explain the type B dips. The first scenario attributes the type B dips to partial covering of an extended X-ray emitting region by material that is totally opaque. The second one is that the source experiences a momentary decrease in its luminosity during a type B dip, due, for instance, to a sudden decrease in the mass accretion rate. The third possibility is that the type B dips originate in pure Thompson scattering, implying that the clumps must be almost fully ionized.

### 3.3 Orbital Variation in Column Density

To investigate the global orbital modulation of X-ray emission previously seen (e.g., Wen et al. 1999), we binned the observation into 10 orbital phases. For each orbital phase, we extracted a spectrum of the source, as before, but used data only from the first xenon layer of each PCU to minimize calibration uncertainties, since it is the most accurately calibrated among all xenon layers. The trade-off is that we lost spectral coverage at energies above about 25 keV. We can afford such a loss here, because we are only interested in spectral variation at the lowest energies due to photoelectric absorption.

For each orbital phase, we fit the spectrum with the same model described in the previous section (broken power law plus Gaussian function) and the results are satisfactory. Figure 8 shows the inferred column density at different orbital phases. Although the column density seems to be higher near superior conjunction, large uncertainties prevent us from drawing a definitive conclusion, especially in the presence of other high points, e.g., at phase 1.3, which may or may not be real. The large uncertainty in determining the column density arises mostly from the fact that the sensitivity of the PCA drops precipitously below about 3 keV. Bałucińska-Church et al. (2000) computed variation in column density through the stellar wind due to orbital motion, based on a simple wind geometry, for two extreme cases with and without the suppression of radiatively-driven wind. Our results seem to lie between those two cases (comparing Figure 8 with Figure 7 of Bałucińska-Church et al. 2000).

## 4 Discussion

As mentioned in the introduction, the phenomenon of X-ray intensity dips in Cyg X-1 is well known and well studied. Over the years, a rich database has been built up for this famous black hole candidate. The database is essential for a systematic study of dips, whose occurrence is somewhat unpredictable and is thus easy to miss in a brief pointed observation. The distribution of the dips over a binary orbital cycle has been reliably established recently (Bałucińska-Church et al. 2000): most of the dips occur near superior conjunction of the X-ray source, but dipping activity has also been seen at orbital phases far from superior conjunction. Our results are consistent with that. We went a step further in this study. We carried out a long observation of Cyg X-1 over one entire orbital period with the large-area detectors aboard RXTE. The data allowed us to conduct more detailed spectral analyses of the dipping phenomenon as well as phase-resolved spectroscopy to quantify the variation in the column density over an orbit.

The main result of this investigation is the recognition of the existence of two types of dips. The spectrum of type A dips shows an energy-dependent reduction at low energies that is characteristic of photoelectric absorption. Detailed spectral modeling shows that the only difference between the spectrum of type A dips and the average non-dip spectrum seems to be the presence of additional column density during a dip. We stress, however, that our data is not of sufficient quality to warrant more complicated modeling, as has been done previously (e.g., Ebisawa et al. 1996). The type A dips seem to occur preferentially around the times of superior conjunction of the black hole. The spectrum of type B dips, on the other hand, shows almost an energy-independent reduction at low energies. Therefore, these dips cannot possibly be due to photoelectric absorption. Moreover, the type B dips appear to distribute randomly over the binary orbit.

Type A dips are much more common than type B dips. They are probably produced by density enhancement in an inhomogeneous wind from the companion star, since the column density changes only moderately during such a dip. They occur more frequently around superior conjunction probably because the line-of-sight follows a longer path through the denser part of the wind (however, see Blondin & Woo 1995). The inhomogeneities in the wind can be caused by a variety of physical processes (see Bałucińska-Church et al. 2000 for a few examples). It remains to be seen whether the observed global orbital modulation of X-ray emission is entirely due to the presence of such dips.

Type B dips might be caused by partial covering of an extended X-ray emission region by an opaque “screen”. In the context of Comptonization models (e.g., review by Liang 1998 and references therein), the presence of an extended emitting region may not be unreasonable for Cyg X-1. The question is what could physically serve as an opaque “screen” to produce the observed type B dips. The answer may lie in the fact that Cyg X-1 is a high-mass X-ray binary (HMXB), in which the mass accretion process is probably mediated by stellar wind from the massive companion star. Numerical simulations of mass accretion in HMXBs show that a tidal stream develops when the companion star is close to filling its Roche lobe and that the stream trails behind the compact object due to the Coriolis force (Blondin et al. 1991). The density in the stream is shown to be as high as 20–30 times the ambient wind density. Since the obscuration of X-rays occurs close to the accretion disk in this case, it is not expected to be very sensitive to the binary motion. Now, the question is whether this scenario can be applied to Cyg X-1, whose orbit is only moderately inclined, with a probable inclination angle of 30°– 40° (Bolton 1975; Guinan et al. 1979; Daniel 1981). A full three-dimensional simulation would be required to show the scale height of the tidal streams. We note that the inclination angle remains to be poorly determined for Cyg X-1, which is one of the largest contributors to the uncertainty in determining the mass of the compact object. For instance, the optical polarization measurements imply a wide range of 25°– 70°(e.g., Long et al. 1980).

Type B dips might also be due to a sudden decrease in the mass accretion rate, causing a decrease in the X-ray luminosity of the source itself. This could, in principle, occur, since, after all, the intrinsic variability of the source is generally thought to be due to fluctuation in the mass accretion rate. However, for Cyg X-1 (and black hole candidates in general), a change in X-ray flux is generally accompanied by a change in spectral hardness. In fact, it has been shown that in the hard state the flux and spectral hardness are, to some degree, anti-correlated (e.g., Wen et al. 2001). However, the correlation is quite weak, so our results cannot rule out this possibility.

It is also possible that type B dips are produced by (nearly) pure Thomson scattering of X-rays in highly ionized dense clumps in the wind. Given that the cross section for photoelectric absorption is, on average, more than two orders of magnitude larger than the Thomson cross section, the density of a clump that produces a type B dip is required to be more than two orders of magnitude larger than that of a clump that produces a type A dip of comparable depth (assuming that the clumps are comparable in size). To abtain an ionization parameter that is, for instance, a factor of 100 larger for type B clumps than for type A clumps, we would then require that the former are, on average, more than 100 times closer to the X-ray source. If so, this scenario could also explain why the type B dips are much less dependent of orbital phases.

Finally, we made an attempt to quantify how the column density varies with orbital motion by carrying out orbital-phase-resolved spectroscopy. The results seem to indicate that the column density is higher at superior conjunction, although error bars are too large for us to be definitive about it. The issue can be resolved by a similar observation with detectors of much improved low-energy response and spectral resolution, such as those aboard Chandra and XMM. Much improved spectral resolution of these detectors will also allow the detection of absorption edges, which provides means of determining column density in a manner that is independent of continuum modeling (Schulz et al. 2001). Furthermore, the highly eccentric orbits of these satellites allow long un-interrupted observations, which is very much important for studying the distribution of dips over a binary orbit.

We wish to thank Dr. Douglas Gies, the referee, for many helpful comments. In particular, he suggested that type B dips might be caused by a change in the luminosity of the source and that the asymmetry in the dip distribution might be due to the compression of the stellar wind by orbital motion. This work was supported in part by NASA through grants NAG5-9098 and NAG5-9998.