text:
# Sulphur chemistry in the envelopes of massive young stars

## 1 Introduction

Spectral line surveys at submillimeter wavelengths have revealed considerable chemical differences between star-forming regions (e.g., Blake et al. 1987; Sutton et al. 1995; Schilke et al. 1997; Helmich & van Dishoeck 1997; Hatchell et al. 1998a; Nummelin et al. 2000; see van Dishoeck 2001 for a complete list). While these differences indicate activity, the dependence of molecular abundances on evolutionary state and physical parameters is poorly understood. Better insight into this relation would be valuable for probing the earliest, deeply embedded phases of star formation where diagnostics at optical and near-infrared wavelengths are unavailable. This is especially true for the formation of high-mass stars, for which the order in which phenomena occur is much less well understood than in the low-mass case, and for which the embedded phase is a significant fraction ($``$10%) of the total lifetime of $``$106 yr.

Sulphur-bearing molecules are attractive as candidate tracers of early protostellar evolution (Hatchell et al., 1998b; Buckle & Fuller, 2003). In models of “hot cores” (Charnley, 1997), where the chemistry is driven by the evaporation of icy grain mantles due to protostellar heating, the abundances of sulphur-bearing molecules exhibit a characteristic variation with time. After its release from the grain mantles, H2S is converted into SO and SO2 in $``$103 yr, and further into CS, H2CS and OCS after $``$105 yr. This type of model has had some success in reproducing the chemistry of star-forming regions (Langer et al., 2000).

Besides thermal heating, shocks are potentially relevant to the chemistry in star-forming regions. Young stars drive outflows which expose their surroundings to the effects of shocks. The submillimeter lines of SO2 in Orion are very broad, and likely arise in the outflow (Schilke et al., 1997). Models of chemistry behind interstellar shocks (Mitchell, 1984; Leen & Graff, 1988; Pineau des Forêts et al., 1993) predict enhancements of H2S, SO and SO2 on time scales of $``$104 yr. Thus, both by ice evaporation and by shock interaction, sulphur could act as a clock on the time scale of $`10^4`$ yr relevant for the embedded phase of star formation.

This paper uses submillimeter observations of sulphur-bearing molecules to test models of chemistry during star formation, and attempts to order these regions chronologically. The sources are nine regions of high-mass star formation with luminosities $`1\times 10^42\times 10^5`$ L at distances of $`14`$ kpc (Table 1). Originally selected for their mid-infrared brightness, the sources were found to have large envelope masses (40-1100 M within $`r=0.150.36`$ pc) based on submillimeter data (van der Tak et al., 2000b). Together with their weak radio emission and their strong molecular outflows, these masses indicate very early evolutionary stages, probably preceding the formation of hot cores. Two of our sources may be more evolved than the others: NGC 6334 IRS1, which has a very rich, ‘hot core’-type submillimeter spectrum (McCutcheon et al. 2000; Thorwirth et al., in prep.), and NGC 7538 IRS1, which has strong free-free emission.

Observations of CH3OH lines towards some of these sources (van der Tak et al., 2000a) indicate an increase in the CH3OH abundance by a factor of 100 in the warm inner envelope, which is likely due to evaporation of CH3OH-rich ices. The excitation of CH3OH, as well as that of CO and C2H2, correlates well with other temperature tracers, such as the gas/solid ratios of CO2 and H2O, the 45/100 $`\mu `$m colour and the fraction of heated solid 13CO2 (Boogert et al., 2000; van der Tak et al., 2000b; Boonman et al., 2003). Since these quantities also correlate with the ratio of envelope mass to stellar mass, their variation likely reflects evolution through the dispersal of the protostellar envelopes. These findings make the sample a good testbed for theories of chemical evolution during the earliest deeply embedded phases of high-mass star formation.

## 2 Observations

The data presented in this paper are part of a targeted spectral line survey of embedded massive stars, aimed to track chemical evolution during the earliest stages of high-mass star formation. Fifteen frequency settings have been observed, which also contain many lines of sulphurless molecules. Some of these data have been published by van der Tak et al. (1999, 2000a, 2000b). These observations were performed with the 15-m James Clerk Maxwell Telescope (JCMT)111The James Clerk Maxwell Telescope is operated by the Joint Astronomy Centre, on behalf of the Particle Physics and Astronomy Research Council of the United Kingdom, the Netherlands Organization for Scientific Research and the National Research Council of Canada. on Mauna Kea, Hawaii between 1995 and 1998. The beam size (FWHM) and main beam efficiency of the antenna were $`18^{\prime \prime }`$ and $`6469`$% at $`230`$ GHz and $`14^{\prime \prime }`$ and $`5863`$% at $`345`$ GHz. The frontends were the receivers A2, B3 and B3i; the backend was the Digital Autocorrelation Spectrometer, covering $`500`$ MHz instantaneous bandwidth. Pointing was checked every 2 hours during the observing and was always found to be within $`5^{\prime \prime }`$. To subtract the atmospheric and instrumental background, the chopping secondary was used with $`180^{\prime \prime }`$ offsets. Total integration times are 30–40 minutes for each frequency setting. Figure 1 shows examples of the data towards all but one of our sources; the data on W3 IRS5 were presented by Helmich & van Dishoeck (1997).

Additional observations of H2S and OCS in the 130-180 GHz window were carried out with the 30-m telescope of the Institut de Radio Astronomie Millimétrique (IRAM)222IRAM is an international institute for research in millimeter astronomy, cofunded by the Centre National de la Recherche Scientifique (France), the Max Planck Gesellschaft (Germany), and the Instituto Geografico Nacional (Spain). on Pico Veleta, Spain, in August 2002. The frontend was the facility receiver C150 and the backend the Versatile Spectral Assembly (VESPA) autocorrelator. Integration times were 10-20 minutes using frequency switching with a throw of 7.3 MHz. Data were calibrated onto $`T_{\mathrm{MB}}`$ scale by multiplying by 1.43, the ratio of forward and main beam efficiencies. The beam size is $`15^{\prime \prime }`$ at these frequencies.

The H2S 393 GHz line was observed in May 1995 at the 10.4-m Caltech Submillimeter Observatory (CSO)333The CSO is supported by NSF grant AST 99-80846., with the facility receiver as frontend and the 50 MHz AOS as backend. The telescope has a beam size of $`25^{\prime \prime }`$ and a main beam efficiency of 60% at this frequency.

Reduction was carried out using the CLASS package developed at IRAM. Linear baselines were subtracted and the spectra were smoothed once and calibrated onto $`T_{\mathrm{MB}}`$ scale. Only the frequency switched 30m data required higher order polynomials to fit the baseline. The final spectra have a resolution of $`0.31.5`$ km s-1 and rms noise levels (in $`T_{\mathrm{MB}}`$) of $`2030`$ mK for the 150 and 230 GHz bands, and $`3050`$ mK for the 345 GHz band. Multiplying these values by the line width (Table 4) yields the uncertainties of the line fluxes in Table 3. Although the absolute calibration is only correct to $``$30%, the relative strength of lines within one frequency setting is much more accurate.

Line identification was performed using the JPL catalog (Pickett et al., 1998)444http://spec.jpl.nasa.gov. A matching frequency is the main criterion for assignment, but the upper energy level of the transition and the complexity of the molecule are also considered. Table 2 has details about the observed transitions.

For each detected feature, a satisfying match could be found in the catalog within 0.5 MHz, except for the H2CS $`10_{19}9_{18}`$ line at 348.5 GHz, whose residual is 2 MHz. Laboratory spectroscopy of H2CS only exists up to 250 GHz (Beers et al., 1972), and the analysis in the JPL catalog only includes pure rotation and centrifugal distortion terms. At higher frequencies and for higher $`J`$-values, higher order terms become important, which may be why the catalog prediction is inaccurate for this line. Our data indicate a frequency of 348534.2 MHz.

## 3 Results

### 3.1 Line profiles

We have fitted the detected lines with Gaussian profiles. Table 3 reports the observed fluxes and Table 4 the widths of the lines. While most lines could be satisfactorily fitted with a single Gaussian, the lines of SO, which are the strongest lines that we have observed, required two Gaussians for a good fit to their profiles (Fig. 2; Table 4). One component agrees in position and width with the values measured for C17O and C34S by van der Tak et al. (2000b) and is attributed to a static envelope. Table 5 reports similar measurements for the source MonR2 IRS3 which van der Tak et al. did not study. The other component has $`2.8`$ times as large a width on average. This high-velocity emission may trace outflowing or infalling motions. The contributions from high- and low-velocity gas to the SO line fluxes are about equal, unlike CS, where the envelope carries 90% of the fluxes (van der Tak et al. 2000b; Table 5). The infrared absorption and submillimeter emission lines of CO towards these sources also show multiple velocity components (Mitchell et al., 1991; Mitchell & Hasegawa, 1991; Mitchell et al., 1992). However, the velocities do not match in detail, as expected because of differences in spectral and angular resolution between infrared and submillimeter data on one hand, and excitation differences between CO and SO on the other.

The fit of two Gaussians to the SO lines is excellent, which is not always the case for CO. An alternative method is to use area integrals, which do not attribute any low-velocity gas to the high-velocity component. Our approach may overestimate the high-velocity contribution to the line flux by factors up to 2.

The ratios of the high- to low-velocity contributions to the SO lines are $`\genfrac{}{}{0pt}{}{_<}{^{}}1`$ for the lines in the 230 GHz band, and $`\genfrac{}{}{0pt}{}{_>}{^{}}1`$ for the lines in the 345 GHz band. This difference could be due to the higher temperatures and densities required for excitation of the 345 GHz lines, or to the smaller beams in which they are observed. Since the high-velocity wings on the CS lines in these sources also become more pronounced in smaller beams, independent of transition (van der Tak et al., 2000b), the broad component probably is more compact than the envelope, but not necessarily warmer or denser. The high-velocity gas in these sources thus seems to be confined to small radii. This situation differs from that in molecular clouds, where line width increases with size, and from that in many molecular outflows, where velocities usually increase with radius (e.g., Lee et al. 2002). However, while CO data show arcminute-sized outflows in our sources, higher velocities naturally occur at smaller radii in the case of infall. For example, gas at 1,000 AU (1′′ at 1 kpc) in free fall onto a 10 M object would have $`V_{\mathrm{inf}}=4`$ km s-1, consistent with the observed width of the broad component (Table 4).

The broad component peaks at slightly more blueshifted velocities (by $`0.17`$ km s-1 on average) than the envelope, a situation found before for CS and H2CO (van der Tak et al., 2000b). In the case of outflow, obscuration by $`A_V=10^4`$ mag of dust could explain the absence of redshifted wings. Since some redshifted SO is seen (Fig. 2), the requirement relaxes to several 1000 mag. Of the sources discussed here, only NGC 6334 IRS1 may have such a large amount of dust in its envelope (Table 1); the other sources may need circumstellar disks seen face-on. In the case of infall, a mixture of absorption and emission is expected at redshifted velocities (Choi, 2002), which may partly cancel each other if unresolved, to give the impression of predominantly blueshifted emission. We conclude that both infall and outflow are plausible explanations for the high-velocity gas, and recommend interferometric observations of SO to decide between these options.

For most sources in our sample, the width of the SO2 lines is between the values for the low- and high-velocity components. The only source where the SO2 lines are strong enough to fit two components is W3 IRS5, where the line profiles do not show wings. However, most of these data, taken from Helmich & van Dishoeck (1997), were observed with an AOS as backend, and have lower spectral resolution than our data. It seems therefore plausible that, like for SO, about half of the SO2 emission arises in high-velocity gas.

For the other molecules, the signal to noise ratio is not high enough to disentangle multiple velocity components. This paper assumes that these lines originate in low-velocity gas, which given our results for SO may overestimate column densities and abundances (§§ 3.3, 4) by factors of up to 2.

### 3.2 Excitation temperatures

For SO and SO2, the number of detected lines is large enough to construct rotation diagrams. This method, described in detail by Blake et al. (1987) and Helmich et al. (1994), assumes that the lines are optically thin and that the molecular excitation can be described by a single temperature, the rotation temperature. If radiative decay competes with collisional excitation, the rotation temperature lies below the kinetic temperature. If the lines are optically thick, the rotation diagram underestimates the column density and overestimates the excitation temperature. Another complication is that the lines are measured in beams of unequal size, so that somewhat different volumes of gas are probed. Despite these caveats, the rotation diagram provides a useful first estimate of excitation conditions and molecular column densities, and a stepping stone towards a more sophisticated analysis (§ 4).

Table 6 presents the results. In the case of SO, separate fits were made to the low- and high-velocity components. For NGC 7538 IRS9, the detected SO2 lines do not cover a large enough range of energy levels to constrain $`T_{\mathrm{rot}}`$. The derived temperatures are 30…70 K for SO and 50…190 K for SO2, reflecting the available range of energy levels. Lower limits to $`T_{\mathrm{rot}}`$ are poor fits to the data, probably caused by optically thick lines. Indeed, the SO/34SO and SO2/34SO2 line ratios of $`1.3`$$`\genfrac{}{}{0pt}{}{_>}{^{}}10`$ are below the isotopic abundance ratio of 32S/34S=22, which indicates optical depths of a few in the SO and SO2 lines.

The last column of Table 6 lists rotation temperatures of CH3OH, measured previously. Strictly, these numbers are only lower limits to the kinetic temperatures in a $`15^{\prime \prime }`$ beam, but they show good correlation with $`T_{\mathrm{ex}}`$(C2H2) which does directly trace the kinetic temperature (van der Tak et al., 2000a). The rotation temperatures of SO and SO2 are not clearly correlated with $`T_{\mathrm{rot}}`$ (CH3OH) within the uncertainties (Fig. 3). Therefore we cannot conclude at this point whether these molecules trace the same gas within the envelopes. In any case, differences in beam sizes and optical depth effects in the lines make it hard to draw firm conclusions before a more detailed radiative transfer analysis is presented in § 4.

For the other molecules, for which only few lines are available, the line ratios have been used to constrain $`T_{\mathrm{ex}}`$. Source-averaged line ratios indicate $`T_{\mathrm{ex}}`$$``$50 K for H2CS, NS and HCS+, $`T_{\mathrm{ex}}`$$``$25 K for H2S and $`T_{\mathrm{ex}}`$$``$100 K for OCS. These values are only rough estimates, especially if based on few sources and/or lines.

### 3.3 Column densities

Table 7 presents beam-averaged molecular column densities, estimated from the observed line fluxes using the excitation temperatures found above. The column densities for H2CS and H2S are the sums of the ortho and para species, assuming an ortho to para ratio of 3. After calibration, the adopted $`T_{\mathrm{ex}}`$ is the main source of uncertainty in the column densities, especially if only one or two lines have been observed. For example, decreasing $`T_{\mathrm{ex}}`$ from 50 to 25 K increases column density estimates by factors of $``$2, while increasing $`T_{\mathrm{ex}}`$ from 50 to 75 K decreases column density estimates by $`\genfrac{}{}{0pt}{}{_<}{^{}}`$20%. In addition, column densities may have been underestimated by factors of a few due to nonnegligible optical depth.

To search for trends of column densities with temperature, Table 7 lists the sources in order of increasing $`T`$(CH3OH). The only trend seems to be that the warmest source has the largest column density of all molecules, except SO and SO2. In order to investigate whether these trends in the column density reflect chemical differences between sources or are due to optical depth or beam size effects, a more detailed radiative transfer analysis is presented in § 4.

### 3.4 Comparison with infrared data

Keane et al. (2001) observed the $`\nu _3`$ band of SO2 around 7.3 $`\mu `$m in absorption towards five of the sources studied here, and found column densities of SO2 of a few $`10^{16}`$ cm-2, a factor of $``$100 higher than the submillimeter data indicate. In contrast, infrared observations of CO and dust yield 3–5 times lower column densities than submillimeter data, as a result of non-spherical geometry on scales $`\genfrac{}{}{0pt}{}{_<}{^{}}1^{\prime \prime }`$ (van der Tak et al., 2000b). Correcting for optical depth in the SO2 submillimeter lines increases the column density estimate, but only by factors of a few (§ 3.2), which does not explain the discrepancy with the infrared values. More likely to be important are chemical effects on small scales. The sources are centrally concentrated, $`n`$(H2)$`r^\alpha `$ with $`\alpha =1.01.5`$, so that absorption data are more sensitive to warm, dense gas at small radii, while emission data probe more extended, cooler and less dense gas. The infrared excitation temperatures of 225–750 K also suggest that the SO2 absorption arises in warm gas.

The infrared estimates of $`N`$(SO2) imply optically thick submillimeter lines of SO2, but their brightness measured with the JCMT is much lower than the infrared $`T_{\mathrm{ex}}`$. The implied limit on the source size is $`\theta _S=15^{\prime \prime }(T_{\mathrm{ex}}/T_0)^{0.5}`$, with $`T_0`$ the typical $`T_{\mathrm{MB}}`$ of 0.1 K. Inferred values of $`\theta _S`$ range from 0$`\stackrel{}{.}`$17 for AFGL 2591 to 0$`\stackrel{}{.}`$32 for MonR2 IRS3, corresponding to linear radii of $`100250`$ AU. The limited spectral resolution of the infrared data prohibits assignment of the absorbers to either velocity component in emission, which however does not affect our conclusion.

## 4 Radiative transfer analysis

### 4.1 Model description

To estimate molecular abundances we have used the Monte Carlo radiative transfer program by Hogerheijde & van der Tak (2000) 555http://talisker.as.arizona.edu/$``$michiel/ratran.html. Starting from radial profiles of the density and kinetic temperature, this program solves for the molecular excitation as a function of radius. Besides collisional excitation, radiation from the cosmic microwave background and thermal radiation from local dust are taken into account. The result is integrated over the line of sight and convolved with the appropriate telescope beam. Observed and synthetic line fluxes are compared with a $`\chi ^2`$ statistic to find the best-fit abundance. Initially, abundances were assumed to be constant with radius.

Molecular data were mostly taken from the database by Schöier et al. (in prep.)666http://www.strw.leidenuniv.nl/$``$moldata, which paper also describes our extension of the SO2 rate coefficients by Green (1995) to higher energy levels. No rate coefficients are available for NS, so its excitation was assumed to be thermalized at the temperature of each grid point of the model. The population of the $`{}_{}{}^{2}\mathrm{\Pi }_{3/2}^{}`$ state of NS, which lies 225 K above the $`{}_{}{}^{2}\mathrm{\Pi }_{1/2}^{}`$ state, was assumed to be negligible. In the case of H2S, Ball et al. (1999) measured the $`1_{10}1_{01}`$ rate coefficient in collisions with He at $`T=1.3635.3`$ K. Our model uses the rate coefficients by Green et al. (1993) for H2O–He, multiplied by 2.6 based on comparison with the Ball et al. data, and scaled for the different reduced mass of the H2S–H2 system. Turner (1996) lists reasons why H2S rate coefficients are likely to be larger than H2O values.

### 4.2 Physical structure

The temperature and density structures of all our sources but MonR2 IRS3 were derived by van der Tak et al. (2000b). We have followed the same procedure to infer the structure of MonR2 IRS3, assuming a power-law density structure $`n=n_0(r/r_0)^\alpha `$. The temperature profile is derived from the observed luminosity of $`1.3\times 10^4`$ L (Henning et al., 1992), while $`n_0`$ is constrained by the submillimeter photometry of Giannakopoulou et al. (1997). These dust models solve the grain heating and cooling self-consistenly, using opacity model 5 from Ossenkopf & Henning (1994). These “OH5” opacities are the only ones that yield dust masses consistent with C17O measurements of AFGL 2591, where CO depletion is known to be negligible (van der Tak et al., 1999). To constrain $`\alpha `$, the CS and C34S line spectrum of MonR2 IRS3 was modeled with the Monte Carlo program. Besides our own JCMT measurements (Table 5), data from Tafalla et al. (1997) and Choi et al. (2000) were used. The best fit was obtained for $`\alpha =1.25`$ (Table 8), with a reduced $`\chi ^2=2.33`$ over 11 degrees of freedom, although $`\alpha =1.0`$ ($`\chi ^2=3.1`$) and $`\alpha =1.5`$ ($`\chi ^2=3.8`$) also give acceptable fits. Radial profiles of 350 $`\mu `$m dust emission for our sources indicate larger values of $`\alpha `$, by up to $`0.5`$, than CS line emission (Mueller et al., 2002). If this difference holds more generally, it may reflect a reduced CS abundance at small radii, although the models of § 5.2 do not predict that.

### 4.3 Molecular abundances

Table 9 lists the results of the radiative transfer analysis, with the same source ordering as in Table 7. Abundances, if assumed constant, are $``$10-8 for OCS, $``$10-9 for H2S, H2CS, SO and SO2, and $``$10-10 for HCS+ and NS. These trends are the same as found above for beam-averaged column densities (§ 3.3); the Monte Carlo analysis does not change results by factors of more than a few. Except for SO, SO2 and OCS, the source-to-source spread in the abundances is less than a factor of 10, and not related to temperature. These molecules (CS, H2S, H2CS, NS and HCS+) appear to trace the chemically inactive outer envelope ($`T<100`$ K), at least in the observed transitions. However, the high OCS abundances occur all in warm sources.

For SO and SO2, the number of detected lines is large enough to investigate their abundances as a function of radius. Figure 4 shows that the Monte Carlo models fit the SO lines uniformly well over the observed range of energy levels for all sources except W 33 A: there is no indication for changes in the SO abundance between $`T=30`$ and $`90`$ K. In the case of W 33 A, some of the observed values appear to be higher than the models, which may be due to an uncertain high-velocity contribution for these lines. Observations of higher-excitation SO lines ($`N=109`$ near 430 GHz or $`N=1110`$ near 473 GHz) would be valuable to investigate the effect of evaporating ice mantles on the abundance of SO.

Figure 5 shows the results of the models for SO2. Models with SO2/H2 $`10^9`$ fit the lines with $`E_u\genfrac{}{}{0pt}{}{_<}{^{}}100`$ K within factors of a few. However, these models underproduce the lines with $`E_u\genfrac{}{}{0pt}{}{_>}{^{}}100`$ K by factors of $`\genfrac{}{}{0pt}{}{_>}{^{}}`$100. This discrepancy is best seen for W3 IRS5, W33A and AFGL 2591, where the most SO2 lines are detected. For these sources, we have run models where the SO2 abundance increases sharply (‘jumps’) when $`T>100`$ K. These models are shown as triangles in Fig. 5 and quantified in Table 10. The SO2 abundance increases by a factor of 100 in the case of W3 IRS5 and W33A, and by a factor of 1000 for AFGL 2591. These models reproduce all observed SO2 lines to factors $``$10. The residuals are considerably bigger than the expected error bars of factors of 2…3, mainly due to the uncertain high-velocity contribution, which suggests that more detailed modeling is needed. However, these ‘jump’ models fit the data much better than models where the SO2 abundance is constant. The SO2 abundances derived from the high-excitation lines are similar to those based on the infrared data, indicating that the same gas is probed. The data for the other sources are consistent with similar abundance jumps, but not enough high-excitation SO2 lines have been detected to constrain such models. However, the data of AFGL 2136, NGC 7538 IRS1 and MonR2 IRS3 allow estimates of SO2 abundances in their inner envelopes, which Table 10 also lists. Except NGC 6334 IRS1, the jumps in the SO2 abundance occur all in sources with high $`T`$(CH3OH), confirming that the abundance enhancement is related to high temperatures, probably through ice evaporation.

The models discussed so far assumed a static envelope and aimed to describe the low-velocity gas. To investigate whether the high-velocity gas could be due to infalling motions, we have modeled the SO emission from W3 IRS5 under the assumption of gravitational infall. In this model, the infall velocity $`V_{\mathrm{inf}}`$ varies with radius $`R`$ as $`\sqrt{GM/R}`$, where $`G`$ is the gravitational constant. The luminosity of W3 IRS5, if due to a single star, suggests a stellar mass of $``$30 M (Vacca et al., 1996), so that $`V_{\mathrm{inf}}`$ ranges from 0.45 to 10.66 km s-1 in the model. However, the synthetic SO line profiles in 14–18′′ beams have widths of 5–6 km s-1 FWHM, much less than the observed value of 13.5 km s-1. Since the other sources have lower luminosities (hence infall speeds), it seems unlikely that the high-velocity SO and SO2 emission is solely due to infall.

## 5 Discussion

### 5.1 Comparison to other objects

Tables 9 and 10 also list molecular abundances in other regions. Only the case of the low-mass protostar IRAS 16293 has been analyzed to the same level of detail as done here. To reduce uncertainties introduced by methodological differences, the table contains abundances scaled to the value of the commonly observed CS molecule. Using ratios also facilitates comparison with chemical models, because these are independent of the total amount of sulphur in the gas phase, which is not well-understood. Ultraviolet observations of diffuse clouds toward $`\zeta `$ Oph indicate that sulphur is undepleted from its solar abundance of $``$2$`\times 10^5`$ (Savage & Sembach, 1996). However, for dense clouds, it is unclear what fraction of sulphur is locked up in dust grains.

The ratios observed in our sources differ considerably from those in PDRs, dark clouds and shocks, indicating that in those regions, other chemical processes dominate. The ratios in IRAS 16293 are at the low end of the values found here, suggesting that in low-mass objects, depletion of molecules onto grains is more important.

In the case of hot cores, both abundance ratios and absolute abundances are similar to those observed here. The chemical similarity of our sources with hot cores prompts us to compare our results with chemical models invoking ice evaporation. Our assumption is that source-to-source differences in abundances are due to age differences, rather than to differences in initial physical or chemical conditions.

### 5.2 Comparison with chemical models

The sulphur chemistry of hot cores, where ice mantles are evaporating off dust grains and start a high-temperature chemistry, was studied by Charnley (1997). The model assumes that the dominant form of sulphur in the ices is H2S, which after evaporation is destroyed by H and H3O+ to form S and H3S+. Reactions with OH and O2 make first SO and then SO2. Some H2S reacts with C+ into CS, which reacts with OH to make OCS. The sulphur chemistry strongly depends on temperature through its connection with the oxygen chemistry. At $`T\genfrac{}{}{0pt}{}{_>}{^{}}300`$ K, the SO2 abundance remains factors of 10–100 below that of SO because O and OH are driven into H2O. Results also depend on the presence of O2 in the ice mantles. Observational limits on solid and gaseous O2 (Vandenbussche et al., 1999; Goldsmith et al., 2000; Pagani et al., 2003) do not firmly rule out either assumption, but do make the model without O2 ice more plausible. However, this model cannot directly be compared with our observations, which probe a mixture of warm and cold gas.

The chemistry of protostellar envelopes was modeled by Doty et al. (2002). In the inner warm region, a hot core chemistry similar to that of Charnley (1997) is used, whereas in the cold outer parts, a low-temperature chemistry including freeze-out is adopted. The initial conditions of the model depend on temperature to mimic the effect of ice evaporation. Although Doty et al. adopt the specific temperature and density structure of AFGL 2591 as modeled by van der Tak et al. (2000b), their results do not change by factors of more than a few for our other sources, which have temperature and density distributions similar to those of AFGL 2591.

The Doty et al. (2002) model has most gas-phase sulphur initially in S at $`T<100`$ K (S/H2=6$`\times `$10-9) and in H2S at $`T>100`$ K (H2S/H2=1.6$`\times `$10-6) to mimic the effect of ice evaporation. This choice of initial conditions is consistent with the nondetection of the \[S I\] line at $`25.249`$ $`\mu `$m toward our sources with ISO. The upper limit of $`\tau <`$0.01 in absorption implies $`N`$(S) $`<5\times 10^{15}`$ cm-2, so that for $`N`$(H2) $`10^{23}`$ cm-2, S/H2$`<5\times 10^8`$, indicating that not all gas-phase sulphur in the inner envelope can be in atomic form. The inferred abundances imply that SO2 is one of the dominant sulphur-bearing gas-phase molecules in the warm inner envelope. However, the sum of all sulphur-bearing molecules in the gas phase is still much lower than the value derived for diffuse clouds (Savage & Sembach 1996). Unless a major sulphur-bearing gas-phase species has been missed in this survey, $``$90% of the sulphur must be in a solid form which has not yet been detected.

Infrared observations of our sources with ISO-SWS do not support the assumption that H2S is the main sulphur reservoir of the grain mantles. Toward W33A, a source with large ice abundances, the 3.98 $`\mu `$m band of solid H2S is not detected to $`\tau <0.02`$, which using a band strength of $`5\times 10^{18}`$ cm mol-1 gives $`N_s`$(H2S) $`<3.6\times 10^{17}`$ cm-2, or $`N_s`$(H2S)/$`N_s`$(H2O) $`<0.03`$ (W. Schutte, priv. comm.). However, solid OCS was detected with an abundance of $`N_s`$(OCS)/$`N_s`$(H2O) $`=(0.41)\times 10^3`$ or OCS/H2 $`=(414)\times 10^8`$ (Palumbo et al., 1997). Evaporation of solid OCS may explain our measured gas-phase OCS abundances which are otherwise unaccounted for. W 33 A is the only source in our sample for which both gas-phase and solid OCS have been detected. The results imply a gas/solid ratio of $``$0.5, much higher than that found for e.g. H2O and CO2 in this source (Boonman et al., 2003; Boonman & van Dishoeck, 2003). Observations of both gas-phase and solid-state OCS toward other sources are needed to determine the role of grain-mantle evaporation. Observations of OCS lines from a wide range of energy levels may reveal an abundance increase at $`T\genfrac{}{}{0pt}{}{_>}{^{}}100`$ K, as indeed found for IRAS 16293 (Schöier et al., 2002). The low $`T_{\mathrm{ex}}`$ of H2S and the high value of OCS (§ 3.2) support these conclusions.

The effect of evaporating OCS-rich ice mantles on hot core chemistry is explored by Doty et al. (2003). Briefly, the gas-phase abundances of CS, HCS+ and OCS increase about linearly with the fraction of solid OCS, while those of H2S, SO and SO2 decrease about linearly. The H2CS molecule is not affected much. Since the gas-phase chemistry of CS, HCS+ and OCS depends only weakly on time, the use of sulphur molecules as chemical clock diminishes when H2S is only a minor sulphur carrier in the ice mantles. However, the ’jump’ models for SO2 (Table 10) indicate that either H2S or SO2 is present in the grain mantles. The limits on solid SO2 from infrared observations are not very stringent (Boogert et al., 1998).

We have compared our observations to the models of Doty et al. (2002), modified to have OCS in the evaporating ice mantles (S. Doty, priv. comm.). This model reproduces our observed abundances of CS, SO and H2CS in the outer envelopes (Table 9) to within factors of a few for chemical ages of $``$105 yr. The abundances of SO2 in the inner envelopes (Table 10) are also reproduced for the same chemical age. The model abundances of HCS+ and SO2 in the outer envelopes are factors of $``$10 below the observed values, which may be due to the adopted initial conditions, or the chemical network, which is based on the UMIST reaction rate database (Le Teuff et al., 2000)777http://www.rate99.co.uk. The outer envelope abundances of H2S and OCS are underproduced even more (factors of $``$100). However, the small number of observed lines forced us to assume constant abundances for these species, while they may in fact be confined to the inner envelopes.

### 5.3 The role of shock chemistry

As an alternative to ice evaporation models, we consider the effect of non-dissociative shocks on the molecular envelopes surrounding the protostars. Such shocks occur by the interaction of the molecular outflows of protostars with their envelopes. Shocks fast enough to dissociate H2 ($`v\genfrac{}{}{0pt}{}{_>}{^{}}40`$ km s-1) would have a major impact on their molecular composition which our data do not indicate. Although our sources show outflows faster than that, both in infrared absorption (Mitchell et al., 1991) and in submillimeter emission (Mitchell & Hasegawa, 1991; Mitchell et al., 1992) of CO, their filling factor must be small.

Our observed SO, SO2 and H2S abundances in the outer envelopes are consistent with those in postshock gas if sulphur was initially in molecular form (Leen & Graff, 1988). The observed similarity of these abundances is not reproduced by the Pineau des Forêts et al. (1993) model which has most sulphur initially in atomic form. However, the two models are hard to compare because one is a two-fluid treatment, while the other treats neutral and charged particles as one fluid. Before concluding that postshock gas reproduces our data, we recommend that a two-fluid calculation is performed with most sulphur initially in molecular form.

Two other molecules suggest that the envelopes of young high-mass stars may have been processed by shocks, although neither one conclusively. First, ISO-SWS data show low ratios of gas-phase to solid-state CO2 ($`\genfrac{}{}{0pt}{}{_<}{^{}}0.1`$; van Dishoeck 1998) and the CO2 gas phase abundance remains low through the 100–300 K temperature regime. After evaporating off grains, CO2 must be promptly destroyed, which shocks can do in reactions with H (Charnley & Kaufman, 2000) or perhaps H2 (Doty et al., 2002; Talbi & Herbst, 2002). Second, Hatchell & Viti (2002) measured NS/CS ratios of 0.02-0.05 in a sample of six hot cores and interpreted these as evidence for shocks. The main reactions to form NS require SH and NH which are produced from OH. Shocks use OH to form H2O and suppress the production of NS. The values of NS/CS=0.001 – 0.01 measured here may indicate that shocks play a role too. However, for $`t=3\times 10^4`$ yr, the Doty et al. (2002) model also predicts NS/CS $`10^3`$, so this ratio cannot be used to demonstrate the influence of shocks.

## 6 Conclusions

Submillimeter observations of SO, SO2, H2S, H2CS, OCS, NS and HCS+ toward nine embedded massive stars show that:

1. High-velocity gas contributes 10–50% to the CS, SO and SO2 emission in $`1520^{\prime \prime }`$ beams. This gas may either trace the outflows which are also seen in CO, or they may trace infall motions.
2. The SO2 abundance increases from dark cloud levels in the outer envelope ($`T<100`$ K) to levels seen in hot cores and shocks in the inner envelope ($`T>100`$ K).
3. Molecular abundances are consistent with a model of ice evaporation in an envelope with gradients in temperature and density for a chemical age of $``$105 yr.
4. The high observed abundance of OCS, the fact that $`T_{\mathrm{ex}}`$(OCS)$``$$`T_{\mathrm{ex}}`$(H2S), and the detection of solid OCS and limit on solid H2S all suggest that OCS is a major sulphur carrier in grain mantles, rather than H2S.
5. For most other sulphur-bearing molecules, the source-to-source abundance variations by factors of up to 10 do not correlate with previously established evolutionary trends in temperature tracers. These species probe the chemically inactive outer envelope. Our data set does not constrain the abundances of H2S and SO in the inner envelope, which, together with SO2, are required to use sulphur as a clock.

These conclusions could be followed up by, respectively:

1. interferometric observations at $`\genfrac{}{}{0pt}{}{_<}{^{}}1^{\prime \prime }`$ resolution of SO lines from a wide range of energy levels.
2. similar interferometric observations of SO2.
3. updating the model of Leen & Graff (1988) using a two-fluid treatment, and extending it to later times.
4. a multi-line study of OCS, covering energy levels well above and well below 100 K.
5. observations of high-excitation SO and H2S lines, e.g. SO 430 and 473 GHz and H2S 369 GHz.

Items 3 and 4 are within the capabilities of current computational and observational facilities. The CSO and JCMT may be able to carry out item 5 in very good weather, but the future Atacama Pathfinder Experiment (APEX) can benefit from an even better site. Current interferometers are not sensitive enough for items 1 and 2, which require observing several lines and several sources. However, future instruments such as the Combined Array for Research into Millimeter Astronomy (CARMA) and the Atacama Large Millimeter Array (ALMA) will be well-suited for these projects.

###### Acknowledgements.

The authors thank the JCMT and IRAM staffs for their support, Malcolm Walmsley, Holger Müller, Jennifer Hatchell and Xander Tielens for useful discussions, and Willem Schutte, Fred Lahuis and Steve Doty for input. Astrochemistry in Leiden is supported by the Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO) through grant 614-41-003 and a Spinoza award.


format: text

text:
# Non-linear 𝜎-models in noncommutative geometry: fields with values in finite spaces

## 1 Introduction

In  we have constructed some noncommutative analogues of two dimensional non-linear $`\sigma `$-models. Since these models exhibit rich and easily accessible geometrical structures, their noncommutative counterparts are very useful to probe the interplay between noncommutative geometry and field theory. We proposed three classes of models: an analogue of the Ising model which admits instanton solutions, the analogue of the principal chiral model together with its infinite number of conserved currents and the noncommutative Wess-Zumino-Witten model together with its modified conformal invariance.

In this short report we extend the first model by constructing instanton solutions for any value of the topological charge (restricted to be $`1`$ in ) and for an arbitrary complex structure parametrized by a complex number $`\tau `$, $`\mathrm{}\tau >0`$.

## 2 The general construction

Ordinary non-linear $`\sigma `$-models are field theories whose configuration space consists of maps $`X`$ from the source space, a Riemannian manifold $`(\mathrm{\Sigma },g)`$, which we assume to be compact and orientable, to a target space, an other Riemannian manifold $`(,G)`$. The corresponding action functional is given, in local coordinates, by

$$S[X]=\frac{1}{2\pi }_\mathrm{\Sigma }\sqrt{g}g^{\mu \nu }G_{ij}(X)_\mu X^i_\nu X^j,$$
(1)
where as usual $`g=detg_{\mu \nu }`$ and $`g^{\mu \nu }`$ is the inverse of $`g_{\mu \nu }`$. The stationary points of this functional are harmonic maps from $`\mathrm{\Sigma }`$ to $``$ and describe minimal surfaces embedded in $``$. Different choices of the source and target spaces lead to different field theories, some of them playing a major role in physics.

When $`\mathrm{\Sigma }`$ is two dimensional, the action $`S`$ is conformally invariant, that is it is left invariant by any rescaling of the metric $`gge^\sigma `$, where $`\sigma `$ is any map from $`\mathrm{\Sigma }`$ to $``$ . As a consequence, the action only depends on the conformal class of the metric and may be rewritten using a complex structure on $`\mathrm{\Sigma }`$ as

$$S[X]=\frac{i}{\pi }_\mathrm{\Sigma }G_{ij}(X)X^i\overline{}X^j,$$
(2)
where $`=_zdz`$ and $`\overline{}=_{\overline{z}}d\overline{z}`$ and $`z`$ is a suitable local complex coordinate.

In order to A noncommutative generalization is constructed by a dualization and reformulation it in terms of the $``$-algebras $`𝒜`$ and $``$ of complex valued smooth functions defined respectively on $`\mathrm{\Sigma }`$ and $``$. Then, embeddings $`X`$ of $`\mathrm{\Sigma }`$ into $``$ correspond to $``$-algebra morphisms $`\pi _X`$ from $``$ to $`𝒜`$, the correspondence being simply given by the pullback, $`f\pi _X(f)=fX`$.
Now, all this makes perfectly sense also for noncommutative algebras $`𝒜`$ and $``$, and we take as configuration space the space of all $``$-algebra morphisms from $``$ to $`𝒜`$. Both algebras are over $``$ and for simplicity we take them to be unital. The definition of the action functional involves noncommutative generalizations of the conformal and Riemannian geometries. According to Connes , the former can be understood within the framework of positive Hochschild cohomology. Indeed, in the commutative two dimensional situation, the tri-linear map $`\varphi :𝒜^3`$ defined by

$$\varphi (f_0,f_1,f_2)=\frac{i}{\pi }_\mathrm{\Sigma }f_0f_1\overline{}f_2$$
(3)
is an extremal element of the space of positive Hochschild cocycles that belongs to the Hochschild cohomology class of the cyclic cocycle $`\psi `$ defined by

$$\psi (f_0,f_1,f_2)=\frac{i}{2\pi }_\mathrm{\Sigma }f_0𝑑f_1df_2.$$
(4)
Now, expressions (3) and (4) still make perfectly sense for a general noncommutative algebra $`𝒜`$. One can say that $`\psi `$ allows to integrate 2-forms in dimension 2, so that it is a metric independent object, whereas $`\varphi `$ defines a suitable scalar product

$$a_0da_1,b_0db_1=\varphi (b_0^{}a_0,a_1,b_1^{})$$
(5)
on the space of 1-forms and thus depends on the conformal class of the metric. Moreover, this scalar product is positive and invariant with respect to the action of the unitary elements of $`𝒜`$ on 1-forms, and its relation to the cyclic cocycle $`\psi `$ allows to prove various inequalities involving topological quantities .

We can compose such a cocycle $`\varphi `$ with a morphism $`\pi :𝒜`$ to obtain a positive cocycle on $``$ defined by $`\varphi _\pi =\varphi (\pi \pi \pi )`$. In order to build an action functional, which assigns a number to any morphism $`\pi `$, we have to evaluate the cocycle $`\varphi _\pi `$ on a suitably chosen element of $`^3`$. Such an element provides the noncommutative analogue of the metric on the target, and we take it as a positive element $`G=_ib_0^i\delta b_1^i\delta b_2^i`$ of the space of universal 2-forms $`\mathrm{\Omega }^2()`$. Thus, the quantity

$$S[\pi ]=\varphi _\pi (G)$$
(6)
is well defined and positive. We shall consider it to be a noncommutative analogue of the action functional of the non linear $`\sigma `$-model.

Clearly, we consider $`\pi `$ as the dynamical variable (the embedding) whereas $`\varphi `$ (the conformal structure on the source) and $`G`$ (the metric on the target) are background structures that have been fixed. Alternatively, one could take only the metric $`G`$ on the target as a background field and use the morphism $`\pi :𝒜`$ to define the induced metric $`\pi _{}G`$ on the source as $`\pi _{}G=_i\pi (b_0^i)\delta \pi (b_1^i)\delta \pi (b_2^i)`$, which is obviously a positive universal 2-form on $`𝒜`$. To such an object one can associate, by means of a variational problem , a positive Hochschild cocycle that stands for the conformal class of the induced metric.

The critical points of the $`\sigma `$-model corresponding to the action functional (6) are noncommutative generalizations of harmonic maps and describe “minimally embedded surfaces” in the noncommutative space associated with $``$.

## 3 Two points as a target space

The simplest example of a target space is that of a finite space made of two points $`=\{1,2\}`$, like in the Ising model. Now, any continuous map from a connected surface to a discrete space is constant and the resulting (commutative) theory would be trivial. However, this is not the case if the source is a noncommutative space and one has, in general, lots of algebra morphisms. The algebra of functions over $`=\{1,2\}`$ is just $`=^2`$ and any element $`f`$ is a couple of complex numbers $`(f_1,f_2)`$ with $`f_i=f(i)`$, the value of $`f`$ at the point $`i`$. As a vector space $``$ is generated by the function $`e`$ defined by $`e(1)=1,e(2)=0`$. Clearly, $`e`$ is a hermitian projection, $`e^2=e^{}=e`$, and $``$ can be thought of as the unital -algebra generated by such an element $`e`$. As a consequence, any $``$-algebra morphism $`\pi `$ from $``$ to $`𝒜`$ is given by a hermitian projection $`p=\pi (e)`$ in $`𝒜`$. The configuration space is then the collection of all such projections and our noncommutative $`\sigma `$-model provides a dynamics of projections. Choosing the metric $`G=\delta e\delta e`$ on the space $``$ of two points, the action functional (6) simply becomes

$$S[p]=\varphi (1,p,p),$$
(7)
where $`\varphi `$ is a given Hochschild cocycle corresponding to the conformal structure.
As we have already mentioned, from general consideration of positivity in Hochschild cohomology this action is bounded by a topological term . In the following, we shall explicitly prove this fact when taking the noncommutative torus as source space.

## 4 The noncommutative torus as a source space

We recall the very basic aspects of the noncommutative torus that we shall need in the following. The algebra $`𝒜_\theta `$ of smooth functions on the noncommutative torus is the unital $``$-algebra made of power series,

$$a=\underset{m,n^2}{}a_{mn}(U_1)^m(U_2)^n,$$
(8)
with $`a_{mn}`$ a complex-valued Schwarz function on $`^2`$ that is, the sequence of complex numbers $`\{a_{mn},(m,n)^2\}`$ decreases rapidly at ‘infinity’. The two unitary elements $`U_1,U_2`$ satisfy the commutation relations

$$U_2U_1=e^{2\pi i\theta }U_1U_2.$$
(9)
There exist on $`𝒜_\theta `$ a unique normalized positive definite trace, denoted by the integral symbol $`\text{ }\text{ }:𝒜_\theta `$, which is given by

$$\text{ }(\underset{(m,n)^2}{}a_{mn}(U_1)^m(U_2)^n):=a_{00}.$$
(10)
This trace is invariant under the action of the commutative torus $`𝕋^2`$ on $`𝒜_\theta `$ whose infinitesimal form is given by two commuting derivations $`_1,_2`$ acting by

$$_\mu (U_\nu )=2\pi i\delta _\mu ^\nu U_\nu ,\mu ,\nu =1,2$$
(11)
and the invariance means just that $`\text{ }\text{ }_\mu (a)=0,\mu =1,2`$ for any $`a𝒜_\theta `$.
The cyclic 2-cocycle allowing the integration of 2-forms is simply given by

$$\psi (a_0,a_1,a_2)=\frac{1}{2\pi i}\text{ }ϵ_{\mu \nu }a_0_\mu a_1_\nu a_2,$$
(12)
where $`ϵ_{\mu \nu }`$ is the standard antisymmetric tensor. Its normalization ensures that for any hermitian projection $`p𝒜_\theta `$, the quantity $`\psi (p,p,p)`$ is an integer: it is indeed the index of a Fredholm operator .

The conformal class of a general constant metric is parametrized by a complex number $`\tau `$, $`\mathrm{}\tau >0`$. Then, up to a conformal factor, the metric is given by

$$g=(g_{\mu \nu })=\left(\begin{array}{cc}1& \mathrm{}\tau \\ \mathrm{}\tau & |\tau |^2\end{array}\right),$$
(13)
with inverse given by $`g^1=(g^{\mu \nu })=\frac{1}{(\mathrm{}\tau )^2}\left(\begin{array}{cc}\left|\tau \right|^2& \mathrm{}\tau \\ \mathrm{}\tau & 1\end{array}\right)`$ and $`\sqrt{detg}=\mathrm{}\tau `$.
By using the two derivations $`_1,_2`$ defined in (11) we still think of ‘the complex torus’ $`𝕋^2`$ as acting on the noncommutative torus $`𝒜_\theta `$ and construct two associated derivations of $`𝒜_\theta `$ given by $`_{(\tau )}=\frac{1}{(\tau \overline{\tau })}(\overline{\tau }_1+_2)`$ and $`\overline{}_{(\tau )}=\frac{1}{(\tau \overline{\tau })}(\tau _1_2)`$. One easily finds that $`_{(\tau )}\overline{}_{(\tau )}=\overline{}_{(\tau )}_{(\tau )}=\frac{1}{4}g^{\mu \nu }_\mu _\nu =\frac{1}{4}\mathrm{\Delta }`$, and the operator $`\mathrm{\Delta }=g^{\mu \nu }_\mu _\nu `$ is just the Laplacian of the metric (13).

By working with the metric (13), the positive Hochschild cocycle $`\varphi `$ associated with the cyclic one (12) will be given by

$$\varphi (a_0,a_1,a_2)=\frac{2}{\pi }\text{ }a_0_{(\tau )}a_1\overline{}_{(\tau )}a_2.$$
(14)
A construction of the cocycle (14) as the conformal class of a general constant metric on the noncommutative torus can be found in .

### 4.1 The action and the field equations

With $`𝒫_\theta =Proj(𝒜_\theta )`$ denoting the collection of all projections in the algebra $`𝒜_\theta `$, we construct an action functional $`S:𝒫_\theta ^+`$ by

$$S_{(\tau )}(p)=\varphi (1,p,p)=\frac{2}{\pi }\text{ }\sqrt{detg}_{(\tau )}p\overline{}_{(\tau )}p.$$
(15)
The action functional can also be written as

$$S_{(\tau )}(p)=\frac{1}{2\pi }\text{ }\sqrt{detg}g^{\mu \nu }_\mu p_\nu p=\frac{1}{\pi }\text{ }\sqrt{detg}g^{\mu \nu }p_\mu p_\nu p.$$
(16)
Here the two derivations $`_\mu `$ are the ones defined in (11) while the metric $`g`$ is the one in (13) which carries also the dependence on the complex parameter $`\tau `$.

By taking into account the nonlinear nature of the space $`𝒫_\theta `$, one finds that the most general infinitesimal variation of its elements (the tangent vectors) is of the form $`\delta p=(1p)zp+pz^{}(1p)`$, with $`z`$ an arbitrary elements in $`𝒜_\theta `$. Then, simple algebraic manipulations give the equations of motion,

$$p\mathrm{\Delta }(p)\mathrm{\Delta }(p)p=0,$$
(17)
where $`\mathrm{\Delta }=g^{\mu \nu }_\mu _\nu `$ is the Laplacian. These are non linear second order ‘differential’ equations and it is not simple to give their solutions in a closed form. We shall show that the absolute minima of (16) in a given connected component of $`𝒫_\theta `$ actually fulfill first order equations which are easier to solve.

### 4.2 Topological charges and self-duality equations

The cyclic 2-cocycle (12) assigns to any projection $`p𝒫_\theta `$ a ‘topological charge’ (the first Chern number)

$$\psi (p):=\frac{1}{2\pi i}\text{ }p\left[_1(p)_2(p)_2(p)_1(p)\right].$$
(18)
Then, the following inequality holds

$$S_{(\tau )}(p)2|\psi (p)|.$$
(19)
Indeed, due to positivity of the trace $`\text{ }\text{ }`$ and its cyclic properties, we have that

$$0\text{ }\sqrt{detg}g^{\mu \nu }\left[_\mu (p)p\pm iϵ_\mu ^\alpha _\alpha (p)p\right]^{}\left[_\nu (p)p\pm iϵ_\nu ^\beta _\beta (p)p\right].$$
(20)
By expanding the LHS and comparing with (16) and (18) we get inequality (19). In (20) the symbol $`ϵ_{\mu ,\nu }`$ stands for the volume form of the metric $`g_{\mu ,\nu }`$. The inequality (19), which gives a lower bound for the action, is the analogue of the one for ordinary $`\sigma `$-models . A similar bound for a model on the fuzzy $`2`$-sphere was obtained in .

It is clear that the equality in (19) occurs when the projection $`p`$ satisfies the following self-dual or anti-self duality equations

$$\left[_\mu p\pm iϵ_\mu ^\alpha _\alpha p\right]p=0\mathrm{and}/\mathrm{or}p\left[_\mu piϵ_\mu ^\alpha _\alpha p\right]=0$$
(21)
The self duality equations (21) can be written as

$$\overline{}_{(\tau )}(p)p=0\mathrm{and}/\mathrm{or}p_{(\tau )}(p)=0,$$
(22)
while the anti-self duality one is

$$_{(\tau )}(p)p=0\mathrm{and}/\mathrm{or}p\overline{}_{(\tau )}(p)=0.$$
(23)
It is straightforward to check that either of the equations (22) and (23) implies the field equations (17) as it should.

## 5 The instantons

The connected components of $`𝒫_\theta `$ are parametrized by two integers $`r`$ and $`q`$ such that $`r+q\theta >0`$.  When $`\theta `$ is irrational, the corresponding projections have trace $`r+q\theta `$ and the topological charge $`\psi (p)`$ appearing in (18) is just $`q`$.

Thus our task is to find projections that belongs to the previous homotopy classes and satisfy (say) the self-duality equation $`(p)p=0`$ or, equivalently, $`pp=0`$. These equations are still non linear and to solve them our next step will be a reduction to a linear problem. The key point is to identify the algebra $`𝒜_\theta `$ as the endomorphism algebra of a suitable bundle and to think of any projection in it as an operator on such a bundle. The bundle in question will be a projective module of finite type on a different copy $`𝒜_\alpha `$ of the noncommutative torus, the two algebras $`𝒜_\theta `$ and $`𝒜_\alpha `$ being related by Morita equivalence.

### 5.1 The modules

Let us then consider another copy $`𝒜_\alpha `$ of the noncommutative torus with generators $`Z_1,Z_2`$ obeying the relation $`Z_2Z_1=e^{2\pi i\alpha }Z_1Z_2`$. When $`\alpha `$ is not rational, every finitely generated projective module over the algebra $`𝒜_\alpha `$ which is not free is isomorphic to a Heisenberg module . Any such a module $`_{r,q}`$ is characterized by two integers $`r,q`$ which can be taken to be relatively coprime with $`q>0`$, or $`r=0`$ and $`q=1`$. We shall briefly describe them. As a vector space

$$_{r,q}=𝒮(\times _q)𝒮()^q,$$
(24)
the space of Schwarz functions of one continuous variable $`s`$ and a discrete one $`k_q`$ (we shall implicitly understand that such a variable is defined modulo $`q`$). By denoting $`\epsilon =r/q\alpha `$, the space $`_{r,q}`$ is made into a right module over $`𝒜_\alpha `$ by

$$(\xi Z_1)(s,k):=\xi (s\epsilon ,kr),(\xi Z_2)(s,k):=e^{2\pi i(sk/q)}\xi (s,k),$$
(25)
with $`\xi _{r,q}`$; the relations $`Z_2Z_1=e^{2\pi i\alpha }Z_1Z_2`$ for the torus $`𝒜_\alpha `$ are easily verified.
On $`_{r,q}`$ one defines an $`𝒜_\alpha `$-valued hermitian structure $`,_\alpha :_{r,q}\times _{r,q}𝒜_\alpha `$, which is antilinear in the first factor.

It is proven in  that the endomorphism algebra $`End_{𝒜_\alpha }(_{r,q})`$, which acts on the left on $`_{r,q}`$, can be identified with another copy of the noncommutative torus $`𝒜_\theta `$ where the parameter $`\theta `$ is ‘uniquely’ determined by $`\alpha `$ in the following way. Since $`r`$ and $`q`$ are coprime, there exist integer numbers $`a,b`$ such that $`ar+bq=1`$. Then, the transformed parameter is given by $`\theta =(a\alpha +b)/(q\alpha +r)`$. Notice that, given any two other integers $`a^{},b^{}`$ such that $`a^{}r+b^{}q=1`$, one would find that $`\theta ^{}\theta `$ so that $`𝒜_\theta ^{}𝒜_\theta `$. Thus, we are saying that the algebra $`End_{𝒜_\alpha }(_{r,q})`$ is generated by two operators $`U_1,U_2`$ acting on the left on $`_{r,q}`$ by

$$(U_1\xi )(s,k):=\xi (s1/q,k1),(U_2\xi )(s,k):=e^{2\pi i(s/\epsilon ak)/q}\xi (s,k),$$
(26)
and one verifies that $`U_2U_1=e^{2\pi i\theta }U_1U_2`$, the defining relations of the algebra $`𝒜_\theta `$.

The crucial fact is that the $`𝒜_\theta `$-$`𝒜_\alpha `$-bimodule $`_{r,q}`$ is a Morita equivalence between the two algebras $`𝒜_\theta `$ and $`𝒜_\alpha `$: there exists also a $`𝒜_\theta `$-valued hermitian structure on $`_{r,q}`$, $`,_\theta :_{r,q}\times _{r,q}𝒜_\theta `$, which is compatible with the $`𝒜_\alpha `$-valued one $`,_\alpha `$,

$$\xi ,\eta _\theta \zeta =\xi \eta ,\zeta _\alpha ,$$
(27)
for all $`\xi ,\eta ,\zeta _{r,q}`$. The hermitian structure $`,_\theta `$ is antilinear in the second factor.

### 5.2 The constant curvature connection

A gauge connection on the right $`𝒜_\alpha `$-module $`_{r,q}`$ is given by two covariant derivative operators $`_\mu :_{r,q}_{r,q}`$, $`\mu =1,2`$ which satisfy a right Leibniz rule

$$_\mu (\xi a)=(_\mu \xi )a+\xi (_\mu a),\mu =1,2.$$
(28)
One also requires compatibility with the $`𝒜_\alpha `$-valued hermitian structure

$$_\mu (\xi ,\eta _\alpha )=_\mu \xi ,\eta _\alpha +\xi ,_\mu \eta _\alpha ,\mu =1,2.$$
(29)
A particular connection on the right $`𝒜_\alpha `$-module $`_{r,q}`$ is given by the operators

$$(_1\xi )(s,k):=\frac{2\pi i}{\epsilon }s\xi (s,k),(_2\xi )(s,k):=\frac{d\xi }{ds}(s,k)$$
(30)
(the discrete index $`k`$ is not touched); this connection is of constant curvature,

$$F_{1,2}:=[_1,_2]_{[_1,_2]}=\frac{2\pi i}{\epsilon }𝕀_{_{r,q}},$$
(31)
with $`𝕀_{_{r,q}}`$ the identity operator on $`_{r,q}`$.

Given any connection $`_\mu `$ on $`_{r,q}`$, one can define derivations on the endomorphism algebra $`End_{𝒜_\alpha }(_{r,q})`$ by commutators: $`\widehat{\delta }_\mu (T):=_\mu TT_\mu ,\mu =1,2`$, for any $`TEnd_{𝒜_\alpha }(_{r,q})`$. Then, by remembering that $`End_{𝒜_\alpha }(_{r,q})𝒜_\theta `$, one finds that the derivations $`\widehat{\delta }_\mu `$ on $`End_{𝒜_\alpha }(_{r,q})`$ determined by the particular connection (30) are proportional to the generators of the infinitesimal action of the torus $`𝕋^2`$ on $`𝒜_\theta `$, that is the canonical derivations (11): $`\widehat{\delta }_\mu (U_\nu )=\frac{2\pi i}{q\epsilon }\delta _\mu ^\nu U_\nu ,\mu ,\nu =1,2`$.

The holomorphic and anti-holomorphic connections $`_{(\tau )},\overline{}_{(\tau )}`$ will be the lift of the derivations $`_{(\tau )},\overline{}_{(\tau )}`$ with respect to the connection (30).

### 5.3 Instantons from Gaussians

We shall look for solutions of the self-dual equations (22) of the form

$$p_\psi :=|\psi \psi ,\psi ^1\psi |,$$
(32)
with $`|\psi `$ a ‘section of a suitable vector bundle’ over the noncommutative torus $`𝒜_\theta `$ and $`\psi ,\psi `$ an invertible element in another noncommutative torus which cannot be $`𝒜_\theta `$ itself but rather is Morita equivalent to it, being indeed $`𝒜_\alpha `$. We shall take $`|\psi `$ to be an element of the Schwarz space $`_{r,q}`$ on which $`𝒜_\theta `$ acts on the left as the endomorphism algebra. Thus, $`_{r,q}`$ will be though of as a right module over the algebra $`𝒜_\alpha `$ and $`\psi ,\psi =\psi ,\psi _\alpha `$ will be required to be an invertible element of $`𝒜_\alpha `$.

Than, let us suppose that $`|\psi _{r,q}`$ be such that $`\psi ,\psi _\alpha `$ is invertible. A simple computation shows that the projection $`p_\psi :=|\psi (\psi ,\psi _\alpha )^1\psi |`$ is a solution of the self-duality equation (22) if and only if the element $`|\psi _{r,q}`$ obeys the equation

$$\overline{}\psi \psi \lambda =0,$$
(33)
with $`\lambda `$ a suitable element in $`𝒜_\alpha `$ and $`\overline{}`$ the anti-holomorphic connection.

For $`\lambda `$ there are simple solutions to (33) given by generalized Gaussians

$$\psi _\lambda (s,k)=A_ke^{i\tau \pi s^2/\epsilon +\lambda (\overline{\tau }\tau )s}.$$
(34)
Any such a function is an element of $`_{r,q}`$ since $`\mathrm{}\tau >0`$. The vector $`A=(A_1,\mathrm{},A_q)^q`$ can be taken in the complex projective space $`^{q1}`$ by removing an inessential normalization. Restrictions on the possible values of the constant parameter $`\lambda `$ will be discussed presently.

A projection corresponding to a Gaussian was constructed in  for the lowest value of the charge $`q=1`$ and for $`\tau =i`$ and the invertibility of the corresponding element $`\psi _\lambda ,\psi _\lambda _\alpha `$ (now $`\alpha =1/\theta `$) was proved for any values of $`0<\theta <1`$ in . By using the methods of  one should be able to prove invertibility of $`\psi _\lambda ,\psi _\lambda _\alpha 𝒜_\alpha `$ for the most general situation.

Gauge transformations (in fact complexfied ones) are provided by invertible elements in $`𝒜_\alpha `$ acting on the right on $`_{r,q}`$,

$$_{r,q}|\psi |\psi ^g=|\psi g_{r,q},gGL(𝒜_\alpha ).$$
(35)
It is clear that projections of the form (32) are invariant under gauge transformations. Notice that we do not require $`g`$ to be unitary.

Now, let $`|\psi `$ be a solution of (33), $`\overline{}\psi \psi \lambda =0`$; and let $`gGL(𝒜_\alpha )`$. Then a simple computation shows that the gauge transformed vector $`|\psi ^g`$ will be a solution of an equation of the form (33): $`\overline{}\psi ^g\psi ^g\lambda _g=0`$ with $`\lambda _g`$ given by

$$\lambda _g=g^1\lambda g+g^1\overline{}_{(\tau )}g.$$
(36)
The moduli space of Gaussians is simply described. It turns out that two Gaussians $`\psi _\lambda `$ and $`\psi _\lambda ^{}`$ are related by a gauge transformation if an only if there exist integers $`(m,n)^2`$ such that

$$\lambda ^{}=\lambda _g,\mathrm{with}g=(Z_1)^m(Z_2)^n.$$
(37)
Furthermore,

$$\lambda _g\lambda =\frac{2\pi i\tau }{\tau \overline{\tau }}\left(m\frac{1}{\tau }n\right)$$
(38)
Thus, gauge nonequivalent constant parameters $`\lambda `$ form a complex ordinary torus $`𝕋_\tau ^2`$ and the moduli space of Gaussians (34) is $`^{q1}\times 𝕋_\tau ^2`$.

It is an open problem to prove whether it is possible to gauge to a Gaussian any solution of the self-duality equation (33). Equivalently, this problem can be stated as follows: given any $`\rho 𝒜_\alpha `$, there exists an element $`gGL(𝒜_\alpha )`$ such that $`\rho =\lambda +g^1\overline{}_{(\tau )}g`$, with $`\lambda `$. Such a gauge transformation can be found if the deformation parameter $`\theta `$ is small enough.

As a final remark, me mention that for special values of the deformation parameter, $`\theta =1/N`$, with $`N`$, the Gaussian (Boca) projection can be mapped, in the large torus limit , to the basic GMS soliton .

## Acknowledgments

The work of L.D. and G.L. was supported in part by the Progetto INdAM 2003 di ricerca interdisciplinare GNFM-GNSAGA.


text:
# New constraints on varying 𝛼

## I Introduction

Cosmology and astrophysics provide a laboratory with extreme conditions in which to test fundamental physics and search for new paradigms. Currently preferred unification theories Polchinski (1998) predict the existence of additional space-time dimensions, which have a number of possibly observable consequences, including modifications in the gravitational laws on very large (or very small) scales and space-time variations of the fundamental constants of nature Uzan (2002); Martins (2002). Recent evidence of a time variation of fundamental constants Webb et al. (2001, 2003); Ivanchik et al. (2003) offers an important opportunity to test such fundamental physics models. It should be noted that the issue is not if such theories predict such variations, but at what level they do so, and hence if there is any hope of detecting them in the near future.

The most promising case is that of the fine-structure constant $`\alpha `$, for which some evidence of time variation at redshifts $`z23`$ already exists Webb et al. (2001, 2003). Since one expects $`\alpha `$ to be a non-decreasing function of time Damour and Nordtvedt (1993); Barrow et al. (2002), it is particularly important to try to constrain it at earlier epochs, where any variations relative to the present-day value should be larger. The cosmic microwave background (CMB) anisotropies provide such a probe, being mostly sensitive to the epoch of decoupling, $`z1100`$. Avelino et al. (2000, 2001); Martins et al. (2002, 2003)

The reason the CMB is a good probe of variations of the fine-structure constant is that these alter the ionisation history of the universe Hannestad (1999); Kaplinghat et al. (1999); Avelino et al. (2000). The dominant effect is a change in the redshift of recombination, due to a shift in the energy levels (and, in particular, the binding energy) of Hydrogen. The Thomson scattering cross-section is also changed for all particles, being proportional to $`\alpha ^2`$. A smaller effect (which has so far been neglected) is expected to come from a change in the Helium abundance. As is well known, CMB fluctuations are typically described in terms of spherical harmonics, $`T(\theta ,\varphi )=_\mathrm{}ma_\mathrm{}mY_\mathrm{}m(\theta ,\varphi )`$ from whose coefficients one defines the angular power spectrum $`C_{\mathrm{}}=<|a_\mathrm{}m|^2>`$. Increasing $`\alpha `$ increases the redshift of last-scattering, which corresponds to a smaller sound horizon. Since the position of the first Doppler peak ($`\mathrm{}_{peak}`$) is inversely proportional to the sound horizon at last scattering, increasing $`\alpha `$ will produce a larger $`\mathrm{}_{peak}`$ Avelino et al. (2000). This larger redshift of last scattering also has the additional effect of producing a larger early ISW effect, and hence a larger amplitude of the first Doppler peak Hannestad (1999); Kaplinghat et al. (1999). Finally, an increase in $`\alpha `$ decreases the high-$`\mathrm{}`$ diffusion damping (which is essentially due to the finite thickness of the last-scattering surface), and thus increases the power on very small scales. These effects have been implemented in a modified CMBFAST algorithm which allows a varying $`\alpha `$ parameter Avelino et al. (2000, 2001). The changes were introduced in the subroutine RECFAST Seager et al. (1999) according to the extensive description given in Hannestad (1999); Kaplinghat et al. (1999).

## II Results and discussion

In Martins et al. (2003) we presented a up-to-date constraints on the value of the fine-structure constant at the epoch of decoupling from the recent observations made by the Wilkinson Microwave Anisotropy Probe (WMAP) satellite. In the framework of models considered, a positive (negative) variation of the value of $`\alpha `$ at decoupling with respect to the present-day value is now bounded to be smaller than $`2\%`$ ($`6\%`$) at $`95\%`$ C.L..

The likelihood distribution function for $`\alpha _{\text{dec}}/\alpha _0`$, obtained after marginalization over the remaining parameters, is plotted in Figure 1. We found, at $`95\%`$ C.L. that $`0.94<\alpha _{\text{dec}}/\alpha _0<1.01`$, improving previous bounds, (see Martins et al. (2002)) based on CMB and complementary datasets.

WMAP satellite data tightens the CMB constraints on the value of the fine-structure constant at the epoch of decoupling. As in other previous works Avelino et al. (2000, 2001); Martins et al. (2002), the current data is consistent with no variation,though the likelihood is skewed towards smaller values at the epoch of decoupling.

These previous works results were somewhat weakened by the existence of various important degeneracies in the data. This issue has been analysed by means of a Fisher Matrix Analysis (FMA) Martins et al. (2003); Rocha et al. (2003). Following Rocha et al. (2003) we present the precision with which cosmological parameters can be reconstructed using both CMB temperature and E-polarization measurements. We consider the WMAP experiment, the planned Planck satellite and an ideal experiment which would measure both temperature and polarization to the cosmic variance limit (in the following, ’CVL experiment‘). Cosmological models are characterized by the 8 dimensional parameter set

$$𝚯=(\mathrm{\Omega }_bh^2,\mathrm{\Omega }_mh^2,\mathrm{\Omega }_\mathrm{\Lambda }h^2,,n_s,Q,\tau ,\alpha ),$$
(1)
We assume purely adiabatic initial conditions and we do not allow for a tensor contribution. Our maximum likelihood model has parameters $`\omega _b=0.0200`$, $`\omega _m=0.1310`$, $`\omega _\mathrm{\Lambda }=0.2957`$ (and $`h=0.65`$), $`=0.9815`$, $`n_s=1.00`$,$`Q=1.00`$, $`\tau =0.20`$ and $`\alpha /\alpha _0=1.00`$.

The experimental parameters used for the Planck analysis are in Table 1, and we use the first 3 channels of the Planck High Frequency Instrument (HFI) only. For the cosmic variance limited (CVL) experiment, we set the experimental noise to zero, and we use a total sky coverage $`f_{\mathrm{sky}}=1.00`$. Although this is never to be achieved in practice, the CVL experiment illustrates the precision which can be obtained in principle from CMB temperature and E-polarization measurements.

Fig. 2 illustrates the effect of $`\alpha `$ and $`\tau `$ on the CMB temperature and polarization power spectra—see Rocha et al. (2003) for a more detailed discussion.

For the CMB temperature, reionization simply changes the amplitude of the acoustic peaks, without affecting their position and spacing (top left panel); a different value of $`\alpha `$ at the last scattering, on the other hand, changes both the amplitude and the position of the peaks (top right panel). The outstanding effect of reionization is to introduce a bump in the polarization spectrum at large angular scales (lower left panel). This bump is produced well after decoupling at much lower redshifts, when $`\alpha `$, if varying, is much closer to the present day’s value. If the value of $`\alpha `$ at low redshift is different from that at decoupling, the peaks in the polarization power spectrum at small angular scales will be shifted sideways, while the reionization bump on large angular scales won’t (lower right panel). It follows that by measuring the separation between the normal peaks and the bump, one can measure both $`\alpha `$ and $`\tau `$. Thus we expect that the existence of an early reionization epoch will, when more accurate cosmic microwave background polarization data is available, lead to considerably tighter constraints on $`\alpha `$.

Table 2 and Fig. 3 summarize the forecasts for the precision in determining $`\tau `$ and $`\alpha `$ (relative to the present day value) with Planck and the CVL experiment. We consider the use of temperature information alone (TT), E-polarization alone (EE) and both channels (EE+TT) jointly. Note that one could use the temperature-polarization cross correlation (ET) instead of the E-polarization, with the same results. As it is apparent from Fig. 3, TT and EE suffer from degeneracies in different directions, because of the reasons explained above. Thus combination of high-precision temperature and polarization measurements can constrain in the most effective ways both variations of $`\alpha `$ and $`\tau `$. Planck will be essentially cosmic variance limited for temperature but there will still be considerable room for improvement in polarization (Table 2). This therefore argues for a post-Planck polarization experiment, not least because polarization is, in itself, better at determining cosmological parameters than temperature. We conclude that Planck alone will be able to constrain variations of $`\alpha `$ at the epoch of decoupling within $`0.34\%`$ ($`1\sigma `$, all other parameters marginalized), which corresponds to approximately a factor 5 improvement on the current upper bound. On the other hand, the CMB alone can only constrain variations of $`\alpha `$ up to $`𝒪(10^3)`$ at $`z1100`$. Going beyond this limit will require additional (non-CMB) priors on some of the other cosmological parameters. This result is to be contrasted with the variation measured in quasar absorption systems by Ref.Webb et al. (2001), $`\delta \alpha /\alpha _0=𝒪(10^5)`$ at $`z2`$. Nevertheless, there are models where deviations from the present value could be detected using the CMB.

Finally Table 3 and Fig. 4 summarize the FMA results for all parameters for WMAP, Planck and a CVL experiment see Rocha et al. (2003) for further details.

## III Conclusions

We presented up-to-date constraints on the value of the fine-structure constant at the epoch of decoupling, using the WMAP satellite data. Within the set of models considered, a variation of the value of $`\alpha `$ at decoupling with respect to the present-day value is now bounded to be smaller than $`2\%`$ ($`6\%`$) at $`95\%`$ confidence level. We have proposed a way of using the existence of an early reionization epoch as suggested by WMAP, to improve these constraints. We have shown that CMB data alone will be able to constrain $`\alpha `$ up to the $`0.1\%`$ level. Tighter constraints than this will require invoking further (non-CMB) priors. These points are discussed in more detail in Rocha et al. (2003).


text:
# Deflation of Eigenvalues for Iterative Methods in Lattice QCD

## 1 INTRODUCTION

Removal of the smallest eigenvalues of a matrix can significantly improve the performance of iterative methods. This is important for lattice QCD where stochastic methods for disconnected diagrams and new matter formulations such as overlap fermions are computationally intensive. In addition, there is often a need for these systems to be shifted to take account of multiple quark masses or other parameters. We think that it is important for any iterative method we develop to be applicable to general non-Hermitian systems, for this is often the case for lattice systems.

GMRES  is an iterative method that can be used in lattice studies with stochastic noise methods. Convergence can be improved by deflating eigenvalues. Two years ago it was shown that deflating eigenvalues  can be especially helpful for lattice problems with multiple right-hand sides . Here we report on generalizations of these techniques. First we show how deflated restarted GMRES can accomodate multiple shifts. We give numerical results for a Wilson matrix with very small eigenvalues. We then consider the case of combining multiple shifts with multiple right-hand sides. Finally we also discuss and give numerical results on the possibility of combining deflated GMRES for the first system with a BiCGStab algorithm which uses those deflated eigenvalues for the subsequent solutions.

## 2 DEFLATION FOR MULTIPLY SHIFTED SYSTEMS

The shifted systems are expressed by

$`(As(i)I)x^i=b,`$ (1)
where the s(i) are the numerical shifts. The method proposed in  is based upon keeping the residual vectors for the shifted systems parallel to one another. That is, one uses only a single Krylov subspace to solve all the systems. This is possible for GMRES only by forcing the residuals to be parallel to one another after a restart and thus the minimum residual condition is applied to only the unshifted (most difficult) system. We find that this also results in an almost optimal solution for the shifted systems.

Our method for deflating GMRES (called GMRES-DR ) simulataneously solves linear equations and computes eigenvalues. This is possible because it has a Krylov subspace for solving the linear equations as well as Krylov subspaces designed for solving the eigenvalue equation. These approximate eigenvectors are refined and improved as the GMRES-DR cycles proceed. We have determined that GMRES-DR can be efficiently applied to multiply shifted systems because of its Krylov properties.

Fig.1 shows the result of a shifted run of GMRES-DR(50,30), i.e., the subspace is 50-dimensional and includes 30 approximate eigenvectors. It is compared to the residual norm for regular GMRES(50). One can see the deflation start to take place around 200 matrix-vector products (the approximate eigenvectors are accurate enough by then).

## 3 MULTIPLE SHIFTS AND ADDITIONAL RIGHT-HAND SIDES

Typically, many right-hand sides are used to do stochastic estimates of matrix elements. It would be beneficial if the deflated eigenvalues and eigenvectors calculated for the first right-hand side could be reused. Then deflation could begin from the start instead of having to wait for the required eigenvectors to be recalculated, with more dramatic results than in Fig.1. A method to accomplish this called GMRES-Proj was given in . GMRES-Proj alternates cycles of GMRES with projections over the approximate eigenvectors. Next we consider using this for the case of multiple shifts for each of the multiple right-hand sides.

GMRES-DR can deflate eigenvalues and handle multiple shifts at the same time for the first right-hand side (the right-hand sides for each shifted system can be kept parallel). However, for the second and subsequent right-hand sides there is no way to deflate eigenvalues by projecting over the previously computed approximate eigenvectors and still keep the right-hand sides for the different shifted systems parallel. If the eigenvectors were exact, this would be possible. So one approach is to use only fairly accurate approximate eigenvectors and do a correcting iteration at the end for each shifted system separately.

We now propose a better fix for this problem which does not require accurate eigenvectors. GMRES-DR produces approximate eigenvectors that are represented in the form of an Arnoldi-like recurrence relation:

$`AV_k=V_{k+1}\overline{H}_k,`$ (2)
where the columns of $`V_k`$ span the space of approximate eigenvectors, $`V_{k+1}`$ has one extra column added to $`V_k`$, and $`\overline{H}_k`$ is a $`(k+1)`$ by $`k`$ matrix. Using this, we can perform the projection over the approximate eigenvectors so that the right-hand sides for the different shifted systems are all multiples of each other, except for error in the direction of $`v_{k+1}`$, the last column of $`V_{k+1}`$. This error can be easily corrected at the end, if we are willing to solve one additional right-hand side, namely $`v_{k+1}`$.

There is also potential to develop a deflated version of BiCGStab for multiply shifted systems using this idea. This will be reported later, but for now we will look at deflating BiCGStab for multiple right-hand sides in the non-shifted case.

## 4 DEFLATED BICGSTAB

The deflation on the additional right-hand sides can employ other methods after the initial eigenvalue generation by GMRES-DR. We have looked at BiCGStab, a non restarted method popular in the lattice community, in this context. We again combine a projection over approximate eigenvectors with the iterative method. However, the projection can only be performed once before the BiCGStab iteration begins. Therefore this projection must reduce the critical eigencomponents to a low enough level that no further reduction is needed during the BiCGStab iteration. The minimal residual projection is the best possible one in the sense that it minimizes the norm of the residual vector. However, it can be shown that even if one knows the exact eigenvectors, the corresponding component in the residual vector will not generally be zeroed out (unless the matrix is Hermitian). To deal with this problem, we will give a projection that is better for this task. It uses both left and right approximate eigenvectors.

Left-right Projection for deflating eigenvectors:

1. Let the current approximate solution be $`x_0`$ and the current system of equations be $`A(xx_0)=r_0`$. Let $`V`$ be orthonormal with columns spanning the subspace of approximate right eigenvectors. Let $`U`$ similarly be orthonormal with columns from approximate left eigenvectors.

2. Solve $`U^{}AVd=c`$, where $`c=U^{}r_0`$.

3. The new approximate solution is $`x_k=x_0+Vd`$.

This projection can reduce the small eigenvector components better than minimum residual and in the case of exact left and right eigenvectors, the components can be zeroed out. Of course, it is necessary to compute left eigenvectors. However, for many lattice problems, including the Wilson-Dirac matrix, there is a simple relationship between the left and right eigenvectors.

Fig.2 above shows a run on a smaller $`4^4`$ lattice (using MATLAB) of the deflated BiCGStab algorithm using eigenvectors calculated from an initial GMRES-DR solution of the first right-hand side. The deflated solutions take a smaller number of iterations than the nondeflated implementation. Notice that deflation does not care about the nature of the right-hand sides used, so these may be noise vectors or units vectors, as one desires. This applies to all techniques used.

## 5 ACKNOWLEDGEMENTS

This material is based upon work supported by the National Science Foundation under Grants No. 0070836 (Theoretical Physics) and 0310573 (Computational Mathematics) and NCSA and utilized the SGI Origin 2000 System at the University of Illinois. RM acknowledges the Baylor University Summer Sabbatical Program.


text:
# Mass spectrum of N* and source optimization

## 1 INTRODUCTION

Lattice QCD calculations of baryon masses are based on the analysis of 2-point correlation functions, defined by

$$C_{ij}(t,t_0)=\underset{𝐱}{}0\left|B_i(𝐱,t)\overline{B}_j(0,t_0)\right|0,$$
(1)
where $`B_i`$ are color-singlet baryon source fields. Summation over $`𝐱`$ projects to zero total momentum. For large $`t`$ the contribution from the lowest mass eigenstate dominates. An effective mass is calculated from $`M_{\mathrm{eff}}=ln\left(\frac{C(t,t_0)}{C(t+1,t_0)}\right)`$.

Ideal sources maximally overlap the states of interest: $`n\left|\overline{B}_i(t_0)\right|0\delta _{ni}`$, where $`n`$ is an intermediate state. By using the orthogonality relation of basis vectors spanning irreducible representations (IRs) of the cubic group, one can attempt to realize this relation by taking $`B_i`$ to be a linear combination of source operators that have a good overlap with state $`n`$.

There are 48 discrete spatial rotations that conserve the cubic symmetry of spin $`1/2`$ objects, and they form a group, called “double covered octahedral group”, or $`{}_{}{}^{2}𝒪`$. The reduction of the group $`{}_{}{}^{2}𝒪`$ falls into three spinorial IRs: $`G_1,G_2`$, and $`H`$ with respective dimensions 2, 2, and 4. Table 1 lists the inverse relation of the reduction of SU(2) to $`{}_{}{}^{2}𝒪`$ .

Fermionic operators on a cubic lattice can be written in terms of basis operators that span IRs of $`{}_{}{}^{2}𝒪`$. Since lattice calculations preserve cubic symmetry, different IRs do not admix. However, each IR is an admixture of different total angular momenta. Total angular momentum is not a “good” quantum number on the lattice because continuous rotational symmetry is broken. It is expected that higher spin states have higher masses, therefore within a given IR the lowest spin state dominates at large $`t`$.

## 2 NONLOCAL SOURCES

A general expression of three-quark baryon operator would be written as

$`B_\mu (𝐱)={\displaystyle A(𝐱,𝐥_1,𝐥_2,𝐥_3,\mu ,\alpha ,\beta ,\gamma )ϵ^{abc}}`$
$`U_1^{a,a^{}}q_\alpha ^a^{}(𝐥_1)U_2^{b,b^{}}q_\beta ^b^{}(𝐥_2)U_3^{c,c^{}}q_\gamma ^c^{}(𝐥_3),`$ (2)
where $`\alpha ,\beta `$, and $`\gamma `$ are Dirac indices, $`a`$, $`b`$, $`c`$, $`a^{}`$, $`b^{}`$, and $`c^{}`$ are color indices, and $`𝐥_1,𝐥_2`$, and $`𝐥_3`$ are spatial displacement vectors. $`U_1^{a,a^{}}`$ is the gauge link connecting the spatial path from $`𝐱`$ to $`𝐥_1`$. The coefficients tie together a set of spatial displacement vectors $`𝐥_1,𝐥_2`$, and $`𝐥_3`$, and Dirac indices $`\alpha ,\beta `$, and $`\gamma `$ in such a way as to yield gauge invariant operators that transform according to IRs of $`{}_{}{}^{2}𝒪`$ and definite isospin.

In order to have sources that are commensurate with the size of baryons, it is appropriate to smear the sources, using for example, Gaussian smearing,

$$F_0=\left(1\frac{\sigma ^2D^2}{4N}\right)^N,$$
(3)
where $`\sigma `$ is the smearing width, $`N`$ is an integer, and $`D^2`$ is the gauge-covariant three-dimensional Laplacian operator. For different radial distributions, additional powers of Laplacian operators are used, e.g., $`F_2=D^2F_0`$, $`F_4=D^2F_2`$, and so on.

The simplest nonlocal operator has $`𝐥_1=𝐥_2=𝐱`$, $`𝐥_3=𝐱+a\widehat{𝐞}_1`$. Applying spatial rotations of the cubic group to this operator leads to $`6\times 8=48`$ dimensional representations, where six is for the directions of $`a\widehat{𝐞}_1`$ and eight for the three spin degrees of freedom. Reduction of these representations results in 6 $`G_1`$’s, 2 $`G_2`$’s, and 8 $`H`$’s. From these “one-link” constructions, we used both positive- and negative-parity nucleon sources from $`G_1`$ basis in the simulation.

More complicated sources have been constructed by choosing the initial operator as, for instance, $`𝐥_1=𝐥_2=𝐱`$, $`𝐥_3=𝐱+a\widehat{𝐞}_1+a\widehat{𝐞}_2`$ in Eq. (2).

## 3 SIMULATION DETAILS AND PRELIMINARY RESULTS

For a test of these sources, we formed two $`3\times 3`$ matrices of correlation functions with positive-parity and negative-parity, by using projection operators $`1\pm \gamma _4`$. All operators transform according to $`G_1`$ IRs. Isospin is selected to be $`I=I_z=1/2`$. The operators used in the simulation are shown below,

$`s_{1/2}^{(\pm )}\left[d(C\gamma _5)u\right](1\pm \gamma _4)(1+\gamma _3\gamma _5)u,`$ (4)
$`D^2s_{1/2}^{(\pm )}\left[d(C\gamma _5)u\right](1\pm \gamma _4)(1+\gamma _3\gamma _5)D^2u,`$ (5)
$`p_{1/2}^{(\pm )}\left[d(C\gamma _5)u\right](1\gamma _4)(1+\gamma _3\gamma _5)D_+u`$
$`+\left[d(C\gamma _5)u\right](1\gamma _4)(1\gamma _3\gamma _5)D_zu,`$ (6)
where $`D_i`$ is a covariant derivative and $`D_+=D_x+iD_y`$. Color and Dirac indices are implicit here. The upper sign is for positive-parity and the lower sign is for negative-parity. All quark fields are smeared using Eq. (3). Equation (4) is the widely used $`G_1`$ irreducible source, which mainly couples to the ground state of the nucleon: we call it “$`s_{1/2}^{(\pm )}`$” even though from Table 1 one sees that there can be admixtures with spin 7/2 or higher. Equation (5) is s-wave with an additional Laplacian operator acting on the third quark. Equation (6) is constructed from one-link sources as described above. The derivative operators and spin-projection operators in Eq. (6) that act on the third quark can be rewritten as $`D_+u_{1/2}+D_zu_{1/2}`$. This corresponds to the Clebsch-Gordan formula for a $`p_{1/2}`$ state, i.e.,

$$p_{1/2}=\sqrt{\frac{1}{3}}Y_{11}\chi _{1/2}+\sqrt{\frac{2}{3}}Y_{10}\chi _{1/2},$$
(7)
where $`Y_{lm}`$ represents a discretized version of a spherical harmonic, and $`\chi _{\pm 1/2}`$ represent up and down two-component spinors. Because of this correspondence, Eq. (6) is called the “$`p_{1/2}^{(\pm )}`$” source.

We use $`16^3\times 32`$ isotropic lattice with periodic boundary conditions. The smearing width that is used in Eq. (3) is $`\sigma =4.5`$, with $`N=50`$. We use 221 configurations111Configurations are from the NERSC lattice archive. generated from the isotropic Wilson gauge action at $`\beta =6.0`$ and the Wilson fermion action with $`\kappa =0.1550`$, which corresponds to pion mass about 542 MeV.

We adopt a variational method to extract effective masses . Correlation matrices are diagonalized at a time slice $`t_v`$. This produces three orthogonal eigenvectors, each of the form $`𝐯=(v_1,v_2,v_3)`$, where $`v_i`$’s are the coefficients of $`s_{1/2}^{(\pm )}`$, $`D^2s_{1/2}^{(\pm )}`$, and $`p_{1/2}^{(\pm )}`$, respectively. The resulting eigenvectors are then used to calculate “diagonal” correlators of the form $`𝐯^TC(t)𝐯`$ at other times. For positive-parity the eigenvectors obtained at $`t_v=10`$ give the most clear mass splitting, which is shown in Fig. (1). Filled diamonds have $`𝐯(0.44,0.57,0.68)`$, i.e., each source is almost equally weighted. This state reaches a plateau for $`5t10`$, suggesting that the ground state of the nucleon has been obtained. Filled circles and open boxes are for other orthogonal eigenvectors. The effective masses based on the other orthogonal eigenvectors are larger than the ground state mass, indicating that excited states of positive-parity are being resolved.

For negative-parity, the eigenvectors at $`t_v=6`$ give the most clear mass splitting, which is shown in Fig. (2). Filled diamonds have $`𝐯=(0.37,0.56,0.71)`$, where the contribution from the $`p_{1/2}^{()}`$ state is significant. The effective mass of this state is about 40% larger than the positive-parity ground state; the corresponding $`N^{}(1/2^{})`$ physical state is about 63% heavier than the nucleon. Filled circles and open boxes have larger error bars, but they indicate larger masses than diamonds.

The present results are preliminary. Work is under way to improve the chiral properties of the action and to use a large set of sources with anisotropic lattices.

This work was supported by the U.S. National Science Foundation under Awards PHY-0099450 and PHY-0300065, and by the U.S. Department of Energy under contract DE-AC05-84ER40150 and DE-FG02-93ER-40762.