format: text

text:
# Surface imaging of HD~199178 (V1794 Cygni) Based on observations made with the 2 m telescope of the National Astronomical Observatory, Rozhen, Bulgaria, and the Nordic Optical Telescope, operated on the island of La Palma jointly by Denmark, Finland, Iceland, Norway, and Sweden, in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofisica de Canarias.

## 1 Introduction

Doppler or surface imaging is, at present, the most powerful tool to map spots on rapidly rotating late-type stars. Even so, the Doppler imaging maps should always be studied critically. The method suffers from difficulties in modeling spectral lines in late-type stars, as well as from uncertainties in stellar and spectral parameters. These problems have been amply illustrated by conflicting results. For instance, the temperature maps obtained with some Doppler imaging methods most often show very large cool spots covering the poles (cf. Vogt vogt1 (1988); Strassmeier hd106225 (1994); Hatzes & Vogt hv (1992)), while the results with the methods developed by Piskunov (pisk (1991)) and Berdyugina (occam (1998)) rarely show polar spots (cf. Piskunov et al. si1 (1990); Strassmeier et al. eieri (1991); Korhonen et al. heidi2 (2000)). Since partly the same stars, and once even the same observations (Strassmeier et al. eieri (1991)), have been modeled, it is evident that the results depend either on the applied methods or the selected stellar and spectral parameters.

There are several studies on how line and stellar parameters influence Doppler imaging (cf. Unruh & Collier Cameron unruh (1995)). But the influence of one important phenomenon, namely surface differential rotation, has usually been neglected. Many active late-type stars show clear variations in their photometric rotation periods (Strassmeier et al. klasu1 (1989); Cutispoto cuti (1998)). Variable periods have also been discovered in the Ca ii H&K emission (Donahue et al. donahue (1996)). These periodicity changes are usually explained as a signature of surface differential rotation combined with latitudinal activity migration. Apart from distorting the spot pattern, differential rotation would alter the shape of the line profile in the same way as a large cool or hot polar spot. This emphasizes the necessity of implementing surface differential rotation into the Doppler imaging techniques.

The generally accepted empirical result is that differential rotation is weaker in more rapidly rotating stars (Hall hall (1991)). Hence the effects of differential rotation have usually been bypassed in Doppler imaging. Although antisolar differential rotation has previously been noted to alter the line profiles in a similar way as a large cool polar spot, Hatzes et al. (polspot (1996)) concluded that this would require an unreasonably strong differential rotation.

HD~199178 (V1794 Cygni) is a rapidly rotating single G5 subgiant (Herbig herbig (1958)), which was classified as a FK Comae-type star by Bopp & Rucinski (bopp1 (1981)). Jetsu et al. (jetsu1 (1990)) derived a photometric rotation period of $`P_{\mathrm{phot}}=3\stackrel{d}{.}337484`$ $`\pm 0\stackrel{d}{.}000043`$. Furthermore, they reported a 9.1 year period in the mean brightness level and indications of a shorter cycle in the light curve amplitude. In a recent analysis of photometry between 1975 and 1995, Jetsu et al. (1999b ) showed that the photometric period varied and interpreted this as a signature of differential rotation.

We present new temperature maps of HD~199178 for the years 1994 and 1995. The maps were calculated with the surface imaging technique developed by Piskunov (pisk (1991)). Previous Doppler images of this star display a large cool polar spot and some smaller high and low latitude features (Vogt vogt1 (1988); Strassmeier et al. klasu2 (1999)). Unlike these earlier studies of HD~199178, we take into account the possible effect of differential rotation on the rotational broadening of photospheric absorption lines. We also demonstrate how different values for the differential rotation coefficient alter the retrieved surface imaging map.

## 2 Stellar parameters of HD 199178

HD~199178 has an apparent visual magnitude of $`V7\stackrel{m}{.}2`$. The Hipparcos-satellite parallax measurement of $`\pi 10.68\pm 0.73`$ mas gives a distance of 94 pc (Perryman hipp (1997)). The long-term mean colour indices $`BV0.79`$ and $`VR0.68`$ determined by Jetsu et al. (jetsu1 (1990)) suggested an average “background” stellar surface temperature of $`T_{\mathrm{eff}}`$ 5375 K. The above results agree with the G5iii-iv spectral classification by Herbig (herbig (1958)).

The new Hipparcos distance measurement allows us to estimate the absolute magnitude and radius of the star. We use $`7\stackrel{m}{.}09`$ as the apparent visual magnitude of the unspotted star (Strassmeier et al. klasu2 (1999)). HD~199178 lies in the direction of the North America and Pelican Nebulae complex, but in the foreground. Assuming an interstellar absorption of $`0\stackrel{m}{.}2`$ (Straiz̆ys et al. ext (1989): Fig. 2) yields an absolute visual magnitude of $`2\stackrel{m}{.}02`$. A bolometric correction of $`BC=0.184`$ (Strassmeier et al. klasu2 (1999)) gives the bolometric magnitude $`M_{\mathrm{bol}}1\stackrel{m}{.}84`$. Finally, if the unspotted temperature of $`T_{\mathrm{eff}}=5400`$ K is combined to the solar $`T_{}=5800`$ K and $`M_{\mathrm{bol}_{}}=4\stackrel{m}{.}64`$, the Stefan-Boltzmann relation gives $`R4.2R_{}`$.

The Barnes-Evans relation is based on the same principle as the above, except that the value of the bolometric correction is not required. We use the formulation by Lacy (lacy (1977)):

$$\mathrm{lg}(R/R_{})=7.47240.2V_02F_v+\mathrm{lg}d,$$
(1)
where $`F_v=3.9770.429(VR)_0`$ and $`[d]=`$ pc. The combination $`V_0=7\stackrel{m}{.}0`$, $`VR=0.68`$ and $`E_{VR}=0.05`$ gives $`F_v3.707`$ and $`R(4.3\pm 0.7)R_{}`$. (Here $`V_0`$ is the mean apparent visual magnitude without interstellar absorption.)

The rotation period $`P3\stackrel{d}{.}3`$ and projected rotational velocity $`v\mathrm{sin}i`$ $``$ 70 km s-1 give a lower limit $`R\mathrm{sin}i4.6R_{}`$ for the radius. There is thus a slight discrepancy between the different radius determinations. Anyway, the above estimates and the lower limit $`R\mathrm{sin}i`$ suggest that $`4.6R_{}R5.0R_{}`$. This is equivalent to an inclination of $`67^{}i90^{}`$. The exceptional nature of HD~199178 (rapid rotation, spots, possibly a coalesced binary) may mean that the equatorial radius will be underestimated. In our calculations we will thus allow for an inclination $`i50^{}`$.

A variable photometric rotation period was detected in the long-term photometry by Jetsu et al. (1999b ). They analyzed 114 sets of photometric data obtained between 1975 and 1995. Each set contained at least 8 nights of observations and had a length of about 30 nights (Jetsu et al. 1999a ). The regular yearly change rate of the photometric period was 3.3 %. In addition, irregular seasonal changes of 7.5 % were found. The varying photometric period was explained by surface differential rotation quantified by $`|\alpha |0.075`$, where

$$\alpha =(\mathrm{\Omega }_{\mathrm{equator}}\mathrm{\Omega }_{\mathrm{pole}})/\mathrm{\Omega }_{\mathrm{equator}}.$$
(2)
Even this lower limit for the differential rotation is considerably larger than that expected from the results by Hall (hall (1991)). If the latitudinal range of the spot activity does not extend all the way from the equator to the pole, then $`|\alpha |`$ could even reach the solar value of 0.2. We thus assume that the realistic differential rotation coefficient values are within the limits $`0.20\alpha 0.20`$.

Differential rotation would alter the shape of the line profile. The case of $`\alpha >0`$ induces a “sharper” profile, while $`\alpha <0`$ causes a more “flat bottomed” profile (cf. Bruning bruning (1981)). Thus a negative $`\alpha `$ modifies the line profile much in the same way as a cool polar spot.

Since there is no unique photometric period, in principle different periods could be used for different seasons. But this would make it harder to compare the surface images. Therefore we chose to use the mean period $`P=3\stackrel{d}{.}3250`$ calculated for the years 1994 an 1995. Our ephemeris for the rotational phases was

$$\mathrm{HJD}_{\mathrm{min}}=(2449511.150\pm 0.070)+(3.3250\pm 0.0012)E,$$
(3)
where the $`\mathrm{HJD}_{\mathrm{min}}`$ is the epoch of the light curve minimum (Jetsu et al. 1999b ).

## 3 The surface imaging technique

Spots on the surface of a rapidly rotating star will cause distortions (“bumps”) in the observed spectral line profiles. As the star rotates, these bumps will move across the absorption line profiles. The main idea of surface imaging is to trace these distortions and retrieve a surface map of the star. We use the surface imaging technique developed by Piskunov (Piskunov et al. si1 (1990), Eqs. 3–6 and 8; Piskunov pisk (1991), Eqs. 5–8), but also include Johnson $`V`$-photometry in the inversion.

Local line profiles and $`V`$-band fluxes are calculated for a set of temperatures and limb angles using stellar model atmospheres. This table is then used for the disk integration of a given surface temperature distribution. The surface imaging problem is solved by searching for a surface temperature distribution that minimizes the discprepancy between the observations and the calculated line profiles and light curve.

Since we always have a confined temperature range of model atmospheres, the solution should be limited to temperatures between the minimum and maximum of the models ($`T_{\mathrm{min}}`$ and $`T_{\mathrm{max}}`$). We do this by adding a “penalty function” to the quantity to be minimized. We use an integral of a double-sided Heaviside step function, which is nearly constant between $`T_{\mathrm{min}}`$ and $`T_{\mathrm{max}}`$, but grows nearly linearly when the temperature is below $`T_{\mathrm{min}}`$ or exceeds $`T_{\mathrm{max}}`$. Thus the penalty function merely prevents extreme temperatures in the solution.

The crucial point of surface imaging is how the line profiles for a given temperature distribution are calculated. Errors in the calculated line profiles will introduce artifacts in the resulting map and prevent the solution from reaching the observational signal-to-noise level. Errors in the line profiles can arise from numerous problems, e.g. wrong line parameters, missing lines, non-LTE effects, wrong limb darkening or inadequacies of the model atmospheres. The fact that the lines used for Doppler imaging are strong absorption lines, whose cores are formed high in the photosphere, puts extra demands on the line modeling.

The errors in the line profiles are not the only source of systematic errors. For example, an erroneous $`v\mathrm{sin}i`$ will manifest itself as a cool or hot belt in the image (cf. Berdyugina occam (1998)). Errors in the microturbulence velocity will mainly change the average temperature of the solution. Doppler images are not sensitive to errors in the macroturbulence velocity (cf. Unruh & Collier Cameron unruh (1995)). A wrong inclination will cause a latitudinal displacement and a slight distortion of the spot. Combined with surface differential rotation, the inclination will also influence the line profile. Decreasing the inclination will increase the effect of differential rotation on the line profile (cf. Bruning bruning (1981)).

Differential rotation will, of course, also change the spot distribution, since different latitudes have different rotational periods. In this study we do not take this effect into account. Firstly, differential rotation becomes important in this sense only if the locations of spots extend over a large latitude range. Secondly, we assume that the magnetic field strength is sufficient to resist the shearing effect of differential rotation within the spot(s).

Since it is not (yet) possible to calculate perfect line profiles, one can always expect some artifacts in Doppler images. Fortunately, most of these artifacts can easily be identified, since certain shapes of artifacts are related to particular error sources. The same error in the overall shape of all calculated line profiles may manifest itself as a spot band at a certain latitude (or a polar spot). However, if the rotational phase coverage is inadequate, this artifact might not be axisymmetric. Minor errors in the calculations (e.g. missing lines in blends) or the data reduction (e.g. remnants of cosmic spikes or uncertainty in the continuum normalization) may cause thin high contrast features, like vertical stripes, arches or even ovals. For example, the remains of a cosmic spike in one spectrum may cause a thin high contrast artifact at a stellar longitude corresponding to the wavelength of the spike and rotation phase of the observation. This in turn might cause additional hot and cold features to cancel out the effect of the artifact in the remaining observed phases. A missing line, on the other hand, may lead to an artifact repeated at each observed phase. Noise in the spectra might be interpreted as a series of small “bumps”. Because such bumps will appear at arbitrary wavelengths, the resulting artifacts cannot be latitudinally confined, and will therefore be seen as vertical stripes in the image.

It is clear, that the retrieved longitude for the spots will usually be more reliable than the latitude. In some cases the longitude of the spot can be estimated even with a single spectrum, while the determination of the latitude requires a good phase coverage. Moreover, the latitude is sensitive to errors in the inclination and the general shape of the rotationally broadened line profile. The influence of all error sources is greatly reduced by using many spectral lines and simultaneous photometry in the inversion.

## 4 Observations

High resolution spectra of HD~199178 were measured in July 1994, August 1994, November 1994 and July 1995. The July 1994 data set was obtained with the 2 m Ritchey-Chretién telescope coudé-spectrograph of the National Astronomical Observatory (Rozhen), Bulgaria. The rest of the spectroscopic observations were made using the SOFIN high resolution échelle spectrograph at the 2.5 m Nordic Optical Telescope (NOT), La Palma, Spain. The observed spectral regions were $`64096414`$ Å (August 1994 and July 1995), $`64276441`$ Å (all seasons) and $`75097514`$ Å (August 1994, November 1994 and July 1995). The typical signal-to-noise ratio was between 200 and 300. The spectral observations are summarized in Table 1.

The spectral observations were reduced with the 3A software system (Ilyin ilyin (2000)). The reduction included bias, cosmic ray, flat field and scattered light correction, wavelength calibration and normalization. Special care was taken during the last two steps, where errors can cause artifacts in the surface images.

The calibration of the wavelength scale was obtained from Th-Ar comparison spectra. For the échelle spectra, we used all available spectral orders of the comparison frames to calculate a two-dimensional ($`k\lambda `$) wavelength solution, taking into account distortions of the line positions (Ilyin ilyin (2000)). Small shifts between the stellar and comparison spectra (e.g. due to flexure of the spectrograph) were corrected using atmospheric lines. The spectra were rebinned to a common wavelength grid.

The photospheric absorption lines of HD~199178 are broad and it is hard to determine the real continuum in the spectrum. Therefore the continuum normalization had to be done in two steps. A “quasi-continuum” consisting of points near to the real continuum level was used as a first approximation. The final normalization was done by utilizing a rotationally broadened synthetic spectrum.

Simultaneous or nearly simultaneous Johnson $`V`$-photometry for July 1994, August 1994 and July 1995 was obtained with the 60 cm telescope at the Mount Maidanak Observatory and the APT Phoenix-10 telescope (subsets SET=108, 109 and 111, Jetsu et al. 1999a ).

## 5 Temperature maps

The following three spectral regions were used for surface imaging: 6409.0 – 6413.5 Å (strongest line Fe i $`\lambda 6411.649`$ Å), 6427.5 – 6441.0 Å (Fe i $`\lambda 6430.846`$ Å; Fe ii $`\lambda 6432.680`$ Å; Ca i $`\lambda 6439.075`$ Å) and 7509.0 – 7513.5 Å (Fe i $`\lambda 7511.020`$ Å). The local line profiles were calculated for 20 limb angles and temperatures 3500 – 6500 K (with a step of 250 K) using a code written by Berdyugina (sloc (1991)) and atmospheric models by Kurucz (kurucz (1993)) with solar element abundances. Spectral line parameters were obtained from the Vienna Atomic Line Database (Kupka et al. vald2 (1999)). All significant atomic lines and the most important molecular lines (CN and TiO) were included in the calculations. The final adjustment of the $`gf`$ values was done after analyzing the wavelength dependence of the deviation between a synthetic spectrum of HD~199178 and the mean of the observations. Of the 111 lines used in the calculations, the $`gf`$ value of 36 lines was changed. The changes in the strongest lines, Fe i $`\lambda 6411.649`$ Å, Fe i $`\lambda 6430.846`$ Å, Fe ii $`\lambda 6432.680`$ Å, Ca i $`\lambda 6439.075`$ Å and Fe i $`\lambda 7511.020`$ Å, were $`0.1<\mathrm{\Delta }\mathrm{lg}(gf)<0.02`$. It should be emphasized that the adopted $`gf`$ values may be biased by uncertainties in element abundances. A list of all atomic lines that we used, can be found in Table 2.

The temperature maps were derived using the INVERS7-9PD code, originally written by N.E. Piskunov, with some changes made by T. Hackman. We first tested different values of rotational velocity (68 km s-1$``$ $`v\mathrm{sin}i`$$``$ 74 km s-1), microturbulence (1.2 km s-1$`v_{\mathrm{micro}}`$ 1.8 km s-1), inclination ($`50^{}i85^{}`$) and surface differential rotation ($`0.20\alpha 0.20`$). Differential rotation was implemented by adjusting the angular rotational velocity at each latitude $`b`$ using the Maunder rotation formula

$$\mathrm{\Omega }(b)=\mathrm{\Omega }_{\mathrm{equator}}(1\alpha \mathrm{sin}^2b).$$
(4)
The best solution (the smallest deviation between the spectroscopic observations and the calculated line profiles) was obtained with the stellar parameters given in Table 3.

Before combining all available spectral regions, separate images were calculated for each line blend in order to check for inconsistencies. The individual maps for each spectral region gave the same main features, but contained more artifacts (small high contrast features, especially at low latitudes) than the maps obtained from the combination of all spectral regions.

All surface images show a large cool spot remaining at approximately the same latitude (Figs. 14 and Table 4). The mean temperature of the visible surface is 5300 K. The polar region is about 4900 K. The lower latitude features in the images are probably artifacts. In all maps there are symmetry effects typical for surface imaging of stars with relatively high inclination ($`i60^o`$), i.e. “shadows” of the high latitude spot below the equator. The bright features, reaching a maximum of $`T_{\mathrm{eff}}`$ 5800 K, are probably also artifacts.

The July 1994 image shows a large high latitude cool spot (3700 K) with a longitudinal appendage (3900 K). In the August 1994 image, the coolest part of the spot (3700 K) is more concentrated and the appendage (4200 K) is weaker. In November 1994, the spot appears to be partly disintegrated. The coolest part of the spot (4100 K) seems to be at the position of the previous appendage. No photometric data was available for the November 1994 image, but there were clear changes in the light curves for July – September 1994 (subsets SET=108–110, Jetsu et al. 1999b ), which support the spot evolution. The July 1995 image resembles the July 1994 image, but with an additional new appendage stretched “southwards”. The coolest part of the spot (3700 K) has regained its strength.

The best surface imaging solution was achieved with the differential rotation coefficient $`\alpha =0.17`$. This implies strong antisolar differential rotation in HD~199178, i.e. that the pole would rotate faster than the equator. It seems that current models for global circulation in outer stellar convection zones predict only solar-like differential rotation (cf. Kitchatinov & Rüdiger kitch (1999)). Therefore we also wanted to test the possibility of no differential rotation and solar-like differential rotation. Our Fig. 5 displays the July 1995 surface imaging solutions obtained with $`\alpha =0`$ (rigid body rotation) and $`\alpha =0.15`$ (solar-like differential rotation). The best fits with $`\alpha =0`$ and $`\alpha =0.15`$ were achieved by increasing $`v\mathrm{sin}i`$ to 72 km s-1 and 74 km s-1, respectively. The mean deviation ($`d`$) of the calculated line profiles from the observations was slightly higher for $`\alpha =0`$ and $`\alpha =0.15`$ than for $`\alpha =0.17`$, i.e. the antisolar case provided a better fit to the data (Table 5).

The images in Fig. 5 demonstrate the difficulties in separating the effect of a polar spot from that of differential rotation. With no differential rotation the polar region is $``$ 800 K cooler than the mean surface temperature, while $`\alpha =0.15`$ gives a cool polar spot with $`T_{\mathrm{eff}}`$ 3900 K, much like the results by Strassmeier et al. (klasu2 (1999)). The same general effect was verified for all other seasonal images. Except for the polar cap, the images in Fig. 5 are very similar. The coolest part of the spot remains at approximately the same position. Similar artifacts can also be seen in all three images, in particular some bright features and smaller spots near the equator.

## 6 Differential rotation or polar spot?

The surface imaging solution with antisolar surface differential rotation might raise the question of whether this is a real phenomenon or just a trick to eliminate a cool polar spot. Thus, what we need is a way to determine which one, the polar spot or differential rotation, gives a better solution. When analyzing the T Tauri star Sz 68, Johns-Krull & Hatzes (sz68 (1997)) tackled this problem by calculating separate temperature maps for different spectral lines. They argued for solid body rotation and a polar spot, because this alternative gave a larger correlation between the individual maps.

We tested whether this approach works in our case. First we simulated synthetic spectroscopic “observations” using as input a stellar surface with a single high latitude spot, but no polar spot, and stellar parameters $`v\mathrm{sin}i`$=70 km s-1, $`\alpha =0.2`$, $`i=60^{}`$ (CASE i). To imitate real observations, we added normally distributed noise ($`\sigma =0.30\%`$) to the synthetic observations. These synthetic spectra were then used to calculate temperature maps with two sets of stellar parameters: $`v\mathrm{sin}i`$=70 km s-1, $`\alpha =0.2`$, $`i=60^{}`$ and $`v\mathrm{sin}i`$=72 km s-1, $`\alpha =0.0`$, $`i=60^{}`$. Maps were calculated both using all spectral regions and separately for each spectral region (6411Å, 6431 Å, 6439 Å and 7511 Å). With the correct stellar parameters, we retrieved a surface imaging solution consistent with the input stellar surface. As expected, neglecting the differential rotation, i.e. using $`\alpha =0.0`$, resulted in a solution with a cool polar spot in addition to the high latitude spot.

For both sets of stellar parameters we calculated the correlations between the maps retrieved for the different spectral lines. The test results showed that the correlation was higher when the differential rotation was wrongly neglected (Tables 6 and 7). The reason for this is clear. Neglecting a substantial differential rotation, one always obtains the same artifact, a polar spot, regardless of the spectral line used for the image. The presence of the same polar spot in all images will naturally increase the correlation between the images. The correlation only measures the consistency between the calculated images. It will not reveal the inconsistency between the original data (no polar spot) and the solution (polar spot).

An alternative approach to check the reliability of Doppler imaging maps, is to estimate the systematic errors in the calculated line profiles. As opposed to noise, the systematic errors of neighbouring wavelength points will correlate. In other words, systematic errors will cause the neighbouring data points to tend to be on the same side of the calculated curve. Not surprisingly, the linear correlation coefficient $`r`$ of neighbouring residuals, i.e. the systematic error, was significantly larger for the solution with the (false) solid body rotation than for the one with the (correct) differential rotation (Table 8: 1st and 2nd line). The systematic error could also be seen in the slightly larger deviation $`d`$.

To show that this criterion (lowest correlation of residuals) does not favour any particular spot distribution, we also simulated synthetic observations with a polar spot and no differential rotation (CASE II). Again the simulated data was used to retrieve temperature maps with and without differential rotation. The correct stellar parameters gave a solution consistent with the input map. In the solution with artificial differential rotation the polar spot was reduced, as expected. But also in this case, an erroneous differential rotation coefficient caused correlation in the residuals and increased the deviation (Table 8: 3rd and 4th line).

We conclude, that aiming for a maximum correlation between images obtained from different spectral lines will favour solutions with polar spots. Furthermore, systematical errors in the calculated lines profiles increase, if the value of $`\alpha `$ is erroneous.

In the simulated tests we correlated the errors of neighbouring points. But in our real data of HD~199178 neighbouring residuals will anyway be correlated, because the spectra have been rebinned to a common wavelength scale. This effect can be avoided by comparing each point in the spectra to the point two steps ahead. The correlation of next-to-neighbouring residuals was checked for the July 1995 surface imaging solutions with the differential coefficients $`\alpha =0.17`$, $`\alpha =0.0`$ and $`\alpha =0.15`$ (Fig. 5).

The results in Table 5 showed that the solution with $`\alpha =0.17`$ gave both the best fit (lowest $`d`$) and smallest correlation of the residuals (lowest r). However, we note that the correlation was highly significant even for the best stellar parameters ($`v\mathrm{sin}i`$=70 km s-1, $`\alpha =0.17`$, $`i=60^{}`$). Apparently there still remained systematic errors in the solution. This could also be seen from the fact that the deviation was larger than the observational noise. The most likely cause for these systematic errors were inaccuracies in the calculated line profiles (missing lines, errors in spectral parameters or element abundances, non-LTE effects, etc.). Even though all systematic errors could not be eliminated, we consider negative surface differential rotation more likely than the polar spot, because the best quality of the surface imaging solution was achieved with $`\alpha =0.17`$. Still, we emphasize that the $`\alpha `$-value was used to adjust the calculated profiles to fit the observations. Therefore our estimate for $`\alpha `$ may be biased by errors in the line calculations.

The presence of differential rotation is supported by the variation of the photometric period, although photometry alone can only give a lower limit for $`|\alpha |`$. Inserting the means of the spot latitudes and photometric periods of Table 4 into Eq. 4, the value of $`\alpha =0.17`$ would give the rotation periods at the equator $`P_{\mathrm{equator}}3\stackrel{d}{.}81`$ and the pole $`P_{\mathrm{pole}}3\stackrel{d}{.}25`$. The reliable photometric periods of HD~199178 determined by Jetsu et al. (1999b ) were within the range $`3\stackrel{d}{.}191\pm 0.017P_{\mathrm{phot}}3\stackrel{d}{.}559\pm 0.017`$. Since HD~199178 is a rapid rotator, one would expected it to have only high latitude spots (Schüssler & Solanki schuss (1992)). Therefore it is natural, that we do not find the equatorial rotation period from photometry. But the fact that our $`P_{\mathrm{pole}}`$ is larger than the smallest measured $`P_{\mathrm{phot}}`$ indicates a slight discrepancy.

By combining the time series analysis of photometry with individual surface images we can study if the photometric rotation period depends on the latitude of the main spot. Antisolar differential rotation would be seen as a shorter period for higher latitude spots. We note that this is indeed the case for the photometric periods and latitudes of the spot in Table 4. But since we have only three measurements, this tendency is far from statistically significant. Furthermore, the differences in the spot latitudes obtained from the surface images are too small to draw reliable conclusions.

There are other methods to measure the surface differential rotation. The approach most commonly used in Doppler imaging, is to cross-correlate constant latitude slices of images and calculate the shift in longitude as a function of the latitude (cf. Donati & Collier Cameron abdor (1997)). This approach requires that there really are long-lived spots at different latitudes. The method may also fail if there are too many artifacts in the images. Artifacts will often appear at nearly the same longitude as the real features or as vertical stripes in the images. A cross-correlation analysis of such images will give a nearly rigid body rotation. In our images, there is only one major spot, and the low latitude features are probably artifacts. For this reason, we do not attempt to cross-correlate our images.

Antisolar differential rotation has previously been found in some active late-type stars. Johns-Krull (ttau (1996)) analyzed the rotationally broadened spectral line profiles of three T Tauri stars and got the best fits with antisolar differential rotation. But he also noted that a polar spot could produce this effect. Weak antisolar differential rotation has been found in Doppler images of the RS CVn stars UX Ari (Vogt & Hatzes uxari (1996)), HR 1099 (Vogt et al. hr1099 (1999)) and HD 106225 (Strassmeier hd106225 (1994); Hatzes hatzes (1998)).

## 7 Conclusions

The surface images of HD~199178 presented here resemble the results recently obtained for FK Comae Berenices (Korhonen et al. heidi (1999) & heidi2 (2000)). The dominating feature is a long-lived high latitude spot (or spot group), with slightly changing position, shape and strength. Such changes could already be seen in the July 1994 and August 1994 images. In the November 1994 image, the spot had weakened and partly disintegrated. In July 1995, the spot had recovered its strength. The changes in the photometric light curves of HD~199178 during 1994 and 1995 can thus be explained by evolution of a single large spot rather than by abrupt changes of the whole spot configuration.

The long-lived spot structure suggests there is a stable active longitude in HD~199178. The high latitude of the spot could be explained by that the Coriolis force dominates over the buoyancy force in the magnetic flux tubes erupting from the deep parts of the convection zone (Schüssler & Solanki schuss (1992)). The changes in the spot latitude could then be connected to variations in the magnetic field strength. The changes in the latitude will in turn lead to drifts in the longitude if there is surface differential rotation.

In previous temperature maps of HD~199178, the “flat bottomed” absorption line profiles have been interpreted as a signature of a large and very cool polar spot (Vogt vogt1 (1988); Strassmeier et al. klasu2 (1999)). We achieved the best fit to the spectral observations by introducing an antisolar surface differential rotation of $`\alpha =0.17`$. With $`\alpha 0`$, a polar spot appeared in the images. We considered the solution with antisolar differential rotation to be more likely than the polar spot, because it reduced both the discrepancy between the observed and calculated line profiles, as well as the correlation between residuals of neighbouring wavelength points (i.e. systematical errors in the calculated line profiles).

During 1994 and 1995 the photometric rotation period varied by 3 %. Jetsu et al. (1999b ) interpreted these variations as a signature of surface differential rotation. Our images show that the photometric minimum of the light curves was apparently caused by the same persistent high latitude spot (group). Therefore the changes in the photometric period could not be explained by a changing pattern of spots. We also note that in our surface images the changes in the spot latitude were less than 10 degrees. The differential rotation would have to be strong in order to cause the observed changes in the photometric rotation period.

There are certainly other alternative explanations for the flat bottomed line profile shape than antisolar differential rotation. The very rapid rotation of FK Comae-type stars, perhaps being recently coalesced binaries, suggests there could be significant deviations from a spherical shape (cf. Weltey & Ramsey nonsph (1994)). Our understanding of absorption line formation in the photosphere of these exceptional stars may also be insufficient. In particular, the strong absorption lines may be subject to non-LTE effects. On the other hand, Bruls et al. (bruls (1998)) showed that chromospheric activity and non-LTE effects could only explain the box-like shape of some (e.g. Fe i $`\lambda `$ 6430 Å), but not all (e.g. Fe i $`\lambda `$ 6411 Å and Ca i $`\lambda `$ 6439 Å) of the lines used in our analysis.

Weak antisolar differential rotation has previously been detected in some RS CVn stars (cf. Vogt et al. hr1099 (1999)). These detections were made by comparing Doppler images having a time difference of up to several months. This approach is, however, only possible for stars with a very small differential rotation, since the largest sensible phase difference is 0.5. Secondly, it is hard to confirm that the same particular spots are really present in different images taken months apart.

Strassmeier et al. (klasu2 (1999)) cross-correlated constant latitude slices of images to determine the differential rotation of HD~199178. They could not see any effects of differential rotation but stated that this “non-detection does not necessarily mean there is no differential surface rotation”. We stress that the cross-correlation will provide reliable results only if there are long-lived spots on different latitudes. That is not the case in our images of HD~199178. Furthermore, the result may be biased by artifacts in the images.

Differential rotation can also be detected by combining Doppler images and time series analysis of photometry. If there is a clear correlation between the latitude of the main spot and the photometric period, one could determine the law of the differential rotation. But this kind of analysis would require more data than was available for the present study.

###### Acknowledgements.

The work of T.H. was supported by a grant from the Jenny and Antti Wihuri foundation, and by Helsinki University research funding for the project “Time series analysis in astronomy” (No. 974/62/98). I. Ilyin reduced the July 1994 spectroscopic data, participated in the observations and gave useful comments on the manuscript. The molecular line data were kindly supplied by S. Berdyugina. We also thank the anonymous referee, whose suggestions helped to improve the paper. The surface imaging temperature maps were calculated using the Cray C94/128 supercomputer at the Centre for Scientific Computing, Espoo, Finland. This research has made use of the Simbad-database operated at CDS, Strasbourg, France.


text:
# Sub-milliarcsec-scale structure of the gravitational lens B1600+434

## 1 Introduction

The gravitational lens B1600+434 was discovered by Jackson et al. (1995). The radio source is a quasar at $`z=1.59`$ and has two images separated by 1.4 arcsec, which are denoted A and B. The main lensing galaxy is identified as an edge-on spiral at $`z=0.41`$, with an additional lensing contribution from a nearby companion galaxy (Jaunsen & Hjorth 1997, Fassnacht & Cohen 1998, Koopmans, de Bruyn & Jackson 1998). The background quasar is variable in both optical and radio bands, thus allowing time delay measurements from monitoring observations (Jaunsen & Hjorth 1997, Koopmans et al. 2000). A time delay of $`47_9^{+12}`$ days has been determined from the radio monitoring (Koopmans et al. 2000), and a delay of $`51\pm 4`$ days from the optical observations (Burud et al. 2000).

The radio observations of Koopmans et al. (2000) showed that the brighter image, A, varied more rapidly than the weaker image, B. The rapid variability in A was not found to be correlated with corresponding rapid variability in B, although differential variability between A and B on longer time scales allowed the measurement of the time delay mentioned above. Correlated variability is one of the fundamental characteristics of gravitational lensing. Additional variability in A can only be explained by invoking propagation effects which affect image A only (Koopmans & de Bruyn 2000; hereinafter KdB). Such effects include a combination of scattering in the lensing galaxy and scintillation in our own galaxy, or microlensing in the halo of the lensing galaxy. KdB considered the possibility that B might be scatter-broadened in the lens and thereby would scintillate less than A. They argued however that the time scale and the frequency dependence of the variability, taken in conjunction with the high scattering measure required in the lens galaxy, did not favour scintillation as a possible cause. Instead, they proposed microlensing of image A in the halo of the lensing galaxy as a preferred alternative if the radio source has a superluminal component.

One of the ways to test these scenarios is to image the structure in the lensed images at high spatial resolution, and to apply the constraints imposed by gravitational lensing. In particular, the surface brightness of gravitationally lensed images must be equal and any significant departure would suggest that one of the images has been affected by scattering more than the other. Since the microlensing hypothesis requires the background quasar be a superluminal source, the images, therefore, might be expected to show evidence of this motion in their morphology, such as a core-jet structure typical of such quasars. In addition, the polarization properties of the images can differ if one of the ray paths is affected either by scattering or microlensing.

With the above aims in mind, we observed B1600+434 at 15 GHz using the VLBA with a resolution of 0.5 milliarcsec (mas), to detemine the fine-scale structure and to measure image sizes for A and B. The observations are described in Section 2 and the results are discussed in Section 3.

## 2 Observations and data analysis

We observed B1600+434 on 2000 September 9/10 (1900 UT to 0500 UT) at 15 GHz using all ten telescopes of the VLBA. The data were recorded in eight 8-MHz channels in dual circular polarization, yielding a total recorded bandwidth of 32-MHz per polarization. Two-bit sampling quantization was used. A single correlation centre was selected at $`\alpha =16^h01^m40.47^s,\delta =43^d16^m47.2^s`$ (J2000). Since the image separation of 1.4 arcsec is large compared with the resolution of 0.5 mas, it is necessary to sample the data finely both in frequency and time to avoid smearing of the visibility function and the associated image distortions in the field of view. The data were therefore sampled at 0.5 MHz in frequency and 0.5 sec in time.

We observed 3C345 (J1640+3946), 1611+343 (J1613+3412) and 1749+096 (J1751+0939) as calibrators. The total integration time on the target source B1600+434 was 195 min. We used 1749+096 to calibrate the instrumental polarization, as well as to determine the bandpass function.

The method of data analysis follows closely that described by Patnaik et al. (1999) as used for the gravitational lens B1422+231. The data were analysed using the software package AIPS, maintained by the National Radio Astronomy Observatory (NRAO). After correcting for the parallactic angle variation of the telescopes, amplitudes were calibrated using the system temperatures and telescope gain factors as a function of elevation. The phase slopes and offsets across each 8-MHz band were aligned using the pulse-cal information from each telescope. Standard fringe-fitting was performed on the calibrators. After careful editing, the bandpass function was determined using the source 1749+096. We followed the usual procedure for calibrating instrumental polarization, using 1749+096 as the polarization calibrator. Typical instrumental polarization amplitude was measured to be $``$1.6%. Since B1600+434 has two images with similar flux densities at a large angular separation, we had to iterate the fringe-fitting process several times, each time using the model of the source determined from the previous solution.

The final map was made with two sub-fields, centred on each lensed image. We also made polarization maps of the source in Stokes $`I,Q`$ and $`U`$.

## 3 Results and Discussion

The maps of the lensed images, made using uniform weighting, are shown in Fig.1. The maps have been restored using a circular gaussian with a FWHM of 0.5 mas. We did not detect any polarized emission from either component A or B, yielding a 5$`\sigma `$ limit to the percentage linear polarization of 2.6% for image A, and 3.3% for image B. The rms noise in the polarization map is 0.09 mJy/beam.

Table 1 lists the results of fitting a single gaussian component to each of the images. Two clear results emerge from these observations. First, both the images are resolved, and A is more compact than B. Secondly, the structure is elongated on either side of the peak in both of the images. In order to quantify the errors of the size measurements we divided the data into two sets by removing alternate scans so as to preserve similar uv-coverage. The resulting maps were then fitted with a single gaussian. The results show that the errors for image A are $`0.01`$ mas in each axis. However, the errors are larger for image B by about 50%. Note that the error estimates derived in this way fold in most systematic calibration or imaging errors which may affect the component size measurement, excluding deviation of the component from non-Gaussian shape. A weak extension is seen to the south of image B which appears to be genuine as it is present in both the datasets mentioned earlier. We are not entirely certain if it is real as there is no counterpart in image A. This can be confirmed by another independent dataset. The measured size of image B, however, is not affected by including this feature.

The implication of the first result is that the surface brightnesses of the two images differ significantly. The flux density ratio of A/B is 1.27$`\pm `$0.05 whereas the ratio of image sizes is 0.80$`\pm `$0.15. Secondly, the image morphology as measured is in apparent conflict with a typical superluminal radio source in the sense that such radio sources are core-dominated and must have one-sided jets due to relativistic motion Doppler boosting of the knot components.

In order to understand the higher modulation index (fractional rms variability) of image A (2.8% at 8.4 GHz), KdB considered the possibility that the size of B (modulation index of 1.6% at 8.4 GHz) might be larger through scatter-broadening by the ionised interstellar medium in the lens galaxy. The fact that the surface brightnesses do not agree, and in particular, that the weaker image B is more extended, suggests that B has been affected by propagation through the lensing galaxy more than A. This would support the idea that B is somewhat broadened by scattering in the lensing galaxy. We note that its ray path passes close to the bulge of this edge-on spiral galaxy whereas the ray path of A goes through the halo of the lens. For the two images to scintillate in our Galaxy, the required sizes are 62 $`\mu `$arcsec and 108 $`\mu `$arcsec for A and B respectively (KdB). Our observations cannot resolve such sizes and our measured sizes are considerably larger than these. We note, however, that the observed images have been parametrised as single gaussian components. This does not exclude any sub-structure in the source. Therefore, our result that image B is larger in size than image A is in general agreement with the expectation from scintillation taking place in our Galaxy.

KdB measured a modulation index $``$1.2% at 1.4 GHz and 3.7% at 5 GHz. If the variability was due to scintillation, then the modulation index should have increased going to lower frequencies. However, if there is frequency-dependent structure then the modulation index can be smaller at lower frequencies. A flat spectrum source likely consists of more than one self-absorbed component (Cotton et al. 1980). It is then possible that the component which is compact and scintillates at 5.0 GHz, which may for example be the core, does not contribute much to the flux density at 1.4 GHz due to self-absorption. Some other component in the jet might be brighter at 1.4 GHz but this component would be more extended and therefore unlikely to scintillate. Therefore, a smaller modulation index at a lower frequency does not necessarily rule out scintillation occurring at a higher frequency. Such frequency-dependent structure has been noted in the gravitational lens B0218+357 (Patnaik & Porcas 1999).

The alternative explanation proposed by KdB is that A is microlensed in the halo of the lensing galaxy. However, for the microlensing to match the time-scale of variability, the model used by KdB requires that the background source be a superluminal source with a knot the size of a few $`\mu `$arcsec containing $``$ 5-11% of the flux density, moving at an apparent speed in the range $``$ 9-26c. The morphology of a superluminal source will of necessity be that of a bright core and one-sided jet, due to Doppler boosting of a relativistically moving knots in the jet. The detection of extensions on either side of the peak in both the lensed images here is contrary to the above expectation. It would, therefore, appear that superluminal motion in the background source is not favoured and hence the microlensing caustic crossing time will be much longer than the observed variability. However, it is also possible that the peak we observe is not necessarily the core; it might be a bright knot and thus still allow the possibility that the source could be intrinsically one-sided. We do not have sufficient resolution to distinguish this level of detailed structure. If we assume that the peak corresponds to the core, which may be a reasonable assumption given the fact that these observations are at 15 GHz and with 0.5 mas resolution where cores usually brighter than knots, then it would impliy that microlensing combined with superluminal motion in the background source may not be a viable explanation for the rapid variability of the lensed image A.

As an example, the two brightest images in the lens system B1422+231 have elongation produced by lensing (Patnaik & Porcas 1998, Patnaik et al. 1999). In this lens system the background radio source is also a flat spectrum object and at milliarcsec-scales does not have the standard core-jet structure. It has been suggested that the magnification in this source might be very high and that the background object could possibly be a radio-quiet object. A similar scenario is possible for B1600+434. Certainly, in both B1422+231 and B1600+434, no jet has been detected and both sources show two-sided emission. If the extension of the images we detect is due to stretching of the extended background emission by gravitational lensing, then we cannot be sure of the direction of the jet. However, should there be superluminal motion (which can be as much as 0.15 milliarcsec per year for the bright quasars), we would detect changes in image size after a year or two.

## 4 Conclusions

We have observed the gravitational lens system B1600+434 at 15 GHz using the VLBA at a resolution of 0.5 milliarcsec. We find that the two lensed images are resolved and have elongation on either side of the peak, and that the surface brightnesses of the images are significantly different. The image shapes do not appear to favour superluminal motion in the background radio source and hence the interpretation that A is microlensed is not favoured by these observations at this epoch. Our results are consistent with the view that image B is perhaps scatter-broadened by the interstellar medium of the lens galaxy, and therefore scintillates less than image A in our Galaxy. However, VLBI observations of the source at a range of other observing frequencies are required to provide further insight into this question.

###### Acknowledgements.

We thank Peter Schneider for helpful discussion and J. Schmid-Burgk and Karl Menten for a careful reading of the manuscript. We would like to thank the referee A.G. de Bruyn for critical comments and helpful advice. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.


text:
# 1 Introduction

## 1 Introduction

Multi-jet production in deep-inelastic scattering (DIS) has been successfully used at HERA to test the predictions of perturbative QCD (pQCD) over a large range of four-momentum transfer squared $`Q^2`$ . Recently the H1 collaboration has determined the strong coupling constant $`\alpha _s`$ and the gluon density in the proton  from the inclusive jet and the dijet cross sections measured in the Breit frame. While these cross sections are directly sensitive to QCD effects of order $`𝒪(\alpha _s)`$, the three-jet cross section in DIS is already proportional to $`\alpha _s^2`$ in leading order in pQCD. The higher sensitivity to $`\alpha _s`$ and the greater number of degrees of freedom of the three-jet final state thus allow the QCD predictions to be tested in more detail in three-jet production. In this paper we present for the first time differential measurements of the three-jet cross section in neutral current DIS and measurements of shapes of angular jet distributions which are sensitive to dynamic effects of the interaction. Similar studies of three-jet production in reactions with initial state hadrons have been carried out previously in hadron-hadron collisions at the SPS , the ISR  and at the TEVATRON  as well as in photoproduction at HERA . The present analysis includes the first comparison of three-jet distributions measured in hadron induced reactions with a perturbative QCD calculation in next-to-leading order $`\alpha _s`$ .

In neutral current DIS the lepton interacts with a parton in the proton via the exchange of a boson ($`\gamma `$, $`Z`$). At a fixed center-of-mass energy the kinematics of the lepton inclusive reaction (for unpolarized lepton and proton beams) is given by two variables which are here chosen to be the four-momentum transfer squared $`Q^2`$ and the Bjorken scaling variable $`x_{\mathrm{Bj}}`$. The subprocess $`1+23+4+5`$ in which three massless jets emerge from the boson-parton reaction is fully described by six further variables which can be constructed from the energies $`E_i`$ and the momenta $`\stackrel{}{p}_i`$ of the jets in the three-jet center-of-mass (CM) frame. It is convenient to label the three jets ($`i=3,4,5`$) in the order of decreasing energies in the three-jet CM frame. These variables are conventionally chosen  to be the azimuthal orientation of the three-jet system, the invariant mass of the three-jet system $`M_{3\mathrm{j}\mathrm{e}\mathrm{t}}`$, the jet energy fractions111Note that from energy and momentum conservation $`X_3+X_4+X_5=2`$ with $`X_5=(2E_5)/M_{3\mathrm{j}\mathrm{e}\mathrm{t}}`$. $`X_3`$, $`X_4`$

$$X_3\frac{2E_3}{M_{3\mathrm{j}\mathrm{e}\mathrm{t}}},X_4\frac{2E_4}{M_{3\mathrm{j}\mathrm{e}\mathrm{t}}},$$
(1)
and two angles $`\theta _3`$ and $`\psi _3`$ that specify the relative orientation of the jets,

$$\mathrm{cos}\theta _3\frac{\stackrel{}{p}_B\stackrel{}{p}_3}{|\stackrel{}{p}_B||\stackrel{}{p}_3|},\mathrm{cos}\psi _3\frac{(\stackrel{}{p}_3\times \stackrel{}{p}_B)(\stackrel{}{p}_4\times \stackrel{}{p}_5)}{|\stackrel{}{p}_3\times \stackrel{}{p}_B||\stackrel{}{p}_4\times \stackrel{}{p}_5|},$$
(2)
where $`\stackrel{}{p}_B`$ denotes the direction of the proton beam. As indicated in Fig. 1, $`\theta _3`$ is the angle of the highest energetic jet with respect to the proton beam direction and $`\psi _3`$ is the angle between the plane spanned by the highest energy jet and the proton beam and the plane containing the three jets. The angle $`\psi _3`$ indicates whether the third jet (i.e. the lowest energy jet) is radiated within ($`\psi _30`$ or $`\psi _3\pi `$) or up to perpendicular to ($`\psi _3\pi /2`$) the plane containing the highest energy jet and the proton beam. In dijet production, in the dijet center-of-mass frame both jets carry half of the available energy and are scattered back-to-back. The presence of a third jet, however, allows the $`\mathrm{cos}\theta _3`$ distribution to be asymmetric in the three-jet CM frame and the energies of the two leading jets to be smaller than half of the total available energy ($`X_3,X_4<1`$).

The observable $`R_{3/2}`$, defined by the ratio of the inclusive three-jet cross section and the inclusive two-jet cross section, is of interest especially for quantitative studies, since in this ratio both experimental and some theoretical uncertainties cancel to a large extent.

This paper presents measurements of the inclusive three-jet cross section in DIS as a function of $`Q^2`$, $`x_{\mathrm{Bj}}`$ and $`M_{3\mathrm{j}\mathrm{e}\mathrm{t}}`$. Distributions of three-jet events are measured in the variables $`X_3`$, $`X_4`$, $`\mathrm{cos}\theta _3`$ and $`\psi _3`$ and are normalized to the integrated three-jet cross section. The ratio $`R_{3/2}`$ is measured as a function of $`Q^2`$. The kinematic range of the analysis covers four-momentum transfers squared $`Q^2`$ between $`5\mathrm{GeV}^2`$ and $`5000\mathrm{GeV}^2`$ and invariant three-jet masses $`M_{3\mathrm{j}\mathrm{e}\mathrm{t}}`$ in the range from $`25\mathrm{GeV}`$ to $`140\mathrm{GeV}`$.

## 2 Event Selection

The analysis is based on data taken in positron-proton collisions with the H1 detector at HERA in the years 1995–1997 with a positron beam energy of $`E_e=27.5\mathrm{GeV}`$ and a proton beam energy of $`E_p=820\mathrm{GeV}`$, leading to a center-of-mass energy $`\sqrt{s}`$ of $`300\mathrm{GeV}`$. A detailed description of the H1 detector can be found in . The main detector components relevant for the present analysis are the liquid argon (LAr) calorimeter , the backward lead-fiber calorimeter (SpaCal)  and the tracking chamber system. In the polar angular222The polar angle $`\theta `$ is defined with respect to the positive $`z`$-axis which is given by the proton beam direction in all reference frames. range $`4^{}<\theta <154^{}`$ ($`153^{}<\theta <177^{}`$) the electromagnetic and hadronic energies are measured by the LAr calorimeter (SpaCal) with full azimuthal coverage. Charged particle tracks are measured in two concentric drift chamber modules ($`25^{}<\theta <165^{}`$) and by a forward tracking detector ($`7^{}<\theta <25^{}`$).

The experimental procedure is similar to the one used in a previous measurement of the inclusive jet and the dijet cross section . Here we briefly summarize only the salient features. The identification and triggering of neutral current DIS events is based on the reconstruction of the event vertex and the detection of the scattered positron as a compact electromagnetic cluster in either the SpaCal (“low $`Q^2`$”) or the LAr calorimeter (“high $`Q^2`$”). The two event samples correspond to integrated luminosities of $`_{\mathrm{int}}=21.1\mathrm{pb}^1`$ (low $`Q^2`$) and $`_{\mathrm{int}}=32.9\mathrm{pb}^1`$ (high $`Q^2`$), respectively. The trigger efficiencies for the final jet event samples are above 98%. The hadronic final state is reconstructed from a combination of tracks with low transverse momentum ($`p_T<2\mathrm{GeV}`$) and energy deposits in the LAr calorimeter and the SpaCal. The kinematic variables $`Q^2`$, $`x_{Bj}`$ and $`y=Q^2/(sx_{\mathrm{Bj}})`$ are reconstructed using the $`e\mathrm{\Sigma }`$-method . The kinematic range of the analysis is specified by

low $`Q^2`$: $`5<Q^2<100\mathrm{GeV}^2\text{and}\text{ }\theta _{\mathrm{positron}}>156^{},`$
high $`Q^2`$: $`150<Q^2<\mathrm{\hspace{0.17em}5000}\mathrm{GeV}^2,`$ (3)
low $`Q^2`$ and high $`Q^2`$: $`0.2<y<\mathrm{\hspace{0.33em}\hspace{0.33em}0.6}.`$
The jet selection is carried out in the Breit frame in which $`\stackrel{}{q}+2x_{\mathrm{Bj}}\stackrel{}{P}=0`$, where $`\stackrel{}{q}`$ and $`\stackrel{}{P}`$ are the momenta of the exchanged boson and of the incoming proton, respectively. Jets are defined by the inclusive $`k_{}`$ clustering algorithm  which is applied to the final state particles, excluding the scattered positron. The parameter $`R_0`$, which defines the minimal separation of jets in pseudorapidity333The pseudorapidity $`\eta `$ is related to the polar angle $`\theta `$ by $`\eta =\mathrm{ln}(\mathrm{tan}\theta /2)`$. and azimuth space, is set to $`R_0=1`$. The clustering of particles is done in the $`E_T`$ recombination scheme in which the resulting jets are massless. A detailed description of the procedure is given in . The jet phase space is defined by cuts on the jet pseudorapidity $`\eta _{\mathrm{jet},\mathrm{lab}}`$ in the laboratory frame and on the transverse jet energy $`E_T`$ in the Breit frame

$$1<\eta _{\mathrm{jet},\mathrm{lab}}<2.5\mathrm{and}E_T>5\mathrm{GeV}.$$
(4)
The inclusive three-jet (dijet) sample consists of all events with three (two) or more jets from which the three (two) jets of highest $`E_T`$ have an invariant mass of

$$M_{3\mathrm{j}\mathrm{e}\mathrm{t}}>\mathrm{\hspace{0.17em}25}\mathrm{GeV}(M_{2\mathrm{j}\mathrm{e}\mathrm{t}}>\mathrm{\hspace{0.17em}25}\mathrm{GeV}).$$
(5)
With these selection criteria applied, the inclusive three-jet event sample at low $`Q^2`$ (high $`Q^2`$) contains $`2903`$ ($`666`$) events, and the corresponding inclusive dijet sample contains $`6746`$ ($`2005`$) events. The jet selection cuts in (4) deplete the phase space regions of $`\mathrm{cos}\theta _3\pm 1`$, $`\psi _30`$ and $`\psi _3\pi `$ . To reduce this influence the following additional cuts are applied for the measurement of the differential distributions of the variables $`X_3`$, $`X_4`$, $`\mathrm{cos}\theta _3`$ and $`\psi _3`$

$`X_3`$ $`<`$ $`0.95,`$
$`|\mathrm{cos}\theta _3|`$ $`<`$ $`0.8,`$ (6)
$`M_{3\mathrm{j}\mathrm{e}\mathrm{t}}`$ $`>`$ $`40\mathrm{GeV}\text{only for }Q^2<100\mathrm{GeV}^2.`$
Due to the limited size of the event sample the last cut is not applied in the high $`Q^2`$ analysis. The cuts in (2) are passed by $`523`$ ($`536`$) three-jet events in the low $`Q^2`$ (high $`Q^2`$) event sample. The size of the photoproduction background has been estimated using two samples of photoproduction events generated by PYTHIA  and PHOJET  and is found to be negligible (i.e. below 2%) in all distributions measured.

## 3 Correction of the Data

The data are corrected for effects of detector resolution and acceptance, as well as for inefficiencies of the selection and for higher order QED effects. The correction functions are determined using event samples generated by the Monte Carlo event generators LEPTO , RAPGAP  and ARIADNE , all interfaced to HERACLES  to take QED corrections into account. For each generator two event samples are generated. The first sample, which includes QED corrections, is subjected to a detailed simulation of the H1 detector based on GEANT . The second event sample is generated under the same physics assumptions, but without QED corrections and without detector simulation. The correction functions are determined bin-wise for each observable as the ratio of its value in the second sample and its value in the first sample.

The applicability of the simulated event samples for the correction procedure has been tested by comparing a variety of their distributions to the data for the three-jet and the dijet event samples . The simulated events are found to give a good description of single jet related quantities such as the surrounding energy flow, their angular distributions and their internal structure. Based on the event simulation the bin widths of all variables have been adjusted to match their resolution. At low $`Q^2`$ (high $`Q^2`$) the correction functions are determined using event samples from RAPGAP and ARIADNE (LEPTO and ARIADNE). In both $`Q^2`$ regions the correction functions agree within typically 10% for the measured cross sections and for the ratio $`R_{3/2}`$ and within typically 5% for the normalized distributions of $`X_3`$, $`X_4`$, $`\mathrm{cos}\theta _3`$ and $`\psi _3`$. The final correction functions are taken to be the average of the two models and half of their difference is used as an estimate of the model dependence. The list of further experimental uncertainties considered in the analysis is identical to the one in the recent measurement of inclusive jet and dijet cross sections . The uncertainty of the measured cross sections are typically 16% and are dominated by systematic effects, the largest contributions being due to the hadronic energy scale of the LAr calorimeter and the model dependence of the correction functions. The uncertainties of the measured rate $`R_{3/2}`$ have equal statistical and systematic contributions while in the normalized distributions the statistical uncertainties dominate.

## 4 Results

The measured three-jet distributions are presented in Figs. 2–5 and in table 1. The inner error bars represent the statistical uncertainties and the outer error bars the quadratic sum of all uncertainties. The data are compared to the pQCD predictions in leading order (LO) and in next-to-leading order (NLO) with and without hadronization corrections. The LO and NLO calculations are carried out in the $`\overline{\mathrm{MS}}`$-scheme for five massless quark flavors using the recent program NLOJET . It was checked that the LO calculations from NLOJET agree with those of MEPJET  and DISENT . Heavy quark mass effects are estimated using the LO calculation MEPJET and are found to lower the three-jet cross section by typically $`5\%`$ at low $`Q^2`$ and $`3\%`$ at high $`Q^2`$. Renormalization and factorization scales ($`\mu _r`$, $`\mu _f`$) are set to the average transverse energy $`\overline{E}_T`$ of the three jets in the Breit frame. The parton density functions of the proton are taken from the parameterization CTEQ5M1 . The strong coupling constant is set to the world average value of $`\alpha _s(M_Z)=0.118`$  and is evolved according to the two-loop solution of the renormalization group equation. Hadronization corrections $`\delta _{\mathrm{hadr}}`$ are determined using LEPTO as the relative change of an observable before and after hadronization. These corrections are in the range $`22\%<\delta _{\mathrm{hadr}}<18\%`$ for the three-jet cross section and $`18\%<\delta _{\mathrm{hadr}}<10\%`$ for the ratio $`R_{3/2}`$ over the whole range of $`Q^2`$ and they agree with those obtained from HERWIG  to within $`2\%`$. Hadronization corrections are negligible for the normalized distributions (Figs. 4 and 5) since they basically change the size of the three-jet cross section but not the shapes of differential distributions.

The three-jet cross section is presented in Fig. 2 (a) as a function of the four-momentum transfer squared $`Q^2`$. The data are compared to the LO and NLO predictions, the latter with and without hadronization corrections. In the lower part of the figure the ratio of the measured cross section and the NLO prediction (corrected for hadronization effects) is shown. Over the whole range of $`Q^2`$ the NLO prediction (corrected for hadronization effects) gives a good description of the data — at high $`Q^2`$, where NLO corrections are small, but also at low $`Q^2`$ where the NLO prediction is a factor of two above the LO prediction. The theoretical prediction is subject to several uncertainties, the dominant sources being the value of the strong coupling constant, the parton density functions of the proton (especially the gluon density) and the renormalization scale dependence of the NLO calculation. In the lower part of Fig. 2 (a) the size of these uncertainties is displayed by three different bands, indicating variations of the strong coupling constant (by $`\mathrm{\Delta }\alpha _s(M_Z)=\pm 0.006`$), the gluon density444The variations of $`\alpha _s(M_Z)`$ and the gluon density in the proton correspond roughly to the uncertainties within which both have been determined in a previous analysis  from jet cross sections in DIS . (overall by $`\mathrm{\Delta }g(x,\mu _f)=\pm 15\%`$) and the renormalization scale (by $`0.5<(\mu _r/\overline{E}_T)<2`$). While at $`Q^250\mathrm{GeV}^2`$ the variation of $`\alpha _s`$ gives the largest effect, the renormalization scale dependence is the dominant source of uncertainty at lower values of $`Q^2`$, i.e. in the region where NLO corrections are also large. Over the whole $`Q^2`$ range the change of the cross section, induced by the variation of the gluon density is approximately half as large as the change induced by the $`\alpha _s`$ variation.

For the ratio $`R_{3/2}`$ of the inclusive three-jet cross section and the inclusive dijet cross section which is shown in Fig. 2 (b) some experimental and theoretical uncertainties cancel. The data are measured as a function of $`Q^2`$ and compared to LO and NLO calculations which are corrected for hadronization effects. While the LO calculation predicts a stronger $`Q^2`$ dependence of $`R_{3/2}`$ than observed in the data, the NLO calculation gives a good description of the data over the whole range of $`Q^2`$. Uncertainties in the theoretical prediction for the ratio $`R_{3/2}(Q^2)`$ are investigated in the same way as for the three-jet cross section. The NLO corrections for the three-jet and the dijet cross sections are of similar size. At low $`Q^2`$ this leads to a smaller NLO correction and renormalization scale dependence for the ratio $`R_{3/2}`$ than for the cross section. Furthermore, when measured in the same region of $`x_{\mathrm{Bj}}`$ and $`Q^2`$, with the same cut on the invariant multi-jet mass, the three-jet and the dijet cross sections probe the parton density functions of the proton in the same range of proton momentum fractions $`\xi =x_{\mathrm{Bj}}(1+M_{\mathrm{n}\mathrm{jet}}^2/Q^2)`$. Since both jet cross sections are dominated by gluon induced processes, the ratio $`R_{3/2}`$ is almost insensitive to variations of the gluon density in the proton. It is recognized that the ratio $`R_{3/2}`$ is experimentally measured and theoretically calculated with small uncertainties over the entire range of $`Q^2`$. For the central value of $`\alpha _s(M_Z)0.118`$ used in the calculations, the theoretical predictions are consistent with the data for the three-jet cross section and the ratio $`R_{3/2}`$.

The three-jet cross section is displayed in Fig. 3 as a function of the Bjorken scaling variable $`x_{\mathrm{Bj}}`$ and the invariant three-jet mass $`M_{3\mathrm{j}\mathrm{e}\mathrm{t}}`$ over the range of $`10^4<x_{\mathrm{Bj}}<0.2`$ and $`25<M_{3\mathrm{j}\mathrm{e}\mathrm{t}}<140\mathrm{GeV}`$ in two regions of $`Q^2`$. In the low $`Q^2`$ region the LO calculation underestimates the cross section at small $`x_{\mathrm{Bj}}`$ and small $`M_{3\mathrm{j}\mathrm{e}\mathrm{t}}`$ by a factor of two. The NLO calculation, however, gives a good description of the data.

In Fig. 4 the distributions of the energy fractions $`X_3`$ and $`X_4`$ of the two higher energetic jets are shown. They are normalized to the total three-jet cross section in the same kinematic range. The data are well described by the QCD prediction in NLO. For these normalized distributions NLO corrections are negligible and the LO calculation (not shown) can not be distinguished from the NLO curves. In addition the figures also include a three-jet distribution generated with a uniform population of the available three-body phase space (labeled “Phase Space”)555The phase space distribution of the three jets is generated according to the observed invariant mass spectrum using the jet separation and selection criteria as applied in the data. To apply the pseudorapidity cuts on the jets in the laboratory frame also the distributions of the variables $`Q^2`$ and $`x_{\mathrm{Bj}}`$ (which define the boost vector) are fitted to the data.. The phase space prediction is, however, similar to the NLO prediction, except for a small shift towards larger values of $`X_3`$ and $`X_4`$, indicating the Bremsstrahlung nature of the process.

The normalized distributions of the angular variables $`\mathrm{cos}\theta _3`$ and $`\psi _3`$ are shown in Fig. 5 for the data at low $`Q^2`$ (left) and at high $`Q^2`$ (right). The $`\mathrm{cos}\theta _3`$ distributions in both $`Q^2`$ regions are peaked at the cut value of $`\mathrm{cos}\theta _3=\pm 0.8`$, corresponding to angles close to the proton and the photon direction. The phase space prediction shows the opposite behavior, peaking at zero and falling towards $`\mathrm{cos}\theta _3\pm 0.8`$. The QCD calculation in NLO shows an asymmetry around zero and gives a reasonable description of the effects seen in the data. The latter is also the case for the measured $`\psi _3`$ distributions which are relatively flat, while the underlying phase space is strongly peaked at $`\psi _3=\pi /2`$. Although at low $`Q^2`$ the NLO prediction is almost a factor of two higher than the LO prediction, the shapes of the angular jet distributions are almost unaffected by the NLO correction. Differences between the QCD calculation and the phase space prediction indicate that due to the Bremsstrahlung nature of the process configurations are preferred in which the plane containing the two less energetic jets coincides with the plane spanned by the proton beam and the highest energetic jet, corresponding to values of $`\psi _30`$ and $`\psi _3\pi `$. These results are in qualitative agreement with measurements in $`\overline{p}p`$ collisions  and in $`\gamma p`$ collisions . Differences in the observed shapes of the two variables at low and at high $`Q^2`$ are, to some extent, due to the phase space resulting from the different cuts in $`M_{3\mathrm{j}\mathrm{e}\mathrm{t}}`$.

## 5 Summary

A measurement of the three-jet cross section in deep-inelastic scattering has been presented. At a positron-proton center-of-mass energy $`\sqrt{s}`$ of $`300\mathrm{GeV}`$ the production rates and angular distributions of three-jet events, selected in the Breit frame with $`E_T>5\mathrm{GeV}`$, have been studied over a large range of four-momentum transfer squared, $`5<Q^2<\mathrm{5\hspace{0.17em}000}\mathrm{GeV}^2`$.

The inclusive three-jet cross section has been measured as a function of $`Q^2`$, Bjorken-$`x`$ and the invariant three-jet mass for invariant three-jet masses above $`25\mathrm{GeV}`$. The ratio $`R_{3/2}`$ of the inclusive three-jet and the inclusive dijet cross section has been measured as a function of $`Q^2`$. The predictions of perturbative QCD in next-to-leading order give a good description of the three-jet cross section and the ratio $`R_{3/2}`$ over the whole range of $`Q^2`$ for values of the strong coupling constant close to the current world average of $`\alpha _s(M_Z)0.118`$. Angular jet distributions and jet energy fractions have been measured in the three-jet center-of-mass frame. The angular orientation of the three-jet system follows the radiation pattern expected from perturbative QCD. While the angular distributions are not consistent with a uniform population of the available phase space, they are reasonably well described by the QCD predictions.

## Acknowledgments

We thank Erwin Mirkes, Zoltan Nagy, Mike H. Seymour and Dieter Zeppenfeld for many helpful discussions. We are grateful to the HERA machine group whose outstanding efforts have made and continue to make this experiment possible. We thank the engineers and technicians for their work in constructing and now maintaining the H1 detector, our funding agencies for financial support, the DESY technical staff for continual assistance and the DESY directorate for the hospitality which they extend to the non-DESY members of the collaboration.


text:
# Constraining the window on sterile neutrinos as warm dark matter

## 1 Introduction

Astrophysics provides an increasing amount of independent indications that the dark matter of the universe is warm, so that the small-scale fluctuations are damped out by free streaming. This is most easily achieved by giving a keV mass to the DM particle, in which case the preferred candidate is the sterile neutrino. Support for warm dark matter (WDM) comes from simulations of the number of satellite galaxies (Colin, Avila-Reese & Valenzuela 2000) and of disk galaxy formation without the need for stellar feedback (Sommer-Larsen & Dolgov 2000), which both find that a DM particle mass of about 1 keV is optimal: a significantly larger mass has little impact on galaxy formation, and a significantly smaller mass would lead to the well known difficulties faced by hot dark matter. A quantitative lower limit on the candidate WDM particle mass is inferred from the existence of a massive black hole at large redshift (Barkana, Haiman & Ostriker 2001) and the requirement of sufficiently early galaxy formation to account for reionization of the universe and the observed Ly-$`\alpha `$ forest properties (Narayanan et al. 2000), constraining the DM mass to be larger than 0.75 keV. A recent discussion of x-ray emission from decays of sterile neutrinos (Abazajian, Fuller & Tucker 2001) has imposed an upper limit of about 5 keV on the neutrino mass. Here we discuss a reinterpretation of these bounds on neutrino mass, and demonstrate that the proper inclusion of the neutrino momentum, arising from the specific production temperature, reduces the allowed sterile neutrino WDM mass to be in the range $`2.6`$ keV $`<m<5`$ keV.

## 2 WDM particle decoupling

All of these studies (Colin et al. 2000; Sommer-Larsen & Dolgov 2000; Barkana et al. 2001; Narayanan et al. 2000) are based on the mass-dependent cut-off on small scales, produced by free-streaming. In the previously cited studies, a “conventional” WDM model was considered for the underlying particle physics (see Bode, Ostriker & Turok (2001) and Sommer-Larsen & Dolgov (2000) for recent overviews of such particle models). In such conventional WDM models, the particles decouple in the early Universe at higher temperatures than do massless neutrinos. Therefore they do not share the entropy release from the successive particle annihilations. Since they were relativistic at decoupling, their distribution function in momentum space is subsequently that of a massless fermion, but with a temperature, $`T_W`$, which is given today by

$$T_{W_0}=T_{\nu _0}\left(\mathrm{\Omega }_Wh^2\frac{94\text{eV}}{m_W}\right)^{1/3},$$
(1)
where $`T_{\nu _0}1.946`$ K, $`H_0=100h\text{km}\text{s}^1\text{Mpc}^1`$, $`\mathrm{\Omega }_W`$ is the present energy density of WDM in units of the critical density, and $`m_W`$ is the WDM mass. The needed entropy release in these conventional models is much bigger than allowed in the standard model, and such WDM candidates should therefore have decoupled before a larger gauge group breaks down.

Now a very natural candidate for the WDM particle is a massive sterile neutrino mixed with an ordinary neutrino (Dodelson & Widrow 1994; Colombi, Dodelson & Widrow 1996; Shi & Fuller 1999; Dolgov & Hansen 2001; Abazajian, Fuller & Patel 2001). Since the mixing angle is temperature dependent (Nötzold & Raffelt 1988) (and for small vacuum mixing angle, $`\mathrm{sin}^22\theta 10^7`$), only a small amount of these heavy neutrinos, relative to ordinary active neutrinos, can be produced at high temperatures. The distribution function of sterile neutrinos is, to a fair approximation, characterized by the temperature of massless neutrinos, but smaller by a factor $`𝒳`$. For a specific choice of $`m_W`$ and $`\mathrm{\Omega }_W`$ today the value of $`𝒳`$ can be found from

$$𝒳=\mathrm{\Omega }_Wh^2\left(\frac{94\text{eV}}{m_W}\right)10^2,$$
(2)
therefore the two models produce the same contribution to $`\rho _{tot}`$ of WDM particles today if

$$𝒳=\left(\frac{T_{W_0}}{T_{\nu _0}}\right)^3.$$
(3)
However, WDM particles with $`m_W`$ have a different distribution function in these two models, and their free-streaming effect is not equivalent. The effect on large scale structure was first discussed in detail by Colombi et al. (1996). We will use conventional WDM (cWDM) when referring to the first case, and sterile neutrino WDM (sWDM) for the second one.

## 3 Lower bounds

The two neutrino models are easily included in a Boltzman code, in order to compute the present matter power spectrum $`P(k)`$. Using the code cmbfast (Seljak & Zaldarriaga 1996), we have found analytical fits for the transfer functions, $`T(k)`$, relating the power spectrum in the WDM to the CDM scenario

$$T^2(k)=\frac{P^W(k)}{P^{CDM}(k)}\text{for}\mathrm{\Omega }_W=\mathrm{\Omega }_{CDM},$$
(4)
where $`P^W`$ is the power spectrum for cWDM, and a similar expression with $`P^\nu (k)`$ for the sWDM model. These transfer functions, which essentially reflect the free streaming cut-off, have the form

$$T(k)=\left[1+\left(\alpha k\right)^{2\nu }\right]^{5/\nu },$$
(5)
where $`k`$ is the wavenumber in units $`h\text{Mpc}^1`$, $`\nu =1.12`$, and $`\alpha `$ depends on the cosmological parameters as

$$\alpha =A\left(\frac{\mathrm{\Omega }_W}{0.3}\right)^b\left(\frac{h}{0.65}\right)^c\left(\frac{m_W}{500\text{eV}}\right)^d.$$
(6)
Numerically we find for cWDM, $`A=1.07`$, $`b=0.11`$, $`c=1.20`$ and $`d=1.11`$ in good agreement with Bode et al. (2001). For the sWDM one can derive similar numbers by noting that the mass in the cWDM case differs by $`(T_{W_0}/T_{\nu _0})`$. This means that if we have a dependence $`\mathrm{\Omega }_W^bh^cm_W^d`$ for the sWDM case in eq. (6), and a dependence $`\mathrm{\Omega }_W^b^{}h^c^{}m_W^d^{}`$ for the cWDM case, then one finds

$$\mathrm{\Omega }_W^bh^cm_W^d=\mathrm{\Omega }_W^{b^{}d^{}/3}h^{c^{}2d^{}/3}m_W^{d^{}+d^{}/3},$$
(7)
which is solved by $`b^{}=b+d/4`$, $`d^{}=3d/4`$ and $`c^{}=c+d/2`$, in good agreement with what we found numerically, by explicitly changing the massive neutrino phase-space distribution function as done by Lesgourgues & Pastor (1999). To be very explicit, this means that for a given cut-off scale of the power spectrum one can find the corresponding mass of the elementary particle, and the mass will differ in the two cases. E.g. if for cWDM one finds $`m_W=0.75`$ keV, then this corresponds in the sWDM case to $`m_W=(\mathrm{\Omega }_Wh^2\mathrm{\hspace{0.17em}94}\text{eV}/0.75\text{keV})^{1/3}0.75\text{keV}2.6`$ keV, when using $`h=0.7`$ and $`\mathrm{\Omega }_W=0.4`$. In other words, if one believes that sterile neutrinos indeed constitute the dark matter, and they are produced as described in (Dodelson & Widrow 1994; Colombi, Dodelson & Widrow 1996; Dolgov & Hansen 2001; Abazajian, Fuller & Patel 2001), then the bounds obtained by Barkana et al. (2001) and Narayanan et al. (2000), $`m_W>0.75\text{keV}`$, should really be multiplied by a factor $`3.4`$, and the lower bound on sterile neutrinos as dark matter is thus about $`2.6`$ keV.

This lower bound may be subject to minor corrections. First, the temperature of the sterile neutrinos is really slightly lower than the active neutrino temperature, since the sterile neutrinos are being produced while the muons are still present in the Universe, $`T130`$ MeV (Langacker 1989; Kainulainen 1990; Barbieri & Dolgov 1990, 1991). Similarly, for the neutrino states being produced above the QCD phase transition, one must also take into account the quark degrees of freedom (Abazajian et al. 2001a). Another effect arises from the fact that the produced sterile neutrino spectrum is not exactly thermal, but slightly warmer in the sense that the higher momentum part is more populated than the lower momentum part (Dolgov & Hansen 2001; Abazajian et al. 2001b). Furthermore, the factor of 3-4 found above depends on the specific values of $`\mathrm{\Omega }_W`$ and $`h`$, and can therefore change slightly. It is also worth noting, that if the sterile neutrinos are produced resonantly (Shi & Fuller 1999), through a pre-existing lepton asymmetry, then the upper limit on 5 keV may weaken substantially (Abazajian et al. 2001b).

## 4 Upper bounds and discussion

Sterile neutrinos in the keV mass range have a decay time that is of cosmological interest. Recently a very interesting paper appeared (Abazajian et al. 2001b), where the signature from decaying sterile neutrinos in galaxies and clusters of galaxies was studied in detail 111One could imagine a slight change of strategy in the analysis of (Abazajian et al. 2001b), namely to consider regions with large concentration of dark matter, but with little baryonic matter. Such “dark blobs” may have been observed by inverting the matter distribution in clusters of galaxies from weak lensing, and comparing with the baryonic matter inferred from optical observations of the cluster (Clowe et al. 2000), but the significance of such blobs is far from being established.. A bound $`m<5`$ keV was derived, using the relation between mass and mixing angle obtained in (Abazajian et al. 2001a). We thus conclude that a sterile neutrino as WDM must lie in the mass range $`2.6`$ keV $`<m<5`$ keV, and it is therefore more interesting than ever before to search for a spectral line with energy $`E=m/2`$ from the decay $`\nu _s\nu _\alpha +\gamma `$.

If future searches for a spectral line from the sterile neutrino decay should give a negative result, then one must find new WDM candidates. An interesting possibility is an active neutrino, which may never reach thermal equilibrium if the reheat temperature at the end of inflation is low enough (Giudice et al. 2001). Such a scenario demands a reheat temperature of a few MeV, and the resulting neutrino distribution function is also warmer than a conventional WDM candidate. Such a solution does, however, not appear too natural in view of the recent neutrino data (Ahmad et al. 2001; Toshito et al. 2001), indicating that all the active neutrino masses are sub eV. Several WDM candidates have different bounds from what is presented here, e.g. gravitinos (Pagels & Primack 1982; Borgani, Masiero & Yamaguchi 1996; Kawasaki, Sugiyama & Yanagida 1997) produced in the very early Universe; or sterile neutrinos if there should exist an initial lepton asymmetry (Shi & Fuller 1999). Even more traditional DM candidates (like a neutralino) can disguise themselves as WDM, namely if they have scattering cross section with photons or neutrinos. Specifically, for a 100 GeV DM particle with scattering cross section about $`\sigma _{DM\gamma }10^{30}\text{cm}^2`$ one obtains a reduction of the matter power spectrum corresponding to a conventional WDM candidate with mass about 1 keV (Bœhm et al. 2001).

In conclusion, we have shown that the proper inclusion of the neutrino momentum distribution changes the lower bound allowed for the simplest models of sterile neutrinos as WDM, and the resulting allowed region thus becomes $`2.6`$ keV $`<m<5`$ keV.

## Acknowledgement

It is a pleasure for SHH and SP to thank Kev Abazajian, George Fuller, Kimmo Kainulainen and Dima Semikoz for discussions and comments. SHH and SP are supported by Marie Curie Fellowships of the European Commission under contracts HPMFCT-2000-00607 and HPMFCT-2000-00445. They acknowledge a visit to LAPTH, supported by CNRS, during which part of this work was carried out.