format: text

text:
# Bright CO ro-vibrational emission lines in the class I source GSS 30 IRS1

## 1 Introduction

The innermost $`50\mathrm{AU}`$ of the circumstellar environments around embedded young stellar objects is poorly constrained observationally due to the large extinction through the embedding material ($`A_V>20\mathrm{mag}`$) and the small angular size ($`<1\mathrm{}`$) of the region. The usual molecular probes in the millimeter-submillimeter region are not effective for the temperatures ($`T>200\mathrm{K}`$) and densities ($`n_{\mathrm{H}_2}>10^7\mathrm{cm}^3`$) thought to be present. However, an understanding of the processes taking place in this regime is essential to obtain a complete picture of the process of accretion and the driving of outflows from low mass protostars as well as the early chemical and physical evolution of circumstellar disks.

One of the most effective probes of warm dense gas is through emission in molecular ro-vibrational transitions. The most common bands readily available from ground-based facilities are the CO overtones and $`\mathrm{H}_2`$ fundamental bands around $`2.2\mu \mathrm{m}`$ (e.g. Reipurth & Aspin 1997) and the fundamental transitions of CO in the $`M`$-band around $`4.7\mu \mathrm{m}`$. Also emission from very hot water gas ($`2000\mathrm{K}`$) near $`2.29\mu \mathrm{m}`$ has been reported toward a few sources (e.g. Najita et al. 2000). Since the upper levels of these transitions lie at temperatures of up to a few thousand Kelvin, they probe hot gas with temperatures between 100 and $`1000\mathrm{K}`$, making them ideal to study the region of interaction between disk, protostar, accretion and outflow. With the new generation of sensitive ground-based high resolution spectrometers for the $`35\mu \mathrm{m}`$ region an efficient window has been opened for the detailed study of the CO ro-vibrational lines toward embedded young stellar objects (YSO). The past generation of instruments was suitable to either observe a low resolution spectrum with a fairly wide spectral range (e.g. Teixeira et al. 1998) or a high resolution echelle spectrum with a very narrow spectral range (e.g. Carr et al. 2001). The main exception are the Fourier Transform Spectroscopy (FTS) observations of Mitchell et al. (1988, 1990), who observed a number of high-mass stars in the entire M band at high spectral resolution ($`R>10^6`$), but such studies are limited to the brightest objects. VLT-ISAAC has a large instantaneous spectral range in the $`M`$-band ($`0.237\mu m`$), which combined with a medium resolution of $`R=\mathrm{5\hspace{0.17em}000}\mathrm{10\hspace{0.17em}000}`$ and a limiting magnitude of $`M10`$ allows the entire fundamental band of gaseous CO of a wide range of young stars to be observed in a short time, including low mass stars down to a few tenths of a solar mass in the nearest star-forming clouds.

The embedded stars studied so far in CO ro-vibrational bands have showed mostly lines in absorption (Mitchell et al. 1990; Boogert et al. 2000, 2002) implying that the warm CO gas is seen in front of a bright infrared continuum, produced by hot dust close to the central object. The fundamental CO lines are usually only seen in emission towards more evolved sources characterized by a class II type spectrum where a circumstellar disk is directly visible (Carr et al. 2001; Blake et al. 2002). CO overtone bandhead emission at $`2\mu \mathrm{m}`$ has been observed in emission toward a few intermediate mass pre-main sequence stars (Thompson 1989; Najita et al. 1996) and T Tauri stars (Carr 1989), and has been associated with hot gas ($`T\mathrm{1\hspace{0.17em}500}\mathrm{5\hspace{0.17em}000}\mathrm{K}`$) located within a fraction of an AU in a Keplerian disk.

We present here the peculiar $`4.54.8\mu \mathrm{m}`$ spectrum of the illuminating source IRS1=Elias 21 of the bipolar reflection nebula GSS 30 located in the core of the $`\rho `$ Ophiuchus molecular cloud at a distance of $`160\mathrm{pc}`$. It has been classified as a low mass class I YSO from its spectral energy distribution (SED) (Grasdalen et al. 1973; Elias 1978; Wilking et al. 1989) and low bolometric luminosity ($`\mathrm{L}_{\mathrm{bol}}=2126\mathrm{L}_{\mathrm{}}`$, Greene et al. 1994; Bontemps et al. 2001). Extensive polarimetric studies in the H and K band of the reflection nebula have shown that the source is surrounded by a large ($`\mathrm{2\hspace{0.17em}000}\mathrm{AU}`$) disk-like envelope and a smaller circumstellar disk of $`150\mathrm{AU}`$, which are inclined about $`25\mathrm{°}`$ away from the plane of the sky, i.e. close to edge-on (Chrysostomou et al. 1996, 1997). The high degree of linear polarization (up to 50%) as well as the presence of circular polarization imply that the light coming from the reflection nebula must have been multiply scattered, before heading into the line of sight.

Two other sources (IRS2 and IRS3) are present toward the $`K`$-band reflection nebula (within 20″of IRS1). IRS2 has a class III SED and is probably a more evolved star (André & Montmerle 1994). IRS3 has a class I SED, but is much fainter than IRS1 in the near-infrared. It has a bolometric luminosity of $`0.13L_{}`$, (Bontemps et al. 2001). IRS3 shows strong $`6\mathrm{cm}`$ emission (see Leous et al. 1991, where IRS3 is designated LFAM 1).

The presence of a molecular outflow has not been firmly established. High velocity red- and blue-shifted CO millimeter emission to the south of GSS 30 has previously been reported by Tamura et al. (1990), but since both lobes are located to the south of the infrared source, the gas is likely to be associated with the VLA 1623 jet, which is passing only $`30\mathrm{}`$ to the SW of IRS1. Using millimeter interferometric line data, Zhang et al. (1997) find evidence for a spherical expansion of the core surrounding the three sources in GSS 30. In addition, the inner region of the core seems to be cleared, which is indicative of a young outflow. This is additionally supported by the presence of variable unbound water maser emission from within $`0\stackrel{}{.}3`$ of IRS1 (Claussen et al. 1996).

## 2 Observations

The observations of GSS 30 IRS1 were carried out using the long wavelength (LW) medium resolution mode on the Infrared Spectrometer And Array Camera (ISAAC) mounted on the Very Large Telescope (VLT-UT1) at the Paranal Observatory on September 3, 2001. The detector for $`35\mu \mathrm{m}`$ observations is a $`1\mathrm{K}\times 1\mathrm{K}`$ Aladdin InSb array, which allows a spectral coverage of $`0.237\mu \mathrm{m}`$ per setting in the $`M`$-band with a spectral resolution of up to $`R=\mathrm{10\hspace{0.17em}000}`$. The detector resolution in the spatial direction in the spectroscopic modes is $`0\stackrel{}{.}148/\mathrm{pixel}`$. The spectrum of this source was taken as a part of a large program to study ices and gas around young low mass stars (van Dishoeck et al. 2001).

The data were obtained under excellent photometric conditions in service mode with $``$ 0$`\stackrel{}{.}`$3 infrared seeing and $`10\%`$ humidity which gave a very stable atmosphere. The frames were both chopped and nodded along the slit with a $`15\mathrm{}`$ chop throw. The 0$`\stackrel{}{.}`$6 slit was used in two settings for 20 minutes per setting yielding a final signal to noise ratio on the continuum of the extracted spectrum ranging from 50 to more than 100 over a spectral range from $`4.5\mu \mathrm{m}`$ to $`4.82\mu \mathrm{m}`$ and with a spectral resolution of $`R=\mathrm{5\hspace{0.17em}000}`$.

The standard star BS6084 (B1III) was observed at air masses between 1.2 and 1.5 immediately before GSS 30 for removal of telluric features, while the source was observed at air masses between 1.5 and 2.1. Because the source was already descending at the time of observation, it was not possible to observe source and standard at the same airmass as is otherwise of importance in order to obtain a good telluric subtraction. However, due to the exceptional stability of the atmosphere through most of the night we were able to correct for the airmass difference using a simple Lambert-Beer law, which states that the depth of the atmospheric lines scales exponentially with the airmass. The final corrected spectrum turned out to be of acceptable quality in spite of the airmass difference. Our experience is that it is crucial to use a B star or an early A star as the standard because many narrow intrinsic stellar lines become visible in the $`35\mu `$m region for spectral resolutions better than $`1000`$ if an F or later type star is used.

The spectrum was reduced with IDL routines using standard methods appropriate for chopped and nodded infrared spectroscopy. Bad pixels were removed, the frames were distortion corrected using the telluric features in a star-trace map as reference, and AB nodded pairs were subtracted to remove the long time scale differences which the chopping is not capable of removing. The spectra were extracted by integrating over a spatial width determined by the positions along the slit where the signal drops below the $`3\sigma `$ level. The raw spectrum was then ratioed by the standard after applying an optimized shift and airmass correction. In the final spectrum, the points where telluric residuals are apparent have been removed, including the region short-wards of $`4.55\mu \mathrm{m}`$, where atmospheric $`\mathrm{CO}_2`$ lines render the spectrum useless. The final result is shown in Fig. 1.

The CO ro-vibrational lines in the spectrum are red-shifted with $`20\mathrm{km}\mathrm{s}^1`$ due to the systemic velocity of the Earth relative to the source, which means that there is a significant overlap between intrinsic CO lines and telluric absorption by CO. In spite of this, the telluric division is very good, especially when taking the unfavorable airmasses of source and standard into account. A few residuals from telluric $`{}_{}{}^{12}\mathrm{CO}`$ absorption are visible and are marked in the final spectrum. Since the resolution of our spectra is $`60\mathrm{km}\mathrm{s}^1`$, none of the features in the source are unaffected by the telluric lines. At higher resolution ($`R>\mathrm{20\hspace{0.17em}000}`$) the telluric lines may saturate, thereby forcing a division by zero when ratioing with the standard and resulting in a loss of information in the affected pixels. At sufficiently high resolution, this will not be a problem since the narrow telluric lines will often be completely offset from the intrinsic lines. However, at a resolution of 5 000, the typical telluric line has a depth which is only a fraction of the intrinsic continuum allowing a significant signal to be present in the standard spectrum even in the center of a telluric line and only very few lines are partially or completely lost due to saturated telluric lines.

By measuring the strength of the small residuals from telluric lines of species not intrinsically present in the source spectrum, the systematic uncertainty in the source lines is estimated to be less than 15%. The statistical uncertainty for the brighest lines is typically a few %, while the faintest lines are $`3\sigma `$ detections. The spectrum was wavelength calibrated relative to the telluric lines and flux calibrated relative to the standard star. The wavelength calibration of high resolution M band spectra is very accurate due to the large number of telluric lines, and a conservative estimate for the precision is $`5\mathrm{km}\mathrm{s}^1`$ over the entire spectrum. Indeed, the standard deviation of the line center velocities given in Table 1 is $`4.2\mathrm{km}\mathrm{s}^1`$. Due to the difference in airmass between standard and source, we estimate the systematic uncertainty in absolute flux calibration to be $`30`$%. Since the bright lines are typically most severely affected by telluric residual, all line fluxes have total uncertainties with values between 15 and 50%. The exact uncertainty for each line is difficult to estimate and we consequently adopt a value of 30% for all lines.

## 3 Results

### 3.1 Signatures of hot CO gas

It is evident that the $`M`$-band spectrum of GSS 30 is dominated by strong emission lines from the ro-vibrational ($`v`$=1–0) transitions of $`{}_{}{}^{12}\mathrm{CO}`$ gas. R- and P-branch lines within the observed spectral range are detected up to $`J`$=15 and $`J`$=18, respectively. The lines are all unresolved at a resolution of $`R=\mathrm{5\hspace{0.17em}000}`$, which gives an upper limit to the intrinsic line width of $`30\mathrm{km}\mathrm{s}^1`$. No apparent decrease in line intensity is seen towards higher rotational quantum numbers, so observing the source in the entire $`M`$-band will clearly detect more lines. The brightest lines reach peak fluxes of more than 200% of that of the continuum. However, since none of the lines are resolved, this is only a lower limit to the peak flux. Fainter detected features include $`{}_{}{}^{13}\mathrm{CO}v`$=1–0 lines as well as $`{}_{}{}^{12}\mathrm{CO}v`$=2–1 lines. The fact that the observed gas has a population excited to the second vibrational level already indicates that the gas is at least warmer than a few hundred K. The $`{}_{}{}^{12}\mathrm{CO}v`$=1–0 lines have intensities of order $`2\times 10^{13}\mathrm{erg}\mathrm{s}^1\mathrm{cm}^2`$, the $`{}_{}{}^{13}\mathrm{CO}v`$=1–0 and $`{}_{}{}^{12}\mathrm{CO}v`$=2–1 lines being 1-2 orders of magnitude fainter. The line fluxes of all unblended lines were extracted by fitting a gaussian as well as a first order polynomial to the local continuum. After subtraction of the continuum, the fluxes were derived by integrating the spectrum in a $`3\sigma `$ region around the line.

The extracted fluxes of all unblended lines are given in Table 1. Many of the faint lines are completely blended with the bright lines, but it was still possible to extract 12 $`{}_{}{}^{13}\mathrm{CO}v`$=1–0 lines and 4 $`{}_{}{}^{12}\mathrm{CO}v`$=2–1 lines, which is sufficient to constrain the most important parameters of the emitting gas.

In addition to the CO emission lines, a faint 4.67 $`\mu \mathrm{m}`$ CO ice feature is visible. Since the feature is partly filled by the bright $`{}_{}{}^{12}\mathrm{CO}\mathrm{P}(1)`$ line it is not possible to obtain a well-defined ice profile. For the same reason the affected emission line is not included in the analysis and modeling. There is also a very broad ($`\mathrm{FWHM}500\mathrm{km}\mathrm{s}^1`$), $`4.65\mu \mathrm{m}`$ Pfund $`\beta `$ hydrogen recombination line visible, although with the bright CO lines superposed.

Finally the $`\mathrm{H}_200\mathrm{S}(9)`$ line is detected with an integrated line flux of $`2\pm 1\times 10^{14}\mathrm{erg}\mathrm{s}^1\mathrm{cm}^2`$, signifying the presence of hot ($`>1000\mathrm{K}`$) molecular hydrogen or possibly shock-excited gas.

### 3.2 Size of the emitting region

In order to determine the parameters of the line emitting gas unambiguously, the size of the emitting region must be known. The FWHM of the continuum in the spatial direction is about $`0\stackrel{}{.}5`$, which is consistent with the seeing throughout the observation. However, the lines appear to be spatially resolved, although most of the line flux is centered on the source. The wings of the spatial profile of each line extend up to $`2\mathrm{}`$ away from the central source along the slit. This is illustrated in Fig. 2, which shows the normalized average of the spatial profile of the lines, with the normalized average of the continuum profile subtracted. It is evident that the line emission is resolved and detected to a distance of more than $`320\mathrm{AU}`$ from the central source. There is also an indication that the northern part of line emission is stronger than the southern part, which agrees with the morphology seen in the $`L`$-band image of IRS1 shown in Fig. 3.

This extension suggests that the line emission is associated with the scattering nebula and is not coming directly from the central source. Since the continuum is less spatially extended than the line emission, the line emitting region appears to be spatially distinct from the continuum emitting region. Possible interpretations are that the line emission can be coming from gas present in the bipolar lobes (e.g. an outflow), from resonance scattering of continuum emission or from scattering of lines emitted closer to the star. Each of these possibilities will be discussed in detail below as the size of the emitting region depends strongly on the scenario.

## 4 Models

### 4.1 Optically thin LTE models

As a simple first approximation, the extracted line intensities were fitted with a one-temperature optically thin LTE model. The molecular data were taken from the HITRAN database (Rothman et al. 1998) and the relative line intensities were calculated directly from a Boltzmann level population distribution. The two parameters of the model, the gas kinetic temperature and the CO column density, should be well constrained in the optically thin limit. If only the $`{}_{}{}^{12}\mathrm{CO}v`$=1–0 transitions are used, an optically thin model is to some extent degenerate in the two parameters. This degeneracy is due to the fact that the populations in the vibrationally excited levels are much more sensitive to temperature changes in the $`400800\mathrm{K}`$ range than the rotational level populations. Since lines from two species are available, the degeneracy is broken. It is reasonable to assume that the emission from both species is coming from the same mass of gas, and that the $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$-ratio is close to 60 (e.g. Bensch et al. 2001).

This model is not able to provide a good fit to the entire line spectrum. However, from Fig. 4 it is evident that the faint $`{}_{}{}^{13}\mathrm{CO}v`$=1–0 and $`{}_{}{}^{12}\mathrm{CO}v`$=2–1 lines can be simultaneously fitted with a single temperature gas in the optically thin limit. The model lines of $`{}_{}{}^{12}\mathrm{CO}v`$=1–0 are however an order of magnitude too bright. Excluding the $`{}_{}{}^{12}\mathrm{CO}v`$=1–0 lines from the fit yields a well-constrained temperature of $`515\pm 25\mathrm{K}`$ and a total $`{}_{}{}^{12}\mathrm{CO}`$ column density of $`4\pm 1\times 10^{14}\mathrm{cm}^2`$. Since the fit includes lines from both CO isotopic species, we conclude that it is unlikely that the apparent discrepancy between the $`{}_{}{}^{12}\mathrm{CO}v`$=1–0 and $`{}_{}{}^{13}\mathrm{CO}v`$=1–0 lines is due to a deviation from the normal interstellar ratio. The simplest explanation is that the $`{}_{}{}^{12}\mathrm{CO}v`$=1–0 lines are optically thick, while the rest remain optically thin. This would conveniently prevent us from seeing the whole column of $`{}_{}{}^{12}\mathrm{CO}`$ and could possibly lower the line intensities enough to explain the data. However, the fundamental ro-vibrational transitions of $`{}_{}{}^{12}\mathrm{CO}`$ do not produce optically thick lines for column densities less than $`10^{17}\mathrm{cm}^2`$, which is still 250 times more than the simple model dictates. Thus the optically thin LTE model is internally inconsistent.

### 4.2 Results of the full radiative transfer model

In order to quantify the discrepancy more accurately, we computed the full radiative transfer for a single temperature static envelope in LTE. For simplicity a spherical geometry was assumed since the results are not very sensitive to geometry. The size and density of the emitting region cannot be constrained by the model, so the quantities fitted were the single gas temperature and the column density of the gas, by keeping the size fixed at the observed 4″.

Furthermore, it is assumed that the line intensity is independent of the underlying continuum. This last assumption may not be completely valid, since the CO lines could be present in absorption in the continuum before the lines from the emitting gas are superposed. The effect of this would be to lower the line peak flux with maximally the continuum flux level. Since the lines are completely unresolved, the true peak flux of the lines is higher than that observed. The observed peak flux of the $`{}_{}{}^{12}\mathrm{CO}v`$=1–0 lines is 2–3 times higher than the continuum flux, so we conclude that this effect can only effectively lower the intensity of these lines with a small fraction at most. The fainter lines can be significantly suppressed by absorption in the continuum if they are optically thick, but this would only increase the discrepancy between the bright and faint lines. Also, if the absorption lines are significantly broader than the superposed emission lines, the observed line intensities could be sufficiently suppressed. However, since we have no possibility of constraining this with the present data, we assume that the possible absorption lines and the emission lines have similar intrinsic widths, i.e. the contribution from absorption in the continuum is negligible. Higher resolution spectra may clarify the matter.

The $`{}_{}{}^{13}\mathrm{CO}v`$=1–0 lines were first fitted independently of the $`{}_{}{}^{12}\mathrm{CO}v`$=1–0 and $`{}_{}{}^{12}\mathrm{CO}v`$=2–1 lines. The model confirms that the two sets of lines cannot be fitted simultaneously, although each set is consistent with a single temperature gas. There is a certain degeneracy in the two free parameters, especially in the $`{}_{}{}^{13}\mathrm{CO}`$ lines, and the model gives a one-dimensional family of good fits for each species. However, the two families do not overlap for the observed line intensities. Fitting the $`{}_{}{}^{12}\mathrm{CO}v`$=1–0 and the $`{}_{}{}^{12}\mathrm{CO}v`$=2–1 lines simultaneously yields poorer fits than fitting them separately. All the lines are in the extreme optically thin limit with typical values of $`\tau 10^3`$ for the $`{}_{}{}^{12}\mathrm{CO}v`$=1–0 lines. The discrepancy between the $`{}_{}{}^{12}\mathrm{CO}v`$=1–0 and $`{}_{}{}^{12}\mathrm{CO}v`$=2–1 lines may be due to non-LTE effects, reflecting that the rotational temperature $`T_{\mathrm{rot}}`$ and the vibrational temperature $`T_{\mathrm{vib}}`$ are not equal. However, non-LTE effects cannot explain why the $`{}_{}{}^{12}\mathrm{CO}v`$=1–0 to $`{}_{}{}^{13}\mathrm{CO}v`$=1–0 line ratios differ with almost an order of magnitude from the predictions of the simplest LTE model, since the excitation structures of different isotopic species must be almost identical.

To explain the data, we propose that all the observed line emission has been scattered by dust in the bipolar cavity seen in the near-infrared. The scattering efficiency is very low at $`4.7\mu \mathrm{m}`$, and thus only a small fraction of the original line emission will be scattered into the line of sight. This means that the line intensities from the emitting region and consequently the column densities can reach much larger values than those previously suggested by the model. To test this hypothesis, we adopted a factor, $`ϵ=S\times (\pi R^2)/\mathrm{\Omega }`$, which accounts for our lack of information about the true solid angle of the emitting region $`\mathrm{\Omega }`$, and the ratio of intrinsic line flux to line flux scattered into the line of sight, $`S`$. $`R`$ is the observed angular radius of the emitting region, which is taken to be $`R=2\mathrm{}`$. Absolute upper and lower bounds to the size of the emitting region are given by the natural conditions $`S>1`$ and $`\mathrm{\Omega }<\pi R^2`$. Decreasing the size of the emitting region will increase the column densities, but since line emission is observed far from the source, we cannot explain the entire discrepancy with a smaller emitting size. To recover the original flux it is necessary to multiply the observed values with this factor. Multiplying all observed line intensities with a sufficiently large factor will allow the $`{}_{}{}^{12}\mathrm{CO}v`$=1–0 lines to become optically thick. We are able to constrain the scattering factor in the model by demanding a simultaneous fit to all the lines.

The best fitting value of the scattering factor is $`ϵ=1.4\times 10^3`$ with a reduced $`\chi ^2`$ of 1.29, which gives a well-constrained single temperature LTE gas at $`515\pm 5\mathrm{K}`$ and $`\mathrm{N}(\mathrm{CO})=1.5\pm 0.5\times 10^{18}\mathrm{cm}^2`$, as indicated in Table 2. With an extent of the emitting region of $`4\mathrm{}`$ an upper limit to the total gas mass of $`200\mathrm{M}_{}`$ is found. Since the total gas mass scales with $`\mathrm{\Omega }^{0.5}`$ an absolute lower limit to the mass is $`1.3\mathrm{M}_{}`$ with a diameter of the emitting region of $`0.19\mathrm{}=30\mathrm{AU}`$, assuming that $`S=1`$ and ignoring that line emission is observed to much larger distances from the central source. Since $`S`$ must be significantly less than unity, the total mass is most likely significantly larger than the lower limit. It is not necessary to invoke any non-LTE effects to obtain a nearly perfect fit to all the observed lines and the assumption of LTE holds such that $`T_{\mathrm{rot}}=T_{\mathrm{vib}}`$. The $`\chi ^2`$ maps of the best-fitting models are shown in Fig. 5.

It is possible to make an estimate for an absolute lower limit to the grain size by assuming $`S=1`$. In this case the scattering efficiency, $`Q_{\mathrm{sca}}>ϵ`$, since it is not possible to scatter more light into the line of sight than out of it. This happens for silicate grain sizes larger than $`0.20.3\mu \mathrm{m}`$ (Draine & Lee 1984; Draine 1985). As a consequence the grains must be large, perhaps as large as $`1\mu \mathrm{m}`$, which is expected for grain growth in dense media (Pendleton et al. 1990). On the other hand, the lower limit on the size of the scattering grains does not fit with the polarization maps of the GSS 30 nebula. Chrysostomou et al. (1996) finds that grains larger than $`0.35\mu \mathrm{m}`$ are inconsistent with the variation of linear polarization with wavelength. However, the polarization was mapped on a much larger angular scale than the size of the nebula at $`4.6\mu \mathrm{m}`$, so the discrepancy may be explained by grain growth within a few hundred AU of the central source.

### 4.3 Other scenarios and pumping mechanisms

One way of avoiding the scattering assumption entirely may be to apply a two-temperature model to the data. In this scenario a small part of the gas close to the central source with a kinetic temperature of $`1000\mathrm{K}`$ would be responsible for the $`{}_{}{}^{12}\mathrm{CO}v`$=2–1 lines as well as part of the high $`J`$ lines of the $`v`$=1–0 lines, while a massive outer component at $`200\mathrm{K}`$ would dominate the emission in the low $`J`$ lines. The scattering assumption would not be needed since the $`{}_{}{}^{12}\mathrm{CO}v`$=1–0 lines require a much higher column density to create the observed line intensity at $`200\mathrm{K}`$ compared to $`500\mathrm{K}`$ and will thus more readily become optically thick. We applied a two-temperature model by requiring that the total column density must be the same as that determined by the single-temperature model. This is reasonable because the optical depth of the lines is given by the ratio between the $`{}_{}{}^{12}\mathrm{CO}`$ and $`{}_{}{}^{13}\mathrm{CO}`$ lines independently of $`ϵ`$ and therefore must be conserved. As in the single-temperature model three parameters are varied, namely the two temperatures and the ratio of the column density between the components. In other words, $`ϵ`$ is exchanged as a free parameter with a second temperature, thus keeping the number of free parameters constant.

The $`{}_{}{}^{13}\mathrm{CO}`$ lines can be fitted well within the two-temperature model, although the parameters are somewhat degenerate in the same way as in the single temperature model. However, it was not possible to fit the $`{}_{}{}^{12}\mathrm{CO}`$ lines within this scenario, with the best reduced $`\chi ^2`$ being of order 10. Furthermore there is a preference towards a single low temperature in the $`{}_{}{}^{12}\mathrm{CO}`$ lines which is inconsistent with the $`{}_{}{}^{13}\mathrm{CO}`$ lines. Inclusion of all lines gives no meaningful fits within the parameter space. The explanation is that it is not possible to combine the $`{}_{}{}^{12}\mathrm{CO}v`$=2–1 and the $`{}_{}{}^{12}\mathrm{CO}v`$=1–0 lines, since the high temperature component easily produces too much $`v`$=1–0 emission in order to fit the $`v`$=2–1 emission. We conclude that a two-temperature model is clearly inferior to the scattering model in explaining the observed line emission.

Radiative pumping by a strong infrared continuum is another possibility that needs to be considered as an alternative excitation mechanism and a way to avoid the scattering scenario by introducing non-LTE effects. This would produce a single vibrational temperature (Scoville et al. 1980). However it does not explain why the rotational and vibrational temperatures seem to be identical. Also, luminosities larger than $`10^3\mathrm{L}_{}`$ are required to effectively pump CO molecules to the second vibrational level or higher at distances from the central source beyond $`1\mathrm{AU}`$ (Scoville et al. 1980). Using the 1-dimensional Monte-Carlo radiative transfer code from Schöier (2000) illustrates the effects of radiative excitation on the line populations. The simplest possible model for radiative excitation assumes a blackbody with a luminosity $`L=25\mathrm{L}_{}`$ at the center of a single density spherical and static envelope . The blackbody temperature is taken to be the colour temperature of the continuum around $`4.7\mu \mathrm{m}`$, determined from the $`M`$-band spectrum to $`420\mathrm{K}`$. The radius used is the observed $`320\mathrm{AU}`$, the gas kinetic temperature is assumed to be a power law beginning at the surface of the central source. The molecular constants are taken from Schöier et al. (2002). This leaves the gas density as the only free parameter, which basically determines the optical depth of the lines and consequently the importance of line trapping effects. The model calculations show that a density of at least $`n_{\mathrm{H}_2}10^6\mathrm{cm}^3`$ is needed to reproduce the approximate flux level of the $`{}_{}{}^{12}\mathrm{CO}v`$=1–0 lines. The vibrational excitation temperature is largely determined by the continuum flux level at $`4.7\mu \mathrm{m}`$ and therefore by the luminosity and to a smaller extent the temperature of the radiating source. Since the dominating mechanism for populating the $`v`$=2 level is via the $`v`$=1 level, and since both the $`v`$=2–1 and $`v`$=1–0 transitions are placed in the same narrow wavelength region, the radiative excitation rates are strongly dependent on local photon density in the $`4.7\mu \mathrm{m}`$ region. Consequently, as long as the direct route through the $`v`$=2–0 overtone transitions around $`2.35\mu \mathrm{m}`$ is not important, the vibrational level population depends on the temperature of the input radiation field only through the flux level at $`4.7\mu \mathrm{m}`$. At the adopted luminosity, the $`{}_{}{}^{12}\mathrm{CO}v`$=2–1 line fluxes are two orders of magnitude too low, reflecting an excitation temperature of $`400\mathrm{K}`$. To significantly populate the second vibrational levels, a luminosity of at least $`10^3\mathrm{L}_{}`$ is required, in agreement with Scoville et al. (1980). Finally, the population of the rotational levels shows significant departures from a Boltzmann distribution, especially when line trapping becomes important. Altogether, this indicates that radiative pumping as an excitation mechanism in the case of GSS 30 IRS1 is unlikely.

Finally, the possibility that the lines are continuum emission which is resonance scattered in CO transitions should be considered. It is known that this takes place in the case of asymptotic giant branch stars, e.g. in the case of Mira (Ryde & Schöier 2001). However, since this mechanism probes the scattering gas rather than the emitting medium, it will not solve the discrepancy. Also it would require that the gas at a distance of $`320\mathrm{AU}`$ from the central source has a single temperature of $`500\mathrm{K}`$, which seems highly unlikely.

### 4.4 Hydrogen recombination lines and mass loss rate

There are a number of hydrogen recombination lines common to YSOs ($`\mathrm{Pf}\beta `$, $`\mathrm{Br}\alpha `$, $`\mathrm{Pf}\delta `$) present in the $`35\mu \mathrm{m}`$ region, which are usually taken as an indicator of accretion activity. This is particularly true in the case of low mass young stars, where the hard radiation from the boundary layer between disk and star is needed to produce the ionizing radiation, since the central star is not hot enough. Because these lines have broad wings ($`3500\mathrm{km}\mathrm{s}^1`$), they are often interpreted as coming from a strong stellar wind, which may be the engine of an outflow. The possibility that the CO gas could be associated with a wind or molecular outflow should therefore be explored. If associated with a stellar wind, the lines could be formed either as cooling emission from the gas as molecules are reformed or as emission from a shock as the wind collides with the circumstellar envelope or disk. It is known that no large scale outflow is present which could be produced by any of the three stars in the cluster, but there is some evidence of an expanding core surrounding the cluster from interferometric observations of rotational $`{}_{}{}^{13}\mathrm{CO}`$ and $`\mathrm{C}^{18}\mathrm{O}`$ lines (Zhang et al. 1997).

Hydrogen recombination lines emitted from a plasma moving at supersonic speeds can be modeled in the Sobolev approximation (Sobolev 1960), which is valid when the thermal line broadening is much smaller than the wind speed. This formalism has been treated for spherical ionized winds in LTE from evolved stars in Castor (1970) and Krolik & Smith (1981), and for Herbig AeBe stars in Nisini et al. (1995). Using this formalism an ionized mass loss rate can be estimated for the star if a wind origin for the hydrogen lines is assumed. The extinction corrected line intensities for the hydrogen recombination lines observed in GSS 30 IRS1 are tabulated in Table 3 along with the mass loss rates derived from each line. The line ratios and the mass loss rate are not sensitive to the temperature, so a constant electron temperature of $`10^4\mathrm{K}`$ has been assumed. Finally, it is assumed that the central star does not shadow a significant part of the radiation from the wind. Changing the model parameters, i.e. the temperature and the size of the ionized region, with a factor of two changes the mass loss rate with a factor of two at most. The velocity structure is such that the wind starts at the surface of the star with $`20\mathrm{km}\mathrm{s}^1`$ and rises quickly to $`450\mathrm{km}\mathrm{s}^1`$ to fit the observed FWHM of the $`\mathrm{Pf}\beta `$ line. The exact shape of the velocity structure is otherwise not significant for the derived mass loss rate.

We find that an ionized mass loss rate of $`1.7\times 10^6\mathrm{M}_{}\mathrm{year}^1=0.6\mathrm{M}_{}\mathrm{year}^1`$ is consistent with the two observed Pfund lines. The $`\mathrm{Br}\alpha `$ line is too faint for this mass loss rate, which may indicate deviations from LTE. We can rule out that it is an effect from a poor extinction correction since the $`\mathrm{Br}\alpha `$ line lies between the $`\mathrm{Pf}\beta `$ and $`\mathrm{Pf}\delta `$ lines in wavelength. The most serious problem with this model is likely that very little is known about the structure of winds from low mass embedded YSOs; in particular, the ionization structure will depend directly on energetic processes which are poorly understood, such as accretion activity or interactions of the circumstellar matter with strong magnetic fields. Also, if the wind is predominantly neutral, then the total mass loss rate will be much higher than the ionized mass loss rate and will become inconsistent with the non-detection of an outflow. It is possible that other mechanisms such as accretion flows must be considered for the formation of recombination lines from low mass young stars. The hydrogen lines are therefore not likely to provide any strong constraints on the emitting CO gas and a detailed discussion on the hydrogen emission is outside the scope of this paper.

## 5 Discussion

In summary, to determine the origin of the CO gas emission it must be explained why only a single well-defined temperature is needed. Which mechanism can heat up to $`100\mathrm{M}_{}`$ of gas to a unique temperature of $`500\mathrm{K}`$, yet keep the intrinsic line width less than $`30\mathrm{km}\mathrm{s}^1`$? Finally, why is no other similar embedded source from the literature showing the same strong CO rovibrational emission as GSS 30 IRS1?

### 5.1 Outflow?

If the warm CO gas resides in the wind/outflow component of the circumstellar environment then an estimate of the time needed to create it can be found from the mass loss rate derived in Section 4.4 under the assumption that the wind is predominantly ionized. A mass between 3 and 100$`M_{}`$ corresponds to a production time of minimally 5 years and maximally a few centuries, while the most likely time is about one century. It seems unlikely that hot, ionized material emitted over a longer period should thermalize at a single temperature. Also, an outflow origin of the line emission would produce broad wings in the lines, which is not observed. If the wind is predominantly neutral, the resulting mass loss rate would be so high that a clear outflow should have been detected.

### 5.2 Inner disk?

Another possibility is that the emission is produced by warm thermalized gas in the disk itself. Since more evolved T Tauri disks are known to exhibit similar CO emission, although with smaller intensities, it is conceivable that we are seeing an equivalent process in the case of a younger disk. A typical circumstellar disk around a low mass young star is often observed to have a mass of a few times the minimum solar nebula, i.e. $`\mathrm{M}_{\mathrm{disk}}\mathrm{a}\mathrm{few}0.01\mathrm{M}_{}`$ (André & Montmerle 1994; Osterloh & Beckwith 1995). The disk mass for GSS 30 IRS1 inferred by the 1.3 mm continuum emission is $`M_{\mathrm{disk}}=0.03M_{}`$ (André & Montmerle 1994). In hydrostatic equilibrium, a disk irradiated by the central star has a surface density $`\mathrm{\Sigma }R^{3/2}`$ (Chiang & Goldreich 1997). If viscous dissipation in the disk is taken into account the surface density attains a flatter $`R`$-dependency, $`\mathrm{\Sigma }R^1`$. In both cases the accumulated disk mass reaches the observed $`10100\mathrm{M}_{}`$ between 2 and 10 AU from the central star. A single temperature of $`515\mathrm{K}`$ is however not consistent with the disk models, which prescribe a large range in temperatures with a small 1000-2000 K component within a few stellar radii to an extended 100-200 K component within a few AU of the central star. However, in no models for passive circumstellar disks is it expected that the disk exhibits a single temperature, indeed quite the contrary is the case. Furthermore, Keplerian rotation within a few AU would produce observable line broadening. The narrow lines could only be explained if we are seeing light emitted vertically from the plane of the disk before being scattered into the line of sight, effectively removing the effects of rotation.

### 5.3 Accretion shock?

A final option is that the emission lines are cooling lines from post-shocked dense gas. Vibrational $`\mathrm{H}_2`$ emission is one of the principal tracers of low density ($`n_{\mathrm{H}_2}10^7\mathrm{cm}^3`$) shocked gas. As mentioned in Section 3.1 the (0–0) $`\mathrm{S}(9)`$ line of $`\mathrm{H}_2`$ at $`4.695\mu \mathrm{m}`$ is seen in our $`M`$-band spectrum, while no $`\mathrm{H}_2`$ lines are seen in the $`K`$-band spectrum of GSS 30 IRS1 by Greene & Lada (1996) although no upper limit is given. Assuming that the observed $`\mathrm{H}_2`$ (0–0) $`\mathrm{S}(9)`$ line is thermalized at the CO temperature, the required molecular gas mass is $`5\mathrm{M}_{}`$. This is assuming that the $`\mathrm{H}_2`$ emission is seen directly and is not corrected for extinction, since the extinction is hard to estimate for embedded stars. A typical Ophiuchus extinction of $`A_\mathrm{V}=25\mathrm{mag}`$ (Teixeira & Emerson 1999) will increase the molecular gas mass by a factor of 2. If the emission in the $`\mathrm{H}_2`$ line is scattered in the same way as the CO emission all values must be multiplied by $`ϵ`$. If the $`\mathrm{H}_2`$ emission is thermalized at $`515\mathrm{K}`$ and directly observed then the integrated line flux expected for the $`2.12\mu \mathrm{m}`$ (1–0) S(1) line is $`7.5\times 10^{13}\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$, which should be observable. If the $`\mathrm{H}_2`$ lines are scattered into the line of sight with an efficiency of $`ϵ`$, then no $`\mathrm{H}_2`$ lines should be visible in the $`M`$-band. Consequently, the $`\mathrm{H}_2`$ line does not seem to be directly associated with the same gas emitting the CO lines since it is at least a factor of 10 too bright, but it may be an indication of a warmer component or the line may be pumped. Sensitive observations of $`\mathrm{H}_2`$ lines in the K and L band are needed to unambigously determine the origin of the molecular hydrogen emission.

Neufeld & Hollenbach (1994) show that if a dense gas of $`10^{7.5}<n_{\mathrm{H}_2}<10^{12}\mathrm{cm}^3`$ is shocked, the vibrational emission from $`\mathrm{H}_2`$ is significantly suppressed partly because the shock dissociates the $`\mathrm{H}_2`$ molecules and partly because the reforming molecules can easily be collisionally de-excited in the dense post-shock gas. The main coolants will be rotational and vibrational lines from $`\mathrm{H}_2\mathrm{O}`$, $`\mathrm{OH}`$ and $`\mathrm{CO}`$. This is valid for a wide range of shock velocities $`5<v_\mathrm{s}<100\mathrm{km}\mathrm{s}^1`$ as long as the shock is of jump-type. A J-shock occurs in any case for $`v_\mathrm{s}30\mathrm{km}\mathrm{s}^1`$ and for smaller shock velocities if the magnetic field is weak or if the length scale over which the shock occurs is small compared to the length scale for acceleration of charged particles. Furthermore, the theoretical post-shock temperature structure gives an interesting prediction. A J-shock will dissociate the gas at the shock front. The gas then cools slowly from $`10^5\mathrm{K}`$ principally through the $`\mathrm{Ly}\alpha `$-line until $`\mathrm{CO}`$ forms at around $`7000\mathrm{K}`$ and $`N(\mathrm{HI})10^{21}\mathrm{cm}^2`$. Once the CO molecules are available the gas cools quickly to about $`500\mathrm{K}`$ through the $`\mathrm{CO}`$ vibrational lines, where a constant temperature plateau is maintained due to the release of chemical potential energy as $`\mathrm{H}_2`$ is formed which is balanced against the thermalized vibrational emission from molecules. The plateau occurs at an atomic hydrogen column density of $`10^{21}10^{22}\mathrm{cm}^2`$, which in the model by Neufeld & Hollenbach (1994) corresponds to a $`\mathrm{CO}`$ column density of $`3\times 10^{17}3\times 10^{18}\mathrm{cm}^2`$. It is stressed that the slower C-shock scenario was not treated in such high density models, and it may produce similar cooling lines. Further modelling is needed to exclude a slow shock. It is evident that the model values are remarkably similar to the values observed in GSS 30 IRS1 with temperature $`T_{\mathrm{gas}}=515\mathrm{K}`$ and column density $`N(\mathrm{CO})=2\times 10^{18}\mathrm{cm}^2`$. However, it is unclear why no hot $`50007000\mathrm{K}`$ component in the fundamental CO band is observed, which is also expected from the shock. Also no overtone emission is seen in the $`K`$-band spectrum of Greene & Lada (1996). One plausible explanation is that the accretion activity is highly variable and that the hot component is only seen as the shock is occurring while the warm thermalized gas is seen in the post-shock medium. This would be supported by the variable water maser emission. The single observed $`\mathrm{H}_2`$ line is not inconsistent with an accretion shock although it is too bright to be associated with a dense shocked gas and may only be an indication of a warmer less dense component.

As described in Neufeld & Hollenbach (1994), the most likely site for a dense shock in an embedded YSO is in a circumstellar disk which is accreting matter. When infalling matter is colliding with the disk at supersonic velocities, a J-shock is expected to occur at the surface of the disk. If the accretion velocities are too low to produce a J-shock, the molecules will not be dissociated and the observed temperature plateau will not be present. The main uncertainty is if the high infall velocities required to create a dissociative shock are possible. However, infall velocities of $`10\mathrm{km}\mathrm{s}^1`$ within $`10\mathrm{AU}`$ are reasonable (Cassen & Moosman 1981). Therefore, if the dissociative accretion shock scenario holds, the intrinsic line widths are predicted to be larger than or very close to $`10\mathrm{km}\mathrm{s}^1`$. An accretion shock at a large distance ($`10\mathrm{AU}`$) from the central star is consistent with the observation that the line emitting region is physically separated from the continuum emitting region. A shock will not heat a large amount of dust and the continuum emission, which has a colour temperature of about $`700\mathrm{K}`$, is likely to be produced by irradiated dust within $`1\mathrm{AU}`$ of the star.

Other observable tracers of a shock in a dense medium are the rovibrational lines from the $`\mathrm{H}_2\mathrm{O}\nu _2`$ bending mode around $`6\mu \mathrm{m}`$ as well as the vibrational transitions of $`\mathrm{OH}`$, all of which are unavailable from ground-based instruments. However, they may be detected by the upcoming Stratospheric Observatory For Infrared Astronomy (SOFIA) mission.

The possible dissociation of the molecules has important consequences for the chemistry in the inner disk. If the emission is indeed coming from an accretion shock, and if a significant amount of matter passes through the shock, it will significantly affect the models of chemical evolution of the innermost parts of disks around low mass YSOs, by evaporating ice mantles and inducing high temperature chemistry.

## 6 Summary

In this paper, we have analyzed the $`4.54.8\mu \mathrm{m}`$ spectrum of the embedded young stellar object GSS 30 IRS1 in the $`\rho `$ Ophiuchus core. The observed emission lines from the rovibrational transitions of $`\mathrm{CO}`$ are fitted with a simple full radiative transfer 1D model. It is found that the line emission must be scattered on a bipolar cavity in order to simultaneously account for the size of the observed emitting region, the absolute flux level of the lines and the ratio of $`{}_{}{}^{12}\mathrm{CO}`$ to $`{}_{}{}^{13}\mathrm{CO}`$ lines. In this case the emission is well fitted by a single temperature gas with $`T=515\pm 5\mathrm{K}`$ and a column density of $`N(\mathrm{CO})=2\pm 0.5\times 10^{18}\mathrm{cm}^2`$. Furthermore, assuming a two-temperature distribution of the gas does not yield satisfactory fits.

The observed emission line spectrum can be best explained by the presence of a J-shock in a dense medium, although other possible scenarios cannot be ruled out. Most likely the shock is produced by the accretion of gas onto a dense disk within a few tens of AU from the star as predicted by Neufeld & Hollenbach (1994). Only this scenario explains why a large amount of $`\mathrm{CO}`$ gas is thermalized at a single temperature of $`500\mathrm{K}`$ and why the emission lines are so narrow.

###### Acknowledgements.

The authors wish to thank the VLT staff for all their help in obtaining the observations and in particular Chris Lidman for many helpful comments on the data reduction. We also acknowledge the many constructive suggestions made by an anonymous referee, which helped improve the quality of this paper. This research was supported by the Netherlands Organization for Scientific Research (NWO) grant 614.041.004, the Netherlands Research School for Astronomy (NOVA) and a NWO Spinoza grant.


text:
# 1 The Model

## 1 The Model

To explain the infrared and millimetre emission oscillations from GRS 1915+105, we have adopted the Blandford & Königl (1979, hereafter BK79) model of a relativistic jet, containing power-law electron distributions. Adiabatic expansion gives the jet a conical shape with an electron density and magnetic field strength that decreases with radius. The jet emits partially self-absorbed synchrotron emission over a radial range covering several orders of magnitude, and hence from a wide range of optical depths. Therefore, the total jet emission consists of the summation self-absorbed synchrotron spectra from regions of the jet of decreasing optical depth (or increasing radius), and hence electrons at each radius, from the base to the top of the jet, emit a spectrum that is progressively shifted towards lower frequencies. This summation results in a ‘flat’ region to the spectrum, which has a slope of zero if energy losses due to adiabatic decompression are ignored (as in the BK79 model). However, if some energy losses do occur then this region has a positive, inverted slope, the magnitude of which depends upon the magnitude of the energy loss.

The extent of the flat spectral region, $`\nu _{\mathrm{max}}/\nu _{\mathrm{min}}`$, is solely determined by the emission region size, which in our model is fixed by the jet velocity and time-scale of variability. At the lower frequency end of the flat spectrum there is a smooth transition to an optically thick synchrotron spectrum with $`F_\nu \nu ^{5/2}`$, and the high frequency end smoothly terminates into an optically thin, $`F_\nu \nu ^{5/8}`$, synchrotron spectrum.

Full detail of the derivation of our model may be found in Collins, Kaiser & Cox (2002). The emission spectrum predicted by this model is given in terms of the optical depth function, $`\tau _\nu (r)`$.

$$F_\nu =\frac{8.8\times 10^{18}}{(D_\mathrm{j}/\text{pc})^2}b_0\left(\frac{\nu }{\mathrm{GHz}}\right)^{5/2}_{r_{\mathrm{min}}}^{r_{\mathrm{max}}}\left(\frac{r}{r_0}\right)^{3/2}\left[1\mathrm{e}^{\tau _\nu (r)}\right]dr\text{ mJy},$$

where,

$$\tau _\nu (r)=1.5\times 10^{12}a_0\left(\frac{r}{r_0}\right)^{25/8}\left(\frac{\nu }{\mathrm{GHz}}\right)^{25/8},$$
(1)
$`D_\mathrm{j}`$ is the distance to the jet, $`r_{\mathrm{max}}/r_{\mathrm{min}}`$ defines the extent of the emission region, and all parameters have SI units unless stated otherwise.

The model has two parameters. The first, which we have named $`a_0`$, controls the optical depth as a function of radius, $`\tau _\nu (r)`$, and thus shifts the flat region of the spectrum to higher or lower frequencies. The second model parameter, which we have named $`b_0`$, only effects the flux normalisation. The $`a_0`$ parameter also has an effect upon the flux normalisation, and therefore we first fit the $`a_0`$ parameter to the spectral shape, before fitting $`b_0`$ to the flux normalisation. The model parameter, $`a_0`$, is related to the physical parameters of the system by $`a_0=w(r=r_0)k(r=r_0)B(r=r_0)^{17/8}`$, where $`k`$ is the normalisation value of the electron density distribution, and $`B`$ is the magnetic field strength. The $`b_0`$ parameter is defined as $`b_0=w(r=r_0)B(r=r_0)^{1/2}`$. The half-width of the jet $`w(r=r_0)`$ is fixed by observational constraints, and therefore from $`a_0`$ and $`b_0`$, we may determine $`k(r=r_0)`$ and $`B(r=r_0)`$.

To implement time dependence in this model as a variation in the flux level of the flat spectrum it would be simplest to make $`b_0`$ a function of time with $`a_0`$ fixed. However, physically this would require the electron density to vary by a process that exactly compensates variations in the magnetic field strength. Therefore we choose to make $`a_0`$ a function of time, and fix $`b_0`$, which may be simply interpreted, physically, as a variation to the electron density injected into the jet. Time dependence is then included into equation 1 by defining $`r_{\mathrm{max}}=v_\mathrm{j}t`$, where $`v_\mathrm{j}`$ is the bulk velocity of the jet material, and by increasing the injected value of $`k_0`$ with time according to a Gaussian function.

From this model we can predict three observable parameters to which our model parameters may be fit; the time lag between the peak fluxes at each frequency of observation, the flux ratio of these peak fluxes, and the flux normalisation value. The value of the $`a_0`$ model parameter affects all three of these observables. The flux ratio is determined by the position of the flat spectral region with respect to the two frequencies of observation, which is defined by $`a_0`$. The time lags are also determined by the position of the flat spectral region, because the optically thin end of the flat spectrum is formed at small radii, whereas the optically thick end is formed at much larger radii. Therefore the closer that both frequencies of observation are to the optically thin end, the smaller the time lag.

## 2 Application to GRS 1915+105

Fender & Pooley (2000) observed large amplitude oscillations in the emission from GRS 1915+105 at an infrared wavelength (2.2 $`\mu `$m) and a millimetre wavelength (1.3 mm). The time lag between the peak fluxes from each wavelength for the first oscillation was approximately 25 s, with the dereddened infrared flux reaching 525 mJy and the millimetre flux reaching 340 mJy, giving a peak flux ratio of 1.5. However, the extinction correction applied to the infrared flux is uncertain, such that the peak flux ratio is $`1.5\pm 0.6`$. Such a small flux ratio over a frequency range covering three orders of magnitude is suggestive of the ‘flat’ spectrum produced by relativistic jets.

For the fixed parameters of the model we assume that the distance to GRS 1915+105 is 11 kpc, the opening half-angle of its jet is $`4^{}`$, the base of the emission region is at $`r_{\mathrm{min}}=10^5`$ m, and the jet material has a bulk velocity of $`v_\mathrm{j}=0.6c`$. The value of the $`a_0`$ parameter was fit to the observed time lag, and, for this value of $`a_0`$, $`b_0`$ was fit to the flux normalisation of the peak millimetre flux. From these model parameters the corresponding physical parameters of jet were determined to be $`B(r=r_{\mathrm{min}})=7.8\times 10^7`$ T and $`k(r=r_{\mathrm{min}})=5.4\times 10^{40}`$ m<sup>-3</sup>. The model also predicts a flux ratio of 1.1 which agrees with the observations.

We have found that it is impossible to obtain flux levels of the order of $`100`$ mJy without unphysically high electron densities, for time lags of less than 300 s, as demonstrated by figure 1. This is not just a problem with this time dependent model, but is a fundamental problem of the BK79 model, as can be seen by performing the same fitting process to the observed flux ratio.

To agree with the observed flux levels and the time lag between peaks with realistic electron densities, the model must either predict a flux $`10^4`$ times larger with the fitted value of $`k_0`$, or predict the observed time lag with a value of $`a_0`$ that is $`10^9`$ times larger. Since relativistic effects are negligible in the GRS 1915+105 jet (see Collins et al. 2002), the later case can be immediately ruled out, as it would require the model’s fixed parameters to be adjusted beyond physical limits. It is also difficult to justify the higher flux required by former solution, as Doppler beaming and geometrical effects are not sufficient.

In conclusion it is impossible to explain synchrotron emission of this strength from a BK79 type jet with justifiable electron densities, when restricted to the observed time lags.


text:
# Molecular Bose-Einstein condensate as an amplifier of weak interactions

Abstract

Collisions and chemical reactions of molecules in Bose-Einstein condensates (BECs) are extremely sensitive to weak fields. This sensitivity arises due to the high density of compound resonances and a macroscopic number of molecules with kinetic energy $`E=0`$ (perfect energy resolution). We suggest that chemical reactions in molecular BECs could be used to enhance effects produced by small external perturbations and search for a parity-violating energy difference in chiral molecules.

PACS: 32.80.Ys, 03.75.Fi, 82.30.Cf

Techniques to trap and cool atoms and molecules are rapidly developing. Bose-Einstein condensation (BEC) has been realized in numerous dilute atomic gases . While a pure molecular BEC has not yet been achieved, ultra-cold molecules have been formed  and, just recently, a coherent superposition of an atomic and a molecular BEC was reported .

Collisions of atoms in a BEC can be controlled by the variation of a relatively weak magnetic field; the magnitude and sign of the scattering length can be changed by varying the field near a Feshbach resonance . In this way the expansion and collapse of an atomic BEC (termed “Bosenova”) has been observed . In collisions of molecules in a BEC, a change of the scattering length could be achieved with a magnetic field much weaker than that used to obtain the same effect in atoms. This is because there is an exponential dependence of the density of resonances on the number of “active” particles. This high sensitivity to weak fields could be used to search for a parity violating energy difference between chiral molecules.

It is well-known that biological molecules have a definite chiral structure (for example, there are only naturally occurring left-handed amino acids and right-handed sugars) . There have been numerous attempts to explain this effect by the influence of the parity violating weak interaction, which breaks the energy equivalence of right- and left-handed molecules . That parity violation can discriminate between molecules of different chirality is easily seen: a parity violating nuclear spin-independent electron-nucleus interaction in atoms  creates a spin helix of the electrons which interacts differently with right- and left-handed molecules. However, the parity violating energy difference $`\mathrm{\Delta }E_{PV}`$ is very small ,

$$\mathrm{\Delta }E_{PV}10^{20}Z^5\eta \mathrm{a}.\mathrm{u}.,$$
(1)
where $`Z`$ is the nuclear charge of the heaviest atom, and $`\eta `$ is an asymmetry factor which can be found from molecular structure calculations. This strong dependence on $`Z`$ originates from the weak ($`Z^3`$) and spin-orbit ($`Z^2`$) interactions. It may appear that in molecules with heavy atoms $`\mathrm{\Delta }E_{PV}`$ could become relatively large due to the $`Z^5`$-dependence. However, the geometrical suppresion factor $`\eta `$ remains very small. We should note that the effect may be larger for molecules with two heavy atoms . For recent calculations of $`\mathrm{\Delta }E_{PV}`$ for various molecules, see, e.g.,  and references therein.

Let us consider how a parity violating energy difference could manifest itself in the collision of two molecules in a BEC. Remember that in order to form a chiral molecule there must be at least four atoms involved ; therefore the collision of two diatomic molecules is sufficient. The distance between energy levels in a combined molecule $`D`$ is much smaller than the distance in an initial molecule, $`DD_i`$. Therefore, in an “ordinary” system it would be impossible to resolve compound resonances of the combined system since we would expect $`k_BTD`$, where $`k_B`$ is the Boltzmann constant, $`T`$ is the temperature. However, in a molecular BEC a finite fraction of molecules (up to 100% for $`TT_c`$) have zero kinetic energy, $`E=0`$ . This means that we have a possibility of perfect energy resolution . This allows one to study molecular collisions and chemical reactions in a unique situation where the reaction is dominated by the closest compound resonance of the combined system.

We can express the cross-section for formation of a chiral compound molecular state by the Breit-Wigner formula

$$\sigma =\frac{\pi }{k^2}\frac{\mathrm{\Gamma }_c\mathrm{\Gamma }}{(EE_0)^2+\mathrm{\Gamma }^2/4},$$
(2)
where $`k`$ is the wave vector, $`\mathrm{\Gamma }_c`$ is the capture width, and $`\mathrm{\Gamma }`$ is the total width of the resonance. The parity violating weak interaction in the chiral molecules shifts the resonance energies, for example, let’s consider that for the left-handed structure $`E_0E_L=E_0\mathrm{\Delta }E_{PV}/2`$ while for the right-handed structure $`E_0E_R=E_0+\mathrm{\Delta }E_{PV}/2`$. Therefore, cross-sections for the formation of left and right molecules, $`\sigma _L`$ and $`\sigma _R`$, from achiral components may be different. We can define an asymmetry parameter

$$P=\frac{\sigma _L\sigma _R}{\sigma _L+\sigma _R}.$$
(3)
The maximum value for $`P`$ is reached when $`E=E_0\pm \mathrm{\Gamma }/2`$. At this energy the asymmetry parameter

$$|P_{\mathrm{max}}|=\frac{\mathrm{\Delta }E_{PV}}{\mathrm{\Gamma }}.$$
(4)
The resonances can be shifted close to zero energy (the energy of the molecular collisions) $`E_0\pm \mathrm{\Gamma }/2=0`$ by application of an external electric or magnetic field.

The width of the level $`\mathrm{\Gamma }`$ may be quite small since the capture width $`\mathrm{\Gamma }_c=0`$ for energy $`E=0`$. The radiative width may be smaller than $`10^7\mathrm{Hz}`$. This can be compared to the weak energy shift. The largest values for $`\mathrm{\Delta }E_{PV}`$ that have been considered in molecular calculations are $`10^4\mathrm{Hz}`$ (e.g., for H2Po2 ).

Another possibility to observe parity violation in molecules may be related to the admixture of an $`s`$-wave to a $`p`$-wave compound resonance. For energy $`E=0`$, only $`s`$-wave molecules have a significant interaction cross section. Consider now a $`p`$-wave compound resonance. It seems to be invisible for $`E=0`$. However, the weak interaction $`W`$ mixes states of opposite parity and produces the combined state $`|\psi =|p+\beta |s`$, thus opening the $`s`$-wave reaction amplitude proportional to $`\beta `$. The mixing coefficient $`\beta =p|W|s/(E_sE_p)`$ is enhanced since the energy interval between the opposite parity compound states $`(E_sE_p)`$ is very small due to the high level density in a combined molecule. Note that this mechanism is responsible for the enhancement of weak interaction effects in neutron-nucleus reactions $`10^6`$ times  (for a review of the experiments, see ). Interference of the very small $`p`$-wave amplitude and the weak-induced $`s`$-wave amplitude may possibly also lead to a difference in the production of right and left molecules in the $`p`$-wave resonance.

This work was supported by the Australian Research Council.


text:
# 3D MHD Modeling of the Gaseous Structure of the Galaxy: Setup and Initial Results.

## 1 INTRODUCTION.

Even though spiral structure is one of the most prominent features of disk galaxies, details of the spiral arms in our own Galaxy remain uncertain. Georgelin & Georgelin (1976) traced the spiral structure of the Milky Way using H II regions, and developed a model with four arms. More recent attempts concluded that the Milky Way might actually have a superposition of two and four arm structures, each one with different pitch angles, which might arise from different components of the galactic disk (Drimmel, 2000; Lépine, Mishurov & Dedikov, 2001), suggesting that the stellar and gaseous disks might not be tightly coupled. Similar behavior has been frequently observed in external galaxies (Puerari & Dottori, 1992; Grosbøl & Patsis, 1998, for example).

Roberts (1969) showed that the gas must generate a large scale shock in the presence of a spiral perturbation. It was proposed that the density enhancement induced by this shock might generate a sequence of molecular clouds and star formation downstream from the shock, which itself was associated with the strong dust lane observed in the inner region of the spiral arms in external galaxies. Two dimensional numerical models by Tubbs (1980) and Soukup & Yuan (1981) showed that the gas forms a vertical shock perpendicular to the plane of the galactic disk. The post-shock gas remained close to hydrostatic equilibrium, even with an adiabatic equation of state. Their results did not show vertical motions larger than $`3\text{ km s}^1`$. In fact, the largest downflow they found was due to the pre-shock gas readjusting its vertical structure as it flows into the arm potential. Therefore, when H I observations on face-on galaxies showed extended velocity components with dispersions of the order of $`20\text{ km s}^1`$, they were attributed to other phenomena, such as galactic fountains, a warping of the H I disk, or intermediate velocity clouds (Dickey, Hanson & Helou, 1990; Kamphuis & Briggs, 1992; Kamphuis & Sancisi, 1993).

Since then, we have realized that the ISM is thicker and with a higher pressure than previously thought. The pressure scale height has been found to be larger than the density scale height, and the non-thermal pressures (turbulent, magnetic and cosmic ray) are at least as large as the thermal component (Badhwar & Stephens, 1977; Reynolds, 1989; Boulares & Cox, 1990). Therefore, less compressible gas needs to be considered in order to generate more realistic models of the ISM. Such a medium, with a larger effective $`\gamma `$ (the ratio of the specific heats) would be more likely to display the vertical motions characteristic of a hydraulic jump. With this in mind, Martos & Cox (1998, MC) performed 2D MHD simulations of the flow of the gaseous disk and found diverse structures that differed from the vertical near-hydrostatics found in previous studies. In many cases, the gas moved up ahead of the stellar arm, sped up over it, and fell behind with large bulk velocity. Frequently, there was a downstream shock at higher $`z`$ as this downflow was arrested, sometimes resulting in secondary midplane density maxima.

The goal of our investigation is to extend calculations like those of MC to three dimensions, to a large fraction of the Galaxy, and to look for its possible observational signatures. In this paper, we present the early results of these simulations. In Section 2 we describe the numerical setup and the procedure to achieve the initial hydrostatic equilibrium, in Section 3 we describe the results of the simulations, in Section 4 we present three examples of synthetic observations that can be done with this type of simulation, and in Section 5 we present our conclusions.

## 2 THE NUMERICAL SETUP.

We performed 3D MHD simulations in polar coordinates using the code ZEUS (Stone & Norman, 1992a, b; Stone, Mihalas & Norman, 1992). This code solves the ideal MHD equations for an inviscid fluid with infinite conductivity in a fixed eulerian grid. For our standard case, the grid extends from 0 to 1 kpc in $`z`$, 3 to 11 kpc in $`r`$, and 0 to $`2\pi /N`$ in the azimuthal angle $`\varphi `$ with 50, 80 and $`200/N`$ grid points in each respective direction, $`N`$ being the number of arms in each case. The boundaries in $`r`$ and $`z`$ are reflective, while those in $`\varphi `$ are periodic.

### 2.1 Hydrostatics, Theory.

The azimuthally averaged gravitational potential used is model 2 from Dehnen & Binney (1998). Using this potential, we set up a hydrostatic interstellar medium based on the scheme introduced by York et al. (1982) and described in the appendix of Benjamin (2002). We extended the procedure to include the effects of an azimuthal magnetic field whose magnitude depends only on the density. Given a density profile in the midplane, an equation of state (isothermal, in our case) and a density-magnetic pressure relation, the density and rotation velocity are uniquely defined everywhere.

Define the function $`G(\rho )`$ as:

$$G(\rho )=\frac{d}{d\rho }(p_T+p_B)\frac{d\rho }{\rho },$$
(1)
where $`\rho `$ is the density and $`p_T`$ and $`p_B`$ are the thermal and magnetic pressures, respectively. Vertical hydrostatics, $`(p_T+p_B)/z=\rho \mathrm{\Phi }/z`$ where $`\mathrm{\Phi }`$ is the gravitational potential, reduces via Equation 1 to:

$`{\displaystyle \frac{(G+\mathrm{\Phi })}{z}}`$ $`=`$ $`0`$ (2)
$`G[\rho (r,z)]`$ $`=`$ $`G[\rho (r,0)]+\mathrm{\Phi }(r,0)\mathrm{\Phi }(r,z).`$ (3)
At any $`z`$, the velocity profile is given by the radial balance between the radial potential and pressure gradients, the magnetic tension and the centrifugal force:

$$\frac{v_\varphi ^2(r,z)}{r}=\frac{\mathrm{\Phi }}{r}+\frac{1}{\rho (r,z)}\left[\frac{}{r}(p_T+p_B)+\frac{2p_B(r,z)}{r}\right].$$
(4)
This balance can be reduced to:

$$\frac{(G+\varphi )}{r}=\frac{v_\varphi ^22p_B/\rho }{r}=\frac{v_\varphi ^2v_A^2}{r},$$
(5)
where $`v_A`$ is the Alfvén velocity. Provided that $`^2(G+\mathrm{\Phi })/zr=^2(G+\mathrm{\Phi })/rz`$, along with Equations 2 and 5, we have:

$`{\displaystyle \frac{(v_\varphi ^2v_A^2)}{z}}`$ $`=`$ $`0`$ (6)
$`v_\varphi ^2(r,z)`$ $`=`$ $`v_\varphi ^2(r,0)v_A^2(r,0)+v_A^2(r,z).`$ (7)
Given the midplane density distribution, $`\rho (r,0)`$, the density off the plane is obtained by solving Equation 3 for $`\rho (r,z)`$, while $`v_\varphi (r,0)`$ is obtained from Equation 4 evaluated at $`z=0`$, after which the equilibrium rotation speed at all other $`z`$ follows from Equation 7.

### 2.2 Hydrostatics, Implementation.

Martos & Cox (1998) performed their 2D MHD calculations of spiral arm structure using the vertical density distribution at the solar circle compiled by Boulares & Cox (1990), modified to have a slightly lower vertical scale height for the warm ionized component. We found that the vertical distribution can be reproduced fairly accurately with thermal and magnetic pressures as follows. The thermal component assumes a neutral gas with a constant temperature of $`5700{}_{}{}^{\text{o}}\text{K}`$ and an isothermal equation of state. The magnetic pressure is taken as:

$$p_B=1.75\times 10^{12}\frac{n}{n+n_c}\text{ dyn}\text{ cm}^2,$$
(8)
where $`n_c=0.04\text{ cm}^3`$. The form of the magnetic pressure is such that it has little gradient at high density, and is proportional to $`\rho `$ at low density. The former accomodates a dense thermally supported core near the midplane, while the latter leads to a higher but constant signal speed at low density, far off the plane.

With a helium abundance equal to 10% of the hydrogen abundance by number, the mean atomic mass is $`(14/11)m_H`$.

In our initial work, we wanted to explore a situation that was not so heavily dominated by magnetic pressure. We therefore raised the temperature to $`10^4{}_{}{}^{\text{o}}\text{K}`$ and reduced the magnetic pressure by a factor of 10, keeping $`n_c=0.04\text{ cm}^3`$. The midplane density distribution was taken as exponential, with a radial scale length of 4 kpc and a density of $`1.11\text{ cm}^3`$ at $`r=8\text{ kpc}`$.

When these parameters are introduced into the above formalism and the hydrostatics found, the vertical half disk column density at $`r=8\text{ kpc}`$ is $`0.12\text{ kpc}\text{ cm}^3`$. Figure 1 shows the density, rotation velocity and magnetic field strength versus radius at the midplane, and versus $`z`$ at $`r=8\text{ kpc}`$. The midplane density varies by less than a factor of 10, while vertically the density drops nearly four orders of magnitude between the midplane and $`z=1\text{ kpc}`$, our present maximum height. At smaller radii, the vertical gravity is stronger and the density gradient in $`z`$ even larger, so that both the highest and lowest densities occur at the inner boundary. A curious feature is that at high $`z`$, the density increases with increasing radius: the disk “flares”.

The rotation velocity varies only slightly with radius, by less than $`15\text{ km s}^1`$, and by much less with $`z`$, only about $`1\text{ km s}^1`$ increase (as per Equation 7, higher Alfvén speed requires higher rotation rate at high $`z`$ in this approximation). The magnetic field strength varies only slowly with radius, and appears roughly Gaussian in $`z`$, with a flat region in the inner 200 pc of the thermal core.

We also report below on a case which has no magnetic field, in which the constant value of the temperature was taken as $`2.5\times 10^4{}_{}{}^{\text{o}}\text{K}`$. More precisely, the thermal pressure at a given density was taken as 2.5 times that of the previous run, because the above temperature was used inconsistently with assuming the gas was still neutral. This increase in thermal pressure was made in order to have a density distribution roughly similar to the magnetic case.

### 2.3 The Spiral Perturbation.

In addition to the axisymmetric potential, we used a spiral perturbation of fixed shape that rotates with $`\mathrm{\Omega }_P=12\text{ km s}^1\text{ kpc}^1`$ and has a pitch angle of $`15\mathrm{deg}`$. Details are reported in Cox & Gómez (2002). All the simulation grid is inside corotation. The depth of the perturbation varies slightly in $`r`$, weakens in $`z`$ and has a sinusoidal profile in $`\varphi `$; in the midplane, its corresponding mass density amplitude is $`52\%`$ of the disk component of the axisymmetric model at $`r=8\text{ kpc}`$, which provides a peak to valley potential difference of about $`(30\text{ km s}^1)^2`$. This mass contrast is consistent with K-band observations performed by Rix & Rieke (1993), Rix & Zaritsky (1995) and Kranz, Slyz & Rix (2001), who quote an arm/interarm contrast in density of old stars between 1.8 and 3 for a sample of spiral galaxies.

### 2.4 Numerical Complications.

In early runs, performed with outflow boundary conditions at the inner, outer and upper boundaries, ZEUS soon reported difficulties with “hot zones” in which the timestep became so short that the calculation terminated. This is almost certainly caused by the enormous density contrast of our hydrostatic solution. By making those three boundaries reflecting, so that the material is unable to flow off the grid, this problem was postponed or eliminated, depending on the case run. The two arm magnetized case ran to about 270 Myr, the unmagnetized case and the four arm magnetized case ran the full 400 Myr asked of them.

In addition to changing the boundary conditions, we also changed the perturbation force field near the inner boundary after noticing that the spiral potential was pushing material against the inner boundary, causing reflected waves that propagated outward. In order to avoid splashing against the inner-$`r`$ boundary, this perturbation is not applied in the inner 1 kpc of the grid, while it is smoothly turned on in the subsequent kpc. Thus, the useful computational grid runs from 5 to 11 kpc. Also, in order to diminish initial transient effects, the spiral perturbation is turned on gradually during the first 50 Myr. This short turn-on time undoubtedly creates part of the transient behavior and may have exaggerated some of the early velocity structure. Our intention is to make runs lasting so long that such transients have died out, and to report on that asymptotic behavior in future work.

We are wary of artifacts that might be caused by our boundary conditions and will continue to experiment with alternatives, including cases with lower overall density contrast that might allow open boundaries.

## 3 BEHAVIOR OF THE SIMULATIONS.

Our primary example is the two armed spiral with moderate (reduced from MC) magnetic pressure. Our calculation space was the upper half (in $`z`$) of half a circular disk, with periodic boundaries in $`\varphi `$. The period of time for a mass element to rotate around this half disk, relative to the rotating pattern, is about 100, 200 and 340 Myr ar $`r=5,8`$ and 10 kpc, respectively. We have chosen a fiducial time of 248 Myr, for our initial examination of the structure. At this time, nearly all mass elements have experienced the spiral perturbation once or twice, but conditions are still transient, representative of local interaction with rather than global accomodation to the perturbation. We will compare this early structure with that of the unmagnetized case at the same time, and the 4 arm magnetized case at half that time, which roughly represents the same level of maturity.

Having examined those single early time characteristics, we will present features of the subsequent development of the 2 arm cases and the much more mature 4 arm case, as indicative of features requiring a longer time to appear.

In the remainder of this paper, we will refer to each case using a three element naming scheme, describing the temperature (in units of $`10^4{}_{}{}^{\text{o}}\text{K}`$), the numerical coefficient in the magnetic pressure-density relation (Equation 8, in units of $`10^{12}\text{ dyn}\text{ cm}^2`$) and the number of spiral arms in the perturbation potential. Therefore, the magnetic two arm case will be denoted (1, 0.175, 2), the non-magnetic two arm case (2.5, 0, 2), and the four arm magnetic case (1, 0.175, 4).

### 3.1 Two Arm Magnetized Case.

Figures 2 and 3 show the results for the (1, 0.175, 2) case. The lower panel shows the density and (in Figure 3) velocity field along a surface of constant radius $`r=8\text{ kpc}`$. The upper panel shows column density for $`z0`$, half the total. In the upper panel, rotation is clockwise, and in Figure 3, the lines show the integrated velocity field at the midplane. For clarity in the visualization, we modified the components of the velocity in the lower panel of this Figure, so the arrows representing velocity in the inertial reference frame are parellel to the flowlines, with relative lengths proportional to the total velocity.

The most important feature in Figure 2 and 3 is the presence of a simple grand design spiral. The density concentration in the midplane contains most of the column density. At higher $`z`$, this feature leans forward. The midplane gaseous arm appears slightly after the perturbation potential minimum (outward in radius), shifting to a better alignment farther above the plane (Figure 4). Note also the strong arm to interarm contrast; a significant fraction of the material is located in the arms. The convergences of the velocity flow field in the upper panel of Figure 3 are consistent with this concentration.

Figure 5 shows the negative of the velocity divergence at the same positions as Figure 3. Only the regions with negative divergence are shown, in order to mark the places where strong compression (and shocks) appear. In the midplane, the compression of the material into the dense features is associated with a complex of intersecting shocks, following a diffuse compression of the material falling toward the arm. Comparison of Figures 3 and 5 shows an important shock 200 pc off the plane, preceding the gaseous arm, also slanting upstream at higher $`z`$. This shock accelerates gas upwards around $`\varphi /\pi 0.6`$. Over the arm, the flow is nearly horizontal. Behind the arm, it falls with vertical velocities of the order of $`20\text{ km s}^1`$, forming secondary shocks. This behavior is more evident along the $`\varphi `$-direction, but it is also visible in radial plots (not shown), since the presence of the stellar arms also induces velocities in that direction. Motions like these are similar to hydraulic jumps and were observed in the simulations performed by MC. MC also found midplane gas concentrations which they attributed to downstream bouncing of the flow. In our calculations, such concentrations are (so far) much smaller or absent in the midplane, but do appear in the lower density high $`z`$ gas, almost as if it tries to form another gaseous arm between the stellar arms. This interarm structure has little column density, but it is evident in the bottom panels of Figures 2 and 3, and the right hand panels of Figure 4.

In Figure 6 we present the vertical velocity structure at 310 pc above the plane. The upper panel shows a gray scale of $`v_z`$ with contours at 0 and $`\pm 20\text{ km s}^1`$. The gas moves up and over the gaseous structures at the arm and at the interarm positions, generating twice as many spirals in the upper panel of this Figure. Along the $`\varphi /\pi =0.5`$ direction (lower panel), this behavior is clearer. While the arm is at $`r8\text{ kpc}`$ in this direction, the vertical velocity behaves similarly at 6 and 10 kpc. Notice that in all cases, the transition from downflow to upflow is sharp, while the downturn of the gas is much smoother.

### 3.2 Two Arm Unmagnetized and Four Arm Cases.

As a comparision, Figures 7 and 8 show the (2.5, 0, 2) case at 248 Myr. Here, the temperature is higher in order to have a similar vertical density distribution. In the upper panel of Figure 7, we again observe the crowding of the midplane velocity field into the arms, as in Figure 3. The vertical structure now presents a strong shock at the leading edge of a nearly vertical gaseous arm and a much smoother density distribution behind it, which makes the gaseous arms fuzzier on the downstream side. There is vertical velocity structure, similar to the previous case, but with lower magnitude and only in the upper half of the grid.

At 124 Myr, the four arm case (1, 0.175, 4) in Figure 9 looks very similar to our standard case. Above 150 pc, strong downflows before the arms change suddenly to upflow at the arm position. In general, the column density is always smaller, since the same amount of mass is being distributed over 4 arms, leading to smaller arm to interarm contrast. This case is identical to the one presented in Figure 2, except for having 4 spiral arms perturbing the potential111As discussed in Cox & Gómez (2002), the scale height for the potential perturbation also depends on the number of arms involved. See also Martos et al. (2002). and that the grid covers only a quarter of the disk. We do not get interarm structures probably because there is not enough room between the arms for it.

Another difference between the four and two arm cases is that, for a given depth, the potential is steeper when four arms are present. Martos et al. (2002) developed a self-consistent model of the galactic spiral arms which resulted in narrower arms with a flat interarm region. They performed 1D MHD simulations with the usual sinusoidal potential and their modified perturbation, and found differences in the gaseous structures generated. The full 3D hydrodynamical effects of the details of the implementation of the perturbation potential will require further investigation.

A nice way to examine the phase between gaseous arms and the perturbing potential is by plotting column density in $`\varphi \mathrm{log}(r)`$ space, so that a logarithmic spiral appears as straight lines. Such plot is presented in Figure 10. Solid lines show contours of the underlying potential perturbation at the midplane, with the tick marks indicating the downhill direction. Dotted lines follow logarithmic spirals with pitch angle equal to the perturbation ($`15\mathrm{deg}`$), along the perturbation minima. Gray scale indicates column density in arbitrary units. As our model is of trailing spirals, the gas flows down from the top. In the (1, 0.175, 2) (upper panel), and (2.5, 0, 2) (middle panel) cases, the gaseous arms are slightly downstream from the potential minimum, by a gradually varying amount. The (1, 0.175, 4) case (bottom panel), on the other hand, follows a tighter spiral with a pitch angle of about $`12\mathrm{deg}`$. Also, the gaseous arms do not extend as far into the inner radii. If this difference in pitch turns out to be a robust feature of future simulations, it would lead to differences in the characteristics of the spiral structure of galaxies when comparing observations of Pop II with H I or Pop I tracers (Drimmel, 2000, for example).

### 3.3 Cyclic Variation.

The runs presented here do not yet correspond to an asymptotic state, but it appears that they already have some cyclic behavior in the their evolution. Figure 11 presents the time sequence for the evolution of the (1, 0.175, 2) case. Each panel corresponds to the lower panel in Figure 2, starting at t=0 with 8 Myr spacing. A high density structure above the arm is fully formed and leans upstream at $`120\text{ Myr}`$. It then contracts back as material from the tip falls down ($`t168\text{ Myr}`$), only to rise again from behind the arm, as in a breaking wave. Although it is hard to see in this Figure, such behavior is evident in an animation. The interarm structure (at the left side of the plots) also shows such cycles, with approximately the same phase. The reader must keep in mind that the gas making these structures is constantly moving around the galaxy and, therefore, these motions are really density waves on top of the galactic rotation. These cycles show a period of approximately 80 to 100 Myr. The (1, 0.175, 4) case presents a similar behavior, with a slightly longer period. The (2.5, 0, 2) case also presents such a cyclic variation, although it is less prominent. In all cases, the behavior is superposed on secular evolution in the arm structures.

### 3.4 Later Stages.

Figure 12 show the (1, 0.175, 4) case at 248 Myr. At this time, the gas has encountered the spiral arms twice as many times as in Figure 9. Fragmented arms are present inside $`r=7\text{ kpc}`$, while grand design arms are still present further out. A $`\varphi \text{ vs. }\mathrm{log}(r)`$ plot (Figure 13) shows the gaseous structures actually distributed along the potential maxima, instead of the minima. In the outer edges, the gaseous arms return to their previous position just downstream of the minimum. A time sequence of this case show that the gaseous arms drift downstream and stabilize at the potential crest. This evolution proceeds from the inner radii out. After a small section of the arms has drifted downstream, the arm breaks and the outer tip remains anchored at its original position. The result is that each individual section follows a tighter spiral, and the locus of all the sections, as a set, follow nearly the same spiral as the pertubation potential. Frequently, observations of the spiral structure in galaxies, even grand design spirals, show this type of feathering, in gaseous or Pop I tracers.

After 400 Myr, the (2.5, 0, 2) case also starts developing an interarm structure at small radii (Figure 14), in this case a bridge between arms. As this entire structure is interior to 5 kpc, the smallest radius at which our perturbation forces are fully activated, this structure may be an artifact. The potential for this formation is already evident from the velocity field in the upper panel of Figure 7.

## 4 SYNTHETIC OBSERVATIONS.

When we have refined our models and followed them into asymptotic behavior, we will examine signatures to compare with observational data. Here, we present some preliminary examples of how our simulated galaxy would look to observers from within, and from the outside.

Gas velocities inconsistent with a circular galactic rotation have been routinely observed (see for example Tripp, Sembach & Savage 1993). In order to study the line-of-sight component of the non-circular motions in our modeled galaxy, we assigned a reference position in the midplane of the disk and calculated the radial velocity from it of all other midplane locations. Upper panels on Figure 15 show that radial component, while the lower panels show the velocity difference between it and a purely azimuthal rotation. The two cases are for an observer in an interarm region and one at the inner edge of a gaseous arm. The circular velocity considered was the solution to our hydrostatics equations discussed in Section 2. The velocity differences frequently exceed $`20\text{ km s}^1`$, which would lead to errors in distance determinations of several kpc, even near the observer. Suggestion of problems in the determination of kinematic distances to pulsars have been found by Gómez et al. (2001). After tuning the input parameters in our simulation to a more realistic picture of the Milky Way, we will examine the degree to which those distance inconsistencies can be accounted for by considering non-circular gas motions.

Motions in the ISM are also observed in the form of intermediate velocity gas above the plane. In the simulations presented here, a large layer of gas falls behind the spiral arms (and sometimes, between them), all with similar speed. So, to an observer situated at the right position, this gas would appear as a falling velocity feature, a “cloud,” even without a localized density enhancement. Figure 16 shows this situation. After choosing a particular position in radius and azimuth, we interpolated the vertical distribution to a finer grid in order to generate a smooth artificial spectrum. All the spectra in the Figure correspond to $`r=8\text{ kpc}`$. Comparison with Figure 3 shows that the four azimuthal positions correspond to interarm, just before the arm, at the midplane density peak, and the downflow region after the arm. The two left hand panels show the characteristic upflow extending to about $`5\text{ km s}^1`$, while the two right hand panels (and to a lesser extent, the lower left one) present an extended wing towards larger negative velocities, together with a small peak at about $`18\text{ km s}^1`$. As seen in Figure 6, the asymmetry in the velocity distribution appears because the upflow happens in a narrower region of higher density and lower velocity than the downflow. Therefore, just by chance, it is easier to pick a region in which the gas seems to be falling than one with upflowing gas, and the velocities are likely higher as well.

Our third example shows the appearance of our galaxy model from outside, in velocity resolved spectra of lines originating in the low density material well off the plane. Vertical motion of the gas have been found in H$`\alpha `$ observations of NGC 5427 by Alfaro et al. (2001). These motions are consistent with what would be expected by corrugation in the velocity field of the gaseous disk induced by a hydraulic jump around the spiral arms. Our simulations do not include the ionization structure of the gas (and could not without including the ionizing agents) but it is reasonable to expect that the lowest density regions will be ionized and that their behavior will approximate that of the warm ionized component of the interstellar medium and have a direct relation with H$`\alpha `$ observations. With this in mind, we simulated the study done by Alfaro et al. (2001) in Figure 17 using our (1, 0.175, 2) case. The continous line is the integrated vertical velocity weighted by the square of the density,

$$w_z(r,\varphi )=\frac{_zv_z(r,\varphi ,z)\times [n(r,\varphi ,z)f(n)]^2}{_z[n(r,\varphi ,z)f(n)]^2},$$
(9)
where $`f(n)=1`$ if $`nn_{lim}`$, $`f(n)=(n/n_{lim})^{2.17}`$ otherwise, and $`n_{lim}=0.01\text{ cm}^3`$. The dotted line show the emission measure for these same grid points, defined as

$$EM=\underset{z}{}[n(r,\varphi ,z)f(n)]^2,$$
(10)
in arbitrary units. The midplane gaseous arm is at $`r=8.5\text{ kpc}`$. As observed by Alfaro et al. (2001), the gas moves up as it approaches the arm, falling behind it. A distinctive feature is that the EM peak is upstream from the gaseous arm at the midplane. Notice that the observations of NGC 5427 are for a region outside corotation, and therefore, the gas there moves from the convex to the concave sides of the gaseous arm, in the opposite direction from the gas in our simulation. In addition, the approximate $`30\mathrm{deg}`$ inclination angle for NGC 5427 should also show the imprint of radial streaming motions in the plane of the galaxy in this kind of observation.

Comparison of our Figures 3 and 17 with Figure 2 in Alfaro et al. (2001) shows that the H$`\alpha `$ profile of a spiral arm is more ragged than the model profile in higher density gas. While the column density plot in our Figure 3 is a smooth spiral, the H$`\alpha `$ emission measure in Figure 17 has many wiggles.

## 5 CONCLUSIONS AND FUTURE WORK.

In this work, we present our early results in the 3D MHD modeling of the large scale interaction of the ISM with a spiral potential. The presence of a thicker, more pressurized gaseous disk, together with the extra freedom the gas has in 3D simulations, allows the generation of density and velocity structures that previous work failed to reveal.

We confirmed and extended the work by MC, in which large scale vertical motions of the gas are an intrinsic feature of the response to the spiral perturbation. The downflow occurs along a much broader region and at higher velocities than the upflow. The falling gas can have large regions with a very similar vertical component, which translates into velocity crowding. In the present models, this gas appears as peaks at about $`20\text{ km s}^1`$ in spectra taken directly up from the midplane. These motions are accompanied with rapid flow above the arms and similar “up and over” motions in the radial direction. So far, there are some hints that such motions might be occuring in NGC 5427 (Alfaro et al., 2001).

We also found significant differences in the midplane line-of-sight velocity distribution as compared with a purely circular rotation model. We think that the presence of streaming motions generated by the spiral arms must be considered when estimating the distance to elements of the ISM using their velocity as reference. In the future, when we obtain a more realistic model for the Milky Way, we may be able to provide a reasonable recipe for translating radial velocities and galactic longitude data to distance in a more reliable way.

Our models have a number of numerical simplifications (low resolution, closed boundaries, short run times) and omission of physical processes (heating and cooling of the gas, cosmic rays, self gravity, ionization, star formation or associated energy injection). Improvement on the run times and resolution will allow us to follow the structures to maturity, better examine cyclic features, explore substructure formation such as feathers, bridging and gaseous interarms, follow the magnetic field energy density and geometry to saturation, and to explore radial migrations of material and angular momentum. Addition of a more realistic equation of state to represent heating and cooling of the gas will allow the formation of truly dense regions. The interaction between magnetized flow and these regions may qualitatively alter the general arm structure, the velocity field and the complexity of the magnetic field configuration.

Our results show that failure to consider high $`z`$ and non-circular motions of the ISM associated just with the response to the spiral potential can easily lead to confusion when interpreting observational results. The study of the gaseous structure of the Milky Way and other galaxies require the consideration of three dimensional effects and a more realistic model of the nature of the ISM and its interaction with other dynamical elements of the system.

We thank R. Benjamin, M. Martos and M. Bershady for useful comments and suggestions, to the NASA Astrophysics Theory Program for financial support under the grant NAG 5-8417, and to México’s Consejo Nacional de Ciencia y Tecnología for support to G. C. G.