text:
# ISO–SWS observations of pure rotational H2O absorption lines toward Orion–IRc2 Based on observations with ISO, an ESA project with instruments funded by ESA Member States (especially the PI countries: France, Germany, The Netherlands and the United Kingdom) and with the participation of ISAS and NASA.

## 1 Introduction

Water is one of the prime species for probing the interaction between young stars and their surroundings, both in terms of its abundance and its peculiar excitation. In high temperature gas, appropriate to shocks and hot cores for instance, all of the oxygen not locked up in CO is predicted to be driven into H2O at temperatures above 230 K by the O + H2 $``$ OH + H and OH + H2 $``$ H2O + H reactions, resulting in very bright H2O lines (e.g., Hollenbach & McKee 1979, Neufeld & Melnick 1987, Kaufman & Neufeld 1996, Charnley 1997). In contrast, the H2O abundance may be at least two orders of magnitude lower in surrounding colder gas (e.g., Zmuidzinas et al. 1995). In addition to collisions, the H2O excitation and line profiles can be strongly affected by mid- and far-infrared radiation from warm dust (e.g., Phillips et al. 1978, Takahashi et al. 1983, 1985), providing detailed information on the physical parameters of the gas and its location with respect to the radiation sources.

Interstellar water has been difficult to detect, apart from its presence as an ice mantle on dust grains or through its maser emission, owing to the severe telluric absorption encountered at Earth-based observatories. Nevertheless, over the last twenty years, many searches for lines of gas-phase H2O and its isotopomers have been made from the ground and airborne altitudes, in particular toward Orion (e.g., Waters et al. 1980, Phillips et al. 1978, Jacq et al. 1990, Wannier et al. 1991, Cernicharo et al. 1994, Zmuidzinas et al. 1995, Tauber et al. 1996, Timmermann et al. 1996, Gensheimer et al. 1996, Cernicharo et al. 1999a,b). However, due to the choice of line, wavelength and beam size, different components are often probed, whilst many ground-based observations refer to masing lines, for which sophisticated shock and maser models are required in order to extract physical parameters.

One of the major aims of the Infrared Space Observatory (ISO) mission has been the routine measurement of thermal gas-phase water lines, and their use as diagnostics of the chemical and physical conditions within molecular clouds. Far-infrared pure rotational H2O emission lines in the 50–200 $`\mu `$m wavelength range have been detected with the Long Wavelength Spectrometer (LWS) of ISO in a number of star-forming regions (e.g., Liseau et al. 1996, Ceccarelli et al. 1998), including Sgr B2 (Cernicharo et al. 1997a) and Orion (Cernicharo et al. 1997b, 1999a; Harwit et al. 1998). Additionally, van Dishoeck & Helmich (1996), van Dishoeck (1998) and Dartois et al. (1998) have observed absorption around 6 $`\mu `$m in the $`\nu _2`$=1–0 band toward a number of deeply embedded, massive young stars. Typical H2O abundances of $`10^5`$ with respect to H2 have been derived from these data. Similar observations have recently been reported for Orion BN/IRc2 by van Dishoeck et al. (1998) and González-Alfonso et al. (1998), although in this case emission is also detected. Wright et al. (1997) have however shown that the detection of the corresponding pure rotational lines at $``$ 30–200 $`\mu `$m in most sources observed at 6 $`\mu `$m is still difficult. The observations of Orion-IRc2 presented here form a notable exception.

Many of the earlier searches for H2O lines have been performed toward the BN/IRc2 complex of infrared sources in Orion, because of the extraordinary brightness of many atomic and molecular lines in this region compared with other clouds (e.g., Genzel & Stutzki 1989, Blake 1997). See van Dishoeck et al. (1998) and references therein for a detailed description of the geometry and the diverse range of physical phenomena in this region. In this paper we present the first detection of numerous pure rotational water lines in absorption toward IRc2 in the 25–45 $`\mu `$m interval with the Short Wavelength Spectrometer (SWS) (de Graauw et al. 1996) on board ISO. Some of these lines have been velocity resolved with the Fabry-Perot, enabling direct information on the location of the absorbing gas to be inferred. These data complement the earlier ground-based data, as well as observations of pure rotational lines with the LWS obtained by Cernicharo et al. (1997b, 1999a) and Harwit et al. (1998) in a much larger beam.

Our interpretation of the ISO spectra of Orion IRc2 suggests that the 25–45 $`\mu `$m H2O spectrum originates in a region where the intrinsically strong lines couple efficiently to an intense continuum. In principle, the formation of such a spectrum should be described for a stratified atmosphere in which lines and continuum are treated self-consistently. We show here that the observed features of the spectrum can be described well by a simple “generalised curve-of-growth”, which includes the effects of coupling to a strong continuum. In the following we describe our observations and present our results, followed by a discussion of the location of the absorbing water vapour, its excitation and finally its abundance.

## 2 Observations and data reduction

### 2.1 Observations

A complete grating scan from 2.4–45.3 $`\mu `$m using the maximum spectral resolution SWS06 observing mode was made on September 6 1997 (revolution 660) centered at $`\alpha (2000)=05^h35^m14.2^s`$, $`\delta (2000)=05^{}22^{}31.5^{\prime \prime }`$, which is about 3′′ W and 1′′ S of the IRc2 position quoted by Downes et al. (1981). The long axis of the ISO–SWS aperture was oriented at 172.81 E of N. The aperture size varies from $`14^{\prime \prime }\times 20^{\prime \prime }`$ at 2.4–12 $`\mu `$m, to $`14^{\prime \prime }\times 27^{\prime \prime }`$ at 12–27.5 $`\mu `$m, $`20^{\prime \prime }\times 27^{\prime \prime }`$ at 27.5–29 $`\mu `$m and $`20^{\prime \prime }\times 33^{\prime \prime }`$ at 29–45.2 $`\mu `$m. A further SWS06 scan was made from 26.3 to 45.2 $`\mu `$m in revolution 826 (February 19 1998) with the long axis oriented at 165.79 E of N. The SWS beam includes both IRc2 and BN, but not “peak 1” or “peak 2” of shocked H2 (Beckwith et al. 1978).

The resolving power of the ISO–SWS grating, $`\lambda /\mathrm{\Delta }\lambda `$, varies from about 1000–2500, implying that the observed line profiles are not resolved and that little velocity information can be obtained to help disentangle the various components which are known to exist toward IRc2. For this reason, follow-up observations of a selection of H2O lines were obtained using the Fabry-Perot SWS07 observing mode in revolutions 823 and 831 (February 16 and 24 1998). Five 250 km s-1 scans were made across each line. The resolving power is of order 30 000, corresponding to a velocity resolution of 10 km s-1, and the aperture long axis was oriented at 164.14 and 168.46 E of N on the two dates. The aperture size is $`10^{\prime \prime }\times 39^{\prime \prime }`$ and $`17^{\prime \prime }\times 40^{\prime \prime }`$ for wavelengths below and above 26 $`\mu `$m respectively. During the revolution 831 observation the wavelength interval between 5.3 and 7.0 $`\mu `$m was simultaneously scanned by the SWS grating. Within this region are several H2 lines (0–0 S(7), S(6) and S(5)), and H2O solid state ice and H2O gas phase $`\nu _2`$=1–0 features.

### 2.2 Data reduction

Data reduction was carried out on the Standard Processed Data file from the Off Line Processing system (OLPv6.1.1 for revolution 660 and OLPv6.3.2 for revolutions 823, 826 and 831), using standard routines within the SWS Interactive Analysis package up to the Auto Analysis Result (AAR) stage. The dark current subtraction was performed interactively for the grating data, but not for the FP data. The most up-to-date calibration files available for wavelength, flux and relative spectral responsivity were used. Even so, there is at least one instrumental artefact in our grating spectra, at 33 $`\mu `$m, due to structure in the relative spectral response calibration file (RSRF).

The high continuum flux encountered in band 4 of the SWS grating (29–45 $`\mu `$m), and the finite response time of the detectors, led to significant memory and response effects near the band edges, at positions of important water and OH lines, and so necessitated a non-standard up-down scan correction. For the SWS an “up” (“down”) scan proceeds toward increasing (decreasing) grating positions but decreasing (increasing) wavelengths. Since the up scan precedes the down scan, in high flux cases it is usually the down scan which is deemed to have the more reliable spectral shape, since by that time the instrument has settled down. Indeed, this was the case in the long-wavelength portion (41–45 $`\mu `$m) of our spectra, so that the up scan was corrected to approximately the same shape as the down scan. On the other hand, the short-wavelength portion (29–33 $`\mu `$m) of the down scan of our spectra had a peculiar concave shape, which did not match cleanly (in terms of shape) to band 3e, which in turn matched cleanly with band 3d. This was true for both the revolution 660 and 826 data sets. In this case therefore the down scan was corrected to have a similar shape to the up scan. The up-down correction obviously entailed splitting the band 4 data into two segments (29–37 $`\mu `$m and 37–45 $`\mu `$m), but this produced a more reliable spectral shape at both ends, e.g. allowing the continuum on either side of the 45.1 $`\mu `$m water line to be determined, and the 28.940 and 28.914 $`\mu `$m OH and H2O $`4_{40}3_{13}`$ lines to be partially resolved.

Following the AAR stage, further data processing, such as flatfielding, sigma clipping and rebinning, was performed using software in the SWS IA package. For instance, for both the grating and FP cases the individual detector scans have been brought to a common flux level using an $`n^{\mathrm{th}}`$ degree polynomial ($`n`$ = 0, 1, 2 or 3) fit to a reference spectrum which may either be the median of the down scans (e.g. grating band 3), the median of both the up and down scans (e.g. the up-down corrected grating band 4), or the mean of all scans (e.g. the Fabry-Perot data). The data have subsequently been sigma clipped such that any points lying more than 3 sigma outside of the average of all data within a bin of width equal to the theoretical spectral resolution have been discarded. The final step in producing a spectrum involves rebinning the much over-sampled data using the mean within a bin which is again equal to the theoretical resolution.

In the case of the 25.94 $`\mu `$m FP scans it was observed that glitches and their associated tails, which had not been flagged as bad data by the pipeline, contributed significant noise to the final spectrum. In this instance an interactive glitch and tail correction was performed before flatfielding, sigma clipping and rebinning. As noted by Schaeidt et al. (1996) and Heras (1997) the photometric accuracies of the SWS grating band 4 (where most of our water detections occur) and Fabry-Perot are 30% and 40% respectively. However, our observed line equivalent widths, being ratios of the line area over adjacent continuum, will be more accurate than this. The accuracy of the grating wavelength calibration is 1/10–1/5 of a resolution element, whilst that of the Fabry-Perot is of order 10-4 $`\mu `$m, or about 1 km s-1 (Valentijn et al. 1996; Feuchtgruber et al. 1997).

## 3 Results

In Fig. 1 the SWS grating spectrum in the 25–45 $`\mu `$m range is presented for the revolution 660 data, with the position and identification of features marked. The data consist of spectral segments from bands 3d, 3e and 4, for which “jumps” in the flux were observed due to the different aperture sizes. However, for display purposes small shifts, of order a few thousand Jy (i.e. a few percent), have been applied to produce a continuous spectrum. Clearly marked in Fig. 1 are 19 absorption lines arising from pure rotational levels in both the ortho- (o) and para- (p) species of H2O. In the few cases where the plot scale makes the line difficult to see in Fig. 1, we present close-ups in Fig. 2. The lines in this wavelength region originate from highly-excited states, with excitation energies of the lower levels $`hcE_{\mathrm{}}/k=200750`$ K, and the upper levels up to 1200 K. The strongest absorptions come from levels belonging to the backbone of the energy ladder. Among the highest levels seen is the upper state of the $`7_{25}6_{16}`$ transition at 29.8367 $`\mu `$m, involving the upper level ($`6_{16}`$) of the well-known 22 GHz H2O maser line (Cheung et al. 1969). Indeed, many of the levels from which absorption is seen are also upper levels of either known or suspected maser transitions (Neufeld & Melnick 1991). It is shown below based on velocity considerations that our lines probably arise in the same region where the masers are situated, indicating that our observations may form a useful data set for studying the maser pump mechanism.

A sub-set of 5 of the 19 lines was chosen to be observed with the Fabry-Perot. These spectra are shown in Fig. 3. In one case ($`4_{32}3_{03}`$ 40.69 $`\mu `$m) and possibly a second ($`5_{41}4_{32}`$ 43.89 $`\mu `$m), two velocity components are detected, one in absorption and the other in emission. A modulation in the baselines, due to the tracking of the Fabry-Perot transmission function across the region of the peak of the grating instrumental profile, has been modelled in terms of two parameters, the grating spectral resolution and the shift between the two response functions. The spectra shown in Fig. 3 have been corrected for this tracking pattern. Also shown in Fig. 3 is the grating spectrum from 5.3–7.0 $`\mu `$m obtained from these new data. It is consistent with that already observed by van Dishoeck et al. (1998) and González-Alfonso et al. (1998).

Table 1 presents an overview of our observations for the Revolution 660 data, including the rest wavelength, upper and lower levels, and equivalent widths $`W=cW_\lambda /\lambda `$ in units of Doppler velocity (km s-1), where $`W_\lambda `$ is the line equivalent width in wavelength units. The superscripts on $`W`$ refer to the observed and calculated equivalent widths, the latter derived from a generalised curve-of-growth analysis to be described in Sect. 4.2. Table 1 also contains the results for three pure rotational OH lines detected in absorption in our spectra. All derived quantities were calculated from unshifted spectral segments. Figure 4(a) presents our results graphically, in the form of a standard population diagram (Goldsmith & Langer 1999), where the column density in the lower level has been determined using the formula $`N_{\mathrm{}}=(8\pi cW_\lambda g_{\mathrm{}})/(\lambda ^4A_u\mathrm{}g_u)`$, where $`g_{\mathrm{}}`$ and $`g_u`$ are the statistical weights of the lower and upper levels, $`A_u\mathrm{}`$ is the Einstein $`A`$-coefficient in s-1 and $`W_\lambda `$ is the equivalent width. This formula is only strictly valid for optically thin lines, and under the assumption that the upper levels are not significantly populated. The error bars represent a 30% uncertainty, so that the scatter in the data is clearly larger than any systematic or statistical uncertainties. Although there is a clear trend in the data, this scatter already hints that a more detailed analysis is required to obtain a reliable total H2O column density. Such an analysis, in the form of a generalised curve-of-growth, is presented in Sect. 4.2.

Almost all lines listed in Table 1 at 27.5–45.2 $`\mu `$m are visible in both the Revolution 660 and 826 spectra, and their measured equivalent widths agree within 30% in the majority of cases. See Fig. 4(a), where both data sets are presented graphically. We present only the Revolution 660 data in Table 1 because the data from that revolution are of a better quality, and also only 15 of the 19 lines were observed in Revolution 826 due to the shorter wavelength scan range. The fact that the pure rotational H2O and OH lines observed with the SWS occur in absorption is in strong contrast with the lines seen with the LWS, which are mostly in emission (e.g., Harwit et al. 1998). Since the LWS lines arise in levels with similar excitation energies at wavelengths where the continuum flux is comparable to that in the SWS range, this suggests that most of the mid-infrared continuum and H2O absorption originate in a small region of size comparable to the SWS beam. In the almost order-of-magnitude larger area of the LWS beam, even more of the absorption is probably filled in by large-scale emission from more extended, lower excitation shocked gas. The qualitative change in the character of the water spectrum is thus probably a result of effects of both geometry (beam-filling) and radiative transfer.

The small differences in equivalent widths found from the grating and Fabry-Perot spectra are due to uncertainties in the calibration and varying beam shapes, which may couple differently to the complex source structure. The continuum radiation at 25–45 $`\mu `$m is produced by extended dust emission over a $``$ 30′′ region centered near IRc2, in contrast with the situation at $`6`$ $`\mu `$m where BN dominates (Wynn-Williams et al. 1984). Since H2O emission is observed on an arcmin scale with the LWS (e.g. Cernicharo et al. 1999a), it is likely that the H2O covers all of the continuum, at least for the lower-excitation lines. However, the absorption of the lower-lying lines may be partially filled in by emission within the SWS beam. As discussed below, these effects must be taken into account explicitly in an analysis of the equivalent widths.

## 4 Analysis and discussion

### 4.1 Spatial origin of the water

It is important to establish where the observed water is located along the line-of-sight towards IRc2. Possible candidates are the quiescent ridge, the shocked (low and/or high velocity) plateau gas, the hot core, or a combination of these components. Fortunately, our Fabry-Perot observations can aid in resolving the ambiguity. For the 3 pure absorption cases, the $`v_{\mathrm{LSR}}`$ is $`8`$$`\pm 3`$ km s-1. This blueshift is similar to velocities found for molecules such as CO, C2H2 and OCS of $`3`$ to $`18`$ km s-1 in absorption by Scoville et al. (1983) and Evans et al. (1991), which they interpret as arising from the low-velocity plateau gas, representing a shell expanding about a point close to IRc2. The most likely candidate for the location of our H2O gas is therefore the shocked, low-velocity plateau gas. Our lines are also resolved with observed full-width at half-peak $`\mathrm{\Delta }V`$ 30 km s-1. Assuming that the intrinsic source profile is Gaussian and that the FP instrumental profile is an ideal Airy profile with FWHM of 10 km s-1, the intrinsic line widths are of order 23 km s-1. This broad linewidth excludes most other components, such as the quiescent ridge and hot core, which have $`\mathrm{\Delta }V<20`$ km s-1 (Blake et al. 1987). H2O differs in this respect from the other molecules seen in absorption, which have $`\mathrm{\Delta }V`$ of order only 3–10 km s-1 (Scoville et al. 1983, Evans et al. 1991).

The peak of the emission component of the 40.69 $`\mu `$m line occurs at $`v_{\mathrm{LSR}}`$=+7 to $``$ +12 km s-1. This velocity range is similar to that found by Cernicharo et al. (1994) for the $`3_{1,3}2_{2,0}`$ 183 GHz line, as well as for a selection of H$`{}_{}{}^{18}{}_{2}{}^{}`$O lines by Zmuidzinas et al. (1995), and is in fact close to $`v_{\mathrm{LSR}}9`$ km s-1 appropriate for the dense, quiescent ridge in which IRc2 is embedded. The emission component is also resolved, with an observed linewidth of $`\mathrm{\Delta }V40`$ km s-1, consistent with the intrinsic widths of $``$ 35 km s-1 of the water lines detected in emission at longer wavelengths with the ISO–LWS by Harwit et al. (1998).

In the comparisons presented above the absorbing and emitting components have been considered separately. However, the profile as a whole for the 40.69 $`\mu `$m line, and to some extent also the 43.89 $`\mu `$m line, is very reminiscent of a P Cygni-type profile, and is similar to those observed for several far-infrared OH lines toward Orion by Betz & Boreiko (1989) and Melnick et al. (1990). In fact, the total blue-shift of the absorption component for all 5 lines, i.e. with respect to the rest velocity of the cloud $`v_{\mathrm{LSR}}9`$ km s-1, is very close to the expansion velocity of the “low-velocity flow”, (or plateau, or “expanding doughnut”) of 18$`\pm `$2 km s-1 quoted by Genzel et al. (1981) from proper motion studies of the 22 GHz H2O maser emission. Together with the fact that the emission component of the 40.69 $`\mu `$m line is near the rest velocity of the cloud, this is precisely the definition of a P-Cygni profile and establishes almost without a doubt that the observed water arises from an outflow centered near IRc2, most likely from source “I” (Menten & Reid 1995). As noted by Genzel et al. (1981) this outflow, or low velocity plateau, can be traced from $`2\times 10^{17}`$ cm (30′′) to within 1014 cm from its dynamical center.

We note that Takahashi et al. (1985) predicted that many H2O lines from a linearly expanding spherical cloud would exhibit P-Cygni type profiles, although the only case in common with our observations, the 40.69 $`\mu `$m line, is not P-Cygni shaped in any of their models. However, they do not include shock emission in their calculations. Further, there is a difference between a true P-Cygni profile and the P-Cygni type profiles observed here and predicted by Takahashi et al. A true P-Cygni profile results from an expanding stellar atmosphere with a central luminosity source, whilst the P-Cygni type profiles of far-infrared water lines are caused by the presence of warm dust throughout the expanding cloud. Further details on water line profiles in regions of intense dust emission may be found in Takahashi et al. (1985) or Doty & Neufeld (1997).

### 4.2 Excitation of the absorbing water vapour

#### 4.2.1 Generalised curve-of-growth method

The ISO spectra of Orion highlight the need for a non-traditional analysis of interstellar spectra. At shorter wavelengths (in the visible and ultraviolet), molecular absorption lines arise in gas that is physically very cold in relation to the radiation temperature of the background star, so that the absorption can be described with a classical curve-of-growth analysis and the effects of stimulated emission can be neglected. At longer wavelengths (radio), an isolated molecular region typically exhibits collisionally excited emission lines on top of a weak continuum so that the most significant radiative coupling is to the cosmic background radiation at $`T_{\mathrm{cbr}}=2.728`$ K. In the present case, we observe a molecular spectrum that is formed in the presence of a very strong continuum, i.e., where the radiation temperature of the continuum may be comparable to the excitation energies of the states involved. This implies both that the classical curve-of-growth analysis of cold absorbers detached from the background continuum is invalid and that the standard analysis of collisional excitation in competition with a very cold, dilute background is inadequate. Indeed, it is the near equality of continuum brightness temperature $`T_{\mathrm{rad}}(\nu )`$ and molecular excitation temperature $`T_{\mathrm{ex}}`$ that explains at least partly why the water spectrum of Orion IRc2 is predominantly in absorption at $`\lambda 50\mu `$m and in emission at longer wavelengths. In detail, such spectra should be modelled as an extended atmosphere in which lines and continuum are treated consistently. We will outline briefly how it is possible to extract important information on excitation and abundance through a simpler analysis of the unresolved lines based on a generalised curve-of-growth. Details of this method will be discussed elsewhere (Black, in preparation).

Consider the situation in which a fractional area $`a_c`$ of a continuum source of intensity $`I_c(\nu )`$ is obscured by a column of gas in the foreground. The observed spectrum in the vicinity of a transition $`u\mathrm{}`$ is given by the sum of the unattenuated continuum (first term), the absorbed continuum (second term), and the emission of the cloud itself (third term):

$`I_{\mathrm{obs}}(\nu )`$ $`=`$ $`(1a_c)I_c(\nu )`$ (1)
$`+`$ $`a_cI_c(\nu )e^{\tau _u\mathrm{}(\nu )}`$
$`+`$ $`B_\nu (T_u\mathrm{})\left(1e^{\tau _u\mathrm{}(\nu )}\right),`$
where $`\tau _u\mathrm{}`$ is the optical depth and $`B_\nu (T_u\mathrm{})`$ is the Planck function evaluated for an excitation temperature $`T_u\mathrm{}`$. This temperature is defined by the column densities of molecules in the upper ($`N_u`$) and lower $`(N_{\mathrm{}})`$ states through:

$$\frac{N_u}{N_{\mathrm{}}}=\frac{g_u}{g_{\mathrm{}}}\mathrm{exp}\left(\frac{h\nu }{kT_u\mathrm{}}\right).$$
$`(2)`$
The optical depth function for a single line can be written (e.g. Rybicki & Lightman 1979, Wannier et al. 1991)

$$\tau _u\mathrm{}(\nu )=3.738\times 10^7\frac{A_u\mathrm{}}{\stackrel{~}{\nu }_u\mathrm{}^3}\frac{N_{\mathrm{}}}{\mathrm{\Delta }V}\frac{g_u}{g_{\mathrm{}}}\left(1\mathrm{exp}\left(\frac{h\nu }{kT_u\mathrm{}}\right)\right)\varphi (\nu )$$
$`(3)`$
where $`N_{\mathrm{}}`$ is in cm-2, $`\stackrel{~}{\nu }_u\mathrm{}=\nu _u\mathrm{}/c`$ is the wavenumber of the line in cm-1, and the line shape function is a gaussian of full-width at half-maximum $`\mathrm{\Delta }V`$ in km s-1, normalized so that $`\varphi (\nu _u\mathrm{})=1.0`$ . The term in parentheses is the correction for stimulated emission.

The generalised equivalent width, in frequency units, is expressed as

$$W_\nu ^{\mathrm{obs}}=\frac{I_c(\nu )I_{\mathrm{obs}}(\nu )}{I_c(\nu )}𝑑\nu $$
$`(4a)`$
$$W_\nu ^{\mathrm{obs}}=\left[\left(a_c\frac{B_\nu (T_u\mathrm{})}{I_c(\nu )}\right)\left(1e^{\tau _u\mathrm{}(\nu )}\right)\right]𝑑\nu .$$
$`(4b)`$
In the limit as $`a_c1`$ and $`T_u\mathrm{}<<h\nu /k`$, the equivalent width approaches the classical form of cold absorbers detached from the background continuum. By convention, the equivalent width is positive for net absorption and negative for emission. Emission can arise either when $`T_u\mathrm{}`$ is large enough that $`B_\nu (T_u\mathrm{})/I_c(\nu )>a_c`$ or when there is a population inversion ($`N_u/N_{\mathrm{}}>g_u/g_{\mathrm{}}`$) so that $`\tau _u\mathrm{}`$ is negative, as in a maser.

For a single unresolved line of equivalent width $`W_\nu ^{\mathrm{obs}}`$, there is not a unique solution, because 4 parameters must be constrained: $`a_c`$, $`N_u`$, $`N_{\mathrm{}}`$, and $`\mathrm{\Delta }V`$. In the classic problem, in which a continuum point source of intensity is located behind, and completely obscured by, a column of cold foreground gas, and there is negligible population in the upper level of an absorption transition, then $`a_c=1`$ and $`N_u/N_{\mathrm{}}0`$ respectively. In that case, two lines arising in a common level suffice to determine $`N_{\mathrm{}}`$ and $`\mathrm{\Delta }V`$.

#### 4.2.2 Results of the generalised curve-of-growth analysis

We have analyzed the equivalent widths using Eq. (4b). If we take the continuum intensity to be approximately constant over the line profile, then the term $`a_cB_\nu (T_u\mathrm{})/I_c(\nu )`$ can be taken outside the integral in Eq. (4b). We adopt an empirical continuum intensity $`I_c(\nu )=f_\nu ^{\mathrm{obs}}/\mathrm{\Omega }_{\mathrm{beam}}`$, where $`f_\nu ^{\mathrm{obs}}`$ is the flux density observed in the ISO–SWS aperture of solid angle $`\mathrm{\Omega }_{\mathrm{beam}}`$. When $`I_c(\nu )`$ is fixed in this way, the equivalent width of a line depends on the parameters $`N_{\mathrm{}}`$, $`T_u\mathrm{}`$, $`\mathrm{\Delta }V`$, and $`a_c`$. We adopt an intrinsic width $`\mathrm{\Delta }V=23`$ km s-1 and assume that $`a_c`$ has the same value for all lines in the following analysis. Thus each line is characterized by two parameters $`N_{\mathrm{}}`$ and $`T_u\mathrm{}`$ (or $`N_{\mathrm{}}`$ and $`N_u`$). The solution is still underdetermined because there are 29 distinct energy states involved in the 19 observed transitions. Fortunately, there are several cases in which some of the same states are involved in two or more transitions. As a first approximation, it is assumed that a single value of the excitation temperature $`T_u\mathrm{}=T_{\mathrm{ex},0}`$ applies to all transitions with lower-state energies $`E_{\mathrm{}}275`$ cm-1. The computed equivalent widths (Eq. 4b) agree with the observed values within 33% for all but the weakest one of these lines, when $`a_c=1`$, $`T_{\mathrm{ex},0}=72.5`$ K, and $`N_{\mathrm{total}}=1.5\times 10^{18}`$ cm-2. This solution is not very sensitive to the value of the continuum covering factor: for $`a_c=0.9`$, we obtain $`T_{\mathrm{ex},0}=71.0`$ K and $`N_{\mathrm{total}}=1.6\times 10^{18}`$ cm-2, while for $`a_c=0.8`$, the corresponding results are 69.1 K and $`1.7\times 10^{18}`$ cm-2.

With the populations of the lowest energy states ($`E_{\mathrm{}}275`$ cm-1) fixed at $`T_{\mathrm{ex},0}=72.5`$ K, the column densities in higher states are adjusted until the standard deviation of the mean $`W^{\mathrm{obs}}W^{\mathrm{GCOG}}`$ is minimized for the entire set of lines. In this solution, the observed $`W^{\mathrm{obs}}`$ and calculated $`W^{\mathrm{GCOG}}`$ equivalent widths all agree within 16%, with the single exception of the $`4_{41}4_{14}`$ line at 37.98 $`\mu `$m. The calculated equivalent widths are listed in Table 1. Note that the adjustment also improves the fit for the lines arising in the low-lying states because $`W^{\mathrm{GCOG}}`$ is sensitive to column densities in both the upper and lower states. The result is a population distribution that is indistinguishable from a thermal distribution at 72.5 K for all states with $`hcE_i/k700`$ K. A few of the observed levels at higher energies ($`hcE_i/k>700`$ K) are well constrained because they are the upper states of more than one transition, and their populations start to clearly show effects of subthermal excitation, $`T_u\mathrm{}55`$ K. Even though the adopted linewidth is large, the derived line-center optical depths, listed in Table 1, reach values as high as $`\tau _0=12.8`$ and 9.7 in the 45.11 and 40.69 $`\mu `$m lines, respectively. The smallest values of optical depth are $`\tau _00.06`$. Finally, ortho and para states of H2O are consistent with the same excitation temperature and with an ortho/para ratio of 3.

The results of our generalised curve-of-growth analysis are shown graphically in Fig. 4(b), in the form of a standard population diagram. The much reduced scatter in the data is clearly a marked improvement over the optically thin approach assuming no stimulated emission used in obtaining Fig. 4(a). It is however interesting that the inferred excitation temperatures from the two approaches are similar. The excitation temperature of $``$ 72 K is similar to the colour temperature of the dust between 20 and 100 $`\mu `$m inferred by Werner et al. (1976), which confirms the expectation that pumping by infrared continuum radiation from dust plays an important role in the H2O excitation. In this case, there is no direct constraint on the kinetic temperature of the H2O-containing gas. Subthermal excitation by collisions at densities below the critical densities of the observed levels ($`10^9`$ cm-3) can also result in excitation temperatures of order 50–100 K, but this merely re-affirms that radiative processes probably dominate the excitation. We also note that the excitation temperature inferred from the emission components of the 40.69 and 43.89 $`\mu `$m lines is $``$ 64 K, also implying a large contribution by radiative processes to the populations. Moreover, the analysis suggests that $`T_u\mathrm{}`$ of the absorption lines is so close to the continuum radiation temperature in this wavelength range, that it is easy to see why the H2O spectrum should go into emission at $`\lambda >50\mu `$m as the opacity and radiation temperature of the continuum decrease with increasing wavelength.

Because the lowest state involved in the 25–45 $`\mu `$m spectrum lies at $`hcE_i/k=197`$ K, we cannot determine the column densities in the most populated states directly from these observations. However, if we assume that the observed excitation pattern applies to the unobserved lower states, then we infer a total column density of $`N`$(H2O) = $`1.5\times 10^{18}`$ cm-2 averaged over the ISO–SWS beam toward Orion IRc2.

The net absorption by water in Orion is sensitive to the competition between the molecular excitation and the continuum brightness. We have assumed that the true continuum intensity in Eq. (4) is equal to the mean surface brightness observed in the spectrum (flux averaged over the aperture). In general this need not be true, and the analysis can include an additional correction factor for the coupling to the continuum. Tests show that if the continuum intensity is thus changed by a factor of two in either direction, the derived total column density and mean excitation temperature change by less than 20%, which we take as our maximum uncertainty.

#### 4.2.3 Shock models

We have also compared our observed line equivalent widths with those computed for the Kaufman & Neufeld (1996) C-shock in Orion with a pre-shock density of 105 cm-3 and shock velocity of 37 km s-1 (Kaufman & Neufeld, private communication) (see Table 1). The water abundance (with respect to H2) used in this shock calculation was 3.5$`\times `$10-4, based on the O and C abundances of Cardelli et al. (1996), with all the carbon locked up in CO and H2O accounting for all the remaining oxygen. Although this model refers to “peak 1” rather than the IRc2 region, it reproduces the equivalent widths quite well within a factor of a few, except for some of the higher-lying lines (e.g., 27.03, 39.37, 40.76, 44.19 $`\mu `$m). Most of these lines are likely to be strongly affected by radiative excitation, which was not included in their model.

We note that the shock velocity used is 37 km s-1, different from the expansion velocity inferred from our FP observations of 18 km s-1. Indeed, our modelling of our H2 emission, using the Kaufman & Neufeld (1996) shock model, indicates that at least 2 shock components are required. One component at 15–20 km s-1 matches well the observed column densities in the $`v`$=0, $`J`$=3 to 7 levels, whilst a second component, at 35 km s-1, is required to match higher $`J`$ levels and the ro-vibrational lines (Wright 1999). Harwit et al. (1998) used a shock velocity of 37 km s-1 to match their 8 LWS water line emission detections from 72 to 125 $`\mu `$m. These data were taken at a slightly different position, centred on BN rather than IRc2 and offset 8.2′′ to the north. The SWS results may therefore be more sensitive to the low-velocity plateau gas identified by several authors (e.g. Genzel & Stutzki 1989), whereas the larger LWS beam will include a significant contribution from the high-velocity gas. We also note that water line emission in the Kaufman & Neufeld (1996) shock models is for many lines insensitive to the shock velocity, especially when the velocity is $``$ 15–20 km s-1. This is presumably because in such cases the water abundance has reached its maximum value and that the gas temperature is well above the upper levels involved. Therefore, the emission depends more on density and not temperature (which in turn is highly sensitive to the shock velocity). Even so, future attempts to match both the Harwit et al. lines and our detections with a slower shock, but including radiative excitation, will be worthwhile.

A long standing problem in shock research is the contribution of water to the total gas cooling behind the shock, and it is pertinent to mention it here. Unfortunately, the situation in Orion IRc2, i.e. a high water abundance (see Sect. 4.3) coupled with an intense mid- and far-infrared radiation field, has not previously been considered in shock models. However, the detection of all our transitions in absorption raises the interesting possibility that at least these transitions are not cooling the gas, but rather heating it. Such heating results from the water molecules absorbing the radiation and thereby being raised to excited states. The excess energy may then be imparted to the gas through collisional de-excitation with other molecules, e.g. H2, as in the scenario proposed by Takahashi et al. (1983). This would of course be in competition with radiative de-excitation (i.e. cooling), and is therefore density dependent. Further, Takahashi et al. (1983) and Neufeld et al. (1995) state that a necessary condition for H2O heating is that the dust temperature is greater than the gas temperature. A full scale shock model, taking into account all transitions, would be required to determine if there was net heating or cooling. We merely mention it here as an interesting aside, and direct the reader to the paper by Harwit et al. (1998) for a full discussion of the water cooling contribution in Orion. Briefly though, from their larger beam LWS observations they find that H2O and CO contribute similar amounts to the overall cooling, but both are only about a tenth of the total H2 cooling. On the other hand, Saraceno et al. (1999), in a study of a sample of shock sources, find that CO cooling typically dominates over H2O by a factor of several.

### 4.3 H2O abundance in Orion

In order to determine an H2O abundance directly from our data, information on the H2 column density is needed. As noted above, from the H2 lines detected elsewhere in our spectra (van Dishoeck et al. 1998), the shock has been characterized using the Kaufman & Neufeld (1996) C-shock models, in a manner similar to that applied by Wright et al. (1996) to Cepheus A. Further details can be found in Wright (1999); in short, the H2 pure rotational lines originating from $`J`$=3 to $`J`$=7 indicate a column $`N`$(H2)=$`(1.2\pm 0.2)\times 10^{21}`$ cm-2 of warm ($`T_{ex}`$=700 K) gas, and the inferred shock velocity of 15–20 km s-1 is similar to the 18 km s-1 expansion velocity deduced from our FP observations, but lower than that obtained from modeling the H2 vibration-rotation lines (35 km s-1). However, supporting evidence for the co-location of the H2O and shocked H2 comes from Scoville et al. (1982), who find that at and around the position of IRc2 the peak emission of the 1–0 S(1) line is at $`v_{\mathrm{LSR}}`$=$`6`$ to $`16`$ km s -1, similar to the velocity of maximum H2O absorption in our spectra.

Using the above H2 column density, and that of H2O calculated from the generalised curve-of-growth analysis, the inferred water abundance would be $`1.25\times 10^3`$. However, the inferred column density from the high-$`J`$ (i.e. 3–7) lines is a lower limit on the total column density, since the bulk of the H2 is likely to be cooler than $`T_{ex}`$=700 K. Such cooler gas would be probed by the 28.2 $`\mu `$m H2 0–0 S(0) line, but which unfortunately has not been detected toward IRc2, due to the extremely high continuum emission. However, a lower limit implies that $`T_{ex}110`$ K, calculated from the populations in $`J`$=3 and 2, inferred from the observed 0–0 S(1) and upper limit S(0) lines respectively. On the other hand, the S(0) line is detected at about the 3$`\sigma `$ level toward the shock “peak 2”, $``$ 30′′ to the south-east of IRc2, and in this case the excitation temperature between the $`J`$=3 and $`J`$=2 levels is $``$ 150 K (Wright 1999).

This temperature of $``$150 K is likely to be indicative of the kinetic temperature, $`T_{kin}`$, at this position, since the density of $``$105 cm-3 (e.g. Genzel & Stutzki 1989) is above the critical density of the H2 0–0 S(0) and S(1) transitions (e.g. Le Bourlot et al. 1999). Although this estimate of $`T_{kin}`$ is at a different position than where the water absorption line detections were made, Wright (1999) shows that in at least 3 shock sources, including Orion, the excitation temperature determined from pure rotational H2 lines ($`J_{up}`$=3–7) is quite invariant across the sources, as is the column density. Other estimates of $`T_{kin}`$ in the low-velocity plateau, tabulated by Genzel & Stutzki (1989), range from 100–500 K.

If we assume that $`T_{kin}`$ is between 100–500 K, then Table 2 presents the expected total H2 column density based on the observed surface brightness of the 0–0 S(1) line toward IRc2 of 1.5$`\times 10^{10}`$ W cm-2 sr-1. We believe that 100 K is too low to be applicable to the shocked plateau gas, since it is lower than $`T_{kin}`$ we find at “peak 2”, and in any case would imply a column density an order of magnitude larger than the value of 1023 cm-2 quoted by Genzel & Stutzki (1989). Also shown in Table 2 is the inferred water abundance, $`N(\mathrm{H}_2\mathrm{O})/N(\mathrm{H}_2)`$, using the observed water column density of 1.5$`\times `$1018 cm-2, and assuming the water and molecular hydrogen to be co-located. Taking the maximum possible water abundance to be $`5\times 10^4`$ (anything higher would violate the constraint imposed by the oxygen abundance itself of \[O\]/\[H\]=3.19$`\times 10^4`$ from Meyer et al. 1998), then the water abundance is in the range of 2.6–50$`\times 10^5`$.

Let us instead assume a priori that our water detections arise only in that portion of the shock where the temperature is high enough for the water creation reactions to rapidly proceed. Several authors (e.g. Hollenbach & McKee 1979, Neufeld et al. 1995, Ceccarelli et al. 1997 and Charnley 1997) have shown that water begins to be quickly formed when the temperature exceeds about 230 K, and has reached its maximum abundance by about 320 K. In that case our inferred water abundance is in the range 2–5$`\times 10^4`$, in remarkable agreement with the result of Harwit et al. (1998) of 3.5–5$`\times `$10-4 considering the different methods used. Further, it is interesting to note that Watson et al. (1985) find a H2 column density of $`3\times 10^{21}`$ cm-2, derived from several observed far-infrared CO lines in a 44′′ beam, implying a H2O abundance of 5$`\times `$10-4. Watson et al. however infer a kinetic temperature of 750 K, close to the excitation temperature derived from our H2 data (Wright 1999), whereas from Table 2 their column density would instead imply a temperature of only $``$330 K.

Alternatively, if we simply use a H2 column density of 1023 cm-2 for the plateau gas, as given by Genzel & Stutzki (1989), then the water abundance is $`1.5\times 10^5`$, in better agreement with that of other shocks (e.g. Liseau et al. 1996, Ceccarelli et al. 1998, Spinoglio et al. 1999), and suggesting that a substantial fraction of the gas-phase oxygen may still be in another form, most likely gas-phase atomic O i. However, this is probably a strict lower limit, since such a large column density probably includes gas at a temperature lower than 230 K, as well as gas behind the infrared continuum source.

An H2O abundance at the high end of the range quoted above, i.e. 2–5$`\times `$10-4, is favored if we assume that water is only efficiently created in regions where the temperature is greater than about 230 K. This water abundance is consistent with all of the gas-phase oxygen not locked up in CO being incorporated into H2O. However, H2O ice evaporates at temperatures above about 90 K, so lower temperatures, and hence a lower abundance, cannot be excluded. The major uncertainty on our abundance value is the total hydrogen column density. As mentioned above, our result is in good agreement with Harwit et al. (1998), but also in agreement with Cernicharo et al. (1999a,b) from LWS and ground based observations respectively. From radiative transfer modelling of their water detections, they find toward the Orion plateau a water abundance of order 1–2$`\times `$10-4, at most a factor of a few lower than what we find. Besides Orion, there is only one other source so far observed by ISO with such a high water abundance, namely L1448-mm by Nisini et al. (1999), with an abundance of order $`5\times 10^4`$.

A similar treatment for the OH lines detected in our spectrum yields a column density of approximately $`5\times 10^{16}`$ cm-2, after fixing Tex to be the same as that for water. The implied OH abundance is then $`(8.8156)\times 10^7`$ for Tkin=150–320 K, or $`(6.215.6)\times 10^6`$ for 230–320 K.

## 5 Conclusions

Through use of the ISO–SWS in its grating mode, 19 pure rotational absorption lines of water have been detected toward Orion IRc2 for the first time. Fabry-Perot spectra of 5 lines reveal that the water is located in an outflow expanding at a velocity of 18 km s-1. The strong mid-infrared continuum toward IRc2 plays a dominant role in the excitation of the molecule and the line formation, which can be modeled using a simple, generalised curve-of-growth technique. This yields a total water column density of order 1.5$`\times 10^{18}`$ cm-2 and excitation temperature of 72 K, similar to the dust continuum colour temperature. Both derived quantities have a maximum uncertainty of about 20%. The data provide support for large abundances of H2O in the outflows of massive stars. Simultaneous analysis of the complete ISO–SWS and LWS data set on H2O, OH and CO may provide further information on the abundance and excitation of these molecules in the various physical components within the complex Orion environment.

###### Acknowledgements.

The authors are grateful to the SWS instrument teams in Groningen and Garching and to the SIDT in Vilspa for making these observations possible. They are indebted to G.A. Blake, A. Boonman, F. van der Tak, G.J. Melnick, A.G.G.M. Tielens and R. Timmermann for useful discussions. They are especially grateful to D.A. Neufeld and M.J. Kaufman for providing them with the equivalent widths of the H2O lines in their models, and to W.F. Thi for developing the code for the up-down scan correction. This work was supported by NWO grant 614.41.003 and by the Spanish DGES grant PB96-0883 and PNIE grant ESP97-1618-E. During the final phases of this work CMW was supported by an ARC Australian Postdoctoral Research Fellowship.


format: text

text:
# The fully frustrated XY model with next nearest neighbor interaction

## Abstract

We introduce a fully frustrated XY model with nearest neighbor (nn) and next nearest neighbor (nnn) couplings which can be realized in Josephson junction arrays. We study the phase diagram for $`0x1`$ ($`x`$ is the ratio between nnn and nn couplings). When $`x<1/\sqrt{2}`$ an Ising and a Berezinskii-Kosterlitz-Thouless transitions are present. Both critical temperatures decrease with increasing $`x`$. For $`x>1/\sqrt{2}`$ the array undergoes a sequence of two transitions. On raising the temperature first the two sublattices decouple from each other and then, at higher temperatures, each sublattice becomes disorderd.

A variety of two dimensional systems undergo a phase transition without a rigorous symmetry breaking, the Berezinskii-Kosterlitz-Thouless (BKT) transition . The transition is driven by the thermally excited vortices which form a two dimensional Coulomb gas . Josephson junction arrays are experimental realization of the two-dimensional XY model where the array’s parameters can be modified in a controlled way. In the last decade there has been a great amount of work on the various aspects of the BKT transition in Josephson arrays . Experimental studies are based on electrical resistance , two-coil inductance , and SQUID  measurements.

A magnetic field applied perpendicularly to the array leads to frustration . If the flux piercing the elementary plaquette is half of the flux quantum $`\mathrm{\Phi }_0=hc/2e`$, the system is called fully frustrated (FF) and undergoes two phase transitions related to the $`Z_2`$ and $`U(1)`$ symmetries. The existence of the two critical temperatures $`T_c^{Z_2}`$ and $`T_c^{U(1)}`$ respectively has been extensively investigated both by analytical methods  and Monte Carlo (MC) simulations . The complete scenario is not fully understood yet. There are numerical evidences either supporting the existence of two very close critical temperatures ($`T_c^{Z_2}>T_c^{U(1)}`$) with critical behavior typical of Ising and BKT transitions respectively  or the existence of a single transition with novel critical behavior .

Untill recently only the Josephson arrays with nearest neighbor (nn) couplings were studied. Recent theoretical  and experimental  works on infinite range array opened a new field of investigation in these systems. In this Letter we study the properties of a two dimensional FF Josephson array with both nn and next-to-nearest neighbor (nnn) couplings . Proximity-junction arrays may be good candidates to experimentally probe the effects discussed in this work. They consist of superconducting islands in good electric contact with a metallic substrate. Due to the proximity effect there is a leakage of Cooper pairs into the normal substrate which extends over a temperature-dependent coherence length $`\xi _N`$ ($`\xi _N=\mathrm{}v_F/k_BT`$ or $`\xi _N=\sqrt{\mathrm{}D/2\pi k_BT}`$ for the ballistic and for the diffusive case respectively, $`v_F`$ is the Fermi velocity, $`D`$ the diffusion constant and $`T`$ the temperature). When $`\xi _N`$ becomes comparable with the lattice constant of the array the nnn coupling becomes comparable with the nn coupling. Since $`\xi _N`$ is strongly dependent on the temperature, the nnn Josephson coupling may be observed cooling down the sample. The main results of this Letter are summarized in the phase diagram of Fig. 5.

The system is defined by the Hamiltonian

$`H={\displaystyle \underset{<<i,j>>}{}}J_{ij}\mathrm{cos}(\theta _i\theta _jA_{ij})`$ (1)
where $`J_{ij}=J>0`$ for nn and $`J_{ij}=xJ`$ for nnn ($`x0`$), the symbol $`<<\mathrm{}>>`$ refers to the sum over nn and nnn. For later convenience we introduce the gauge invariant phase difference $`\varphi _{ij}=\theta _i\theta _jA_{ij}`$. The variables $`\theta _i`$ are the phases of the superconducting order parameter of the i-th island. The magnetic field enters through $`A_{ij}=(2\pi /\mathrm{\Phi }_0)_i^j𝐀𝑑𝐥`$ ($`𝐀`$ is the vector potential). The relevant parameter which describes the magnetic frustration is $`f=A_{ij}/(2\pi )`$, where the summation runs over the perimeter of the elementary plaquette. We study the case $`f=1/2`$ on a square lattice.

Ground States \- The model of Eq.( 1) combines the characteristics of both the FF and unfrustrated XY models. While the elementary square plaquette is FF, the square plaquettes formed by nnn couplings are not frustrated. For $`x<x_01/\sqrt{2}`$ we find that the ground state of this system is exactly the same as in the FF model without nnn couplings and is characterized by $`\varphi _{ij}=\pm \pi /4`$ for all pairs of nn sites. This state combines continuous \[$`U(1)`$\] and discrete ($`Z_2`$) degeneracies. Its energy per site $`E=\sqrt{2}J`$ is independent of nnn coupling $`xJ`$. Moreover if one considers a straight domain wall separating the two ground states with opposite orientations of chiralities it turns out that neither the form of such state nor its energy change with addition of nnn interactions.

For $`x>x_0`$ the ground state is the same as in the absence of nn coupling when the system splits into two unfrustrated $`XY`$-models. The relative phase shift between the two sublattices can be arbitrary. The energy of this state $`E=2xJ`$ depends only on nnn coupling $`xJ`$ and does not depend on nn coupling.

At the special point $`x=x_0`$ the energies of both the above ground states coincide. Morover they can be transformed into each other by a continuous transformation without increasing the energy of the system. Therefore the manifold of the ground states also includes an additional set of eight-sublattice ”intermediate” states which can be parametrised by a rotation angle $`\chi `$ ($`\chi =\pi `$ corresponding to low-$`x`$ ground state and $`\chi =0`$ corresponding to high-$`x`$ ground state with a particular relative phase shift between the sublattices) as it is shown in Fig. 1.

Phase Diagram \- We studied the finite temperature behaviour of the model by means of the (low temperature) spin-wave free energy analysis and by Monte Carlo simulations.

For $`x<x_0`$ neither the spectrum of spin waves (in the long wavelength limit) nor the domain wall energy depend on $`x`$. Therefore we can expect only a weak dependence of $`T_c^{Z_2}`$ and $`T_c^{U(1)}`$ on $`x`$, due to the change of the effective interaction between the different types of fluctuations.

For $`x>x_0`$ the system (at finite temperatures) turnes out to be equivalent to two coupled XY models, the effective forth-order coupling between the two sublattices provided by the free energy of spin waves. Although this coupling is weak (always much smaller than the temperature) at low temperatures it is relevant and imposes the presence of a transition at $`T=T_D`$. This transition separates the phases with coupled and decoupled sublattices. In the low temperature phase, where the two sublattices are locked, the spin-wave contribution to the free energy imposes a relative phase shift of $`\pm \pi /4`$ (or equivalently $`\pm 3\pi /4`$) between the two sublattices. At $`T>T_D`$ a second phase transition of the BKT type takes place in each of the decoupled sublattices. For $`xx_0`$ the temperature of this transition depends only on nnn coupling and it is proportional to $`x`$.

The spin wave spectrum remains rigid down to $`x=x_0`$. This indicates (and it is confirmed by the MC simulations) that the critical temperature of this BKT transition remain finite when $`xx_0^+`$. Below this temperature there is a transition between the FF low-$`x`$ and the unfrustrated high-$`x`$ phases. We evaluated numerically the spin wave free energy of the intermediated ground state as a function of $`\chi `$. This dependence is described by a convex function implying a first order phase trasition line separating the low-$`x`$ and high-$`x`$ quasi-ordered phases.

The critical properties were investigated evaluating the helicity modulus $`\mathrm{\Gamma }`$ and the staggered chiral magnetization $`M`$ by means of standard MC simulations for different $`x`$. The order parameter $`M`$, which controls the Ising-like transition, is defined as

$$M=\frac{1}{L^2}\left|\underset{i}{}(1)^{i_x+i_y}m_i\right|$$
(2)
where $`\stackrel{}{r}_i/a=(i_x,i_y)`$ is the position vector (in unit of lattice step $`a`$) of the site $`i`$ and $`m_i=\frac{1}{\sqrt{8}}(\mathrm{sin}\varphi _{i,j_1}+\mathrm{sin}\varphi _{j_1,j_2}+\mathrm{sin}\varphi _{j_2,j_3}+\mathrm{sin}\varphi _{j_3,i})`$ is the chirality of the plaquette with center in $`(i_x+1/2,i_y+1/2)`$ and with site indexes $`i`$,$`j_1`$,$`j_2`$,$`j_3`$ (in clockwise order).

The helicity modulus $`\mathrm{\Gamma }=^2/\delta ^2`$, used to signal the existence of a BKT transition, is defined through the increase of the free energy $``$ due to a phase twist $`\delta `$ imposed in one direction . In order to obtain a precise determination of $`T_c^{Z_2}`$ and of the critical exponent $`\nu `$ associated to the divergence of the correlation length we have calculated the Binder’s cumulant  of the staggered chiral magnetization $`M`$

$$U_M=1\frac{M^4}{3M^2^2}.$$
(3)
Since $`U(T_c^{Z_2},L)`$ does not depend on lattice size $`L`$ for large systems, $`T_c^{Z_2}`$ can be identified without making any assumption on the critical exponents. Once a satisfactory estimation of $`T_c^{Z_2}`$ is obtained the critical exponent $`\nu `$ is estimated through a data collapsing with $`\nu `$ left as the only free parameter. Estimation of $`U_M`$ have been obtained averaging, at least, $`10^7L^2`$ MC configurations by using a standard Metropolis algorithm. The largest lattice studied is $`L=72`$. The result for $`x=0.5`$ are shown in Figs.2,3a. We estimate $`k_BT_c^{Z_2}/J=0.403\pm 0.003`$. The data collapsing, shown in the inset of Fig.2, gives an estimate of $`1/\nu =1.0\pm 0.1`$. The critical temperature $`T_c^{Z_2}`$ decreases with increasing $`x`$ for $`x<x_0`$. For $`xx_0`$, there is no sign of a Ising-like transition (see Fig.3b).

Following the procedure proposed in Refs. , the critical temperature $`T_c^{U(1)}`$ is estimated by using the following ansatz for the size dependence of $`\mathrm{\Gamma }`$

$$\frac{\pi \mathrm{\Gamma }}{2T_c}=\gamma \left[1+\frac{1}{2(\mathrm{ln}L\mathrm{ln}l_0)}\right]$$
(4)
where $`l_0`$ is a fit parameter. This critical scaling is based on the mapping between a neutral Coulomb gas and a XY model. Therefore Eq.(4) can be used both as a test for the existence of a BKT transition and for a precise evaluation of the critical temperature. A very good scaling was obtained with $`\gamma =1`$ (the ordinary BKT transition) in the low-$`x`$ phase and $`\gamma =2`$ (corresponding to BKT transition on each of the two sublattices with lattice constant $`\sqrt{2}`$) in the high-$`x`$ phase. In Fig.4 we show this analysis for the cases $`x=0.5`$ and $`x=1`$. $`T_c^{U(1)}(x)`$, as well as $`T_c^{Z_2}(x)`$, decreases with increasing $`x`$ up to $`xx_0`$. Our results cannot discriminate between the $`T_c^{Z_2}=T_c^{U(1)}`$ and the $`T_c^{Z_2}>T_c^{U(1)}`$ hypothesis since the two temperatures are compatible within the numerical precision (the mean value of $`T_c^{Z_2}(x)`$ remains always above the corresponding mean value of $`T_c^{U(1)}(x)`$). For $`x0.8`$, instead, $`T_c^{U(1)}(x)`$ increases, quickly tending towards the value expected for $`x\mathrm{}`$ i.e. $`k_BT_c^{U(1)}/(xJ)=0.89`$.

We finally discuss the transition related to the decoupling of sublattices in the high-$`x`$ phase. The order parameter $`S=_{j,\alpha }s_{j,j+e_\alpha }`$ can be defined on the bonds of the lattice and can be chosen in the following gauge-invariant form:

$`s_{j,j+e_x}`$ $`=`$ $`(1)^{j_x}\mathrm{exp}\left[i(1)^{j_x+j_y}\varphi _{j,j+e_x}\right]`$ (5)
$`s_{j,j+e_y}`$ $`=`$ $`i(1)^{j_y}\mathrm{exp}\left[i(1)^{j_x+j_y}\varphi _{j,j+e_y}\right]`$ (6)
This form of the order parameter is chosen in such way that for any of the high-$`x`$ ground states the value of $`s`$ will be the same for all the bonds. At low temperatures $`S`$ manifests a true long-range order. The MC simulations confirm the $`\pi /4`$ relative phase shift anticipated by the spin-wave analysis. By increasing the temperature the long-range order in $`S`$ can be expected to disappear as a separate phase transition whereas the unbinding of vortex pairs in each of the sublattices has to occur at still higher temperatures. We performed MC simulations to evaluate the transition temperature $`T_D`$. The results of this computation (not reported here) show that $`T_D`$ is very close, but lower, than the transition to the disordered phase $`T_c^{U(1)}(x)`$. More extensive simulations are needed to determine the critical behaviour of decoupling transition. The dotted line in Fig.5 shows the qualitative behaviour of the decoupling transition as a function of $`x`$.

In conclusion we have introduced a frustrated XY model with nnn interaction. The model can be experimentally realized in Josephson junction arrays in a transverse magnetic field. Signatures of nnn Josephson couplings might already have been seen in specially designed setups . The analysis presented here leads to the phase diagram shown in Fig.5. For $`0<x<x_0`$ the critical temperatures $`T_c^{Z_2}`$ and $`T_c^{U(1)}(x)`$ decrease with increasing $`x`$. For $`x>x_0`$, there is no sign of a Ising-like transition and the system behaves like the unfrustrated XY model. At $`xx_0`$ in the low $`T`$ region there is a first order phase transition between the low and high-$`x`$ phases. Finally for $`x>x_0`$ the array undergoes a sequence of two transitions. On raising the temperature, first the two sublattices become decoupled and then, at higher temperatures, each sublattice becomes disordered.

We thank H. Courtois and B. Pannetier for fruitful discussions. We acknowledge the financial support of INFM (PRA-QTMD) and the European Community (Contract FMRX-CT-97-0143). G.F. is grateful to A. Mastellone and J. Siewert for hospitality in Catania.


text:
# Bose-Einstein condensation in inhomogeneous Josephson arrays

## Abstract

We show that spatial Bose-Einstein condensation of non-interacting bosons occurs in dimension $`d<2`$ over discrete structures with inhomogeneous topology and with no need of external confining potentials. Josephson junction arrays provide a physical realization of this mechanism. The topological origin of the phenomenon may open the way to the engineering of quantum devices based on Bose-Einstein condensation. The comb array, which embodies all the relevant features of this effect, is studied in detail.

today

The recent impressive experimental demonstration of Bose-Einstein Condensation (BEC)  has stimulated a new wealth of theoretical work aimed to better understanding its basic mechanisms  and, possibly, to exploit its consequences for the engineering of quantum devices.

It is well known  that for an ideal gas of Bose particles BEC does not occur in dimension $`d2`$, and an ′′ad hoc$`^{\prime \prime }`$ external confining potential is needed to reach the required density of states. The same is true for free bosons living on regular periodic lattices, while the result cannot be extended to more general discrete structures lacking translational invariance.

In the following we shall prove that even for $`d<2`$  non-interacting bosons may lead to Bose-Einstein condensation into a single non-degenerate state, provided one resorts to a suitable discrete non-homogeneous support structure: indeed, when the bosonic kinetic degrees of freedom do not depend on metric features only, the particles may feel a sort of effective interaction due to topology. The proposed mechanism for BEC in lower dimensional systems is then a pure effect of the structure of the ambient space and avoids as well the need of resorting to external random potentials as the ones investigated by Huang in ; this is a very desirable feature in view of engineering real quantum devices.

In practice, the behavior of free bosons over generic discrete structures is made experimentally accessible through the realization of suitable arrays of Josephson junctions. The latter are devices that can be engineered in such a way as to realize a variety of non-homogeneous patterns. We shall show indeed that classical Josephson junction arrays arranged in a non-homogeneous geometry - not even necessarily planar - provide an example of the proposed mechanism for BEC, leading to a single state spatial condensation.

Theoretical studies of Josephson junction arrays are based on the short-range Bose-Hubbard model, since the phase diagram of Josephson junction arrays may be derived  from an Hamiltonian describing bosons with repulsive interactions over a lattice. In $`d=1`$ the phase diagram has been studied by analytical  and quantum Monte Carlo methods ; experimentally, Josephson junction arrays are used to study interacting bosons in one dimension. For a generic array the corresponding Hamiltonian is given by

$`H^{BH}=U{\displaystyle \underset{i}{}}n_i^2+{\displaystyle \underset{ij}{}}A_{ij}\left(Vn_in_jJ(a_i^{}a_j+a_j^{}a_i)\right),`$ (1)
where $`A_{ij}`$ is the adjacency matrix: $`A_{ij}=1`$ if the sites $`i`$ and $`j`$ are nearest neighbors and $`A_{ij}=0`$ otherwise; $`a_i^{}`$ creates a boson at site $`i`$ and $`n_ia_i^{}a_i`$. The phase diagram structure reflects the competition between the boson kinetic energy (hopping, favouring boson mobility) and repulsive interaction (Coulomb, working so as to suppress dynamics). In a realistic experimental setup , the parameters $`U`$ and $`V`$ depend on the ratio between the intergrain capacitance $`C`$ and the gate capacitance $`C_0`$, while the parameter $`J`$ describes Cooper pair hopping. Josephson junction arrays allow for a good experimental control of $`C/C_0`$ and $`J`$, which can also be varied over a wide range. For $`UJ,V`$, and for bipartite arrays, Hamiltonian $`H^{BH}`$ maps onto the quantum spin-$`\frac{1}{2}`$ XXZ model . On the other hand, in the weak coupling limit (′′classical$`^{\prime \prime }`$ Josephson junctions) $`U,VJ`$, realizable when $`C/C_00`$, the hopping term dominates the physical behaviour of the system, which is then described by the tight-binding Hamiltonian

$$H=t\underset{ij}{}A_{ij}a_i^{}a_j,$$
(2)
where $`t`$ is an effective hopping parameter which accounts for a renormalization of the Josephson coupling $`J`$. For a non translation-invariant geometry of the array the tight-binding model, which describes ′′free$`^{\prime \prime }`$ bosons over a regular lattice, cannot any longer be interpreted as representing non-interacting particles, just due to the ambient graph topology. We shall show the dramatic effect of topology already on the simple graph referred to as square comb . This provides an explicit and remarkable example of topology-induced mechanism leading to a spatial Bose-Einstein condensation in low dimension.

The ′′square comb$`^{\prime \prime }`$ is the graph made of $`N`$ ′′fingers$`^{\prime \prime }`$ of $`N`$ sites represented in Fig. 1. whose total number of sites is $`N^2`$. In the following the generic vertex $`i`$ is labelled with the ′′coordinate$`^{\prime \prime }`$ indices $`(x_i,y_i),i𝐙_N`$ (where the latter requirement is imposed to guarantee periodic boundary conditions).

Fig. 1. The comb graph: the dots represent Josephson junctions and the links describe the topology of the array connections, with no reference to their embedding in real euclidean space.

The tight-binding model on the comb graph, is given by the Hamiltonian (2) with the adjacency matrix $`A_{ij}=\left(\delta _{x_i,x_i+1}+\delta _{x_i,x_i1}\right)\delta _{y_i,0}\delta _{y_j,0}+\left(\delta _{y_i,y_i+1}+\delta _{y_i,y_i1}\right)\delta _{x_i,x_j}`$. By exploiting the comb translation invariance in the direction of the backbone one can perform a Fourier Transform along the $`x`$ direction obtaining a new Hamiltonian in the variables $`k_x`$:

$$H=t\underset{y,y^{}}{}\underset{k_x}{}(\delta _{y,0}\delta _{y,y^{}}\mathrm{cos}k_x+\overline{A}_{yy^{}})a_{k_x,y}^{}a_{k_x,y^{}},$$
(3)
where $`k_x=2\pi n/N`$ and $`n=0,1,,N1`$ and $`\overline{A}_{yy^{}}`$ in the adjacency matrix for a linear chain (i.e. for each comb finger). The operator $`a_{k_x,y}`$ is given by: $`a_{k_x,y}_{x=1}^Ne^{ik_xx}a_{(x,y)}`$. Notice that (3) is the sum of $`N`$ commuting Hamiltonians representing an one-dimensional tight-binding model with a local potential at site $`0`$ of value: $`t\mathrm{cos}k_x`$.

Each of the one-dimensional Hamiltonians appearing in eq. (3) can be diagonalized, in the thermodynamic limit. To do this one uses the property that for $`y2`$ and $`y2`$ the eigenvectors of $`(\delta _{y,0}\delta _{y,y^{}}\mathrm{cos}k_x+\overline{A}_{yy^{}})`$ are those of the Hamiltonian describing a free particle on the linear chain. Such eigenvectors are: $`e^{\pm ik_yy}`$ with eigenvalue $`2\mathrm{cos}k_y`$, $`e^{\pm k_yy}`$ with eigenvalue $`2\mathrm{cosh}k_y`$ and $`(1)^ye^{\pm k_yy}`$ with eigenvalue $`2\mathrm{cosh}k_y`$. If one requires that the eigenvalue equations hold also at sites $`1`$, $`0`$,$`1`$ and that the eigenvectors are normalizable, one finds that the spectrum of the one dimensional problem is given by the isolated point $`E=2t\mathrm{sgn}(\mathrm{cos}k_x)\sqrt{1+\mathrm{cos}^2k_x}`$ and by a continuous part $`R_0=\{E||E|<2t\}`$ with a density of states given by $`\rho _0(E)=(1/\pi )(N1)(4t^2E^2)^{1/2}`$. From the spectra of the one dimensional problems corresponding to different values of $`k_x`$ one obtains in the thermodynamic limit the density of states for the tight-binding Hamiltonian on the comb-graph (Fig. 2).

Fig. 2. The density of states of the Hamiltonian (2). The solid line indicates the continuous part of the spectrum $`\rho _0`$, which is normalized to $`N(N1)`$. The dot lines denote the sets of zero measure, the densities $`\rho _+`$ and $`\rho _{}`$ are normalized to $`N`$. The $`x`$-axis scale is in units of $`t`$.

The spectrum is made of three parts $`R_\alpha ,\alpha =0,\pm `$ : for $`ER_0`$ the density of states is given by:

$`\rho _0(E)={\displaystyle \frac{1}{\pi }}N(N1)(4t^2E^2)^{1/2}.`$ (4)
Since $`_{R_0}\rho _0(E)𝑑E=N(N1)`$, in the thermodynamic limit, almost all the states, i.e. all the states apart a set of measure zero, belong to this region and $`lim_N\mathrm{}(_{R_0}\rho _0(E)𝑑E)/N^2=1`$. In the other two regions, $`R_{}=\{E|\sqrt{8}tE<2t\}`$ and $`R_+=\{E|2t<E\sqrt{8}t\}`$, the density of states is given by:

$`\rho _{}(E)=\rho _+(E)\rho _\pm (E)=hN{\displaystyle \frac{|E|}{\sqrt{8t^2E^2}\sqrt{E^24t^2}}},`$ (5)
where $`h`$ must be chosen so that $`_{R_{}R_+}\rho _\pm (E)𝑑E=N`$. There follows that in the thermodynamic limit only a subset of states of measure zero belongs to these regions of the spectrum, in that one has $`lim_N\mathrm{}N^2_{R_\pm }\rho _\pm (E)𝑑E=0`$. The states with $`ER_{}`$ play a fundamental role in the study of bosonic particles on comb structures. The lowest energy eigenstate of (2) (corresponding to $`E_0=\sqrt{8}t`$) is represented in Fig. 3. It is constant along the $`x`$ direction for any fixed $`y`$, while along $`y`$ it decreases exponentially with the distance from the backbone.

Fig. 3. The eigenvector corresponding to the lowest energy state. It is constant along the backbone direction and it decreases exponentially as $`\mathrm{exp}(\mathrm{arcsh}(2)|x|)`$ along the fingers.

If one introduces a finite bosonic filling $`f`$, i.e. if one fills the comb with $`fN^2`$ non interacting bosons, fixing the number of particles in the gran canonical partition function amounts to choosing the fugacity $`z(N,\beta ,f)`$ as

$$f=N^2\underset{\alpha =,0,+}{}_{R_\alpha }\frac{1}{z^1e^{\beta E}1}\rho _\alpha (E)𝑑E,$$
(6)
with $`0<ze^{\beta t\sqrt{8}}`$. Since $`_{R_+}(z^1e^{\beta E}1)^1\rho _+(E)𝑑E<cN`$, ($`c`$ is a number independent of $`N`$ and $`z`$), the third term of (6) vanishes in the thermodynamic limit. On the other hand, the first term can be positive and finite in the thermodynamic limit if $`z(N,\beta ,f)e^{\beta t\sqrt{8}}`$ when $`N\mathrm{}`$. Denoting by $`n_0`$ the fraction of particles with energy smaller than $`2t`$, $`n_0=lim_N\mathrm{}f^1N^2_R(z^1(N)e^{\beta E}1)^1\rho _{}(E)𝑑E`$, one has that $`n_0`$ is a finite fraction of particle condensed in a subset of states of measure zero. In the thermodynamic limit, if $`\beta _c`$ is the inverse temperature for which:

$`1={\displaystyle \frac{1}{\pi f}}{\displaystyle _{R_0}}{\displaystyle \frac{dE}{(e^{\beta _c(E+t\sqrt{8})}1)\sqrt{4t^2E^2}}},`$ (7)
one has that if $`\beta \beta _c`$ there always exixts a real and positive $`z`$ solution of the equation

$`f={\displaystyle \frac{1}{\pi }}{\displaystyle _{R_0}}{\displaystyle \frac{dE}{(z^1e^{\beta E}1)\sqrt{4t^2E^2}}},`$ (8)
and $`n_0=0`$: there is no condensation. For $`\beta >\beta _c`$ one has $`z=e^{\beta t\sqrt{8}}`$ and $`n_0`$ is given by

$`n_0=1{\displaystyle \frac{1}{f\pi }}{\displaystyle _{R_0}}{\displaystyle \frac{dE}{(e^{\beta (E+t\sqrt{8})}1)\sqrt{4t^2E^2}}}.`$ (9)
For $`T<T_c(\beta _c)^1`$ there is Bose-Einstein condensation and in Fig. 4 we plot the fraction $`n_0`$ of particles in the condensate as a function of $`T`$ for several values of the filling ($`f=0.5,1,2`$). The points where the curves intersect the T-axis are the critical temperatures for the different fillings.

Fig. 4. The fraction $`n_0`$ of condensed particles as a function of the temperature $`T`$ for different fillings $`f=0.5,1,2`$.

For $`TT_c`$ the order parameter depends linearly on $`T`$ since $`n_0(TT_c)`$, and it is analytic in $`T`$ for $`T`$ close to $`T_c`$. The behaviour for $`T0`$ is given by $`n_0\sqrt{T}\mathrm{exp}((2\sqrt{2}2)t/KT)`$. This is different from the customary Bose-Einstein condensation in a 3-$`d`$ box, where $`n_0=1(T/T_c)^{3/2}`$. The critical temperature $`T_c`$ exhibits the dependence on the filling $`f`$ shown in Fig. 5, asymptotically ($`f>1`$) linear and of the form $`T_c\left(\mathrm{ln}f\right)^1`$ for $`f1`$. Both behaviours characterize a gapped system. The gap $`(2\sqrt{2}2)t`$ measures the difference between the ground state energy and the bottom of the spectral region $`R_0`$

Fig. 5. The critical temperature $`T_c`$ as a function of the filling $`f`$.

The average energy per particle

$`<E>=\sqrt{8}tn_0+{\displaystyle \frac{1}{f\pi }}{\displaystyle _{R_0}}{\displaystyle \frac{EdE}{(z^1e^{\beta E}1)\sqrt{4t^2E^2}}},`$ (10)
shows that all the particles in the condensate have energy $`t\sqrt{8}`$, namely they are all in the ground state. Summarizing, in the thermodynamic limit and for $`T>T_c`$ almost all particles have energies between $`2t`$ and $`2t`$ with the distribution $`(z^1e^{\beta E}1)^1\rho _0(E)`$, while for $`T<T_c`$ a finite fraction $`n_0`$ of particles is condensed in the state of lowest energy $`E_0=\sqrt{8}t`$. These particles occupy the inhomogeneous state described in Fig. 3, in which the sites closer to the backbone have larger filling than the farther ones.

In conclusion, we exhibited an explicit example (comb array) of Bose-Einstein condensation into a single state for non-interacting bosons, induced in dimension $`d<2`$ by inhomogeneities without disorder and with no confining external potential in the discrete geometrical structure. The ensuing condensate shows deconfinement in one direction, i.e. along the comb backbone, and localization along the orthogonal direction; this is expected to lead to detectable singularities in the response functions.

The model Hamiltonian used is physically implementable by classical Josephson junction arrays, which it is possible to engineer in any desired geometric setting. Thus Bose-Einstein condensates arise as an intrinsic device feature, without need of fine-tuning any external control parameter.

Furthermore, it has been evidenced that a comb graph may emerge as a relevant geometrical structure if one consider a chain of classical Josephson junctions aligned along the direction of the backbone and interacting with a suitable external environment .

Finally, the same devices, as recently proposed , might lend themselves to be used for the realization of BEC-based encoding and manipulation of quantum information. The idea here is that because of the spontaneous symmetry breaking that characterizes it, a Bose-Einstein condensate should be quite naturally described by a non-linear quantum mechanics. Such non-linearity can be thought of as due just to the effective interaction of the bosons composing the condensate. In such a case the scenario recently described by Abrams and Lloyd whereby non-linear quantum mechanics in the sense of Weinberg implies polynomial time efficiency in dealing with NP-complete and $`\mathrm{}`$P complex computations would hold. It is intriguing that for the system presented here the relevant non-linearity, indeed present, is due to geometry and topology rather than to physical interactions.